Lab 10

Due Thursday April 9 at 11:59 PM

This lab provides extra practice with the kind of maximum likelihood problem that you will encounter on the last problem set and the final exam.

Task 0 (the TA does it)

Let’s see a worked example first. Imagine we observe data

\[ X_1,\,X_2,\,...,\,X_n\overset{\text{iid}}{\sim}\text{N}(\theta, 1), \]

where the true value of the mean \(\theta>0\) is unknown. We will use the method of maximum likelihood to propose an estimator for \(\theta\). Spoiler alert: In this case, it’s the sample mean \(\bar{X}_{n}\)! See below for the derivation.

We consider data

\[ X_1,\,X_2,\,...,\,X_n\overset{\text{iid}}{\sim}\text{N}(\theta, 1). \] Recall that the associated density is given by \[f(x;\ \theta) = \frac{1}{\sqrt{2\pi}}e^{-\frac{(x - \theta)^{2}}{2}}, \quad x \in \mathbb{R}.\] So, the joint likelihood function is given by \[ \begin{aligned} L\left(\theta;\ X_{1:n}\right) & = \prod_{i = 1}^{n}{f\left(X_{i};\ \theta\right)} \\ & = \prod_{i = 1}^{n}{\frac{1}{\sqrt{2\pi}}e^{-\frac{\left(X_{i} - \theta\right)^{2}}{2}}} \\ & = (2\pi)^{-\frac{n}{2}}\prod_{i = 1}^{n}{e^{-\frac{\left(X_{i} - \theta\right)^{2}}{2}}} \\ & = (2\pi)^{-\frac{n}{2}}e^{-\frac{1}{2}\sum_{i = 1}^{n}{\left(X_{i} - \theta\right)^{2}}}. \end{aligned} \] The log-likelihood is then given by \[ \begin{aligned} \ell\left(\theta;\ X_{1:n}\right) & = \ln L\left(\theta; \ X_{1:n}\right) \\ & = \ln\left[(2\pi)^{-\frac{n}{2}}e^{-\frac{1}{2}\sum_{i = 1}^{n}{\left(X_{i} - \theta\right)^{2}}}\right] \\ & = \ln\left[(2\pi)^{-\frac{n}{2}}\right] + \ln\left[e^{-\frac{1}{2}\sum_{i = 1}^{n}{\left(X_{i} - \theta\right)^{2}}}\right] \\ & = -\frac{n}{2}\ln(2\pi) - \frac{1}{2}\sum_{i = 1}^{n}{\left(X_{i} - \theta\right)^{2}}. \end{aligned} \] We want to find \[\hat{\theta}_{n} = \underset{\theta}{\arg\max}{\ L\left(\theta;\ X_{1:n}\right)} = \underset{\theta}{\arg\max}{\ \ell\left(\theta;\ X_{1:n}\right)}.\] We optimize by taking the derivative with respect to \(\theta\), setting it equal to zero, and solving: \[ \begin{aligned} \frac{\text{d}\ell}{\text{d}\theta} & = \frac{\text{d}}{\text{d}\theta}\left[-\frac{n}{2}\ln(2\pi) - \frac{1}{2}\sum_{i = 1}^{n}{\left(X_{i} - \theta\right)^{2}}\right] \\ & = -\frac{1}{2}\sum_{i = 1}^{n}{\frac{\text{d}}{\text{d}\theta}\left[\left(X_{i} - \theta\right)^{2}\right]} \\ & = -\frac{1}{2}\sum_{i = 1}^{n}{-2\left(X_{i} - \theta\right)} \\ & = \sum_{i = 1}^{n}{X_{i}} - n\theta \\ \Rightarrow 0 & = \sum_{i = 1}^{n}{X_{i}} - n\hat{\theta}_{n} \\ \Rightarrow n\hat{\theta}_{n} & = \sum_{i = 1}^{n}{X_{i}} \\ \Rightarrow \hat{\theta}_{n} & = \frac{1}{n}\sum_{i = 1}^{n}{X_{i}} \\ \Rightarrow \hat{\theta}_{n} & = \bar{X}_{n}. \end{aligned} \] For the sake of completeness, we can check the second derivative to verify that this is, in fact, a maximum: \[ \begin{aligned} \frac{\text{d}^{2}\ell}{\text{d}\theta^{2}} & = \frac{\text{d}}{\text{d}\theta}\left[\frac{\text{d}\ell}{\text{d}\theta}\right] \\ & = \frac{\text{d}}{\text{d}\theta}\left[\sum_{i = 1}^{n}{X_{i}} - n\theta\right] \\ & = -n. \end{aligned} \] Since \(\frac{\text{d}^{2}\ell}{\text{d}\theta^{2}} = -n < 0\), we know that \(\hat{\theta}_{n} = \bar{X}_{n}\) is the maximum likelihood estimator.

