#
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# problem on epage 55
# acf plot on epage 40

# fix the random.seed, so I'll always get the same answer:
set.seed(10131985)

Nit = c(0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,6,6,6)
AOB = c(4.26,4.15,4.68,6.08,5.87,6.92,6.87,6.25,6.84,6.34,6.56,6.52,7.39,7.38,7.74,7.76,8.14,7.22)

AOBm = tapply(AOB,Nit,mean) # the means 
Nitm = tapply(Nit,Nit,mean)

postscript("../../WriteUp/Graphics/Chapter1/ex_1_22_b.eps", onefile=FALSE, horizontal=FALSE)

plot(Nit,AOB,xlim=c(0,6),ylim=c(min(AOB),max(AOB)),pch=19) # plot the raw data 
library(splines)

fitAOB = lm(AOBm~ns(Nitm,df=2)) # a natural spline fit
xmin=min(Nit); xmax=max(Nit)

# save the residuals
resids = fitAOB$residuals

nBoot = 100
B = array(0,dim=c(nBoot,length(fitAOB$coefficients)))
for( ii in 1:nBoot ){
  ystar  = AOBm + sample( resids, replace=T )
  bModel = lm(ystar~ns(Nitm,df=2))
  B[ii,]  = as.vector(bModel$coefficients)
  # lets plot each of these samples
  lines(seq(xmin,xmax,0.5), predict(bModel,data.frame(Nitm=seq(xmin,xmax,0.5))),col="gray") # plot the spline fit
}

lines(seq(xmin,xmax,0.5), predict(fitAOB,data.frame(Nitm=seq(xmin,xmax,0.5))),col="red") # plot the original spline fit on the top 

dev.off()