Skip to content

Instantly share code, notes, and snippets.

@spikar spikar/central_limit.R
Last active Sep 3, 2019

Embed
What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.