TipFun fact

For any known variance \(\sigma^{2} > 0\), the MLE of the mean for iid \(\text{N}\left(\theta, \sigma^{2}\right)\) data is \(\hat{\theta}_{n} = \bar{X}_{n}\). Try that out for yourself!

Task 1 (you do it)

We now pivot to our distribution of choice for this lab. If a non-negative random variable \(X\) has the Rayleigh distribution with parameter \(\theta>0\), then its CDF is

\[ F(x;\,\theta)=1-\exp\left(-\frac{x^2}{2\theta}\right)\quad x\geq 0, \]

and we say \(X\sim\text{Rayleigh}(\theta)\). What is the PDF in this parametric family?

Task 2 (we do it)

Imagine we observe data

\[ X_1,\,X_2,\,...,\,X_n\overset{\text{iid}}{\sim}\text{Rayleigh}(\theta), \]

where the true value of the parameter \(\theta>0\) is unknown. Use the method of maximum likelihood to propose an estimator for \(\theta\).

The likelihood function is

\[ \begin{aligned} L(\theta;\,X_{1:n}) &= \prod_{i=1}^n f(X_i;\,\theta) \\ &= \prod_{i=1}^n \frac{X_i}{\theta} \exp\left(-\frac{X_i^2}{2\theta}\right) \\ &= \theta^{-n}\left(\prod_{i=1}^nX_i\right)\left[\prod_{i=1}^n\exp\left(-\frac{X_i^2}{2\theta}\right)\right] \\ &= \theta^{-n}\left(\prod_{i=1}^nX_i\right)\exp\left(-\frac{1}{2\theta}\sum\limits_{i=1}^nX_i^2\right). \end{aligned} \]

The log-likelihood function is

\[ \begin{aligned} \ell(\theta;\,X_{1:n}) &= \ln L(\theta;\,X_{1:n}) \\ &= \ln\left[\theta^{-n}\left(\prod_{i=1}^nX_i\right)\exp\left(-\frac{1}{2\theta}\sum\limits_{i=1}^nX_i^2\right)\right] \\ &= \ln(\theta^{-n})+\ln\left(\prod_{i=1}^nX_i\right)+\ln\left[\exp\left(-\frac{1}{2\theta}\sum\limits_{i=1}^nX_i^2\right)\right] \\ &= -n\ln \theta+\sum\limits_{i=1}^n\ln X_i-\frac{1}{2\theta}\sum\limits_{i=1}^nX_i^2 . \end{aligned} \]

To maximize this, we take the second derivative, set it equal to zero, and solve:

\[ \begin{aligned} \frac{\text{d}}{\text{d}\theta}\ell(\theta;\,X_{1:n}) = - \frac{n}{\theta} + \frac{1}{2\theta^2} \sum\limits_{i=1}^nX_i^2 &= 0 \\ \frac{1}{2\theta^2} \sum\limits_{i=1}^nX_i^2 &= \frac{n}{\theta} \\ \frac{1}{2n} \sum\limits_{i=1}^nX_i^2 &= \theta. \end{aligned} \]

So the maximum likelihood estimator for the Rayleigh distribution is

\[ \hat{\theta}_n=\frac{1}{2n} \sum\limits_{i=1}^nX_i^2. \]

Task 3 (we do it)

Your estimator in the previous part is a random variable that inherits its randomness from the \(X_i\). Derive the full sampling distribution of your estimator.

For these problems, I work “from the inside out.” What’s the distribution of \(X_i\)? What’s the distribution of \(X_i^2\)? What’s the distribution when you add those up? What’s the distribution when you scale the sum by \(1/2n\)? To answer those questions, you need tools like the change-of-variables formulas and our MGF updating rules for linear combinations.

We already know that \(X_1\sim\text{Rayleigh}(\theta)\), so we derive the distribution of \(Y_1=g(X_1)=X_1^2\) using change of variables. Since \(\text{Range}(X_1)=(0,\,\infty)\), the transformation is strictly increasing, so we apply the change of variables formula:

