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.

**Contents**hide

#### Best Answer

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{alphatheta^{alpha}}{left(x+thetaright) ^{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))) `

### Similar Posts:

- Solved – Calculate the uncertainty of a MLE
- Solved – When should I *not* use R’s nlm function for MLE
- Solved – Variance-covariance matrix of the parameter estimates wrongly calculated
- Solved – Weibull MLE: what is the method/algorithm used to perform the optimization
- Solved – Defining gradient function argument in optim function-R