Central Limit Theorem Simulation

What the theorem says

Let X₁, …, X_n be independent draws from a distribution with mean μ and finite variance σ², and let X̄ be their average. Two facts hold for every n, with no approximation:

The central limit theorem adds the shape: as n grows, the distribution of (X̄ − μ)/(σ/√n) approaches the standard normal. The lower chart overlays N(μ, σ²/n) on the histogram of means, and the readout means within μ ± 1.96σ/√n should approach 95%.

The upper chart shows the source itself. The ticks along its bottom are the latest sample and the blue line is that sample’s mean — one more count in the lower histogram. Use Step to draw samples one at a time.

How fast the means become normal

The theorem says nothing about how large n must be. Skewness is a good yardstick, because for a mean of n independent draws it is exactly the source’s skewness divided by √n:

• Uniform, die, bimodal — symmetric, skewness 0. Uniform means look normal by n = 5. The bimodal source is more instructive: means of two draws have three peaks, at −2, 0 and 2, before they merge.
• Exponential — skewness 2. At n = 30, the size often quoted as “large enough”, the means still have skewness 0.37; getting it down to 0.1 takes n = 400.
• Bernoulli, p = 0.1 — skewness (1 − 2p)/√(p(1 − p)) ≈ 2.67, and discrete: every mean is a multiple of 1/n, so the histogram bins are aligned to that grid. Even at n = 100 the skewness is 0.27.

The same σ/√n rule sets the accuracy of every Monte Carlo estimate; see how accurate a Monte Carlo estimate is.

How the simulation works

1. Draw n independent values from the source and average them. Repeat m times.
2. Compute the mean, standard deviation (dividing by m − 1) and skewness (third central moment over the second to the power 1.5) of the m means, and the share within μ ± 1.96σ/√n.
3. Histogram the means over μ ± 4σ/√n. For the die and Bernoulli sources every mean is a whole number divided by n, so bin edges fall halfway between those values; otherwise some bins would catch more lattice points than others and the histogram would show false spikes.

Example result

With seed 42, 10,000 means of n = 5 exponential draws averaged 0.9970 (μ = 1.0000) with a standard deviation of 0.4408 against σ/√n = 0.4472. Their skewness was 0.877, close to the theoretical 0.894, so the histogram still has a longer right tail than the normal curve; 95.7% of the means fell inside μ ± 1.96σ/√n.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Draw m samples of size n from the source and keep each sample's mean.
// The CLT says the means are approximately N(μ, σ²/n) for large n.
export function sampleMean(source, n, rng, values) {
  let sum = 0;
  for (let j = 0; j < n; j++) {
    const x = draw(source, rng);
    if (values) values[j] = x;
    sum += x;
  }
  return sum / n;
}

// Mean, standard deviation and skewness of the simulated means. Skewness is
// the third central moment over the 1.5th power of the second; the skewness
// of a mean of n independent draws is exactly (source skewness)/√n.
export function moments(xs, count) {
  let sum = 0;
  for (let i = 0; i < count; i++) sum += xs[i];
  const mean = sum / count;
  let m2 = 0, m3 = 0;
  for (let i = 0; i < count; i++) {
    const d = xs[i] - mean;
    m2 += d * d;
    m3 += d * d * d;
  }
  const sd = count > 1 ? Math.sqrt(m2 / (count - 1)) : null;
  m2 /= count;
  m3 /= count;
  return { mean, sd, skew: m2 > 0 ? m3 / (m2 * Math.sqrt(m2)) : null };
}

def sample_mean(source, n, rng):
    """Mean of n independent draws from the source distribution."""
    total = 0
    for _ in range(n):
        total += draw(source, rng)
    return total / n


def moments(xs):
    """Mean, sample standard deviation and skewness m3 / m2**1.5 of xs."""
    count = len(xs)
    mean = sum_in_order(xs) / count
    m2 = m3 = 0
    for x in xs:
        d = x - mean
        m2 += d * d
        m3 += d * d * d
    sd = math.sqrt(m2 / (count - 1)) if count > 1 else None
    m2 /= count
    m3 /= count
    return mean, sd, (m3 / (m2 * math.sqrt(m2)) if m2 > 0 else None)

Run the full model from a terminal:

node simulation.js --source exponential --n 5 --seed 42
python simulation.py --source exponential --n 5 --seed 42

Download the model

Assumptions & limitations

• Independent draws from a fixed distribution. The theorem needs independence (or weak dependence) and finite variance. Every source here qualifies; a Cauchy distribution does not, and means of Cauchy draws never settle down.
• Display ranges. The means histogram spans μ ± 4σ/√n and the exponential source is drawn up to 6; values beyond are rare and still count in every statistic.
• Skewness is noisy. The sample skewness of the means has its own sampling error, so small differences from the theoretical value are expected, especially with few means.
• About means, not data. The theorem is about the distribution of the average. It does not make the underlying data normal.

Related simulations

• OS-010Galton boardBalls bouncing left or right through rows of pegs pile up into a binomial bell curve.
• OS-007Dice rollRoll k dice with s sides and check the histogram of sums against the exact distribution.
• OS-006Coin flipFlip a fair or biased coin: the share of heads settles down, and long streaks keep happening.
• OS-001Monte Carlo πEstimate π from random points in a unit square and watch the error shrink like 1/√n.

Background: How Accurate Is a Monte Carlo Estimate?

References

1. Grinstead, C. M. & Snell, J. L. (1997). Introduction to Probability, 2nd ed. American Mathematical Society. math.dartmouth.edu/~prob/prob/prob.pdf
2. Feller, W. (1968). An Introduction to Probability Theory and Its Applications, Vol. 1, 3rd ed. Wiley.
3. Owen, A. B. (2013). Monte Carlo theory, methods and examples. Stanford University. artowen.su.domains/mc