Confidence interval simulation: cola example

Here is a simulation where we take samples of size 6 from a normal distribution with mean 355 (pretend this is the unknown) and standard deviation 2, which we assume is a known standard deviation. We take 100 intervals, keep track of how many of the 100 successfully capture the true mean of 355 (hits) and then graph the 100 intervals to see a graphical display of how confidence intervals work. Remember, the probability is in the method.

sig <- 2  # we pretend this is the known sd
n <- 6   
mu <- 355  # we pretend this is the unknown mean
z <- qnorm(.975)  # we simulate 95% confidence intervals
hits <- 0  # hits keeps track of times the CI captures mu
ucl <- numeric(100) # ucl=upper confidence limit for i-th interval
lcl <- numeric(100) # lcl=lower confidence limit for i-th interval
x <- numeric(100)  # x collects sample means

for (i in 1:100) {
    x[i] <- mean(rnorm(n,mu,sig))  # i-th sample mean
    lcl[i] <- mean(x[i]) - z*sig/sqrt(n)  # lower confidence limit
    ucl[i] <- mean(x[i]) + z*sig/sqrt(n)  # upper confidence limit
    hits <- hits + (mu <= ucl[i] & lcl[i] <= mu)
    }

plot(1,1,xlim=c(350,360),ylim=c(0,101),type='n')
for (i in 1:100) {
     # we'll plot CI schematics from bottom up
     lines(c(lcl[i],ucl[i]),c(i,i))
     points(x[i],i,pch=16,cex=.5,)
     }
abline(v=mu,lty=2)