Spritz - 1 year ago 60

R Question

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

Answer Source

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