# Monte Carlo comparison of performances of different confidence intervals
# for the estimation setting pertaining to Table 13.3 on p. 175 of E&T.

library(bootstrap)
set.seed(321)

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

  z <- rnorm(16)

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

  theta <- function(x, xdata, ...) { exp( mean( x ) ) }
  ci.result <- boott( z, theta, VS=TRUE, v.nbootg=100, v.nbootsd=100, v.nboott=999, perc=c(0.025,0.975) )
  if ( ci.result$confpoints[1] > 1 ) { 
    results[6,2] <- results[6,2] + 1 
   } 
  else {
    if ( ci.result$confpoints[2] < 1 ) { results[6,1] <- results[6,1] + 1 } 
   }

  bss <- matrix( sample( z, 16*999, replace=T ), nrow=999 )
  mean.rep <- apply( bss, 1, mean )
  theta.rep <- exp( mean.rep )
  sorted.theta.rep <- sort( theta.rep )
  if ( sorted.theta.rep[25] > 1 ) { 
    results[3,2] <- results[3,2] + 1 
   } 
  else {
    if ( sorted.theta.rep[975] < 1 ) { results[3,1] <- results[3,1] + 1 } 
   }

  theta.hat <- exp( mean( z ) )
  se.hat <- sd( theta.rep )
  if ( theta.hat - 1.96*se.hat > 1 ) { 
    results[2,2] <- results[2,2] + 1 
   } 
  else {
    if ( theta.hat + 1.96*se.hat < 1 ) { results[2,1] <- results[2,1] + 1 } 
   }
  
  # Note: I'll make use of the transformation respecting property for the BCA and ABC intervals.
  bca.res <- bcanon( z, 999, mean, alpha=c( 0.025, 0.975 ) )
    if ( bca.res$confpoint[1,2] > 0 ) { 
        results[4,2] <- results[4,2] + 1 
       } 
      else {
        if ( bca.res$confpoint[2,2] < 0 ) { results[4,1] <- results[4,1] + 1 } 
       }

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

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