## Chapter 4
##
## Code to support Video Vignette
##
## Videos and supporting code are not a complete portrayal
## of Chapter content.

##
## a: Central Limit Theorem
##

## exponentials are not very Gaussian
y <- seq(0, 5, length=1000)
plot(y, dexp(y), type="l", lwd=2)

## E(Y) = 1 using defaults
## Var(Y) = 1 too, so CLT says N(1, 1/sqrt)

## means from samples of size n=2
y2 <- colMeans(matrix(rexp(2*1000), ncol=1000))
hist(y2, main="", freq=FALSE, xlab="averages of 1000 samples of size 2")
lines(y, dnorm(y, 1, 1/sqrt(2)), col="red")
legend("topright", "AN", lty=1, col=2, bty="n")

## means from samples of size n=5
y2 <- colMeans(matrix(rexp(5*1000), ncol=1000))
hist(y2, main="", freq=FALSE, xlab="averages of 1000 samples of size 5")
lines(y, dnorm(y, 1, 1/sqrt(5)), col="red")
legend("topright", "AN", lty=1, col=2, bty="n")

## means from samples of size n=10
y2 <- colMeans(matrix(rexp(10*1000), ncol=1000))
hist(y2, main="", freq=FALSE, xlab="averages of 1000 samples of size 10")
lines(y, dnorm(y, 1, 1/sqrt(10)), col="red")
legend("topright", "AN", lty=1, col=2, bty="n")

## means from samples of size n=100
y2 <- colMeans(matrix(rexp(100*1000), ncol=1000))
hist(y2, main="", freq=FALSE, xlab="averages of 1000 samples of size 100")
lines(y, dnorm(y, 1, 1/sqrt(100)), col="red")
legend("topright", "AN", lty=1, col=2, bty="n")

## means from samples of size n=1000
y2 <- colMeans(matrix(rexp(1000*1000), ncol=1000))
hist(y2, main="", freq=FALSE, xlab="averages of 1000 samples of size 1000")
lines(y, dnorm(y, 1, 1/sqrt(1000)), col="red")
legend("topright", "AN", lty=1, col=2, bty="n")

##
## Coin flipping example from Ch 1
##

## buddy's coin flips
data <- c("H", "T", "H", "H", "H", "T", "T", "H", "T", "T", "H", "H", "H",
          "T", "T", "H", "H", "H", "H", "T", "H", "H", "H", "T", "H", "H", "T", 
          "H", "T", "T")

## convert to binary
y <- as.numeric(data == "H")
n <- length(y)

## using statistic ybar, but could use thetahat from Ch 3
ybar <- mean(y)

se <- sqrt(ybar*(1 - ybar)/n)               ## that (=ybar) & n in Ch 1&3 
z <- (ybar - 0.5)/se                        ## testing H0: theta = 0.5
p.clt <- 2*pnorm(-abs(z)) 
p.clt                                       ## in Ch1 we got 0.36

##
## discussion about accuracy between CLT and MC simulation from Ch1
##