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)$. The command optim can be used to minimise a function of a vector, starting from a given value. From section 16.4 we know that a good starting value is a simple, consistent estimator, such as 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 command optim can be used to compute this matrix, if we ask for it explicitly with 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 ratio $\alpha/\beta$ is well determined by the data, but each parameter individually is determined only poorly. This is the problem of a nuisance parameter as discussed in lecture 14.

The gamma log-likelihood function has only one peak, so the choice of starting value matters little. 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 needs more than twice as many evaluations. When the likelihood has several local maxima, as in section 5.3, a poor start can end at the wrong peak, since an optimiser only finds a local maximum. Thus we always start from a consistent estimator, and we re-run optim from a different start to see whether the solution changes.

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$, it says nothing about moderate $n$. To check the coverage 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, and at $n=20$ both are slightly above it: the interval is a little too wide rather than too narrow.

The following task is in the style of a practical report. Use your student ID to seed the random number generator, simulate a gamma sample of size $30$ with a shape and rate of your choice, and fit the model by maximum likelihood. For each parameter, report the estimate and its standard error.

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.