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@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.
