2.6 Simulation

Any probability, average or distribution that is hard to work out on paper can be estimated by letting a computer play out the random process many times and counting. The answer gets more precise with the square root of the number of runs, and the same trick checks whether a statistical method really does what it promises.

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from scipy import stats

from statly import style

style.use()
rng = np.random.default_rng(1)

This book has used simulation from the first chapter: the many-small-influences heights, the central limit theorem, the bootstrap. This chapter makes the method itself the subject.

2.6.1 Estimating a probability by counting

Monte Carlo simulation estimates a probability as the share of simulated runs in which the event happens. Take the birthday problem: in a group of 23 people, how likely is it that two share a birthday? Simulate many groups of 23 random birthdays and count the groups with a match.

def shared_birthday(size: int, groups: int) -> np.ndarray:
    """For each simulated group, whether any two of `size` people share a birthday."""
    birthdays = rng.integers(0, 365, (groups, size))
    birthdays.sort(axis=1)
    return (np.diff(birthdays, axis=1) == 0).any(axis=1)


def exact(size: int) -> float:
    """1 minus the probability that all birthdays differ."""
    return 1 - np.prod((365 - np.arange(size)) / 365)


sizes = np.arange(2, 61)
counts = np.unique(np.geomspace(10, 20_000, 40).astype(int))
fig, axes = style.subplots(2, height=3.2, width=10.5)
axes[0].plot(sizes, [exact(n) for n in sizes], color=style.INK, lw=1.5)
axes[0].scatter(
    sizes, [shared_birthday(n, 2000).mean() for n in sizes], s=10, color=style.CORAL, zorder=3
)
axes[0].axvline(23, color=style.GREY, ls="--", lw=1)
axes[0].yaxis.set_major_formatter(style.PERCENT)
axes[0].set(xlabel="people in the group", ylabel="chance of a shared birthday")
for _ in range(20):
    outcomes = shared_birthday(23, counts.max())
    axes[1].plot(
        counts, np.cumsum(outcomes)[counts - 1] / counts, color=style.BLUE, lw=0.7, alpha=0.6
    )
axes[1].axhline(exact(23), color=style.INK, ls="--", lw=1.2)
axes[1].set(xscale="log", xlabel="simulated groups (log scale)", ylabel="estimated chance")
axes[1].yaxis.set_major_formatter(style.PERCENT)
fig.tight_layout()
plt.show()

Figure 2.6.1: Left: the probability that at least two people in a group share a birthday, by group size: exact (line) and estimated from 2,000 simulated groups per size (dots). It passes 50% at 23 people. Right: the estimate for 23 people as the number of simulated groups grows, in 20 independent runs; the estimates close in on the exact value (dashed) and their spread shrinks with the square root of the number of groups.

A Monte Carlo estimate is a share from N independent runs, so its precision follows the binomial (Chapter 2.2):

$$\text{SE}(\hat p) = \sqrt{\frac{p(1-p)}{N}}$$

Ten times more precision costs a hundred times more runs. For p near 50%, 10,000 runs give about ±1 percentage point (95%).

2.6.2 Areas, averages and integrals

The same idea estimates any average, and through it, areas and integrals. Scatter random points in a square; the share that falls inside a shape estimates its area. An expected value $E[g(X)]$ is estimated by the average of $g$ over simulated draws of $X$, which is how complicated expectations are computed in practice.

points = rng.uniform(-1, 1, (2000, 2))
inside = (points**2).sum(axis=1) <= 1
fig, ax = plt.subplots(figsize=(4.4, 4.4))
ax.scatter(*points[inside].T, s=4, color=style.BLUE)
ax.scatter(*points[~inside].T, s=4, color=style.GREY)
ax.add_patch(plt.Circle((0, 0), 1, fill=False, color=style.INK))
ax.set(
    aspect="equal",
    xticks=[],
    yticks=[],
    title=f"π ≈ 4 × {inside.mean():.3f} = {4 * inside.mean():.3f}",
)
plt.show()

Figure 2.6.2: Estimating π by throwing 2,000 random points into a square of side 2. The share landing inside the circle of radius 1 estimates its area divided by the square’s, π/4, so four times that share estimates π.

2.6.3 Propagating uncertainty

A project has five tasks, each with an uncertain duration. Adding up the most likely durations gives one number, and the project almost always takes longer: the durations are skewed (things overrun more than they finish early), and the sum of skewed uncertainties is not the sum of their most likely values. Simulating each task’s duration, adding them, and repeating gives the whole distribution of the total, and with it the answers that matter: how likely is it to finish within the deadline?

MOST_LIKELY = np.array([8, 12, 6, 15, 9])
SPREAD = 0.35  # standard deviation of the log duration


def task_durations(n: int) -> np.ndarray:
    """Log-normal durations whose most likely value (mode) is MOST_LIKELY."""
    mu = np.log(MOST_LIKELY) + SPREAD**2  # the mode of a log-normal is exp(mu − sigma²)
    return np.exp(rng.normal(mu, SPREAD, (n, 5)))


