[R] optim() not finding optimal values

[Ravi Varadhan]

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.