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

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.

• You can get the standard error for the parameters of the MLE via the Hessian but not a standard error for the MLE. – Henrik Aug 31 '14 at 12:29
• @Henrik Thank you. That would be a good start if you wouldn't mind giving some details (maybe as an answer?). In practice I am looking for a confidence interval for the MLE. How would one get that? – felix Aug 31 '14 at 14:01

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}}{x+\theta}$$

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

LLG <- function(pars, X){
a <- pars[]
q <- pars[]
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[]
q <- pars[]
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)))

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.

• We try to avoid link-only answers (which are more of a comment). See this page of the help. If you think a linked answer here completely answers a question, consider whether you should indicate that it is a duplicate. – Glen_b Aug 31 '14 at 21:36
• This is really more of a comment than an answer. When your reputation is >50 (presumably soon), you will be able to comment anywhere, but you should not post comments as answers before then. Can you make this into more of an answer? – gung Aug 31 '14 at 21:36