# R code to obtain confidence intervals for Control group data on p. 11 of E&T.
# Can compare with results given in Ch. 12 of E&T.

# Need to load the bootstrap library in order to make use of the boott function.
# (Note: Some functions pertaining to bootstrapping are in the library bootstrap,
# while other bootstrap-related functions are in the library boot.)
library(bootstrap)

set.seed(321)

control <- c(52,104,146,10,51,30,40,27,46)

# Here is a 90% Student's t confidence interval:
t.test(control, conf.level=0.90)
# This differs from the one given on p. 159 of E&T because E&T use the plug-in
# estimate of the distribution standard deviation, and R uses the common sample
# standard deviation, s.

# Now I'll use the bootstrap-t method.  Because the interval is for the dist'n
# mean, s/sqrt(n) can be used for the estimated standard error and nested boot-
# strapping can be avoided.  
bss <- matrix( sample( control, 9*999, replace=T ), nrow=999 )
mean.rep <- apply( bss, 1, mean )
sd.rep <- apply( bss, 1, sd )
z.rep <- ( mean.rep - mean( control ) )/( sd.rep/sqrt(9) )
sorted <- sort( z.rep )
# upper confidence bound:
mean( control ) - sorted[50]*sd( control )/sqrt( length( control ) )
# lower confidence bound:
mean( control ) - sorted[950]*sd( control )/sqrt( length( control ) )
# These confidence bounds (based on p. 160 of E&T) are close to the ones given
# on p. 161 of E&T.  (They are not expected to exactly match the ones from E&T
# due to different random numbers being used.)

# Now I'll use the boott function from the bootstrap library.  I'm going to reset
# the random number seed to what I initially set it to in order to see if the 
# result is the same as I got above.
set.seed(321)

# Since the estimator the confidence interval is based on is the sample mean, for
# which a good estimate of it's standard error is known, I'll avoid letting boott
# use nested bootstrapping and supply a function to compute the estimated standard
# error needed for each bootstrap replicate of z/t.
se.mean <- function(x, ...) { sqrt( var(x)/length(x) ) }
# Note: One needs the ... in the argument list for the function used for sdfun.

ci.results <- boott( control, theta=mean, sdfun=se.mean, nboott=999, perc=c(0.05,0.95) )
ci.results$confpoints
# Exactly the same results!