# Law of large numbers
# --------------------


# A bimodal distribution
m.1 <- -1
s.1 <- 0.5
m.2 <- 1
s.2 <- 0.7
mu <- 0
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(-5,5), ylim=c(0,1.5), lwd=3, ylab="f(x)", main="Probability density function", cex.main=2, cex.lab=1.5, axes=FALSE)
  axis(1)
  axis(2)
  abline(v=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) {
  if(size==1) {
    title <- "Distribution of 1 sample"
  } else {
    title <- bquote(Distribution~of~mean~of~.(size)~independent~samples)
  }
  hist(sample_mean_bimod(size, ceiling(n/size)),
       main=title,
       probability=TRUE, breaks=seq(-5,5,by=0.1/sqrt(sqrt(size))), xlab="x", ylim=c(0,1.5), cex.main=1.5)
  abline(v=mu, col="red", lwd=3)
}

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