#
# Epage > 137
#
# 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.
# 
#-----

if(is.null(version$language) == FALSE){ 
  require(alr3)
}else{
  library(alr3)
}

if(is.null(version$language) == FALSE) data(lakemary)

attach(lakemary)

LInfinity0     <- 1.05*max(Length) # this is an initial guess for L_{\infty}
minusKTMinusT0 <- log( 1 - Length / LInfinity0 ) # this is the response of the linear fit 

mf <- lm( minusKTMinusT0 ~ Age ) # fit a linear model and extract \beta_{0} and \beta_1
b0 <- coefficients(mf)[1]
b1 <- coefficients(mf)[2]

# here are our initial guesses at the parameters K and t_0: 
K0  <- -b1
t00 <- -b0/b1 

# now we can fit the full nonlinear model: 
n1 <- nls( Length ~ LInfinity * ( 1 - exp(-K*(Age - t0)) ), 
           data=lakemary, start=list(LInfinity=LInfinity0, K=K0, t0=t00) )

# how close was our initial guess: 
print( c(LInfinity0,K0,t00) ) 

# plot the data and the nonlinear fit 
postscript("../../WriteUp/Graphics/Chapter11/prob_2_scatterplot.eps", onefile=FALSE, horizontal=FALSE)
plot( Age, Length,xlab="Age", ylab="Length" )
x <- seq(0,6,0.1)
lines(x,predict(n1,data.frame(Age=x)))
dev.off() 

# bootstrap uses boot.case written originally by Lexin Li, and simplified
# by S. Weisberg for use with the book.  It works for any regression object
# fix the random.seed, so I'll always get the same answer:
if(is.null(version$language) == FALSE){
  set.seed(10131985)
  n1.boot <- boot.case(n1,B=999) 
  library(car)
  
  postscript("../../WriteUp/Graphics/Chapter11/prob_2_param_scatterplotmatrix.eps", onefile=FALSE, horizontal=FALSE)
  scatterplot.matrix(n1.boot,diagonal="histogram",
       col=palette(),#[-1],
       lwd=0.7,pch=".",
       labels=c(expression(L[infty]),expression(K),expression(t[0])),
       ellipse=FALSE,smooth=TRUE,level=c(.90)) 
  dev.off() 
  
  # summary statistics on the bootstrapped samples:
  # 
  n1.boot.summary <- data.frame(rbind(
  apply(n1.boot,2,mean),
  apply(n1.boot,2,sd),
  apply(n1.boot,2,function(x){quantile(x,c(.025,.975))})))
  row.names(n1.boot.summary)<-c("Mean","SD","2.5%","97.5%")
  n1.boot.summary

  # summary statistics on the large sample statistics:
  # 
  n1.sum <- summary(n1)$param
  n1.ls.summary <- data.frame(rbind(
  n1.sum[,1],n1.sum[,2],n1.sum[,1]-1.96*n1.sum[,2],
  n1.sum[,1]+1.96*n1.sum[,2]))
  row.names(n1.ls.summary) <- row.names(n1.boot.summary)
  n1.ls.summary

  # compare the two methods:
  # 
  round(cbind(n1.ls.summary,n1.boot.summary),2)
  
}