set.seed(12343)

alpha <- 0.05 
normval <- qnorm(1 - alpha/2) 
n <- 100
reps <- 10000  
normmat <- matrix(0, nrow = reps, ncol = 2) # Empty matrix to store the sample mean s.e
mu<-10

#y <- 1:reps
#ymat <- rbind(y, y) 

for (i in 1:reps) { 
  samp<-rnorm(n,10,2) 
  sampmean <- mean(samp) 
  sampse <- sqrt(var(samp) / n) 
  normmat[i,1 ] <- sampmean
  normmat[i,2 ] <- sampse 
    } 

# CP Function
coverage <- function(sampmean, sampse, true, level = 1-alpha, df = Inf){ # Estimate, 
  # standard error,
  # true parameter, 
  # confidence level, 
  # and df  
  qtile <- level + (1 - level)/2 # Compute the proper quantile
  lower.bound <- sampmean - qt(qtile, df = df)*sampse # Lower bound
  upper.bound <- sampmean + qt(qtile, df = df)*sampse # Upper bound 
  # Is the true parameter in the confidence interval? (yes = 1)
  true.in.ci <- ifelse(true >= lower.bound & true <= upper.bound, 1, 0)
  cp <- mean(true.in.ci) # The coverage probability
  mc.lower.bound <- cp - 1.96*sqrt((cp*(1 - cp))/length(sampmean)) # Monte Carlo error  
  mc.upper.bound <- cp + 1.96*sqrt((cp*(1 - cp))/length(sampmean))  
  return(list(coverage.probability = cp, # Return results
              true.in.ci = true.in.ci,
              ci = cbind(lower.bound, upper.bound),
              mc.eb = c(mc.lower.bound, mc.upper.bound)))
}

cp.mean <- coverage(normmat[ , 1], normmat[ , 2], mu, df = n-1)
cp.mean[1]

par(mar = c(5, 6, .5, .5))
plot(seq(1, 100, length = 100), seq(9, 11, length = 100), type = "n",
     axes = FALSE, xlab = "", ylab = "")
title(xlab = expression("Sample No."), cex.lab = 1.0)
title(ylab = expression("95% Confidence Interval"), line = 3.75, cex.lab = 1.0)
box()
axis(1, at = seq(0, 100, 10), cex.axis = 1.25)
axis(2, at = seq(9, 11, .5), cex.axis = 1.25, las = 2)
abline(h = mu, lwd = 2)

for (i in 1:100){
  points(i, normmat[i, 1], lwd = 2, col = ifelse(cp.mean$true.in.ci[i] == 1,
                                                 "gray70", "gray20"), pch = 19)
  segments(i, cp.mean$ci[i, 1], i, cp.mean$ci[i, 2], lwd = 2,
           col = ifelse(cp.mean$true.in.ci[i] == 1, "gray70", "gray20"))
}
legend("bottomleft", bty = "n", c(expression("CI includes true"~mu),
                               expression("CI does not include true"~mu)),
       fill = c("gray70", "gray20"), cex = 1.0)