PCW4-2.R

library(boot)

#estimate the mean via bootstrapping
boot.mean <- function(data,index) return(mean(data[index]))

#calculate the CI via t-distribution
t.dist.ci <- function(samp) {
  df <- length(samp) - 1
  factors <- qt(c(0.025, 0.975), df = df)
  samp.mean <- mean(samp) 
  samp.se <- sd(samp)/sqrt(length(samp))
  samp.ci <- c(samp.mean + factors[1]*samp.se, samp.mean + factors[2]*samp.se) 
  return(samp.ci)
}

#calculate the CI via bootstrapping
boot.ci <- function(samp){
  temp.boot <- boot(samp, boot.mean, 1000)
  samp.ci <- quantile(temp.boot$t, c(0.025, 0.975))
  return(samp.ci)
}

# Iterate many times
iterate <- function(iterations) {
  
  # Create storage vectors 
  t.accurate <- rep(NA, iterations)
  boot.accurate <- rep(NA, iterations)
  both.accurate <- rep(NA, iterations)
  for (i in 1:iterations) {
    
    # Sample from the normal
    sample <- rnorm(15)
    
    # Calculate CIs
    t.ci <- t.dist.ci(sample)
    boot.ci <- boot.ci(sample)
    
    # Test whether they are accurate
    t.accurate[i] <- t.ci[1] <= 0 && 0 <= t.ci[2]
    boot.accurate[i] <- boot.ci[1] <= 0 && 0 <= boot.ci[2]
    both.accurate[i] <- t.accurate[i] && boot.accurate[i]
      
    
  }
  return(c(mean(t.accurate), mean(boot.accurate), mean(both.accurate)))
}

accuracy <- iterate(10000)
sprintf("The t-distribution confidence intervals were accurate %.1f%% of the time.", 100*accuracy[1])
sprintf("The bootstrapped confidence intervals were accurate %.1f%% of the time.", 100*accuracy[2])
sprintf("Both confidence intervals were accurate %.1f%% of the time.", 100*accuracy[3])