\[ \begin{aligned} f_{Y_1}(y) &= f\left(g^{-1}(y);\,\theta\right)\left|\frac{\text{d}}{\text{d}y}g^{-1}(y)\right| \\ &= f\left(y^{1/2};\,\theta\right)\left|\frac{\text{d}}{\text{d}y}y^{1/2}\right| \\ &= \frac{y^{1/2}}{\theta} \exp\left(-\frac{(y^{1/2})^2}{2\theta}\right) \left|\frac{1}{2}y^{-1/2}\right| \\ &= \frac{y^{1/2}y^{-1/2}}{2\theta} \exp\left(-\frac{y}{2\theta}\right) \\ &= \frac{1}{2\theta} \exp\left(-\frac{y}{2\theta}\right) ,&&y>0. \end{aligned} \]

This is the density of \(\text{Exponential}(1/2\theta)\), otherwise known as \(\text{Gamma}(1,\,1/2\theta)\). Given this, we have the following:

\[ \begin{aligned} X_i\overset{\text{iid}}{\sim}\text{Rayleigh}(\theta) \quad &\implies \quad X_i^2\overset{\text{iid}}{\sim}\text{Gamma}(1,\,1/2\theta) \\ &\implies \sum\limits_{i=1}^nX_i^2\sim\text{Gamma}(n,\,1/2\theta) \\ &\implies \frac{1}{2n}\sum\limits_{i=1}^nX_i^2\sim\text{Gamma}(n,\,n/\theta). \end{aligned} \] The last implication follows from Problem Set 5 #10d. So, the full sampling distribution of the estimator is

\[ \hat{\theta}_n\sim \text{Gamma}(n,\,n/\theta). \]

Task 4 (we do it)

Now that you know the sampling distribution of the estimator, compute its mean squared error (MSE) using the bias-variance decomposition. Is the estimator unbiased? Is it consistent?

Using the mean and variance formulas for the gamma distribution, we see that

\[ \begin{aligned} E(\hat{\theta}_n)&=\frac{n}{n/\theta}=\theta\\ \text{var}(\hat{\theta}_n)&=\frac{n}{(n/\theta)^2}=\frac{\theta^2}{n}. \end{aligned} \]

Since \(E(\hat{\theta}_n)=\theta\), the bias of the estimator is zero, and so the mean-squared-error is simply \(\text{var}(\hat{\theta}_n)=\theta^2/n\). As \(n\to\infty\), this goes to zero. As such, the estimator is unbiased and consistent.

Task 5 (you do it)

Compute the quantile function of the Rayleigh distribution and use it to generate 100 draws from the distribution that has \(\theta=4\).

Task 6 (you do it)

Write a for loop that simulates the sampling distribution for your estimator. Let \(\theta_{0} = 4\). Here is a schematic of what you are doing:

\[ \begin{matrix} \text{0. True parameter value} &&& \theta_0 && \\ &&& \downarrow && \\ \text{1. True distribution} &&& \text{Rayleigh}(\theta_0) && \\ &\swarrow &\swarrow& \cdots &\searrow&\searrow \\ \text{3. Simulated data}&x_{1:n}^{(1)} &x_{1:n}^{(2)}& \cdots &x_{1:n}^{(M-1)}&x_{1:n}^{(M)} \\ &\downarrow &\downarrow& \cdots &\downarrow&\downarrow \\ \text{4. Different estimates}&\hat{\theta}_n^{(1)} &\hat{\theta}_n^{(2)}& \cdots &\hat{\theta}_n^{(M-1)}&\hat{\theta}_n^{(M)} \\ &\searrow &\searrow& \cdots &\swarrow&\swarrow \\ & && \text{Histogram} && \\ \end{matrix} \]

At the end, add on top of the histogram a curve of the density that you derived in Task 3. They ought to match.

TipTemplate

Here is some template code you can fill in:

set.seed(123)
true_theta <- 4      # true value of parameter
n <- 25                 # sample size
M <- 2500               # how many repetitions
estimates <- numeric(M) # preallocate storage for the sample of estimates

for(m in 1:M){
  # Step 1: simulate a new fake dataset of size n
  # Step 2: apply your formula from Task 2 to compute new estimate
  # Step 3: store it
}

# Step 4: plot histogram of estimates
# Step 5: add density curve of exact sampling distribution