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

library(bootstrap)
set.seed(321)

# I'll initialize the results matrix.
results <- matrix( rep( 0, 18 ), nrow=6 )
rownames(results) <- c("Student's t","percentile","BCA","ABC","bootstrap t","exact")
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 5 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:50) {

  x <- rexp(16)

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

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

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

  resamp.mean <- function( p, x ) { sum( p*x )/sum( p ) }
  abc.res <- abcnon( x, resamp.mean, alpha=c( 0.025, 0.975 ) )
  if ( abc.res$limits[1,2] > 1 ) { 
      results[4,2] <- results[4,2] + 1 
     } 
    else {
      if ( abc.res$limits[2,2] < 1 ) { results[4,1] <- results[4,1] + 1 } 
     }

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

  numerator <- 2*sum(x)
  if ( numerator/qchisq( 0.975, 32 ) > 1 ) { 
     results[6,2] <- results[6,2] + 1 
    } 
   else {
     if ( numerator/qchisq( 0.025, 32 ) < 1 ) { results[6,1] <- results[6,1] + 1 } 
    }
 
 }

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