# R code for Problem 11,  STAT 789,  Summer 2006

library(bootstrap)

set.seed(321)
z <- rnorm(40)

# (a)
t.test( z, conf.level=0.95 )

# (b)
se.mean <- function(x, ...) { sqrt( var(x)/length(x) ) }
ci.results <- boott( z, theta=mean, sdfun=se.mean, nboott=1999, perc=c(0.025,0.975) )
ci.results$confpoints

set.seed(321)
x <- rexp(40)

# (c)
t.test( x, conf.level=0.95 )

# (d)
ci.results <- boott( x, theta=mean, sdfun=se.mean, nboott=1999, perc=c(0.025,0.975) )
ci.results$confpoints

# (e)
theta <- function(x, ...) { var(x)*(length(x) - 1)/length(x) }
ci.results <- boott( x, theta, nbootsd=250, nboott=1999, perc=c(0.025,0.975) )
ci.results$confpoints

# (f)
ci.results <- boott( x, theta, VS=TRUE, perc=c(0.025,0.975) )
ci.results$confpoints

### For just fun I'll override the default argument values in order to hopefully improve accuracy.
ci.results <- boott( x, theta, VS=TRUE, v.nbootg=200, v.nbootsd=250, v.nboott=999, perc=c(0.025,0.975) )
ci.results$confpoints
# Note: Based on a comparison with the exact interval given below, the extra computation did not result
# in an improvement on what was obtained using the default settings.

### For comparison purposes, here are confidence bounds for an exact confidence interval
### obtained using the parametric model.
# upper confidence bound:
( 2*sum(x)/qchisq(0.025, 80) )^2
# lower confidence bound:
( 2*sum(x)/qchisq(0.975, 80) )^2

# Here is a confidence interval for the variance based on an assumption of normality.  
# It shouldn't be trusted, but it's the only one that some people know.
# upper confidence bound:
(length(x) - 1)*var(x)/qchisq(0.025, 39) 
# lower confidence bound:
(length(x) - 1)*var(x)/qchisq(0.975, 39)