# Monte Carlo comparison of performances of different confidence intervals
# for a distribution mean --- normal distribution setting.

library(bootstrap)
set.seed(321)

# I'll initialize the results matrix.
results <- matrix( rep( 0, 12 ), nrow=4 )
rownames(results) <- c("Student's t","percentile","BCA","bootstrap t")
colnames(results) <- c("miss left","miss right","miss total") 
results

# I'll go through a loop 2500 times.  On each pass through the loop I'll
# generate a new data vector, compute 4 c.i's, and for each one determine 
# it it misses the estimand by falling to the left of it or falling to the 
# right of it.
for (i in 1:2500) {

  z <- rnorm(36)

  t.result <- t.test( z, conf.level=0.95 )
  if ( t.result$conf.int[1] > 0 ) { 
    results[1,2] <- results[1,2] + 1 
   } 
  else {
    if ( t.result$conf.int[2] < 0 ) { results[1,1] <- results[1,1] + 1 } 
   }

  bss <- matrix( sample( z, 36*3999, replace=T ), nrow=3999 )
  mean.rep <- apply( bss, 1, mean )
  sorted.mean.rep <- sort( mean.rep )
  if ( sorted.mean.rep[100] > 0 ) { 
    results[2,2] <- results[2,2] + 1 
   } 
  else {
    if ( sorted.mean.rep[3900] < 0 ) { results[2,1] <- results[2,1] + 1 } 
   }

  bca.res <- bcanon( z, 3999, mean, alpha=c( 0.025, 0.975 ) )
  if ( bca.res$confpoint[1,2] > 0 ) { 
      results[3,2] <- results[3,2] + 1 
     } 
    else {
      if ( bca.res$confpoint[2,2] < 0 ) { results[3,1] <- results[3,1] + 1 } 
     }


  sd.rep <- apply( bss, 1, sd )
  mean.z <- mean( z )
  z.rep <- ( mean.rep - mean.z )/( sd.rep/6 )
  sorted.z.rep <- sort( z.rep )
  se.hat <- sd( z )/6
  if ( mean.z - sorted.z.rep[3900]*se.hat > 0 ) { 
    results[4,2] <- results[4,2] + 1 
   } 
  else {
    if ( mean.z - sorted.z.rep[100]*se.hat < 0 ) { results[4,1] <- results[4,1] + 1 } 
   }
 
 }

results[,3] <- results[,1] + results[,2]
results/2500