totals = task_durations(20_000).sum(axis=1)
naive = MOST_LIKELY.sum()
fig, ax = plt.subplots(figsize=(8, 3.2))
style.hist(ax, totals, bins=80)
ax.axvline(naive, color=style.INK, ls="--", lw=1.2)
for q, color in ((50, style.BLUE), (90, style.CORAL)):
    value = np.percentile(totals, q)
    ax.axvline(value, color=color, lw=1.5)
    ax.text(
        value + 0.5,
        ax.get_ylim()[1] * 0.9,
        f"{q}%: {value:.0f} days",
        color=color,
        fontsize=8.5,
    )
style.note(
    ax, f"P(total > {naive}) = {(totals > naive).mean():.0%}", x=0.62, y=0.6, color=style.INK
)
ax.set(xlabel="total project duration (days)", ylabel="simulated projects")
plt.show()

Figure 2.6.3: Five tasks with skewed durations (each a log-normal distribution with its most likely value marked). The sum of the most likely values is 50 days (dashed). In 20,000 simulated projects the total is above 50 days in most runs; the median is about 59 days and there is a 10% chance of exceeding 74.

2.6.4 Checking a method by simulation

Statistical methods come with promises: “this interval covers the truth 95% of the time”, “this test is wrong 5% of the time when nothing is going on”. Simulation checks a promise under conditions where the truth is known: generate many data sets from a known distribution, apply the method to each, and count how often it keeps its promise.

SIZES = np.array([5, 10, 20, 40, 80, 160])
DISTRIBUTIONS = {
    "normal": (stats.norm(), style.BLUE),
    "skewed (log-normal)": (stats.lognorm(1.2), style.CORAL),
}


def coverage(dist: stats.rv_continuous, n: int, runs: int = 10_000) -> float:
    """Share of t-intervals from samples of size n that contain the true mean."""
    samples = dist.rvs(size=(runs, n), random_state=rng)
    means, ses = samples.mean(axis=1), samples.std(axis=1, ddof=1) / np.sqrt(n)
    half = stats.t.ppf(0.975, n - 1) * ses
    return float(np.mean(np.abs(means - dist.mean()) <= half))


fig, ax = plt.subplots(figsize=(8, 3.2))
for name, (dist, color) in DISTRIBUTIONS.items():
    ax.plot(SIZES, [coverage(dist, n) for n in SIZES], "o-", color=color, lw=2, label=name)
ax.axhline(0.95, color=style.INK, ls="--", lw=1)
ax.set(
    xscale="log",
    xlabel="sample size (log scale)",
    ylabel="coverage of the 95% interval",
    ylim=(0.75, 1),
)
ax.set_xticks(SIZES, [str(n) for n in SIZES])
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8)
plt.show()

Figure 2.6.4: How often the usual 95% t-interval for a mean covers the true mean, for samples from a normal distribution (blue) and from a strongly skewed log-normal one (red), by sample size; 10,000 simulated samples per point. For normal data the promise holds at every size. For skewed data, small samples fall well short, and even at 160 observations the interval covers only about 91% of the time.

The same recipe gives the power of any test or design that has no formula: simulate data with the effect you hope to find, run the analysis, and count how often it detects the effect (Chapter 3.10).

2.6.5 Where the random numbers come from

Computers produce pseudo-random numbers: a deterministic sequence that passes every test of randomness. Starting it from the same seed gives the same numbers, which makes a simulation reproducible; every figure in this book is drawn with a fixed seed. Numbers from any distribution can be made from uniform ones by the inverse transform: if U is uniform between 0 and 1 and F is the cumulative curve of the wanted distribution (Chapter 2.1), then $F^{-1}(U)$ has that distribution.

uniform = rng.uniform(size=10_000)
values = stats.expon.ppf(uniform)
grid = np.linspace(0, 6, 300)

fig, axes = style.subplots(2, height=3.1, width=10)
axes[0].plot(grid, stats.expon.cdf(grid), color=style.CORAL, lw=2)
for u in (0.2, 0.6, 0.9):
    x = stats.expon.ppf(u)
    axes[0].annotate(
        "", xy=(x, u), xytext=(0, u), arrowprops=dict(arrowstyle="->", color=style.BLUE)
    )
    axes[0].annotate(
        "", xy=(x, 0), xytext=(x, u), arrowprops=dict(arrowstyle="->", color=style.BLUE)
    )
axes[0].set(xlabel="value", ylabel="uniform number U", xlim=(0, 6))
width = style.hist(axes[1], values, bins=60, range=(0, 6))
axes[1].plot(grid, stats.expon.pdf(grid) * len(values) * width, color=style.CORAL, lw=2)
axes[1].set(xlabel="value", ylabel="draws")
fig.tight_layout()
plt.show()

Figure 2.6.5: The inverse transform. Left: uniform random numbers (on the vertical axis) are mapped through the cumulative curve of an exponential distribution (red) to values on the horizontal axis; three examples are traced. Uniform numbers spread evenly up the curve land densely where the curve is steep. Right: 10,000 values made this way, with the exponential density.

2.6.6 In one picture

• Monte Carlo: play the random process many times and count; $\text{SE} = \sqrt{p(1-p)/N}$.
• Averages over simulated draws estimate expectations, areas and integrals.
• Simulating each uncertain input propagates uncertainty to the result.
• Simulation checks methods: coverage of intervals, error rates and power of tests.
• Fixed seeds make simulations reproducible; the inverse transform turns uniform numbers into any distribution.

Next: Part 3, how to test whether a pattern could be chance.