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
______________________________________________
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.
