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:rvarad...@jhmi.edu] > Sent: Saturday, June 26, 2010 11:52 PM > To: Ravi Varadhan > Cc: Derek Ogle; R (r-help@R-project.org) > 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: 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,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. > > > > > > ______________________________________________ > > > 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.