# Standard error from Hessian matrix when likelihood is used (rather than Ln L)

I understand that at MLE point, the inverse of the Hessian matrix can be used as approximation of V-Cov matrix:

Llikelihood <- function(par, x) {
return(sum(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = TRUE)))
}
o <- optim(par=c(mean=3, sd=4), fn=Llikelihood, x=x, hessian = TRUE, control=list(fnscale=-1))
o$$hessian sqrt(diag(solve(-o$$hessian)))

It works well. I read also that if I use the Likelihood itself (not Ln L), I must divide the hessian by the likelihood:

likelihood <- function(par, x) {
return(prod(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = FALSE)))
}
o <- optim(par=c(mean=3, sd=4), fn=likelihood, x=x, hessian = TRUE, control=list(fnscale=-1))
-o$$hessian/o$$value
sqrt(diag(solve(-o$$hessian/o$$value)))

However the standard errors are not the same:

mean        sd
0.1974196 0.1396388 For Ln L

and

mean        sd
0.4008885 0.3632539  For L

I would be interested to have some advice to understand what's happened in this context.

Thanks

PS I know that Hessian and Standard errors have been discussed several time, but I don't find an answer of my question.

• You appear to have demonstrated that what you read is incorrect. What kind of advice are you looking for beyond that? – whuber Dec 14 '18 at 20:35
• What is the point of your condescending response? If at least you pointed the place that you think I read wrongly! But no, in your answer you do not bring any element. And you ask me to clarify my question! I can answer you in the same way as what you do: just read the initial message! And we are well advanced! If you prefer: How to use the Hessian matrix when likelihood is used for MLE rather than Ln Likelihood? – MarcG Dec 15 '18 at 2:46
• I'm sorry, but could you tell me what the "initial message" might be? BTW, it's my job as moderator to request clarification of posts that might be ambiguous or overly broad. It's not considered condescending and if you read it that way, I apologize. But please read over the material in our help center so you can understand how this site works. – whuber Dec 15 '18 at 16:37
• The "initial" message is the message appearing at the beginning of this post. Better to not take too much time on this post. The answer by Warren Weckesser solve my question. Thanks – MarcG Dec 17 '18 at 5:25

In theory, you are correct, the two computations should produce the same result. Here's a brief explanation.

Define

$$l(x) = \ln L(x)$$

then, using ' for differentiation,

$$l'(x) = \frac{L'(x)}{L(x)}$$

and

$$l''(x) = \frac{L''(x)}{L(x)} - \left(\frac{L'(x)}{L(x)}\right)^2$$

At a critical point $$x_0$$, $$L'(x_0)$$ is 0, so

$$l''(x_0) = \frac{L''(x_0)}{L(x_0)}$$

(The above notation implies $$x$$ is a scalar, but the same idea applies when $$x$$ is a vector.)

So why didn't the two versions of your code produce the same results?

The problem seems to be numerical error in the computations, which can be pretty high. The parameters returned by optim are not exact, and the Hessian is computed using finite differences, which are very susceptible to numerical errors. (Numerical differentiation tends to amplify errors in the inputs.)

I was able to get the results to agree by making several changes to the use of optim:

• use a better estimate for the initial guess;
• use method="CG" along with the control maxit=10000;
• use reltol=1e-12;
• use fnscale=-likelihood(p0, x) when applying optim to likelihood.

I stopped tweaking the code when I got pretty good agreement between the standard errors computed using the two methods. I haven't gone back to figure out which of the changes are most important.

Here's a self-contained R script to demonstrate, followed by its output.

x = c(1, 3, 4, 5, 7.5, 10, 12, 23, 39, 40)

# The *exact* MLE for the normal distribution is
#    mu = mean(x),  sigma = sqrt(var(x)*(length(x) - 1)/length(x))

p0 = c(mean=mean(x), sd=sd(x))

Llikelihood <- function(par, x) {
return(sum(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = TRUE)))
}
o1 <- optim(par=p0, fn=Llikelihood, x=x, hessian=TRUE, method="CG", control=list(fnscale=-1, reltol=1e-12, maxit=10000))

cat("\n*** Using log-likelihood ***\n\n")
cat("Estimated parameters:\n")
o1$$par cat("\nEstimated Hessian:\n") o1$$hessian
cat("\nEstimated standard errors:\n")
sqrt(diag(solve(-o1$hessian))) cat("\n*** Using (plain) likelihood ***\n\n") likelihood <- function(par, x) { return(prod(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = FALSE))) } o2 <- optim(par=p0, fn=likelihood, x=x, hessian=TRUE, method="CG", control=list(fnscale=-likelihood(p0, x), reltol=1e-12, maxit=10000)) cat("Estimated parameters:\n") o2$$par cat("\nEstimated Hessian of neg. log likelihood:\n") o2$$hessian/o2$$value cat("\nEstimated standard errors:\n") sqrt(diag(solve(-o2$$hessian/o2$value)))

