[R] optim() not finding optimal values

Mon Jun 28 15:55:53 CEST 2010

Oops,

I my previous email, the second line in the `SPsse.log' function should have
been:

  par <- exp(par)  # log-transformation

Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Ravi Varadhan
Sent: Monday, June 28, 2010 9:48 AM
To: 'Rubén Roa'; 'Derek Ogle'; 'R'
Subject: Re: [R] optim() not finding optimal values

Ruben,

Transforming the parameters is also a good idea, but the obvious caveat is
that the transformation must be feasible. The log-transformation is only
feasible for positive parameter domain.  This happens to be the case for
OP's problem.  In fact, the log-transform does better than ratio scaling in
terms of improving rate of convergence for this particular model/data
combination.   

SPsse.log <- function(par,B,CPE,SSE.only=TRUE)  {
  n <- length(B)                             # get number of years of data
  par <- tan(par)  # log-transformation
  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)
}

> SPoptim3 <- optim(log(parsbox),SPsse.log,B=d$catch,CPE=d$cpe,
method="BFGS") 
> SPoptim3
$par
        B0          K          q          r 
13.5035475 13.9634098 -8.8142230 -0.9030033 

$value
[1] 1619.481

$counts
function gradient 
      54        6 

$convergence
[1] 0

$message
NULL

> exp(SPoptim3$par)  # transforming to original scale
          B0            K            q            r 
7.320086e+05 1.159396e+06 1.486044e-04 4.053505e-01 
>

I don't think there is a need to set B0=K since convergence is pretty good
with scaling and/or with log-transformation.

Best,
Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Rubén Roa
Sent: Monday, June 28, 2010 2:24 AM
To: Derek Ogle; R (r-help at R-project.org)
Subject: Re: [R] optim() not finding optimal values

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.
> 

______________________________________________
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.