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: rvarad...@jhmi.edu ----- Original Message ----- From: Ravi Varadhan <rvarad...@jhmi.edu> Date: Sunday, June 27, 2010 0:42 am Subject: Re: [R] optim() not finding optimal values To: Derek Ogle <do...@northland.edu> Cc: "R (r-help@R-project.org)" <r-help@r-project.org> > 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: rvarad...@jhmi.edu > > > ----- Original Message ----- > From: Derek Ogle <do...@northland.edu> > Date: Saturday, June 26, 2010 4:28 pm > Subject: [R] optim() not finding optimal values > To: "R (r-help@R-project.org)" <r-help@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@r-project.org mailing list > > > > PLEASE do read the posting guide > > and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ 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.