# Central limit theorem
# ---------------------


# A die
die.space <- c(1,2,3,4,5,6)
die.prob <- rep(1/6, 6)
die.mu <- sum(die.space*die.prob)
die.sd <- sqrt(sum((die.space-die.mu)^2*die.prob))
ddie <- function(x) { # probability mass function
  return(ifelse(x %in% die.space, die.prob[x], 0))
}
rdie <- function(n) {
  return(sample(die.space, n, replace=TRUE, prob=die.prob))
}


{
  plot(die.space, die.prob, type="h", col="blue", xlim=c(-1,8), ylim=c(0,0.3), lwd=3, ylab="p(x)", main=bquote(pmf~of~X), cex.main=2, cex.lab=1.5, axes=FALSE)
  lines(c(-1,8), c(0,0), col="blue", lwd=3)
  points(die.space, 0*die.space, col="white", pch=19)
  points(die.space, die.prob, col="blue", pch=19, cex=1.5)
  axis(1)
  axis(2)
  abline(v=die.mu, col="red", lwd=3)
}



sample_mean_die <- function(size, n=1000000) {
  z <- replicate(size, rdie(n))
  return(rowMeans(z))
}

plot_sample_mean_die_standardized <- function(size, n=1000000) {
  par(pty="s", mfrow=c(1,2))
  
  x <- sample_mean_die(size, ceiling(n/sqrt(size)))
  limdist <- function(z) {dnorm(z, mean=0, sd=die.sd)}
                         
  hist(x, main=bquote(Distribution~of~bar(X)[n]),
       probability=TRUE, breaks=seq(0.5,6.5, by=1/size), xlab="x", xlim=c(-1,8), ylim=c(0,0.3), cex.main=1.5, col=rgb(1,0,0,0.5), border=rgb(1,0,0,0.5))
  legend("topleft", legend=bquote(n==.(size)), cex=1.5, bty="n", x.intersp=0, y.intersp=0)
 
  hist(sqrt(size)*(x-die.mu), main=bquote(Distribution~of~sqrt(n)%.%(bar(X)[n]-mu[X])),
       probability=TRUE, breaks=sqrt(size)*(seq(0.5,6.5, by=2/size)-die.mu), xlab="x", xlim=c(-5.5,5.5), ylim=c(0,0.3), cex.main=1.5, col=rgb(1,0,0,0.5), border=rgb(1,0,0,0.5))
  legend("topleft", legend=bquote(n==.(size)), cex=1.5, bty="n", x.intersp=0, y.intersp=0)
  curve(limdist, col="blue", lwd=3, lty=2, add=TRUE)
  
  par(pty="m", mfrow=c(1,1))
}


plot_sample_mean_die_standardized(1)
plot_sample_mean_die_standardized(2)
plot_sample_mean_die_standardized(3)
plot_sample_mean_die_standardized(5)
plot_sample_mean_die_standardized(10)
plot_sample_mean_die_standardized(100)