Monte Carlo

mu <- 7; sigma <- 2;
final=function(n, no.rep=1000){
l95 <- rep(NA,no.rep)
r95 <- rep(NA,no.rep)
l99 <- rep(NA,no.rep)
r99 <- rep(NA,no.rep)
for(i in 1:no.rep){         #重複1000次
print(i)
set.seed(i)
x <- rnorm(n,mu,sigma)
l95[i] <- mean(x)-qnorm(0.975)*sqrt(sigma^2/n)
r95[i] <- mean(x)+qnorm(0.975)*sqrt(sigma^2/n)
l99[i] <- mean(x)-qnorm(0.995)*sqrt(sigma^2/n)
r99[i] <- mean(x)+qnorm(0.995)*sqrt(sigma^2/n)
}
return(c(mean((l95<=mu) & (mu<=r95)),   # 檢查覆蓋率(coverage)
mean((l99<=mu) & (mu<=r99))))
}
i=1
result=matrix(NA,20,2)
for (n in seq(from=5,to=100,by=5)){
  result[i,]=final(n)
  i=i+1
}