Code from my curve fitting blog post about Ecolog.

exponentialEcolog.R

###### The exponential growth of ecolog
###### Raw data can be found here
###### https://listserv.umd.edu/cgi-bin/wa?A2=ind1108e&L=ecolog-l&P=3401

#####Data from Aug 31. 2011 post
require(date)
subscribers <- c(100,6000,7000,8000,9000,10000,11000,12000,13000)
start.date <- mdy.date(1,1,92)
month <- c("01/01/1992","09/01/2006","11/01/2007","10/01/2008","03/01/2009","04/01/2010","09/01/2010","02/01/2011","09/01/2011")
days <- as.date(month)-start.date
z <- nls(subscribers~100*exp(r*days),start=list(r=.001))
######Plot the data ####
plot(days,subscribers,xlab="Days since Sept 1992",ylab="Number of Subscribers",xlim=c(0,7500),ylim=c(0,18000))
days.v <- 1:8000
fit.z <- 100*exp(coef(z)*days.v)
lines(days.v,fit.z)
#### now add the new data point ####
new.day <- mdy.date(6,23,2012)-start.date
points(days.v[new.day],fit.z[new.day],col=2,pch=19)
points(new.day,15000,pch=19,col=1)


#### so how can we find an estimation of what the starting value should be?
dat <- cbind(days[-1],subscribers[-1])
fp.est <- function(param,dat){
  z <- nls(dat[,2]~param*exp(r*dat[,1]),start=list(r=.001))
  return(deviance(z))
}

tmp <- optim(1,fp.est,dat=dat,method="L-BFGS-B",lower=1,upper=6000)



##### Ok, can we improve with a jacknife style estimate? ######
##### Generate all the possibe combinations
comb.mat <- combn(1:9,8)
dat <- rbind(dat,c(new.day,15000))
est <- vector()
for(i in 1:dim(comb.mat)[2]){
  my.dat <- dat[comb.mat[,i],]
  tmp <- optim(1,fp.est,dat=my.dat,method="L-BFGS-B",lower=1,upper=6000)
  est[i] <- tmp$par
}

tmp <- optim(1,fp.est,dat=dat,method="L-BFGS-B",lower=1,upper=6000)
####Now calculate the bias and variance from the jackknife
bias  <- -(1/8)*sum(8*(tmp$par-est))
var_jack <- (1/(8*9))*sum((8*(tmp$par-est))^2-(9*bias)^2)
sd_jack <- sqrt(var_jack)


#### Make the final plot

z <- nls(dat[,2]~round(tmp$par)*exp(r*dat[,1]),start=list(r=.001))
######Plot the data ####
plot(dat,xlab="Days since Sept 1992",ylab="Number of Subscribers",xlim=c(0,7500),ylim=c(0,18000))
days.v <- 1:8000
fit.z <- round(tmp$par)*exp(coef(z)*days.v)
lines(days.v,fit.z)
#### now add the new data point ####
new.day <- mdy.date(6,23,2012)-start.date
points(days.v[new.day],fit.z[new.day],col=1,pch=19)
points(new.day,15000,pch=19,col=2)