[R] optim() not finding optimal values
Derek Ogle
DOgle at northland.edu
Sun Jun 27 17:08:01 CEST 2010
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:rvaradhan at jhmi.edu]
> Sent: Saturday, June 26, 2010 11:52 PM
> To: Ravi Varadhan
> Cc: Derek Ogle; R (r-help at 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: rvaradhan at jhmi.edu
>
>
> ----- Original Message -----
> From: Ravi Varadhan <rvaradhan at jhmi.edu>
> Date: Sunday, June 27, 2010 0:42 am
> Subject: Re: [R] optim() not finding optimal values
> To: Derek Ogle <DOgle at northland.edu>
> Cc: "R (r-help at R-project.org)" <r-help at 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: rvaradhan at jhmi.edu
> >
> >
> > ----- Original Message -----
> > From: Derek Ogle <DOgle at northland.edu>
> > Date: Saturday, June 26, 2010 4:28 pm
> > Subject: [R] optim() not finding optimal values
> > To: "R (r-help at R-project.org)" <r-help at 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 at r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible
> code.
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> >
> > PLEASE do read the posting guide
> > and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list