i'm trying to write a code in r in order to find the maximum likelihood(not the log-likelihood) values for a univariate normal distribution. I know that there are other methods but a deep understand of numerical optimization is needed for me for further works. When i call the 'optim' function it seems it doesn't iterate at all and return the values that i passed as initial parameters. This doesn't happen if i pass to the optimizer the function that instead computes the log-likelihood. Any idea why? I cannot see where my error is. I could just say that maybe the product of the densities is too much close to the zero and the calculator can't handle it. Here is my code.Thanks a lot!
set.seed(123);
x=rnorm(10, mean = 2, sd = 5)
LikeNormUnivar<-function(param,data){
mu=param[1];
sdev=param[2];
densityvector=dnorm(data, mean = mu, sd = sdev, log = FALSE)
like=prod(densityvector)
return(-like)
}
theta.start = c(2,4)
ans = optim(par=theta.start, fn=LikeNormUnivar, data=x,control=list(trace=TRUE),
method="BFGS")
ans$par
I figured out what was going on by adding the line cat(mu,sdev,like,"\n") at an appropriate place in the function to see what was going on. Basically, on the scale on which the BFGS is estimating derivatives by finite differences, there's not enough variation. Setting retol=1e-16 in the control list works. Better, try minimizing the negative log likelihood ... and fitting on the log-standard-deviation scale is also a good idea, e.g.
LikeNormUnivar <- function(param,data){
mu=param[1]
sdev=exp(param[2])
loglik=dnorm(data, mean = mu, sd = sdev, log = TRUE)
return(-sum(loglik))
}