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

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)))