Lecture 18 Maximum Likelihood in R

The model of the day is $X_{1},\dots,X_{n}$ i.i.d. from the gamma distribution with shape $\alpha$ and rate $\beta$, both unknown. Example 4.9 gave the method of moments estimates $\tilde{\alpha}=\bar{X}^{2}/V$ and $\tilde{\beta}=\bar{X}/V$, with $V$ the plug-in variance, and section 5.3 noted that the likelihood equations have no closed-form solution. We simulate one sample of size $n=40$ with $\alpha=2.5$ and $\beta=0.5$. The command dgamma with log = TRUE can be used to evaluate the logarithm of the density at every observation, and the sum of these values is $\ell(\alpha,\beta)$. Since the R optimisers minimise, the function returns the negative log-likelihood, and since they hand over all parameters as one vector, it takes the vector par, shape first and rate second.

The true parameter values lead to a much smaller negative log-likelihood than the pair $(1,1)$. Here we assume that we know the rate $\beta=0.5$ and that we want to estimate the shape parameter $\alpha$. The command optimise can be used to find the minimum of a function of one variable in a given interval, without requiring a starting value. Before you try the command, make a sketch of the negative log-likelihood as a function of $\alpha$ on the interval from $0.1$ to $20$. How many dips will there be? Will the minimum lie above or below $2.5$?

The estimate $\hat{\alpha}=2.62$ is a little above the true value, and the minimum $101.9$ lies below the value $102.1$ where the truth is located. Since both parameters are now unknown, we have to use optim, which minimises a function of a vector, and which requires a starting value. From section 16.4 we know that a good way to find a starting value is to use a simple, consistent estimator. One possible choice for this is the method of moments estimate. Before you try the code, predict whether the maximum likelihood estimates will agree with the method of moments estimates to one decimal place.

The extra argument data = x is passed on to our function. The result fit is a list: fit$par holds the estimates, $\hat{\alpha}=2.26$ and $\hat{\beta}=0.42$, which agree with the start to one decimal place; fit$convergence is $0$ if the search finished properly, and any other value means the estimates cannot be trusted; fit$counts records that $43$ evaluations were needed. By default optim uses the Nelder–Mead method, which needs no derivatives; appendix B shows the alternative "BFGS".

From theorem 17.1 and the vector version of this result from lecture 17 we know that the MLE is approximately normally distributed with covariance matrix $\mathcal{I}_{n}(\theta_{0})^{-1}$. As in section 17.3, we can replace $\mathcal{I}_{n}(\hat{\theta})$ by the observed information, i.e. the matrix of second derivatives of the negative log-likelihood at the maximum. The optim function can compute this matrix for us if we ask for it explicitly, using hessian = TRUE. The matrix can be inverted using the solve command, the diag command gives the diagonal of this matrix, and the square roots of the diagonal elements are the standard errors. Before you look at the output, predict which of the two standard errors will be larger and by how much.

The table prints the results of a maximum likelihood estimate in an easy to read format. Each row corresponds to one parameter, and shows the parameter estimate, the standard error and the interval $\hat{\theta}\pm 1.96\,\mathop{\mathrm{se}}\nolimits(\hat{\theta})$ from corollary 17.3. Both intervals capture the true parameter values. The standard error for the shape parameter is approximately five times larger than the standard error for the rate parameter, but, since the estimates differ by a similar factor, the two standard errors are comparable relative to the estimates. The final row of the table, computed using the cov2cor function, shows that the two estimates have correlation $0.89$. It is easy to show that, if we increase the shape and the rate simultaneously, while keeping the mean $\alpha/\beta$ fixed, the log-likelihood changes only slightly. Thus the data constrains the ratio $\alpha/\beta$ well, but poorly constrains each parameter individually. This is the problem of a nuisance parameter as discussed in lecture 14.

The gamma log-likelihood function has only one peak, so the method is forgiving. Before you try the code, make a prediction: what happens if we start with very different values, e.g. with shape $30$ and rate $30$? Will optim still find the maximum? Will it take more than $43$ evaluations to find the maximum, when we used the “good” start?

