11

Does this code demonstrate the central limit theorem? This is not a homework assignment! Au contraire, I'm a faculty teaching some methods to non-stats students. 

library(tidyverse)
#Make fake data
population<-rnorm(1000000, mean=100, sd=10)

#Draw 100 samples of size 5
map(1:100, ~sample(population, size=5)) %>% 
  #calculate their mean
  map(., mean) %>% 
  #unlist 
  unlist() %>% 
  #draw histogram of sample means
  hist(, xlim=c(80,120))

#Repeat but with sample size 500
map(1:100, ~sample(population, size=500)) %>% 
  map(., mean) %>% 
  unlist() %>% 
  hist(., xlim=c(80,120))

#Repeat but with sample size 1000
map(1:100, ~sample(population, size=1000)) %>% 
  map(., mean) %>% 
  unlist() %>% 
  hist(., xlim=c(80,120))
spindoctor
  • 489
  • 6
  • 16
  • 11
    You would need to show more than that the distribution gets smaller, which demonstrates convergence (but not necessarily consistency); you would need to show that the distribution becomes more Gaussian. However, since your underlying population is awfully close to a Gaussian as it is, I doubt you'll see any significant deviations from Gaussianity even for a sample of size five. You might want to try an Exponential population distribution, then show how the histograms of the sample means approach a "bell curve" as $n \to \infty$. – jbowman Nov 26 '19 at 19:27
  • 3
    Even an exponential is pretty moderately skew, but it's a good starting place. ... spindoctor: Note that what you can demonstrate is simply the approach of distributions of standardized sample means more closely toward a standard normal distribution for a variety of population distributions. This is not really demonstrating the Central limit theorem, but is a closely related phenomenon. It may also be worth investigating (i) an example where the CLT doesn't apply, and (ii) also an extreme example where the CLT definitely still applies. e.g. (i) the Cauchy, & (ii) lognormal with $σ=4$, say. – Glen_b Nov 27 '19 at 00:44
  • 3
    quibble: "illustrate" might be a better verb that "demonstrate." In a mathematical context, the latter word suggests that it is a proof, but such computer simulations are not adequate for proving the CLT. – John Coleman Nov 27 '19 at 12:33
  • Note that such a test is not a proof of a theorem, it only provides some level of sanity checks. I.e. If I was a PhD advisor and a student claimed a similar result to the Central Limit Theorem, I'd suggest he would use something like this test to build confidence on whether the claim is right or wrong, as the result could be tainted by minor mistakes (though the code could be as well). Only after this sanity check passing I would delve into pages of proof to check its consistency. – Mefitico Nov 28 '19 at 16:05
  • @jbowman My bad. In the example that I use in my class, I actually do the exercise above but from a population that is heavily skewed. – spindoctor Dec 01 '19 at 01:17

3 Answers3

24

Here's a complete study in a few lines.

For a given set of sample sizes n and underlying distribution r, it generates n.sim independent samples of each size from that distribution, standardizes the empirical distribution of their means, plots the histogram, and overplots the standard Normal density in red. The CLT says that when the underlying distribution has finite variance, the red curve more and more closely approximates the histogram.

Figures

The first three rows illustrate the process for sample sizes of $10,20,100,500$ and underlying Normal, Gamma, and Bernoulli distributions. As sample size increases the approximation grows noticeably better. The bottom row uses a Cauchy distribution. Because a key assumption of the CLT (finite variance) does not hold in this case, its conclusion doesn't hold, which is pretty clear.

Execution time is about one second.

f <- function(n, r=rnorm,  n.sim=1e3, name="Normal", ...) {
  sapply(n, function(n) {
    x <- scale(colMeans(matrix(r(n*n.sim, ...), n))) # Sample, take mean, standardize
    hist(x, sub=name, main=n, freq=FALSE, breaks=30) # Plot distribution
    curve(dnorm(x), col="Red", lwd=2, add=TRUE)      # Compare to standard Normal
  })
}
n <- c(5,20,100,500)
mfrow.old <- par(mfrow=c(4,length(n)))
f(n)
f(n, rgamma, shape=1/2, name="Gamma(1/2)")
f(n, function(n) runif(n) < 0.9, name="Bernoulli(9/10)")
f(n, rt, df=1, name="Cauchy")
par(mfrow=mfrow.old)
whuber
  • 281,159
  • 54
  • 637
  • 1,101
5

Here's an example of one of my suggestions from comments. Means of samples of size n=100000 (takes about 20 seconds or so, be patient):

  ln.mean = replicate(1000,mean(rlnorm(100000,0,4)))
  hist(ln.mean,n=100)

histogram of means of sample size 100000

Even at this huge sample size, the distribution of sample means is still really skew -- but the central limit theorem nevertheless applies here - even the "classic" CLT.

Glen_b
  • 257,508
  • 32
  • 553
  • 939
  • +1 Another way to generate the same result uses the code in my answer: `f(10^5, rlnorm,name="Lognormal", sdlog=4)`. (16 seconds.) It's somewhat amusing that even means of sample sizes of $10^7$ are heavily skewed--but at least one can observe the progression towards Normality by then. – whuber Nov 27 '19 at 16:55
4

Maybe use something like the following (simpler, more direct) R code to show that averages of a dozen standard uniform random variables are difficult to distinguish from normal.

set.seed(1126)
a = replicate(5000, mean(runif(12))
shapiro.test(a)

        Shapiro-Wilk normality test

data:  a
W = 0.99965, p-value = 0.565

plot(qqnorm(a))

enter image description here

Then use R code to show that averages of 50, or even 100, standard exponential random variables are easy to distinguish from normal. What is the distribution of $A = \bar X_{100}?$

set.seed(1127)
a = replicate(5000, mean(rexp(100)))
shapiro.test(a)$p.val
 [1] 1.675877e-06

However, averages of 1000 standard exponentials are more difficult to distinguish from normal.

set.seed(1127)
a = replicate(5000, mean(rexp(1000)))
shapiro.test(a)$p.val
[1] 0.2413559
BruceET
  • 47,896
  • 2
  • 28
  • 76