# Another example of the percentile interval, based on the situation 
# described on p. 171.

library(bootstrap)

set.seed( 321 )

x <- rnorm( 10 )

pt.est <- exp( mean( x ) )
pt.est
# point estimate

# I'll get a 95% percentile interval.
bss <- matrix( sample( x, 10*3999, replace=T ), nrow=3999 )
th.rep <- exp( apply( bss, 1, mean ) )
sorted.th.rep <- sort( th.rep )

# upper confidence bound:
sorted.th.rep[3900]
# lower confidence bound:
sorted.th.rep[100]
# A bit different from the percentile interval given by(13.6) on p. 171 of E&T,
# but it needs to be recalled that they started with a different randomly generated 
# sample of size 10.  (This is clear since their point estimate differs from the 
# one computed above.)

# Let's compare with an exact parametric interval, based on an assumption of
# normality, but not that it is a standard normal distribution.
t.result <- t.test( x, conf.level=0.95 )
exact <- exp( t.result$conf.int )
exact

# Now I'll compute a standard interval, of the form given by (12.6) on
# p. 154 of E&T, which is based on an assumption of approximate normality
# for theta hat.
se.hat <- sd(th.rep)
# upper confidence bound:
pt.est + qnorm(0.975)*se.hat
# lower confidence bound:
pt.est - qnorm(0.975)*se.hat

# Let's further compare with a bootstrap-t interval.
theta <- function(x, xdata, ...) { exp( mean( x ) ) }
ci.results <- boott( x, theta, nbootsd=400, nboott=3999, perc=c(0.025,0.975) )
ci.results$confpoints

# Now I'll try the bootstrap-t using automatic variance stabilization.
ci.results <- boott( x, theta, VS=TRUE, v.nbootg=400, v.nbootsd=400, v.nboott=3999, perc=c(0.025,0.975) )
ci.results$confpoints