How to compute (or numerically estimate) the standard error of the MLE

Question

I have a model for which I know the log likelihood function, the gradient of the log likelihood and the Hessian of the log likelihood. For given data I can compute the MLE using a generic optimizer (Nelder-Mead). How do I compute (or estimate) the standard error for the MLE?

If there is existing software that makes this easy that would be even better of course.

Answer 1

If you know the gradient and Hessian of the log-likelihood, you can write quick functions in R similar to the one you need for the LL itself. If you pass the gradient, you can use (L)BFGS in R as opposed to Nelder-Mead, which should converge a bit faster. Regardless, once you have the point of convergence, you can plug the values for the point of convergence into the function for the Hessian, and the sqrt of the diagonals is your estimated error. Here is an example using the Pareto distribution for which: $$ f(x) = \frac{\alpha\theta^{\alpha}}{\left(x+\theta\right) ^{\alpha + 1}} $$

LL <- function(pars, X){
  a <- pars[[1]]
  q <- pars[[2]]
  return(-sum(a * log(q) + log(a) - (a + 1) * log (X + q)))
}

LLG <- function(pars, X){
  a <- pars[[1]]
  q <- pars[[2]]
  ga <- -sum(log(q) + 1 / a - log(X + q))
  gq <- -sum(a / q - (a + 1) / (X + q))
  Z <- c(ga, gq)
  names(Z) <- c('a', 'q')
  return(Z)
}  

LLH <- function(pars, X){
  a <- pars[[1]]
  q <- pars[[2]]
  n <- length(X)
  haa <- n / a ^ 2
  hqq <- n * a / q ^ 2 - sum((a + 1) / (X + q) ^ 2)
  haq <- hqa <- sum(1 / (X + q)) - n / q
  Z <- matrix(c(haa, hqa, haq, hqq), ncol = 2)
  rownames(Z) <- colnames(Z) <- c('a', 'q')
  return(Z)
}

I tend to use nloptr for linear-search optimization, the call for `optim' would be similar. So assuming your data is stored as DATA:

Fit <- nloptr(x0 = c(2, 1e6), eval_f = LL, eval_grad_f = LLG, lb = c(0,0), X = DATA,
              opts = list(algorithm = "NLOPT_LD_LBFGS", maxeval = 1e5))

Your values are in Fit$solution so your Fisher information matrix estimate is the inverse of the Hessian (not negative Hessian since we are minimizing NLL, not maximizing LL) and so the estimate of the standard error can be calculated using:

sqrt(diag(solve(LLH(Fit$solution, DATA))))

and their correlation would be the off-diagonal in:

cov2cor(solve(LLH(Fit$solution, DATA)))

Answer 2

There is a very nice discussion on this topic in CV:

Basic question about Fisher Information matrix and relationship to Hessian and standard errors

You can use the R package numDeriv to approximate the Hessian.