Illustration of the Central Limit Theorem with R

Preliminars

The runif function generates draws from a uniform distribution:
runif(5, min = 2, max = 3)
## [1] 2.686 2.684 2.099 2.413 2.590
The rnorm function generates draws from a normal distribution:
rnorm(5, mean = 0, sd = 1)
## [1] -2.13356  0.20484 -1.59384 -0.07009  1.03166
The density function computes kernel density estimates:
# estimation of the density of the uniform distribution
plot(density(runif(10000, min = 2, max = 3)))
plot of chunk unnamed-chunk-3
# estimation of the density of the uniform distribution
plot(density(rnorm(10000, mean = 0, sd = 1)))
plot of chunk unnamed-chunk-4

Detailed code

The following function performs an animated plot of the mean of \(k\) uniform random variables, with \(k\) varying from \(1\) to \(n\), and using \(N\) individuals. The code was initially written by Bill Venables, updated by Ted Harding, and then by Emmanuel Rousseaux.
clt.illustration <- function(n=30, N=10000, from = NA){
  if(is.na(from)) {
    graphics.off() # clean graphical device
    from <- 1 # we start with 1 random variable
  }
  par(mfrow = c(1,2), pty = "s") # split the screen in 2, for plotting the density
  # and the QQ-plot
  for(k in from:n) {
    # we store our k r.v. in a matrix with N rows and k columns
    k.rv.unif <- matrix(runif(N*k), N, k)
    # we compute the mean for each row to have the estimation of the mean r. v.
    rv.mean <- rowMeans(k.rv.unif)
    # centering and reduction
    rv.mean.centered.reduced <- (rv.mean - 0.5)*sqrt(12*k) # Bienaymé formula
    # we compute an histogram
    h <- hist(
      rv.mean.centered.reduced,
      breaks = "FD", # we use Freedman-Diaconis algorithm to estimate cut-offs (a
      # range-fixed cutting)
      plot=FALSE
    )
    # we plot the histogram
    if(k == 1) {s <- ''} else {s <- 's'}
    plot(h,xlim = c(-4,4),main = paste0(k, ' random variable', s),
         xlab = "Centered and reduced mean R.V. ",
         ylab = "Frequency"
    )
    # then we want to plot a centered reduced normal law
    # we first need to get min and max bounds on x abscissa
    x.min <- par("usr")[1]
    x.max <- par("usr")[2]
    # we generate 500 equidistant points on this range
    x <- seq(x.min, x.max, length.out = 500)
    range.between.2.cut.offs <- h$breaks[2]-h$breaks[1]
    y <- N*dnorm(x)*(range.between.2.cut.offs) #dnorm compute the value of the
    # normal density function with mean = 0 and sd = 1 by default
    lines(
      x, # abscissa
      y, # ordinate
      col = "red"
    )
    # we plot the QQ-plot
    qqnorm(
      rv.mean.centered.reduced,
      ylim = c(-4,4),
      xlim = c(-4,4),
      pch = ".",
      col = "blue"
    )
    abline(0, 1, col = "red")
    Sys.sleep(1) # 1s sleep to have time to watch
  }
}

Illustration

Here are snapshots for some values of \(k\).
plot of chunk unnamed-chunk-6
plot of chunk unnamed-chunk-7
plot of chunk unnamed-chunk-8