central_limit.R

simulate <- function(nsim,nvec){
  simdx <- c()
  for(i in 1:length(nvec)) 
    simdx <- c(simdx,rep(1:nsim,each=nvec[i])+(i-1)*nsim)
  dt <- data.table(sim=simdx)
  bigN <- nrow(dt)
  dt$n <- rep(rep(nvec,nvec),each=nsim)
  dt$one <- 1
  dt$simc <- dt[,cumsum(one),by=sim]$V1
  dt$one <- NULL
  return(dt)
}

## Central limit theorem demo
dt <- simulate(100000,c(1,2,3,5))
bigN <- nrow(dt)
dt$x1 <- runif(bigN,-1,1)
dt$x2 <- runif(bigN,-0.5,0.5)
dt$x3 <- runif(bigN,-2,2)
zstats <- dt[,lapply(.SD,mean),by=sim]
zstats <- zstats[,.(n,x1*sqrt(12*n)/2,x2*sqrt(12*n)/1,x3*sqrt(12*n)/4)]
names(zstats) <- c('n','x1','x2','x3')
zstats

ggplot(zstats[n==1],aes(x=x1)) + geom_histogram(position="identity",binwidth = 0.01)
ggplot(zstats[n==2],aes(x=x1)) + geom_histogram(position="identity",binwidth = 0.01)
ggplot(zstats[n==3],aes(x=x1)) + geom_histogram(position="identity",binwidth = 0.01)
ggplot(zstats[n==5],aes(x=x1)) + geom_histogram(position="identity",binwidth = 0.01)