Raw Message

Message-ID: <5.2.0.9.2.20051111152820.00d0eff8@hsph.harvard.edu>
Date: 2005-11-11T20:59:04Z
From: Laura Forsberg White
Subject: optim not giving correct minima

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