optim() not finding optimal values

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

optim() not finding optimal values

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

______________________________________________
[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.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Nikhil Kaza-2
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[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.
>
> ______________________________________________
> [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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Ravi Varadhan
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
[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: [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[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.
>  
>  ______________________________________________
>  [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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: 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: [hidden email]


----- Original Message -----
From: Ravi Varadhan <[hidden email]>
Date: Sunday, June 27, 2010 0:42 am
Subject: Re: [R] optim() not finding optimal values
To: Derek Ogle <[hidden email]>
Cc: "R ([hidden email])" <[hidden email]>


> 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: [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[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.
>  >  
>  >  ______________________________________________
>  >  [hidden email] mailing list
>  >  
>  >  PLEASE do read the posting guide
>  >  and provide commented, minimal, self-contained, reproducible code.
>  
>  ______________________________________________
>  [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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Derek Ogle
Ravi,

Thank you very much for the pointer to parscale.  This is extremely useful -- in this and some other problems that I am working on.  Thanks again for the valuable help.

> -----Original Message-----
> From: Ravi Varadhan [mailto:[hidden email]]
> Sent: Saturday, June 26, 2010 11:52 PM
> To: Ravi Varadhan
> Cc: Derek Ogle; R ([hidden email])
> Subject: Re: [R] optim() not finding optimal values
>
> 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: [hidden email]
>
>
> ----- Original Message -----
> From: Ravi Varadhan <[hidden email]>
> Date: Sunday, June 27, 2010 0:42 am
> Subject: Re: [R] optim() not finding optimal values
> To: Derek Ogle <[hidden email]>
> Cc: "R ([hidden email])" <[hidden email]>
>
>
> > 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: [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[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,62
> 300,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,11
> 7.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
> >
> >  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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Rubén Roa
In reply to this post by Derek Ogle
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: [hidden email]
> [mailto:[hidden email]] En nombre de Derek Ogle
> Enviado el: sábado, 26 de junio de 2010 22:28
> Para: R ([hidden email])
> 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.
>
> ______________________________________________
> [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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

dave fournier
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Ravi Varadhan
In reply to this post by Rubén Roa
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: [hidden email] [mailto:[hidden email]] On
Behalf Of Rubén Roa
Sent: Monday, June 28, 2010 2:24 AM
To: Derek Ogle; R ([hidden email])
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: [hidden email]
> [mailto:[hidden email]] En nombre de Derek Ogle
> Enviado el: sábado, 26 de junio de 2010 22:28
> Para: R ([hidden email])
> 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.
>
> ______________________________________________
> [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-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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optim() not finding optimal values

Ravi Varadhan
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: [hidden email] [mailto:[hidden email]] 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: [hidden email] [mailto:[hidden email]] On
Behalf Of Rubén Roa
Sent: Monday, June 28, 2010 2:24 AM
To: Derek Ogle; R ([hidden email])
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: [hidden email]
> [mailto:[hidden email]] En nombre de Derek Ogle
> Enviado el: sábado, 26 de junio de 2010 22:28
> Para: R ([hidden email])
> 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.
>
> ______________________________________________
> [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-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-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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.