# ANOVA
# Dennis E. Slice (c) 2012

# Scrub the deck
rm(list = ls())

pause <- function() {
	print("Return...")
	readline()
}

# N(0,1) pdf
x = seq(-4, 4, 0.1)
y = dnorm(x)
ybar_10 = dnorm(x, sd = 1/sqrt(10))
ybar_20 = dnorm(x, sd = 1/sqrt(20))
ybar_50 = dnorm(x, sd = 1/sqrt(50))

yRange = range(c(ybar_10, ybar_20, ybar_50))
xRange = range(x)
plot(x, y, type = "l", ylim = yRange, main = "probability density", 
	ylab = "p", xlab = "x", lwd = 2)
pause()
lines(x, ybar_10, col = "red", lwd = 2)
pause()
lines(x, ybar_20, col = "green", lwd = 2)
pause()
lines(x, ybar_50, col = "blue", lwd = 2)
pause()

# Distribution of means.
nSamples = 1000
n = 50
mean_vec = rep(0, nSamples)
for (i in 1:nSamples) {
	mean_vec[i] = mean(rnorm(n))
}
hist(mean_vec, freq = FALSE, xlim = xRange/2, ylim = yRange, 
	col = "lightblue", main = "1000 sample means, n=50")
lines(x, y, lwd = 2, lt = 3)
lines(x, ybar_10, col = "red", lwd = 2, lt = 3)
lines(x, ybar_20, col = "darkgreen", lwd = 2, lt = 3)
lines(x, ybar_50, col = "blue", lwd = 3)