[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