Demo: Why the Bell Curve Is Everywhere

Roll one die and every face is equally likely, so the distribution is flat. Roll ten dice and add them up.

What shape does the distribution of the sum make, and where does it peak? Sketch your guess before we run anything.

import numpy as np
import matplotlib.pyplot as plt
import scienceplots

plt.style.use(["science", "no-latex"])
TEAL = "#009090"

die = np.ones(6) / 6  # a fair die: each face has probability 1/6

fig, axes = plt.subplots(1, 4, figsize=(9, 2.4))
pmf = die.copy()
for ax, k in zip(axes, [1, 2, 3, 10]):
    while len(pmf) < 6 * k - k + 1:      # convolve until pmf is the sum of k dice
        pmf = np.convolve(pmf, die)
    ax.bar(np.arange(k, 6 * k + 1), pmf, color=TEAL, width=0.7)
    ax.set_xlabel(f"Sum of {k} dice" if k > 1 else "One die")
axes[0].set_xticks(np.arange(1, 7))
axes[0].set_ylabel("Probability")
fig.tight_layout()
plt.show()

Flat, then triangular, then… a bell.

Can we predict which bell? From lecture, the expectation of one die is \(\mathbb{E}[X] = 3.5\) and its variance is \(\mathrm{Var}(X) = \frac{35}{12}\). By linearity of expectation, the sum of \(k\) dice has expectation \(3.5k\). Because the dice are independent, every covariance term in the variance of a sum drops out, so the sum has variance \(\frac{35}{12}k\).

So the claim is that the sum of ten dice follows \(\mathcal{N}(35,\ 350/12)\), with nothing fitted to the bars and nothing used but the two numbers from lecture.

k = 10
pmf = die.copy()
for _ in range(k - 1):
    pmf = np.convolve(pmf, die)
values = np.arange(k, 6 * k + 1)

mu, var = 3.5 * k, (35 / 12) * k
x = np.linspace(values[0], values[-1], 500)
bell = np.exp(-((x - mu) ** 2) / (2 * var)) / np.sqrt(2 * np.pi * var)

fig, ax = plt.subplots(figsize=(7, 3.2))
ax.bar(values, pmf, color=TEAL, width=0.8, label="Sum of 10 dice (exact)")
ax.plot(x, bell, color="black", linewidth=1.6, label=f"$\\mathcal{{N}}({mu:.0f},\\ {var:.1f})$ (predicted)")
ax.set_xlabel("Sum")
ax.set_ylabel("Probability")
ax.legend(frameon=False)
plt.show()

What if the die is loaded?

Maybe the bell only appears because the die was fair? Below is a heavily loaded die. The recipe doesn’t change: compute its expectation and variance, multiply each by \(k\), and predict the bell.

Change the probabilities to anything you like (they just need to sum to 1). Can you design a die that escapes the bell?

loaded = np.array([0.40, 0.10, 0.05, 0.05, 0.10, 0.30])  # try your own!
faces = np.arange(1, 7)

mu1 = (faces * loaded).sum()               # E[X]
var1 = ((faces - mu1) ** 2 * loaded).sum() # Var(X)

k = 30
pmf = loaded.copy()
for _ in range(k - 1):
    pmf = np.convolve(pmf, loaded)
values = np.arange(k, 6 * k + 1)

mu, var = k * mu1, k * var1
x = np.linspace(values[0], values[-1], 500)
bell = np.exp(-((x - mu) ** 2) / (2 * var)) / np.sqrt(2 * np.pi * var)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 2.8))
ax1.bar(faces, loaded, color=TEAL, width=0.6)
ax1.set_xlabel("One loaded die")
ax1.set_ylabel("Probability")
ax2.bar(values, pmf, color=TEAL, width=0.8, label=f"Sum of {k} loaded dice")
ax2.plot(x, bell, color="black", linewidth=1.6, label=f"$\\mathcal{{N}}({mu:.1f},\\ {var:.1f})$")
ax2.set_xlabel("Sum")
ax2.legend(frameon=False, loc="upper right")
fig.tight_layout()
plt.show()

Punchline: sum enough independent effects and the result is Gaussian. The shape of the individual pieces washes out, and only their expectation and variance survive (the central limit theorem).

Measurement noise is exactly such a sum of many small independent effects. That is why, when we derive the mean squared error loss in the Linear Models unit, our noise model will be the Gaussian.