n = 100
mu = 25
sigma = 3

y = rnorm(n, mu, sigma)

plot(density(y))
rug(x)

y.mean = mean(y); y.mean
y.sd = sd(y)

# 95% Confidence Interval 
y.ci = t.test(y)$conf.int[1:2]

# 95% of the population, not a PI

#y.popi = c(mu - 1.96*sigma, mu + 1.96*sigma)

y.popi = c( y.mean + qt(.025,df=n-1)*y.sd, 
          y.mean + qt(.975,df=n-1)*y.sd)

# 95% Prediction Interval

y.pi = c( y.mean + qt(.025,df=n-1)*y.sd*sqrt(1+1/n), 
          y.mean + qt(.975,df=n-1)*y.sd*sqrt(1+1/n))


y.ci
y.popi
y.pi

# Simulation of the 95% confidence level of a 95% prediction interval.
# Two levels of simulation.

B = 10000
Count.ci = 0
Count.popi = 0 
Count.pi = 0

for (i in 1:B){
	y = rnorm(n, mu, sigma)
      y.mean = mean(y)
	y.sd = sd(y)
	y.new = rnorm(1, y.mean, y.sd)
      
      y.ci = t.test(y)$conf.int[1:2]
	if (y.new > y.ci[1] & y.new < y.ci[2]) Count.ci = Count.ci + 1

	y.popi = c( y.mean + qt(.025,df=n-1)*y.sd, y.mean + qt(.975,df=n-1)*y.sd)
	if (y.new > y.popi[1] & y.new < y.popi[2]) Count.popi = Count.popi + 1

	y.pi = c( y.mean + qt(.025,df=n-1)*y.sd*sqrt(1+1/n), 
          y.mean + qt(.975,df=n-1)*y.sd*sqrt(1+1/n))
	if (y.new > y.pi[1] & y.new < y.pi[2]) Count.pi = Count.pi + 1
}

Count.ci/B
Count.popi/B
Count.pi/B