10
$\begingroup$

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))



$\endgroup$
  • $\begingroup$ How could samples of size 5 demonstrate anything about a limiting theorem? Could you explain? $\endgroup$ – whuber Nov 26 '19 at 18:46
  • $\begingroup$ Sorry. I have edited the answer to show that I try to show students the effect of increasing the sample size, i.e. that the resulting distribution of sample means gets smaller $\endgroup$ – spindoctor Nov 26 '19 at 18:58
  • 11
    $\begingroup$ 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$. $\endgroup$ – jbowman Nov 26 '19 at 19:27
  • 3
    $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Nov 27 '19 at 0:44
  • 3
    $\begingroup$ 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. $\endgroup$ – John Coleman Nov 27 '19 at 12:33
23
$\begingroup$

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)
$\endgroup$
  • 1
    $\begingroup$ The return value of par are the original values of the parameters that you have changed so oldpar <- par(…) and later par(oldpar) would be easier and more reliable if the default should have been changed somewhere else. $\endgroup$ – Martin Ueding Nov 27 '19 at 17:10
  • 2
    $\begingroup$ @Martin Thank you; that is better practice. $\endgroup$ – whuber Nov 27 '19 at 17:13
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ +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. $\endgroup$ – whuber Nov 27 '19 at 16:55
4
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.