I use optim() in R to do a lot of MLE. I've used the approach for a lot of problems, but the one I'm working on right now consists of fitting the parameters of the generalized extreme value distribution (GEV). There are three parameters (location scale shape), and I want to get the standard errors from this fit. An example of a paper doing this is Katz et al. 2005, and here is the online code for the article. I am unaffiliated with this publication. I should also point out that the GEV question is just an example of the type of approach I take -- most of my MLE fits are time series models (state-space). I could give more examples if needed.

1) Which of the arguments to optim() are of particular relevance to the statistical validity of using the optim() Hessian to computer standard errors?

  • The first thing that comes to mind is the type of algorithm -- In R, optim() is a wrapper for several types of optimization algorithms (see above link for reference on optim()). I don't even think all of them return a Hessian.
  • I also think that the parameters relating to the reltol and abstol parameters could be important, but I don't fully understand what they do.

2) Is there some pre-treatment of the data (e.g., standardizing data before doing the MLE) that is advisable if I want to use the hessian to compute S.E.'s of the fitted parameters?

3) In general, of what do I need to be careful. I just don't want to incorrectly use this approach.

4) Is there anything you can think of that I definitely should not do?


  • 1
    $\begingroup$ One problem is likely to be when an MLE is on the boundary of the parameter space. $\endgroup$ May 6, 2013 at 12:22

1 Answer 1


Unless you have specified a function for computing the Hessian, optim() will return a numerical approximation which is obtained by taking differences. Depending on your function, this may actually yield a non-invertible Hessian (or other poor approximation), even if you are close to the maximum. You should experiment a bit with different values of 'ndeps' (in the 'control'-structure) to verify that your Hessian (or at least the S.E.s) is reasonably insensitive to the step size. Or use the numDeriv-package after optim() has found a solution, this uses a more elaborate method for computing a numerical Hessian.

The reltol and abstol parameters control how close you get to the maximum. Depending on your function, there may be problems with the Hessian if you are too far away. My experience is that MLEs need quite tight tolerances.

  • $\begingroup$ You guessed correctly that I don't have specified functions for computing the Hessian. Are any of the algorithms in optim() typically better at approximating the Hessian (e.g., Nelder Mead)? I like your tip about numDeriv, though. Seems like a solid idea. $\endgroup$
    – rbatt
    May 7, 2013 at 17:24
  • 2
    $\begingroup$ They all use the same method for computing the Hessian, the optimHess() function. What varies though, is how close to the maximum they get. But this depends on the function. BFGS and L-BFGS-B are generally quite good. For some functions, I've resorted to writing an analytic Hessian, and used a Newton method I wrote myself, but this isn't always feasible. For MLEs, an alternative is to use the Fisher matrix in place of the Hessian. $\endgroup$ May 8, 2013 at 22:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.