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

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.
 A: 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

