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


# A bimodal distribution
m.1 <- -1
s.1 <- 0.5
m.2 <- 1
s.2 <- 0.7
bimod.mu <- 0
bimod.sd <- sqrt(1/2*(s.1^2 + m.1^2 + s.2^2 + m.2^2))
dbimod <- function(x) {
  return((dnorm(x,mean=m.1,sd=s.1) + dnorm(x,mean=m.2,sd=s.2))/2)
}
rbimod <- function(n) {
  x.1 <- rnorm(n, mean=m.1, sd=s.1)
  x.2 <- rnorm(n, mean=m.2, sd=s.2)
  u <- sample(c(0,1), n, replace=TRUE)
  return(u*x.1 + (1-u)*x.2)
}


{
  curve(dbimod, col="blue", xlim=c(-4,4), ylim=c(0,1), lwd=3, ylab="f(x)", main=bquote(pdf~of~X), cex.main=2, cex.lab=1.5, axes=FALSE)
  axis(1)
  axis(2)
  abline(v=bimod.mu, col="red", lwd=3)
}



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

plot_sample_mean_bimod <- function(size, n=1000000) {
  hist(sample_mean_bimod(size, ceiling(n/size)),
       main=bquote(Distribution~of~bar(X)[n]),
       probability=TRUE, breaks="Scott", xlab="x", xlim=c(-4,4), ylim=c(0,1), cex.main=1.5)
  abline(v=bimod.mu, col="red", lwd=3)
  legend("topleft", legend=bquote(n==.(size)), cex=1.5, bty="n", x.intersp=0, y.intersp=0)
}



plot_sample_mean_bimod(1)
plot_sample_mean_bimod(2)
plot_sample_mean_bimod(3)
plot_sample_mean_bimod(5)
plot_sample_mean_bimod(10)
plot_sample_mean_bimod(100)
plot_sample_mean_bimod(1000)


plot_sample_mean_bimod_standardized <- function(size, n=1000000) {
  par(pty="s", mfrow=c(1,2))
  
  x <- sample_mean_bimod(size, ceiling(n/sqrt(size)))
  limdist <- function(z) {dnorm(z, mean=0, sd=bimod.sd)}
                         
  hist(x, main=bquote(Distribution~of~bar(X)[n]),
       probability=TRUE, breaks="Scott", xlab="x", xlim=c(-4,4), ylim=c(0,1), 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, main=bquote(Distribution~of~sqrt(n)%.%bar(X)[n]),
       probability=TRUE, breaks="Scott", xlab="x", xlim=c(-4,4), ylim=c(0,1), 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_bimod_standardized(1)
plot_sample_mean_bimod_standardized(2)
plot_sample_mean_bimod_standardized(3)
plot_sample_mean_bimod_standardized(5)
plot_sample_mean_bimod_standardized(10)
plot_sample_mean_bimod_standardized(100)