###--------------------------------------------------------------------------- ### use animation package ###--------------------------------------------------------------------------- install.packages("animation") library(animation) ###------------------------------------------------------- ###ˇˇLLN: Law of large number lln<-function (FUN, pars=NULL, np = 30, n = ani.options("nmax"),pch = 20,col.poly = "bisque", col.mu = "gray", ...) { dist.name<-deparse(substitute(FUN)) if(dist.name=='rbinom'){FUN<-function(n,pars) rbinom(n,size=pars[1],prob=pars[2]);mu<-pars[2];} if(dist.name=='rpois'){FUN<-function(n,pars) rpois(n,lambda=pars);mu<-pars;} if(dist.name=='rnorm'){FUN<-function(n,pars) rnorm(n,mean=pars[1],sd=pars[2]);mu<-pars[1];} if(dist.name=='rexp'){FUN<-function(n,pars) rexp(n,rate=pars); mu<-1/pars;} if(dist.name=='runif'){FUN<-function(n,pars) runif(n,min=pars[1],max=pars[2]);mu<-sum(pars)/2;} if(dist.name=='rchisq'){FUN<-function(n,pars) rchisq(n,df=pars);mu<-pars;} m = x = NULL for (i in 1:n) { d = colMeans(matrix(replicate(np, FUN(i*100,pars)), i*100)) m = c(m, d) x = rbind(x, range(d)) } rg = range(m) xax = pretty(1:n) for (i in 1:n) { dev.hold() plot(1:n, ylim = rg, type = "n", xlab = paste("n =100*", i), ylab = expression(bar(x)), xaxt = "n",main=dist.name) axis(1, xax[xax <= i]) polygon(c(1:i, i:1), c(x[1:i, 1], x[i:1, 2]), border = NA, col = col.poly) points(rep(1:i, each = np), m[1:(i * np)], pch = pch, ...) abline(h = mu, col = col.mu) ani.pause() } } #LLN for Binomial lln(FUN=rbinom,pars=c(1,0.5)) #LLN for Poisson lln(FUN=rpois,pars=2) #LLN for Uniform lln(FUN=runif,pars=c(0,1)) #LLN for Exponential lln(FUN=rexp,pars=2) ###------------------------------------------------------- ###ˇˇCLT: Central limit theorem clt<-function (obs = 300, FUN =rexp, mu=0,sds=1,nmax = ani.options("nmax"),col = c("bisque", "red", "blue", "black"),xlim, ...) { x = matrix(nrow = nmax, ncol = obs) for (i in 1:nmax) x[i, ] = apply(matrix(replicate(obs, FUN(i)), i), 2, mean) if (missing(xlim)) xlim = quantile(x, c(0.005, 0.995)) for (i in 1:nmax) { dev.hold() hist(x[i, ], freq = FALSE, main = "", xlab = substitute(italic(bar(x)[i]), list(i = i)), col = col[1], xlim = xlim) lines(density(x[i, ]), col = col[2],lwd=2) if(!is.na(mu) && !is.na(sds)) curve(dnorm(x, mu, sds/sqrt(i)), col = col[3], lty = 2, lwd=2, add = TRUE) legend("topright", legend = c("Normal","Est. pdf"),lty=2:1, lwd=2, col=c(col[3],col[2]), bty = "n") ani.pause() } } ani.options(interval = 0.5) par(mar = c(3, 3, 1, 0.5), mgp = c(1.5, 0.5, 0), tcl = -0.3) #Poisson case f<-function(n) rpois(n,lambda=4); clt(FUN = f, mu=4,sds=2) #binomial case f<-function(n) rbinom(n,size=1,prob=0.5) clt(FUN = f, mu=0.5, sds=0.5) #exponential distribution case f<-function(n) rexp(n,rate=2); clt(FUN = f, mu=1/2,sds=1/2) #uniform distribution case f<-function(n,pars) runif(n,min=0,max=1); clt(FUN = f,mu=1/2,sd=1/sqrt(12)) #chi-square distribution f<-function(n) rchisq(n,df=2); clt(FUN = f,mu=2,sd=2) #When CLT does not work f = function(n) rcauchy(n, location = 0, scale = 2) clt(FUN = f, mean = NA, sd = NA) ###--------------------------------------------------------------------------- ###Alternative methods ###--------------------------------------------------------------------------- # # law of large number # set.seed(123) numsim = 100 matx = matrix(rnorm(numsim^2), nrow=numsim) runavg = function(x) { cumsum(x)/(1:length(x)) } # runavg calculates the running average of a vector ramatx = t(apply(matx, 1, runavg)) # apply() we can get the running average of each row. matplot(1:numsim,t(ramatx),pch=20, cex=0.2,type="l") # # sample mean from normal distribution# # set.seed(234) mu=100; sigma=10 n=10 xbar=rep(0,500) for (i in 1:500) {xbar[i]=mean(rnorm(n,mean=mu,sd=sigma))} hist(xbar,prob=TRUE,breaks="FD",xlim=c(80,120)) n=500 xbar=rep(0,500) for (i in 1:500) {xbar[i]=mean(rnorm(n,mean=mu,sd=sigma))} hist(xbar,prob=TRUE,breaks="FD",col=rgb(1,0,0,1/4), xlim=c(80,120),add=T) # #illustration of the Central Limit Theorem # plot_sample_means(runif, n=1, title="Sample means from uniform distribution") plot_sample_means(runif, n=2, title="Sample means from uniform distribution") plot_sample_means(runif, n=10, title="Sample means from uniform distribution") plot_sample_means(rexp, n=6, title="Sample means from the exponential distribution", rate=1) plot_sample_means(rexp, n=12, title="Sample means from the exponential distribution", rate=1) plot_sample_means(rexp, n=48, title="Sample means from the exponential distribution", rate=1) #CLT does not work for Cauchy! plot_sample_means(rcauchy, n=48, title="Sample means from the Cauchy distribution") # # required program # plot_sample_means <- function(f_sample, n, m=300,title="Histogram", ...) { # define a vector to hold our sample means means <- double(0) # generate 300 samples of size n and store their means for(i in 1:m) means[i] = mean(f_sample(n,...)) # scale sample means to plot against standard normal scaled_means <- scale(means) # set up a two panel plot par(mfrow=c(1,2)) par(mar=c(5,2,5,1)+0.1) # plot histogram and density of scaled means hist(scaled_means, prob=T, col="light grey", border="grey", main=NULL, ylim=c(0,0.4)) lines(density(scaled_means)) # overlay the standard normal curve in blue for comparison curve(dnorm(x,0,1), -3, 3, col='blue', add=T) # adjust margins and draw the quantile-quantile plot par(mar=c(5,1,5,2)+0.1) qqnorm(means, main="") # return margins to normal and go back to one panel par(mar=c(5,4,4,2)+0.1) par(mfrow=c(1,1)) # add a title par(omi=c(0,0,0.75,0)) title(paste(title, ", n=", n, sep=""), outer=T) par(omi=c(0,0,0,0)) # return unscaled means (without printing) return(invisible(means)) }