[R] Starting estimates for nls Exponential Fit

[Katharine Mullen]

You used starting values:
   pa <- c(1,2,3)

but both algorithms (port and Gauss-Newton) fail if you use the slightly
different values:
   pa <- c(1,2,3.5)

Scaling does not fix the underlying sensitivity to starting values.
pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both
alg. also fail if there is insufficient spread (e.g., both fail for
pa <- c(1,1.25,1.5) ).

For the record, DE uses (by default at least) a random start and for bad
starts will sometimes fail to give good results; decrease the probability
this happens by increasing itermax from the default itermax=200, as in:

ss <- DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
              control=list(trace=FALSE, itermax=1000))

On Wed, 2 Dec 2009, Prof. John C Nash wrote:

> Kate Mullen showed one approach to this problem by using DEOptim to get
> some good starting values.
>
> However, I believe that the real issue is scaling (Warning: well-ridden
> hobby-horse!).
>
> With appropriate scaling, as below, nls does fine. This isn't saying nls
> is perfect -- I've been trying to figure out how to do a nice job of
> helping folk to scale their problems. Ultimately, it would be nice to
> has an nls version that will do the scaling and also watch for some
> other situations that give trouble.
>
> Cheers, JN
>
>
> ## JN test
> rm(list=ls())
>
> ExponValues <- c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
>                  2468.32,2778.47)
>
> ExponCycles <- c(17,18,19,20,21,22,23,24,25)
>
> mod <- function(x) x[1] + x[2]*x[3]^ExponCycles
>
> modj <- function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles)
>
> fun <- function(x) sum((ExponValues-mod(x))^2)
>
>
>
> pa<-c(1,2,3)
> lo<-c(0,0,0)
> up<-c(20,20,20)
> names(pa) <- c("Y0", "a", "E")
>
> ## fit w/port and GN
> Emodjn<- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)),
>              start=pa, algorithm='port', lower=lo, upper=up,
>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>
> Emodjn1 <- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa,
>              control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE))
>
> ## fit
> matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type="l")
>
> ______________________________________________
> 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.
>