[R] optim not giving correct minima

Laura Forsberg White lforsber at hsph.harvard.edu
Fri Nov 11 21:59:04 CET 2005


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




More information about the R-help mailing list