The search reaches the same maximum, up to two decimal places, but it takes more than twice as long and R complains about several warnings and the presence of NaNs (not a number) in the output: the search tried to use negative parameter values where dgamma is not defined. We can prevent this by adding a line at the top of the function which returns Inf if any(par <= 0). This will prevent the optimiser from trying these values. For a “good” start this is not required. The real problem of a poor start is when the likelihood has several local maxima, as in section 5.3. We illustrate this effect using a toy function which is not a model: the logarithm of a two-component normal mixture density, evaluated at $\theta$, has a local maximum at $\theta=2$ and a lower local maximum at $\theta=-2$ (see figure 18.1). Which of these two peaks will optimise report for the interval from $-6$ to $4$, which contains both maxima? Which optim result will you get if you start at $-3$ and at $3$, respectively? To find the maximum, optim needs to use a different method than the default method, and we can use the BFGS method. Also, we need to tell optimise that we want to find a maximum by using the argument maximum = TRUE.

Both optimisers find a local maximum, but neither of them can determine whether this is the global maximum: optim starts at and ends at the first peak it encounters, and optimise assumes that there is only one peak. There is no command which allows us to fix this problem. The only way to avoid this problem is to use some theory and to follow certain habits: always start with a consistent estimator (as we did for the gamma distribution earlier in this lecture), always plot the log-likelihood function (if the number of parameters is one), and always start the optim procedure again to see whether the solution is different.

While corollary 17.3 shows that the interval $\hat{\theta}\pm 1.96\,\mathop{\mathrm{se}}\nolimits(\hat{\theta})$ covers the truth with probability approaching $0.95$, in section 17.5 we have seen that sometimes, for moderate $n$, this is not the case. To verify the claim for the gamma shape and rate, we use simulation in this section, following the approach to coverage introduced in lecture 12. The complete script for our simulation is given in the following box, with the lines shuffled (except for the bodies of the function and of the loop). Before you read the rest of this section, try to do the following: (1) write down the lines in a way which forms a working script, and (2) predict the coverages for $n=20$: will they be above or below $0.95$? Which of the two, shape or rate, will be worse?

Every name must exist before it is used: the function and the four assignments come first, with coverage after ns; then the seed, then the loop, and the table last. The one new detail is the return value, the two estimates followed by the two standard errors, which replicate stacks as the columns of a matrix; thus res[1:2, ] holds the estimates and res[3:4, ] the standard errors. The script in working order is R/S18-coverage.R; it fits the model $5000$ times for each sample size, takes a few seconds, and prints the following table.

At $n=100$ both coverages are within simulation error of the nominal $0.95$, as expected from corollary 17.3. At $n=20$ both coverages are above the nominal level, around $0.96$: the interval is slightly too wide rather than too narrow. The total hides an imbalance again, and the script also counts for the shape at $n=20$ how many misses are below the truth and how many are above.

Almost three quarters of the misses are below the true shape. Since the standard error increases with the estimate, an underestimated shape will result in a short interval which does not reach the truth, whereas an overestimated shape will result in a long interval which usually still covers the truth. The same mechanism was at work for the exponential rate in lecture 12.

The following task is in the style of a practical report. The published scripts of this lecture can be used as a starting point. Use your student ID to seed the random number generator, and choose your own shape and rate. Simulate one gamma sample of size $n=30$. Use optim to fit the model by maximum likelihood, starting with estimates from the method of moments. For each parameter, report the estimate, the standard error, and the interval $\hat{\theta}\pm 1.96\,\mathop{\mathrm{se}}\nolimits(\hat{\theta})$, in three separate sentences. For each interval, state whether the interval contains the value you chose. Estimate the coverage of both intervals when $n=30$, using $2000$ replicates. Write a short paragraph about whether you would trust the intervals for this sample size. Finally, fit the normal model with unknown mean and variance, using optim, to a normal sample of size $30$. Check that the estimates match the closed forms from example 5.4, and that the standard error of the mean matches $\hat{\sigma}/\sqrt{n}$.

Before the next session, you should type the commands used in this session, and should check that you get the same results as the text.