[R] optim not giving correct minima

[Prof Brian Ripley]

Try the other methods that optim() supplies, and try supplying a 
analytical derviative (which looks easy enough).  On a problem like this I 
would expect to use BFGS with analytical derivatives.

If you don't want to do any of that, at least explore the control 
parameters.  It doesn't look to me as if you have attempted to scale the 
problem as the optim help page suggests you should.  Also, you say a 
`likelihood', but it is usual to maximize a log-likelihood. (Without 
knowing what you are trying to do in detail, I cannot tell if you are in 
fact using a log-likelihood.)

On Fri, 11 Nov 2005, Laura Forsberg White wrote:

> Hello,
>
> I am trying to use optim() on a function involving a summation.  My
> function basically is a thinned poisson likelihood.  I have two parameters
> and in most cases optim() does a fine job of getting the minima.  I am
> simulating my data based on pre specified parameters, so I know what I
> should be getting.  However when my true parameters fall in a particular
> range, optim() gives incorrect results.  I have generated a grid of
> parameter values and calculated my likelihood through those to see that the
> values that optim() gives is clearly not correct.  I can see that my
> likelihood does in fact have a unique maximum.  Any ideas why this might be?
>
> The data given below was generated such that the true parameters should be
> (0.3001, -1.8971).  Here is an example piece of data and the function:
>
> #function to maximize
> likN_alpha <- function(params,N){
>
>     thetas <- exp(params)
>
>     k <- length(thetas)
>     N <- c(rep(0,(k-1)),N)
>     l <- length(N)
>
>     lik <- 0
>
>     for(i in (k):(l-1)){
>         lambda <- thetas%*%N[i:(i-k+1)]
>         lik <- -lambda + N[i+1]*log(lambda) + lik
>     }
>
>     return(-lik)
> }
>
> # data to maximize over
> N <- c( 3, 3, 10, 19, 36, 54, 78,116,177, 265, 388, 598,
> 890,1328,1910,2736,3982,5908,8471,12440,17964,26207,37688,54795,79270,114752,166594,242438,
> 352753,512054)
>
> #optim() command
>
> optim(log(1.5*rep(1/2,2)),likN.alpha,N=N)
>
> # If I use constrained optimization and a slightly different
> parameterization, then the results are fine, at least in this case, but not
> always.
>
> Thanks for any help you might be able to offer!
>
> laura
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595