Output of the script:

*** Using log-likelihood ***

Estimated parameters:
mean       sd
14.45000 13.83194

Estimated Hessian:
mean         sd
mean -0.05226777  0.0000000
sd    0.00000000 -0.1045355

Estimated standard errors:
mean       sd
4.374043 3.092915

*** Using (plain) likelihood ***

Estimated parameters:
mean       sd
14.45000 13.83194

Estimated Hessian of neg. log likelihood:
mean            sd
mean -5.226776e-02  8.106461e-11
sd    8.106461e-11 -1.045355e-01

Estimated standard errors:
mean       sd
4.374043 3.092915

From what you know and as @WarrenWeckesser pointed out, in theory those two approaches are equivalent. From Yogi Berra (and others): "In theory there is no difference between theory and practice. In practice there is."

While the structure of your code clearly shows you know how to do this, there wasn't an MWE: minimal working example. x was missing and we didn't know how it was generated. Although certainly a guess would be from a normal distribution with mean 3 and standard deviation 4 - but we were given no idea of the sample size. And that is critical. Almost all of such issues are associated with round-off error or objective functions being much, much smaller than the default stopping-rule thresholds.

When the sample size is small (say n = 10), then there is good agreement without having to make any "adjustments".

> set.seed(12345)
> x = rnorm(10,3,4)
>
> Llikelihood <- function(par, x) {
+    return(sum(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = TRUE)))
+ }
> o <- optim(par=c(mean=3, sd=4), fn=Llikelihood, x=x, hessian = TRUE, control=list(fnscale=-1))
> sqrt(diag(solve(-o$$hessian))) mean sd 0.9764625 0.6905435 > > likelihood <- function(par, x) { + return(prod(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = FALSE))) + } > o <- optim(par=c(mean=3, sd=4), fn=likelihood, x=x, hessian = TRUE, control=list(fnscale=-1)) > sqrt(diag(solve(-o$$hessian/o$value))) mean sd 0.9768148 0.6911671 But with a sample size of 100 we get two different results: mean sd 0.4436389 0.3136806 mean sd 0.2787466 0.1555169 From the use of the log likelihood we have for the result: > o$par
mean       sd
3.981301 4.436389

$value [1] -290.8821$counts
43       NA

$convergence [1] 0$message
NULL

$hessian mean sd mean -5.080900934 0.001172538 sd 0.001172538 -10.163067714 And from the likelihood we have > o$par
mean   sd
3.4  4.0

$value [1] 5.173568e-128$counts
3       NA

$convergence [1] 0$message
NULL

\$hessian
mean            sd
mean 3.583380e-127 1.085842e-126
sd   1.085842e-126 1.151216e-126

Notice the wild difference in the exponents and that there are only 3 function evaluations for the likelihood function and very round final numbers. Clearly the likelihood approach is stopping too soon because of the smallness of the objective function compared to the thresholds used to decide when to stop.

That can be moderated by multiplying the likelihood (in this case) by 10^128:

> likelihood <- function(par, x) {
+  return(10^128 * prod(dnorm(x, mean=par["mean"], sd=abs(par["sd"]), log = FALSE)))
+ }
> o <- optim(par=c(mean=3, sd=4), fn=likelihood, x=x, hessian = TRUE, control=list(fnscale=-1))
> sqrt(diag(solve(-o$$hessian/o$$value)))
mean        sd
0.4436584 0.3137145

Now we get nearly the identical values for the standard errors of the mean and standard deviation as with the log likelihood approach.

So the morals are: (1) Give MWE's, (2) look at optim's complete output, and (3) question default stopping rules.