Each common distribution is the fingerprint of a mechanism. Know how the values came about, and you know which shape to expect.
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from sklearn import datasets
from scipy import stats
from statly import style
style.use()
rng = np.random.default_rng(1)
engel = (
sm.datasets.engel.load_pandas().data
) # income and food spending, Belgian households (1857)
iris = datasets.load_iris(as_frame=True).frame # flowers of three iris speciesA distribution can be drawn in two ways. The density from Chapter 1.1 shows how common each value is, and the area under it is probability. The cumulative curve adds that area up from the left: at each value it shows the share of values below it.
MEAN, SD = 170, 7 # heights in cm
heights = np.linspace(MEAN - 4 * SD, MEAN + 4 * SD, 400)
below_180 = stats.norm.cdf(180, MEAN, SD)
def plot_density(ax: Axes) -> None:
"""The density, with the area below 180 cm shaded."""
curve = stats.norm.pdf(heights, MEAN, SD)
ax.plot(heights, curve, color=style.BLUE)
shaded = heights <= 180
ax.fill_between(heights[shaded], curve[shaded], color=style.BLUE, alpha=0.25, lw=0)
style.note(ax, f"area below 180 cm:\n{below_180:.0%}", x=0.03)
ax.set(title="Density", xlabel="height (cm)", ylabel="how common", yticks=[])
def plot_cumulative(ax: Axes) -> None:
"""The cumulative curve: the share of people below each height."""
ax.plot(heights, stats.norm.cdf(heights, MEAN, SD), color=style.BLUE)
ax.plot(
[180, 180, heights[0]], [0, below_180, below_180], color=style.CORAL, ls="--", lw=1.2
)
style.note(ax, f"{below_180:.0%} are\nbelow 180 cm", x=0.03)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.set(title="Cumulative", xlabel="height (cm)", ylabel="share below")
fig, axes = style.subplots(2, height=3.1, width=9)
plot_density(axes[0])
plot_cumulative(axes[1])
fig.tight_layout()
plt.show()
Read the cumulative curve backwards and you get percentiles. Go across from 50% and down: that is the median, the height half the people are below. Go across from 90%, and you find the height that 90% are below, the 90th percentile. Later chapters use this curve to compare data with a distribution.
In symbols, with for the density and for the cumulative curve:
The integral sign is a stretched S for “sum”: it adds up the area under between two points. stands for the value of a random case before we see it, a random variable.
A distribution is a model of the whole population: an endless pile of values to draw from. Data are a sample, a handful drawn from that pile. A small sample shows the shape only roughly; the larger the sample, the closer its histogram comes to the curve.
SAMPLE_SIZES = (20, 500, 50_000)
def plot_sample(ax: Axes, n: int) -> None:
"""A histogram of n random heights against the distribution they came from."""
sample = rng.normal(MEAN, SD, n)
ax.hist(sample, bins=np.arange(140, 201, 3), density=True, color=style.BLUE, alpha=0.55)
ax.plot(heights, stats.norm.pdf(heights, MEAN, SD), color=style.CORAL, lw=1.6)
ax.set(title=f"{n:,} people", xlabel="height (cm)", ylabel="how common", yticks=[])
fig, axes = style.subplots(3, height=2.8, width=10.5)
for ax, n in zip(axes, SAMPLE_SIZES):
plot_sample(ax, n)
fig.tight_layout()
plt.show()
Everything that follows rests on this picture: the data we have are one sample, and another sample would have looked a little different. How different is what tests and intervals measure.
The normal distribution has two parameters. The mean sets where the bell stands. The standard deviation sets how wide it is: roughly the typical distance of a value from the mean. With those two numbers, every share is known:
SETOSA = iris.loc[iris.target == 0, "sepal length (cm)"]
VIRGINICA = iris.loc[iris.target == 2, "sepal length (cm)"]
def plot_bands(ax: Axes) -> None:
"""The shares within one, two and three standard deviations of the mean."""
z = np.linspace(-4, 4, 400)
curve = stats.norm.pdf(z)
ax.plot(z, curve, color=style.BLUE)
for k, height, share in ((3, 0.03, "99.7%"), (2, 0.11, "95%"), (1, 0.19, "68%")):
inside = abs(z) <= k
ax.fill_between(z[inside], curve[inside], color=style.BLUE, alpha=0.12, lw=0)
ax.plot([-k, k], [height, height], color=style.INK, lw=1, marker="|", ms=8)
ax.text(0, height + 0.008, share, ha="center", fontsize=9)
ax.set_xticks(range(-3, 4), ["−3", "−2", "−1", "0", "+1", "+2", "+3"])
ax.set(
title="68, 95, 99.7",
xlabel="standard deviations from the mean",
ylabel="how common",
yticks=[],
)
def plot_two_species(ax: Axes) -> None:
"""Two species on the centimetre scale: different centres, different spreads."""
style.hist(ax, SETOSA, bins=12, color=style.BLUE)
style.hist(ax, VIRGINICA, bins=12, color=style.CORAL)
ax.set(title="Centimetres", xlabel="sepal length (cm)", ylabel="flowers")
def plot_z_scores(ax: Axes) -> None:
"""The same two species as z-scores: (value − mean) / standard deviation."""
for values, color in ((SETOSA, style.BLUE), (VIRGINICA, style.CORAL)):
z = (values - values.mean()) / values.std()
style.hist(ax, z, bins=np.linspace(-3, 3, 13), color=color, density=True)
z = np.linspace(-3.2, 3.2, 300)
ax.plot(z, stats.norm.pdf(z), color=style.CORAL)
ax.set(title="z-scores", xlabel="standard deviations from own mean", ylabel="how common")
ax.set_yticks([])
fig, axes = style.subplots(3, height=3, width=10.5)
for ax, plot in zip(axes, [plot_bands, plot_two_species, plot_z_scores]):
plot(ax)
fig.tight_layout()
plt.show()
A z-score says how many standard deviations a value lies from the mean:
It turns any normal distribution into the same one, the standard normal, with mean 0 and standard deviation 1. Its bell has a formula, which is the same for every normal once values are measured as z-scores:
(“mu”) is the mean and (“sigma”) the standard deviation. The exponent falls with the square of the distance, which is why the tails thin out so fast: z = 2 is 7 times rarer than the middle, z = 4 is 3,000 times rarer. The factor in front makes the total area 1. A setosa flower with z = 2 is as unusual for a setosa as a virginica with z = 2 is for a virginica, although their lengths in centimetres differ. Beyond z = 2 lie only 2.5% of values on each side, which is why “two standard deviations” is the usual line for unusual.
Because the tails thin out so fast, a modest difference in the mean becomes a huge difference far out. Take two populations of men with the same spread, 7 cm, and mean heights of 174 cm and 183 cm, roughly those of Spain and the Netherlands.
GROUPS = ((174, style.BLUE, "mean 174 cm"), (183, style.CORAL, "mean 183 cm"))
def plot_tail(ax: Axes, line: float) -> None:
"""Two bells with the area beyond `line` shaded, and the ratio of the two areas."""
x = np.linspace(line - 2, 215, 300) if line > 195 else np.linspace(145, 215, 400)
shares = []
for mean, color, label in GROUPS:
curve = stats.norm.pdf(x, mean, 7)
ax.plot(x, curve, color=color, lw=1.8, label=label)
beyond = x >= line
ax.fill_between(x[beyond], curve[beyond], color=color, alpha=0.3, lw=0)
shares.append(stats.norm.sf(line, mean, 7))
ax.axvline(line, color=style.INK, ls="--", lw=1)
style.note(
ax,
f"{line} cm or taller:\n{shares[0]:.2%} against {shares[1]:.2%}\n"
f"{shares[1] / shares[0]:.0f} times as many",
x=0.03 if line < 195 else 0.4,
y=0.95,
color=style.INK,
)
ax.legend(fontsize=8, loc="center right")
ax.set(xlabel="height (cm)", ylabel="how common", yticks=[], title=f"Beyond {line} cm")
fig, axes = style.subplots(2, height=3.1, width=10)
plot_tail(axes[0], 190)
plot_tail(axes[1], 200)
fig.tight_layout()
plt.show()
The same arithmetic explains why small differences between groups in an average can mean large differences among the extremes, in either direction, and why such claims are very sensitive to the assumed shape of the tails.
The bell needs many small influences. If one influence is large compared with the rest, the distribution becomes a mix of two bells. Sex is such an influence on height: men are on average about 13 cm taller than women.
def plot_mixture(ax: Axes, gap: float) -> None:
"""Two equal groups `gap` cm apart: the groups as lines, the total as bars."""
group = rng.integers(0, 2, 50_000)
values = 163 + gap * group + rng.normal(0, 7, 50_000)
bins = np.arange(120, 220, 2)
ax.hist(values, bins=bins, color=style.GREY, alpha=0.45, label="everyone")
for g, color, label in ((0, style.BLUE, "group 1"), (1, style.CORAL, "group 2")):
counts, edges = np.histogram(values[group == g], bins=bins)
ax.step(edges[:-1], counts, where="post", color=color, lw=1.5, label=label)
ax.legend(fontsize=8)
ax.set(title=f"Groups {gap} cm apart", xlabel="height (cm)", ylabel="people")
fig, axes = style.subplots(2, height=3.1, width=10)
plot_mixture(axes[0], 13)
plot_mixture(axes[1], 25)
fig.tight_layout()
plt.show()
A histogram with two humps, or an unusually flat top, is a hint that the data mix groups that should be looked at separately.
A measurement with a limit cannot follow the bell up to the limit. Say reading ability is normal, and a test has 10 questions, each answered correctly with a chance equal to the student’s ability. If the test fits the students, the scores form a bell. If it is too easy, many students reach the maximum and pile up there: a ceiling effect. A test that is too hard gives a floor effect, a pile at zero.
CENTRES = (
(0.5, style.BLUE, "A test that fits"),
(0.9, style.CORAL, "Too easy: ceiling"),
(0.1, style.GREEN, "Too hard: floor"),
)
def plot_scores(ax: Axes, centre: float, color: str, title: str) -> None:
"""Test scores out of 10 for students whose ability is normal around `centre`."""
ability = np.clip(rng.normal(centre, 0.15, 10_000), 0, 1)
scores = rng.binomial(10, ability)
ax.bar(range(11), np.bincount(scores, minlength=11), color=color, alpha=0.7)
ax.set(title=title, xlabel="score out of 10", ylabel="students", xticks=range(0, 11, 2))
fig, axes = style.subplots(3, height=2.9, width=10.5)
for ax, (centre, color, title) in zip(axes, CENTRES):
plot_scores(ax, centre, color, title)
fig.tight_layout()
plt.show()
A pile at a limit hides real differences: every student at 10 out of 10 gets the same score, however able. Unlike the long tail of the log-normal in Section 2.1.4, no transformation undoes it. Averages and tests also suffer: two groups that both sit at the ceiling look alike even if they are not. Censored-data methods (Chapter 9.4) and models for bounded outcomes treat such data properly.
Chapter 1.2 computed the mean and variance of data. A distribution has them too. The expected value is the balance point of the distribution, every value weighted by how common it is:
and the variance is the expected squared distance from it, . The sample mean estimates ; the more data, the closer it gets.
A few rules do a lot of work later in the book. Shifting values does not change their spread, and stretching them stretches the standard deviation by the same factor:
Adding two random quantities adds their means. Their variances add too, plus a term for how they move together:
For independent quantities the covariance is 0, so variances add, standard deviations do not: they combine like the two sides of a right-angled triangle.
DRAWS = 200_000
waits = rng.exponential(1, DRAWS) # standard deviation 1
errors = rng.uniform(-2, 2, DRAWS) # standard deviation 4 / √12 ≈ 1.15
total = waits + errors
def plot_sum(ax: Axes) -> None:
"""The two parts and their sum, as outlines of histograms."""
edges = np.linspace(-3, 8, 80)
for values, color, label in (
(waits, style.GOLD, "waiting time"),
(errors, style.GREEN, "error"),
(total, style.BLUE, "their sum"),
):
ax.hist(
values,
bins=edges,
histtype="step",
lw=1.8,
color=color,
density=True,
label=f"{label}: sd {values.std():.2f}",
)
ax.legend(fontsize=7.5, loc="upper right")
ax.set(title="A sum", xlabel="value", ylabel="how common", yticks=[])
def plot_triangle(ax: Axes) -> None:
"""Standard deviations of independent parts combine like the sides of a right triangle."""
a, b = waits.std(), errors.std()
ax.plot([0, a, a, 0], [0, 0, b, 0], color=style.INK, lw=0)
ax.plot([0, a], [0, 0], color=style.GOLD, lw=4)
ax.plot([a, a], [0, b], color=style.GREEN, lw=4)
ax.plot([0, a], [0, b], color=style.BLUE, lw=4)
ax.add_patch(plt.Rectangle((a - 0.12, 0), 0.12, 0.12, fill=False, ec=style.INK))
ax.text(a / 2, -0.12, f"{a:.2f}", ha="center", va="top", color=style.GOLD)
ax.text(a + 0.06, b / 2, f"{b:.2f}", va="center", color=style.GREEN)
ax.text(a / 2 - 0.12, b / 2 + 0.1, f"{np.hypot(a, b):.2f}", ha="right", color=style.BLUE)
ax.set(title="Sides of a triangle", aspect="equal", xlim=(-0.3, 1.6), ylim=(-0.4, 1.4))
ax.set_axis_off()
def plot_averages(ax: Axes) -> None:
"""The spread of averages of n values: σ/√n."""
ns = np.array([1, 2, 4, 8, 16, 32, 64])
spreads = [rng.exponential(1, (20_000, n)).mean(axis=1).std() for n in ns]
grid = np.linspace(1, 64, 200)
ax.plot(grid, 1 / np.sqrt(grid), color=style.CORAL)
ax.plot(ns, spreads, "o", ms=6, color=style.BLUE)
ax.set(title="Averages", xlabel="values averaged (n)", ylabel="sd of the average")
fig, axes = style.subplots(3, height=3, width=10.5)
plot_sum(axes[0])
plot_triangle(axes[1])
plot_averages(axes[2])
fig.tight_layout()
plt.show()
The last panel is the reason behind the standard error of Chapter 3.2. An average of n independent values with standard deviation σ has variance , so
Four times as much data halves the uncertainty; a hundred times as much divides it by ten.
In Chapter 1.1, many small effects that add up made a bell. When the effects multiply, the result has a long right tail. A business that grows by 5% a year grows by more in money each year, because 5% of a larger number is larger. Incomes, house prices, city sizes and the time to finish a task all work like this.
growth = rng.uniform(0.94, 1.07, (10_000, 100)) # each year: shrink up to 6%, or grow up to 7%
sizes = 100 * growth.prod(axis=1)
def plot_sizes(ax: Axes) -> None:
"""Sizes after 100 years: most are modest, a few are huge."""
style.hist(ax, sizes[sizes < np.quantile(sizes, 0.99)], bins=50, color=style.GOLD)
ax.set(title="Effects that multiply", xlabel="size after 100 years", ylabel="businesses")
def plot_log_scale(
ax: Axes, values: np.ndarray, ticks: list[int], color: str, title: str, xlabel: str
) -> None:
"""A histogram with log-spaced bins on a log axis: a log-normal shows up as a bell."""
style.hist(ax, values, bins=np.geomspace(values.min(), values.max(), 30), color=color)
ax.set_xscale("log")
ax.set_xticks(ticks, [f"{t:,}" for t in ticks])
ax.minorticks_off()
ax.set(title=title, xlabel=xlabel)
fig, axes = style.subplots(3, height=3, width=10.5)
plot_sizes(axes[0])
plot_log_scale(
axes[1], sizes, [50, 100, 200, 400], style.GOLD, "The same, log scale", "size (log scale)"
)
plot_log_scale(
axes[2],
engel.income.to_numpy(),
[500, 1000, 2000, 4000],
style.BLUE,
"Incomes",
"income (log scale)",
)
axes[1].set_ylabel("businesses")
axes[2].set_ylabel("households")
fig.tight_layout()
plt.show()
A log scale turns multiplying into adding: every step of the same size on the axis is the same factor. On that scale, many multiplied effects become many added ones, and they make a bell. That is the definition of a log-normal distribution: its logarithm is normal. In practice it means that a long right tail is often tamed by taking logs, a trick used again and again in regression.
Customers arrive at a shop at random moments, about one every two minutes. How long until the next one? Most gaps are short, a few are long. That shape is the exponential distribution, and it has one parameter: the average wait.
AVERAGE_GAP = 2.0 # minutes between arrivals, on average
arrivals = np.sort(rng.uniform(0, 60_000, 30_000)) # 1,000 hours of random arrivals
gaps = np.diff(arrivals)
def plot_timeline(ax: Axes) -> None:
"""One hour of arrivals: each line is a customer coming in."""
first_hour = arrivals[arrivals < 60]
ax.vlines(first_hour, 0, 1, color=style.BLUE, lw=1.2)
ax.set(title="One hour", xlabel="minutes", xlim=(0, 60), ylim=(0, 1.4))
style.clean(ax, y=False)
def plot_gaps(ax: Axes) -> None:
"""Gaps between arrivals: an exponential distribution."""
shown = gaps[gaps < 12]
width = style.hist(ax, shown, bins=48)
t = np.linspace(0, 12, 300)
ax.plot(t, stats.expon.pdf(t, scale=AVERAGE_GAP) * len(gaps) * width, color=style.CORAL)
ax.set(title="Gaps: exponential", xlabel="minutes to the next arrival", ylabel="gaps")
def plot_third(ax: Axes, k: int = 3) -> None:
"""The wait for the k-th arrival: k exponential gaps added up, a gamma distribution."""
waits = gaps[: len(gaps) // k * k].reshape(-1, k).sum(axis=1)
width = style.hist(ax, waits[waits < 20], bins=40, color=style.GREEN)
t = np.linspace(0, 20, 300)
curve = stats.gamma.pdf(t, k, scale=AVERAGE_GAP) * len(waits) * width
ax.plot(t, curve, color=style.CORAL)
ax.set(title="Third arrival: gamma", xlabel="minutes to the third arrival", ylabel="waits")
fig, axes = style.subplots(3, height=3, width=10.5)
for ax, plot in zip(axes, [plot_timeline, plot_gaps, plot_third]):
plot(ax)
fig.tight_layout()
plt.show()
The exponential has an odd property: it has no memory. If you have already waited five minutes, the wait from now on is still two minutes on average. Random arrivals do not know how long it has been. That makes it a good model for things that happen at random (calls, failures of electronics, radioactive decay) and a bad one for things that wear out.
Add up several exponential waits and you get a gamma distribution: the wait for the third customer is the sum of three gaps. It has a hump, because three short gaps in a row are rare. The more gaps added, the more it looks like a bell: the central limit theorem again.
Some values cannot go below or above a limit. The simplest case is the uniform distribution: every value in a range equally likely, like the second hand of a clock at a random glance. It looks dull, but it will turn up in testing: when nothing is going on, p-values are uniform (Chapter 3.1).
Shares, rates and proportions lie between 0 and 1. The beta distribution is made for them. Its two parameters pull the hump towards 0 or towards 1, and how strongly decides how narrow it is.
food_share = (engel.foodexp / engel.income).to_numpy()
SHAPES = [(1, 1), (2, 5), (5, 2), (12, 12)] # (α, β) for each curve
def plot_shapes(ax: Axes) -> None:
"""Beta curves for a few pairs of parameters."""
p = np.linspace(0.001, 0.999, 400)
for (a, b), color in zip(SHAPES, style.PALETTE):
ax.plot(p, stats.beta.pdf(p, a, b), color=color, label=f"α={a}, β={b}")
ax.xaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8)
ax.set(title="Beta shapes", xlabel="share", ylabel="how common", yticks=[])
def plot_food_share(ax: Axes) -> None:
"""Each household's food share, with the beta distribution that fits it best."""
width = style.hist(ax, food_share, bins=25, color=style.GREY)
a, b, *_ = stats.beta.fit(food_share, floc=0, fscale=1)
p = np.linspace(0.25, 1, 300)
ax.plot(p, stats.beta.pdf(p, a, b) * len(food_share) * width, color=style.CORAL)
style.note(ax, f"α = {a:.0f}, β = {b:.0f}", x=0.03)
ax.xaxis.set_major_formatter(style.PERCENT)
ax.set(title="Food share", xlabel="share of income spent on food", ylabel="households")
fig, axes = style.subplots(2, height=3.1, width=9)
plot_shapes(axes[0])
plot_food_share(axes[1])
fig.tight_layout()
plt.show()
With α = β = 1 the beta is flat: the uniform is a special case. When α is larger than β the hump moves towards 1, and the larger both are, the narrower it gets. A normal curve would put some probability below 0% or above 100%; the beta never does.
Two more distributions do not describe data at all. They describe numbers we compute from data, and they are the reason statistical tests work.
The t distribution is what an average looks like when we measure it in standard deviations that we had to estimate from the same small sample. The estimate is sometimes too small, which makes some values look extreme, so the t has fatter tails than the normal. With more data the estimate gets better and the t turns into the normal.
The chi-square distribution is what a sum of squared z-scores looks like. Squares are never negative, so it starts at 0 with a long right tail. It is the yardstick for tests on tables (Chapter 3.5).
DRAWS = 50_000
def plot_t_family(ax: Axes) -> None:
"""t distributions for small and larger samples, against the normal."""
z = np.linspace(-5, 5, 400)
ax.plot(z, stats.norm.pdf(z), color=style.GREY, ls="--", label="normal")
for df, color in ((1, style.CORAL), (4, style.GOLD), (30, style.BLUE)):
ax.plot(z, stats.t.pdf(z, df), color=color, lw=1.6, label=f"t, sample of {df + 1}")
ax.legend(fontsize=8)
ax.set(title="The t family", xlabel="standard errors from the mean", ylabel="how common")
ax.set_yticks([])
def plot_t_simulated(ax: Axes, n: int = 4) -> None:
"""Averages of n normal values, divided by their own estimated standard error."""
samples = rng.normal(0, 1, (DRAWS, n))
t = samples.mean(axis=1) / (samples.std(axis=1, ddof=1) / np.sqrt(n))
edges = np.linspace(-6, 6, 61)
style.hist(ax, t, bins=edges, density=True)
z = np.linspace(-6, 6, 400)
ax.plot(z, stats.norm.pdf(z), color=style.GREY, ls="--")
ax.plot(z, stats.t.pdf(z, n - 1), color=style.CORAL)
ax.set(title=f"Samples of {n}", xlabel="standard errors from the mean", ylabel="how common")
ax.set_yticks([])
def plot_chi_square(ax: Axes, k: int = 3) -> None:
"""Sums of k squared z-scores: a chi-square distribution."""
sums = (rng.normal(0, 1, (DRAWS, k)) ** 2).sum(axis=1)
style.hist(ax, sums[sums < 15], bins=50, color=style.GREEN, density=True)
x = np.linspace(0.01, 15, 400)
ax.plot(x, stats.chi2.pdf(x, k), color=style.CORAL)
ax.set(title="Chi-square", xlabel=f"sum of {k} squared z-scores", ylabel="how common")
ax.set_yticks([])
fig, axes = style.subplots(3, height=3, width=10.5)
for ax, plot in zip(axes, [plot_t_family, plot_t_simulated, plot_chi_square]):
plot(ax)
fig.tight_layout()
plt.show()
Both have one parameter, the degrees of freedom: roughly, how many independent pieces of information went into the number. A sample of 4 leaves 3 for the t, because one was used up to estimate the mean.
| Mechanism | Distribution | Mean | Variance | Examples |
|---|---|---|---|---|
| many small effects, added | normal (μ, σ) | μ | σ² | heights, measurement errors, averages |
| many small effects, multiplied | log-normal (μ, σ of the log) | incomes, prices, city sizes | ||
| waiting for one random event | exponential (rate λ) | 1/λ | 1/λ² | time to the next call or failure |
| waiting for k random events | gamma (k, scale θ) | kθ | kθ² | time to the third customer, rainfall |
| everything in a range equally likely | uniform (a, b) | (a + b)/2 | (b − a)²/12 | rounding errors, p-values with no effect |
| a share between 0 and 1 | beta (α, β) | α/(α + β) | conversion rates, food share of income | |
| an average in estimated standard errors | t (ν degrees of freedom) | 0 | ν/(ν − 2) | t-tests, confidence intervals |
| a sum of k squared z-scores | chi-square (k) | k | 2k | tests on tables, variances |
Next: distributions for counts and yes/no outcomes.