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

library(bootstrap)
set.seed(321)

# I'll initialize the results matrix.
results <- matrix( rep( 0, 12 ), nrow=4 )
rownames(results) <- c("exact","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)

  numerator <- 35*var( z )
  if ( numerator/qchisq( 0.975, 35 ) > 1 ) { 
    results[1,2] <- results[1,2] + 1 
   } 
  else {
    if ( numerator/qchisq( 0.025, 35 ) < 1 ) { results[1,1] <- results[1,1] + 1 } 
   }

  bss <- matrix( sample( z, 36*3999, replace=T ), nrow=3999 )

  plug.in.var <- function( x ) { 
    piv <- var( x )*( 1 - 1/length( x ) )
    piv
   }

  var.rep <- apply( bss, 1, plug.in.var )
  sorted.var.rep <- sort( var.rep )
  if ( sorted.var.rep[100] > 1 ) { 
    results[2,2] <- results[2,2] + 1 
   } 
  else {
    if ( sorted.var.rep[3900] < 1 ) { results[2,1] <- results[2,1] + 1 } 
   }

  bca.res <- bcanon( z, 3999, plug.in.var, 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 } 
     }


  se.hat <- function( x ) {
    deviations <- x - mean( x )
    dev.2 <- deviations^2
    dev.4 <- dev.2^2
    m.2 <- sum( dev.2 )/length( x )
    m.4 <- sum( dev.4 )/length( x )
    sehat <- sqrt( ( m.4 - m.2^2 )/( length(x) - 1 ) )
    sehat
   } 
   
  se.hat.rep <- apply( bss, 1, se.hat )
  var.z <- plug.in.var( z )
  z.rep <- ( var.rep*36/35 - var.z )/( se.hat.rep )
  sorted.z.rep <- sort( z.rep )
  est.se <- se.hat( z )
  if ( var.z - sorted.z.rep[3900]*est.se > 1 ) { 
    results[4,2] <- results[4,2] + 1 
   } 
  else {
    if ( var.z - sorted.z.rep[100]*est.se < 1 ) { results[4,1] <- results[4,1] + 1 } 
   }
 
 }

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