[R] Starting estimates for nls Exponential Fit

[dave fournier]

Thanks to Dennis Murphy who pointed out that ExponCycles
is undefined. It is an R gotcha. I had shortened the name but R still
remembered it so the script worked but only on my computer.
This should fix that.

ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9


#keep y1 and yn fixed and get initial value for E
Esp1 <- optim(c(E=E),method ="BFGS",
  function(x)
  {
    E=x[1]
    a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
    Y0=y1-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

E=Esp1[1]

Esp <- optim(c(y1=y1,yn=yn,E=E),method ="BFGS",
  function(x)
  {
    E=x[3]
    a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
    Y0=x[1]-a*E^Expon[1];
    diff=ExponValues-(Y0+a*E^Expon)
    return(1000*sum(diff*diff))
  })$par

y1=Esp[1]
y2=Esp[2]
E=Esp[3]

a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];




-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com