# STAT 824 sp 2023 Lec 08 in class

rm(list=ls())

n <- 20
alpha <- 1.1
beta <- 3
mu <- alpha*beta

S <- 500
B <- 500
z_val <- qnorm(0.975)
cov.asymp <- cov.boot.U <- cov.boot.S <- cov.boot.pctl <-  numeric(S)
for(s in 1:S){
  
  # generate data
  X <- rgamma(n,shape = alpha, scale = beta)
  
  # compute interval based on the assumption of asymptotic Normality
  Xbar <- mean(X) 
  sigma.hat <- sd(X)
  lo.asymp <- Xbar - z_val * sigma.hat / sqrt(n)
  up.asymp <- Xbar + z_val * sigma.hat / sqrt(n)
  cov.asymp[s] <- (lo.asymp < mu) &  (up.asymp > mu)

  # make Monte Carlo draws for the bootstrap
  Xbar.boot <- numeric(B)
  Tn.S <- numeric(B)
  for(b in 1:B){

    X.boot <- X[sample(1:n,n,replace=TRUE)]
    Xbar.boot[b] <- mean(X.boot)
    sigma.hat.boot <- sd(X.boot)
    Tn.S[b] <- sqrt(n) * (Xbar.boot[b] - Xbar)/sigma.hat.boot
    
  }
  
  # make studentized interval
  Tn.S <- sort(Tn.S)
  lo.boot.S <- Xbar - Tn.S[ceiling(0.975*B)] * sigma.hat / sqrt(n)
  up.boot.S <- Xbar - Tn.S[ceiling(0.025*B)] * sigma.hat / sqrt(n)
  cov.boot.S[s] <- (lo.boot.S < mu) &  (up.boot.S > mu)
  
  # make unstudentized interval
  Xbar.boot <- sort(Xbar.boot)
  lo.boot.U <- 2*Xbar - Xbar.boot[ceiling(0.975*B)]
  up.boot.U <- 2*Xbar - Xbar.boot[ceiling(0.025*B)]
  cov.boot.U[s] <- (lo.boot.U < mu) &  (up.boot.U > mu)
  
  # make percentile interval
  lo.boot.pctl <- Xbar.boot[ceiling(0.025*B)]
  up.boot.pctl <- Xbar.boot[ceiling(0.975*B)]
  cov.boot.pctl[s] <- (lo.boot.pctl < mu) &  (up.boot.pctl > mu)
  
}

# compute coverage probabilities of intervals
mean(cov.asymp)
mean(cov.boot.U)
mean(cov.boot.S) # the king
mean(cov.boot.pctl)