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


# A trimodal distribution
m.1 <- -2
s.1 <- 1.5
m.2 <- 1
s.2 <- 0.5
m.3 <- 3
s.3 <- 1

trimod.mu <- 0
trimod.sd <- sqrt(1/4*(2*s.1^2 + 2*m.1^2 + s.2^2 + m.2^2 + s.3^2 + m.3^2))
dtrimod <- function(x) {
  return(dnorm(x,mean=m.1,sd=s.1)/2 + dnorm(x,mean=m.2,sd=s.2)/4 + dnorm(x,mean=m.3,sd=s.3)/4)
}
rtrimod <- function(n) {
  x.1 <- rnorm(n, mean=m.1, sd=s.1)
  x.2 <- rnorm(n, mean=m.2, sd=s.2)
  x.3 <- rnorm(n, mean=m.3, sd=s.3)
  u <- sample(c(1, 1, 2, 3), n, replace=TRUE)
  return(ifelse(u==1, x.1, ifelse(u==2, x.2, x.3)))
}


{
  curve(dtrimod, col="blue", xlim=c(-10,10), ylim=c(0,0.5), 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=trimod.mu, col="red", lwd=3)
}



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

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