# graphical representation of the central limit theorem and the unbiased estimation of the mean # for more R-scripts: www.random-stuff.de pdf("central-limit-theorem.pdf", width=8, height=7) par(bty="n", mar=c(6, 6, 3, 3.5)) z1 <- rnorm(4500, mean=20, sd=5) z2 <- rnorm(3600, mean=28, sd=5) z3 <- rnorm(1300, mean=40, sd=7) z4 <- rnorm(1200, mean=10, sd=7) z5 <- rnorm(1500, mean=50, sd=5) z6 <- rnorm(1500, mean=10, sd=5) z7 <- rnorm(1400, mean=37, sd=3) pop <- c(z1, z2, z3, z4, z5, z6, z7) meansN2 <- numeric() meansN10 <- numeric() meansN40 <- numeric() for (i in 1:1000) { x <- runif(2, min=1, max=15000) x <- round(x) sample <- c(pop[x]) meansN2 <- append(meansN2, mean(sample)) } for (i in 1:1000) { x <- runif(10, min=1, max=15000) x <- round(x) sample <- c(pop[x]) meansN10 <- append(meansN10, mean(sample)) } for (i in 1:1000) { x <- runif(40, min=1, max=15000) x <- round(x) sample <- c(pop[x]) meansN40 <- append(meansN40, mean(sample)) } pop2 <- c(pop) pop2 <- smooth(pop) pop <- table(cut(pop,breaks=seq(from=-10, to=70, length.out=80))) meansN2 <- table(cut(meansN2,breaks=seq(from=-10, to=70, length.out=80))) meansN10 <- table(cut(meansN10,breaks=seq(from=-10, to=70, length.out=80))) meansN40 <- table(cut(meansN40,breaks=seq(from=-10, to=70, length.out=80))) pop <- smooth(pop) meansN2 <- smooth(meansN2) meansN10 <- smooth(meansN10) meansN40 <- smooth(meansN40) plot(pop, type="l", lwd=1.5, xlab="variable", ylab="frequency", ylim=c(0,525)) lines(meansN2, lty="dotted", lwd=1.3) lines(meansN10, lty="dashed", lwd=1.3) lines(meansN40, lty="twodash", lwd=1.3) legend(55, y=467,c("population","n = 40","n = 10", "n = 2"), lty=c("solid","twodash","dashed", "dotted"), bty="n", y.intersp=1.16, lwd=1.2, cex=0.9) box(which = "inner", lty = "solid") lines(x=c(36.3, 36.3),y=c(-7, 515),col="#9E9E9E",lwd=2) text(x= 36.3,y= 535, cex=1.42, label=expression(paste(mu)),col="#9E9E9E") #quartz.save("central-limit-theorem.pdf", type="pdf") # for Mac PCs dev.off() print(paste("Plot was saved in:", getwd()))