kevinushey

case2 <- function(n = 10, mu = 3, sigma = sqrt(5), p = 0.025, rep = 100){
  
  xbar <- vapply(1:rep, FUN.VALUE=numeric(1), function(i) {
    norm <- rnorm(mean = mu, sd = sigma, n = n)
    return(mean(norm))
  })
  
  low <- xbar - qnorm(1-p) * (sigma / sqrt(n))
  up <- xbar + qnorm(1-p) * (sigma / sqrt(n))
  rem <- mu > low & mu < up