According to various sources, the variance of ML estimators can be obtained from the Hessian matrix of the likelihood function. If $H$ is the Hessian of the negative log-likelihood function, then $H^{-1}$ is the variance-covariance matrix of the ML estimators.

Knowing that, we can do a simple test. We can calculate the variance of a ML estimate of a Poisson parameter using this method, and then compare the result with the correct answer that we know because, in this case, it has a closed form solution. In this way, we can see if the Hessian method is sufficiently accurate.

So, here is the code in R:

x <- rpois(500, 5)

log.likelihood <- function(lambda, x) {
  sum(x * log(lambda) - lambda - log(factorial(x)))

mle <- optim(4.5, function(lambda)
             -log.likelihood(if (!lambda > 0) 0.01 else lambda, x),
             method='BFGS', hessian=TRUE)

Then we can check whether $H^{-1}$ is equal to $\frac{1}{n}\text{var}(x)$ which is the actual variance of a ML Poisson parameter estimate.

> 1 / mle$hessian
[1,] 0.01005594
> var(x) / length(x)
[1] 0.01032308

In conclusion, it looks like while the Hessian method provides an approximation to the variance of the parameter, the approximation is somewhat off. Any suggestion about how the accuracy could be improved?


Here is a plot of the errors made by different variance estimation methods for various sample sizes. The errors appear to decrease rapidly as the sample size grows. However, when the sample is small, errors can become substantial.

plot The hs method uses the Hessian produced by the optimization algorithm. The richardson and complex methods use Hessian matrices estimated independently, using the hessian() function of the numDeriv package, and both produce the same result.

  • 1
    $\begingroup$ No, I mean to improve the accuracy for any sample size. Bootstrap can be an alternative in certain circumstances, but it's not feasible when the ML estimation is computationally expensive. E.g., if it takes 30 secs to estimate the model, 10000 repetitions would take over 80 hours! $\endgroup$ – Ernest A Nov 10 '12 at 17:01
  • 1
    $\begingroup$ Use the same method on log-lambda, not lambda, and the asymptotic accuracy should work better. Also note that for this example you can do the whole computation using summary(glm(x~1, family=poisson)). $\endgroup$ – guest Nov 11 '12 at 6:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.