### Coverage Probability of the 100*(1-alpha)% Classical Confidence Intervals for theta.

### In repeated sampling X_1,X_2,... from the (unknown) population X ~ F(theta), 
### 100*(1-alpha)% of the computed CIs will include the parameter of interest theta.

### Coverage Probability of the 100*(1-alpha)% Bootstrap Confidence Intervals for theta.

### Note: Running the bootstrap on a sample from the population, resamples the sample
### B times to produce ONE bootstrap CI.

### In repeated sampling X_1,X_2,... from the (unknown) population X ~ F(theta), 
### 100*(1-alpha)% of the computed bootstrap CIs will include the parameter of interest 
### theta.

### Below is an S-Plus program to compute the coverage probability of the 
### one sample t-ci and one sample boostrap emp-ci and bca-ci.

####################################################################### 

# R is number of samples
R <- 1000
# n is number of observations per sample
n <- 31
# 31 was chosen so that the t-statistic will have 30 df

# We assume that the unknow population X ~ F(theta) is F = Normal
# and theta = (mu,sigma).

mu <- 2    #arbitrary value of mu
sigma <- 3 #arbitrary value of sigma

count.cci <- 0
count.bci.emp <- 0
count.bci.bca <- 0

zstar <- qnorm(0.975) 		#(for a 95% C.I. with normal dist.)
tstar <- qt(0.975,df = 30)	#(for a 95% C.I. with t dist., 30 df)

for(i in 1:R){
	x <- rnorm(n,mu,sigma)
	xbar <- mean(x)
	xsd <- sqrt(var(x))
	UCL <- xbar + tstar * xsd / sqrt(n)
	LCL <- xbar - tstar * xsd / sqrt(n)
	if( mu > LCL && mu < UCL) count.cci <- count.cci + 1
	x.boot <- bootstrap(x,mean,trace=F)
	x.boot.emp <- limits.emp(x.boot)
	if( mu > x.boot.emp[1] && mu < x.boot.emp[4]) count.bci.emp <- count.bci.emp + 1
	x.boot.bca <- limits.bca(x.boot)
	if( mu > x.boot.bca[1] && mu < x.boot.bca[4]) count.bci.bca <- count.bci.bca + 1
}

count.cci 			#this number is near 950
count.bci.emp
count.bci.bca

count.cci/R 		#this number is near 0.95
count.bci.emp/R
count.bci.bca/R