[R] optim() not finding optimal values

[Ravi Varadhan]

A slightly better scaling is the following:

par.scale <- c(1.e06, 1.e06, 1.e-05, 1)  # "q" is scaled differently 

> SPoptim <- optim(pars, SPsse, B=d$catch, CPE=d$cpe, control=list(maxit=1500, parscale=par.scale))
> SPoptim
$par
          B0            K            q            r 
7.320899e+05 1.159939e+06 1.485560e-04 4.051735e-01 

$value
[1] 1619.482

$counts
function gradient 
     585       NA 

$convergence
[1] 0

$message
NULL


Note that the Nelder-Mead converges in half the number of iterations compared to that under previous scaling.

Ravi. 
____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Ravi Varadhan <rvaradhan at jhmi.edu>
Date: Sunday, June 27, 2010 0:42 am
Subject: Re: [R] optim() not finding optimal values
To: Derek Ogle <DOgle at northland.edu>
Cc: "R (r-help at R-project.org)" <r-help at r-project.org>


> Derek,
>  
>  The problem is that your function is poorly scaled.   You can see 
> that the parameters vary over 10 orders of magnitude (from 1e-04 to 
> 1e06).   You can get good convergence once you properly scale your 
> function.  Here is how you do it:
>  
>  par.scale <- c(1.e06, 1.e06, 1.e-06, 1.0)
>   
>  SPoptim <- optim(pars, SPsse, B=d$catch, CPE=d$cpe, 
> control=list(maxit=1500, parscale=par.scale))
>  
>  > SPoptim
>  $par
>            B0            K            q            r 
>  7.329553e+05 1.160097e+06 1.484375e-04 4.050476e-01 
>  
>  $value
>  [1] 1619.487
>  
>  $counts
>  function gradient 
>      1401       NA 
>  
>  $convergence
>  [1] 0
>  
>  $message
>  NULL
>  
>  
>  Hope this helps,
>  Ravi.
>  
>  ____________________________________________________________________
>  
>  Ravi Varadhan, Ph.D.
>  Assistant Professor,
>  Division of Geriatric Medicine and Gerontology
>  School of Medicine
>  Johns Hopkins University
>  
>  Ph. (410) 502-2619
>  email: rvaradhan at jhmi.edu
>  
>  
>  ----- Original Message -----
>  From: Derek Ogle <DOgle at northland.edu>
>  Date: Saturday, June 26, 2010 4:28 pm
>  Subject: [R] optim() not finding optimal values
>  To: "R (r-help at R-project.org)" <r-help at r-project.org>
>  
>  
>  > 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
>  >  
>  >  PLEASE do read the posting guide 
>  >  and provide commented, minimal, self-contained, reproducible code.
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.