2
$\begingroup$

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

Note the - sign in front of the 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.

$\endgroup$
  • $\begingroup$ You appear to have demonstrated that what you read is incorrect. What kind of advice are you looking for beyond that? $\endgroup$ – whuber Dec 14 '18 at 20:35
  • $\begingroup$ 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? $\endgroup$ – MarcG Dec 15 '18 at 2:46
  • $\begingroup$ 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. $\endgroup$ – whuber Dec 15 '18 at 16:37
  • $\begingroup$ 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 $\endgroup$ – MarcG Dec 17 '18 at 5:25
3
$\begingroup$

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
$\endgroup$
0
$\begingroup$

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
function gradient 
      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
function gradient 
       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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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