[R] optim() not finding optimal values
Ravi Varadhan
rvaradhan at jhmi.edu
Mon Jun 28 15:47:33 CEST 2010
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.
More information about the R-help
mailing list