[R] optim() not finding optimal values

[Rubén Roa]

Derek,

As a general strategy, and as an alternative to parscale when using optim, you can just estimate the logarithm of your parameters. So in optim the par argument would contain the logarithm of your parameters, whereas in the model itself you would write exp(par) in the place of the parameter.

The purpose of this is to bring all parameters to a similar scale to aid the numerical algorithm in finding the optimum over several dimensions.

Due to the functional invariance property of maximum likelihood estimates your transformed pars back to the original scale are also the MLEs of the pars in your model. If you were using ADMB you'd get the standard error of the pars in the original scale simply by declaring them sd_report number class. With optim, you would get the standard error of pars in the original scale post-hoc by using Taylor series (a.k.a. Delta method) which in this case is very simple since the transformation is just the exponential.

In relation to your model/data combination, since you have only 17 years of data and just one series of cpue, and this is a rather common case, you may want to give the choice to set B0=K, i.e. equilibrium conditions at the start, in your function, to reduce the dimensionality of your profile likelihood approximation thus helping the optimizer.

HTH

____________________________________________________________________________________ 

Dr. Rubén Roa-Ureta
AZTI - Tecnalia / Marine Research Unit
Txatxarramendi Ugartea z/g
48395 Sukarrieta (Bizkaia)
SPAIN



> -----Mensaje original-----
> De: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] En nombre de Derek Ogle
> Enviado el: sábado, 26 de junio de 2010 22:28
> Para: R (r-help at R-project.org)
> Asunto: [R] optim() not finding optimal values
> 
> I am trying to use optim() to minimize a sum-of-squared 
> deviations function based upon four parameters.  The basic 
> function is defined as ...
> 
> SPsse <- function(par,B,CPE,SSE.only=TRUE)  {
>   n <- length(B)                             # get number of 
> years of data
>   B0 <- par["B0"]                            # isolate B0 parameter
>   K <- par["K"]                              # isolate K parameter
>   q <- par["q"]                              # isolate q parameter
>   r <- par["r"]                              # isolate r parameter
>   predB <- numeric(n)
>   predB[1] <- B0
>   for (i in 2:n) predB[i] <- 
> predB[i-1]+r*predB[i-1]*(1-predB[i-1]/K)-B[i-1]
>   predCPE <- q*predB
>   sse <- sum((CPE-predCPE)^2)
>   if (SSE.only) sse
>     else list(sse=sse,predB=predB,predCPE=predCPE)
> }
> 
> My call to optim() looks like this
> 
> # the data
> d <- data.frame(catch= 
> c(90000,113300,155860,181128,198584,198395,139040,109969,71896
> ,59314,62300,65343,76990,88606,118016,108250,108674), 
> cpe=c(109.1,112.4,110.5,99.1,84.5,95.7,74.1,70.2,63.1,66.4,60.
> 5,89.9,117.0,93.0,116.6,90.0,105.1))
> 
> pars <- c(800000,1000000,0.0001,0.17)                   # put 
> all parameters into one vector
> names(pars) <- c("B0","K","q","r")                      # 
> name the parameters
> ( SPoptim <- optim(pars,SPsse,B=d$catch,CPE=d$cpe) )    # run optim()
> 
> 
> This produces parameter estimates, however, that are not at 
> the minimum value of the SPsse function.  For example, these 
> parameter estimates produce a smaller SPsse,
> 
> parsbox <- c(732506,1160771,0.0001484,0.4049)
> names(parsbox) <- c("B0","K","q","r")        
> ( res2 <- SPsse(parsbox,d$catch,d$cpe,SSE.only=FALSE) )
> 
> Setting the starting values near the parameters shown in 
> parsbox even resulted in a movement away from (to a larger 
> SSE) those parameter values.
> 
> ( SPoptim2 <- optim(parsbox,SPsse,B=d$catch,CPE=d$cpe) )    # 
> run optim()
> 
> 
> This "issue" most likely has to do with my lack of 
> understanding of optimization routines but I'm thinking that 
> it may have to do with the optimization method used, 
> tolerance levels in the optim algorithm, or the shape of the 
> surface being minimized.
> 
> Ultimately I was hoping to provide an alternative method to 
> fisheries biologists who use Excel's solver routine.
> 
> If anyone can offer any help or insight into my problem here 
> I would be greatly appreciative.  Thank you in advance.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>