# optim() not finding optimal values Classic List Threaded 9 messages Open this post in threaded view
|

## 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 <- 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. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: optim() not finding optimal values

 Your function is very irregular, so the optim is likely to return   local minima rather than global minima. Try different methods  (SANN, CG, BFGS) and see if you get the result   you need. As with all numerical optimsation, I would check the   sensitivity of the results to starting values. Nikhil Kaza Asst. Professor, City and Regional Planning University of North Carolina [hidden email] On Jun 26, 2010, at 4:27 PM, Derek Ogle wrote: > 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 <- 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. > > ______________________________________________ > [hidden email] 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. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: optim() not finding optimal values

 In reply to this post by Derek Ogle 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  1619.487 \$counts function gradient     1401       NA \$convergence  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: [hidden email] ----- Original Message ----- From: Derek Ogle <[hidden email]> Date: Saturday, June 26, 2010 4:28 pm Subject: [R] optim() not finding optimal values To: "R ([hidden email])" <[hidden email]> > 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 <- 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. >   >  ______________________________________________ >  [hidden email] mailing list >   >  PLEASE do read the posting guide >  and provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: optim() not finding optimal values

Open this post in threaded view
|

## Re: optim() not finding optimal values

Open this post in threaded view
|

## Re: optim() not finding optimal values

Open this post in threaded view
|

## Re: optim() not finding optimal values

 In reply to this post by Derek Ogle  If you are going to make this program available for general   use you want to take every precaution to make it bulletproof.   This is a fairly informative data set. The model will undoubtedly   be used on far less informative data.  While the model looks   pretty simple it is  very challenging from a numerical point of view.   I took a moment to code it up in AD Model Builder. The true minimum is   1619.480495 So I think Ravi has finally arrived pretty close to the answer.   One way of judging the difficulty of a model is to look at the   eigenvalues of the Hessian at the minimum. They are        3.122884668e-09 1.410866202e-08  1866282.520 1.330233652e+13   so the condition number is around 1.e+21. One begins to see why these   models are challenging.  The model as formulated represents the state   of the art in fisheries models circa 1985.   A lot of progress has been made since that time.   Using B_t for the biomass and C_t for the catch the equation   in the code is             B_{t+1} = B_t + r *B_t*(1-B_t/K) -C_t  (1)   First  notice that   for (1) to make sense the following conditions must be satisfied        B_t > 0 for all t        r > 0        K>0   Strictly speaking it is not necessary that B_t<=K but if B_t>K and r   is large then B_{t+1} could be <0.  So formulation (1) gives   Murphys law a good chance.  How to fix it. Notice that (1) is really   a rough approximation to the solution of a differential equation       B'(t) =  r *B(t)*(1-B(t)/K) -C  (2)   where in (2) C is a constant catch rate.  To fix (1) we use   a semi-implicit differencing scheme. Because it is useful to   allow smaller step sizes than one we denote them by d.        B_{t+d} = B_t + d* r *B_t*(1-B_{t+d}/K) -d*C_t*B_{t+d}/B_t  (1)   The idea is that the quantity  1-x with x>0 will be replaced by   1/(1+x).  Expanding 2 and solving for B_{t+d} yields       B_{t+d} = (1+d*r) B_t / (1+d*r*B_t/K +d*C_t/B_t)  (3)    So long as r>0, K>0 C_t>0 then starting from an initial value    B_0 > 0 ensures that B_t> 0 for all t>0.  We can let    d=1/nsteps where nsteps is the number of steps in the    approximate integration for each year    which can be increased until the solution is judged to be close    enough to the exact solution from (2)    Notice that in (3) as C_t --> infinity  B_{t+d} --> 0    So that you can never catch more fish than you have.    I coded up this version of the model in AD Model Builder and    fit it to the data. It is now much more resistant to bad    starting values for the parameters etc.    If anyone wants the tpl file for the model in ADMB they can    contact me off list. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.