2.3 Fitting and checking · Visual Statistics

2.3.1 Maximum likelihood

A distribution has parameters, and the data have to choose them. The usual way is maximum likelihood: try parameter values, and for each ask how probable the data we actually saw would be. The likelihood of a curve is the product of its heights at every data point. A curve that puts its height where the data are wins.

Code

SD = setosa.std()
CANDIDATES = [(4.85, style.GOLD), (setosa.mean(), style.BLUE), (5.2, style.CORAL)]


def log_likelihood(mean: float) -> float:
    """How probable the 50 lengths are under a normal with this mean (on a log scale)."""
    return float(stats.norm.logpdf(setosa, mean, SD).sum())


def plot_candidates(ax: Axes) -> None:
    """The data and three candidate curves, each with its log-likelihood."""
    x = np.linspace(4, 6, 300)
    for mean, color in CANDIDATES:
        label = f"mean {mean:.2f}: {log_likelihood(mean):.0f}"
        ax.plot(x, stats.norm.pdf(x, mean, SD), color=color, lw=1.8, label=label)
    ax.plot(setosa, np.zeros_like(setosa), "|", ms=12, color=style.INK)
    ax.legend(title="log-likelihood", fontsize=8, title_fontsize=8, loc="upper right")
    ax.set(title="Three candidates", xlabel="sepal length (cm)", ylabel="how common")
    ax.set_yticks([])


def plot_curve(ax: Axes) -> None:
    """The log-likelihood across all means: the best one is the top of the hill."""
    means = np.linspace(4.7, 5.3, 200)
    ax.plot(means, [log_likelihood(m) for m in means], color=style.INK, lw=1.5)
    for mean, color in CANDIDATES:
        ax.plot(mean, log_likelihood(mean), "o", ms=7, color=color)
    ax.axvline(setosa.mean(), color=style.BLUE, ls="--", lw=1)
    ax.set(title="Likelihood hill", xlabel="mean of the normal (cm)", ylabel="log-likelihood")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_candidates(axes[0])
plot_curve(axes[1])
fig.tight_layout()
plt.show()

Figure 2.3.1: Sepal lengths of the 50 setosa flowers (ticks along the bottom). Left: three normal curves with different means; the middle one puts most of its height where the data are. Right: the log-likelihood for every possible mean peaks at the average of the data.

Multiplying 50 small heights gives a tiny number, so we add their logarithms instead: the log-likelihood. Its peak is the best fit. For a normal distribution the peak lands exactly on the average of the data, and the best standard deviation is the spread of the data. For other distributions there is no such shortcut, and a computer climbs the hill. Almost every model in this book, from the t-test to logistic regression, is fitted this way or close to it.

In symbols, for data $x_1, \dots, x_n$ and a density $f$ with parameters θ (“theta”):

$L(\theta) = \prod_{i=1}^{n} f(x_i; \theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i; \theta)$

and the estimate $\hat\theta$ is the θ where $\ell$ is highest. An older and simpler method, the method of moments, sets the distribution’s mean and variance equal to the data’s and solves for the parameters. It is quick, but usually less precise; Chapter 2.2 used it at first for the doctor visits, and the maximum likelihood fit followed the data better.

2.3.2 How sure is the fit?

The likelihood also tells us how precisely the data pin down a parameter. With few data, many values are nearly as likely as the best one and the hill is flat. With more data the hill gets steep, and only values close to the top are believable.

Code

rng = np.random.default_rng(2)
RATES = np.linspace(0.2, 3.6, 600)


def relative_log_likelihood(waits: np.ndarray) -> np.ndarray:
    """Log-likelihood of each rate for exponential data, minus its maximum."""
    ll = len(waits) * np.log(RATES) - RATES * waits.sum()
    return ll - ll.max()


def plot_hills(ax: Axes) -> None:
    """One hill per sample size, with the cut-off for a 95% interval."""
    for n, color in ((10, style.GOLD), (50, style.GREEN), (250, style.BLUE)):
        waits = rng.exponential(1, n)
        ll = relative_log_likelihood(waits)
        inside = RATES[ll >= -1.92]
        ax.plot(
            RATES,
            ll,
            color=color,
            label=f"{n} waits: rate {inside.min():.2f} to {inside.max():.2f}",
        )
    ax.axhline(-1.92, color=style.INK, ls="--", lw=1)
    ax.legend(fontsize=8, loc="lower right")
    ax.set(xlabel="rate of the exponential", ylabel="log-likelihood (top at 0)", ylim=(-8, 0.5))


fig, ax = plt.subplots(figsize=(8, 3.4))
plot_hills(ax)
plt.show()

Figure 2.3.2: Log-likelihood of the rate of an exponential waiting time (true rate 1), for samples of 10, 50 and 250 waits, each shifted so its top is at 0. The more data, the narrower the hill. Rates above the dashed line (1.92 below the top) form an approximate 95% interval.

The steepness at the top, the curvature, gives the standard error: for large samples $\text{SE}(\hat\theta) \approx 1 / \sqrt{-\ell''(\hat\theta)}$ , where $\ell''$ is the curvature. That is how statistics programs compute the standard errors of maximum likelihood estimates, including every regression coefficient in Part 4. Cutting the hill 1.92 below its top (half the 95% point of a chi-square with one degree of freedom, 3.84) gives a likelihood interval, which also works when the hill is lopsided.

2.3.3 Comparing nested distributions: the likelihood ratio test

Sometimes one candidate is a special case of another. The exponential is a gamma with shape 1, so the gamma can never fit worse. Does it fit enough better to be worth its extra parameter? The likelihood ratio test compares the two log-likelihoods:

$\Lambda = 2\,\big(\ell_{\text{larger}} - \ell_{\text{smaller}}\big)$

If the smaller model is true, Λ follows a chi-square distribution with as many degrees of freedom as the extra parameters.

Code

def lr_statistic(values: np.ndarray) -> float:
    """2 × (log-likelihood of the gamma fit − log-likelihood of the exponential fit)."""
    gamma_ll = stats.gamma.logpdf(values, *stats.gamma.fit(values, floc=0)).sum()
    expon_ll = stats.expon.logpdf(values, *stats.expon.fit(values, floc=0)).sum()
    return 2 * (gamma_ll - expon_ll)


simulated = np.array([lr_statistic(rng.exponential(1, 100)) for _ in range(2_000)])


def plot_null(ax: Axes) -> None:
    """Simulated Λ under the smaller model, against the chi-square curve."""
    style.hist(ax, simulated[simulated < 10], bins=40, density=True)
    x = np.linspace(0.05, 10, 300)
    ax.plot(x, stats.chi2.pdf(x, 1), color=style.CORAL)
    ax.axvline(3.84, color=style.INK, ls="--", lw=1)
    style.note(ax, f"beyond 3.84: {(simulated > 3.84).mean():.1%}", x=0.45, color=style.INK)
    ax.set(
        title="If the exponential is true",
        xlabel="Λ",
        ylabel="how common",
        yticks=[],
        ylim=(0, 1.2),
    )


def plot_income(ax: Axes) -> None:
    """Λ for the incomes, on the same kind of scale."""
    value = lr_statistic(income)
    x = np.linspace(0.05, 12, 300)
    ax.plot(x, stats.chi2.pdf(x, 1), color=style.CORAL)
    ax.annotate(
        f"incomes: Λ = {value:.0f} →", xy=(11.8, 0.3), ha="right", fontsize=10, color=style.BLUE
    )
    ax.set(title="Household incomes", xlabel="Λ", yticks=[], ylim=(0, 1.2))


fig, axes = style.subplots(2, height=3, width=9.5)
plot_null(axes[0])
plot_income(axes[1])
fig.tight_layout()
plt.show()

Figure 2.3.3: Left: Λ for 2,000 simulated samples of 100 truly exponential waits, comparing a gamma fit with an exponential fit; it follows a chi-square with one degree of freedom (red), and 5% of samples land beyond 3.84 (dashed). Right: for the 235 household incomes Λ is enormous, far beyond anything chance gives: the gamma’s extra shape parameter is clearly needed.

The same test compares nested regression models (Chapter 4.4) and survival models (Chapter 9.5).

2.3.4 A smooth shape without a name

Sometimes no named distribution fits, and we just want a smooth picture of the shape. A kernel density estimate puts a small bell on every data point and adds them up. The width of those bells plays the role the bin width plays in a histogram.

Code

petal = iris["petal length (cm)"].to_numpy()
GRID = np.linspace(0, 8, 400)


def plot_how(ax: Axes, shown: int = 8, width: float = 0.35) -> None:
    """A few points, each with its own small bell; the estimate is their sum."""
    points = np.sort(np.random.default_rng(4).choice(petal, shown, replace=False))
    bells = stats.norm.pdf(GRID[:, None], points, width) / shown
    ax.plot(GRID, bells, color=style.GREY, lw=0.9, ls="--")
    ax.plot(GRID, bells.sum(axis=1), color=style.BLUE)
    ax.plot(points, np.zeros(shown), "|", ms=12, color=style.INK)
    ax.set(title="Eight flowers", xlabel="petal length (cm)", ylabel="how common", yticks=[])


def plot_width(ax: Axes, factor: float, title: str) -> None:
    """The estimate for all flowers, with bells of a given width."""
    style.hist(ax, petal, bins=30, color=style.GREY, density=True)
    ax.plot(GRID, stats.gaussian_kde(petal, bw_method=factor)(GRID), color=style.BLUE)
    ax.set(title=title, xlabel="petal length (cm)", yticks=[])


fig, axes = style.subplots(4, height=2.8, width=11)
plot_how(axes[0])
plot_width(axes[1], 0.04, "Narrow")
plot_width(axes[2], 0.15, "Sensible")
plot_width(axes[3], 0.6, "Wide")
fig.tight_layout()
plt.show()

Figure 2.3.4: Petal lengths of the 150 iris flowers. Left: how it works, with eight flowers, each carrying a small bell (dashed); their sum is the estimate. Then the full data with narrow bells (wiggly), a sensible width, and wide bells (the two humps melt into one).

The estimate shows the shape and no more: it has no parameters to interpret and says nothing about values beyond the data. It is a picture, not a model, but often the best first picture.

2.3.5 The QQ plot

A histogram with a curve over it hides the tails, where the bars are too short to see. The quantile-quantile plot shows them. Sort the data, and set each value against the value the fitted distribution says should be at that place in the order: the smallest against the distribution’s 1st percentile, the middle against its median, and so on. If the distribution is right, the points lie on the diagonal.

Code

def plot_qq(ax: Axes, values: np.ndarray, title: str, unit: str, color: str) -> None:
    """Each sorted value against where a fitted normal says it should be."""
    values = np.sort(values)
    places = (np.arange(1, len(values) + 1) - 0.5) / len(values)
    expected = stats.norm.ppf(places, values.mean(), values.std())
    ends = [expected.min(), expected.max()]
    ax.plot(ends, ends, color=style.INK, lw=1, ls="--")
    ax.plot(expected, values, "o", ms=3.5, color=color, alpha=0.7)
    ax.set(title=title, xlabel=f"expected if normal ({unit})", ylabel=f"data ({unit})")


fig, axes = style.subplots(3, height=3.1, width=10.5)
plot_qq(axes[0], setosa, "Sepal length: normal", "cm", style.BLUE)
plot_qq(axes[1], income, "Income: not normal", "income", style.CORAL)
plot_qq(axes[2], np.log(income), "Log income: close", "log income", style.GREEN)
fig.tight_layout()
plt.show()

Figure 2.3.5: QQ plots against a fitted normal. Left: setosa sepal lengths lie on the line (the steps come from measuring to 0.1 cm). Middle: incomes bend away at both ends: the top tail is far longer than a normal allows. Right: the logarithm of income is much closer to the line.

How to read the bends: on the right, points above the line mean the top values are larger than the distribution expects, a long right tail. On the left, points above the line mean the bottom values are larger than expected, a short left tail. Together, as for income, the points curve upwards like a smile: right skew. Left skew bends the other way. An S-shape that crosses the line means both tails are longer than normal, or both shorter.

2.3.6 Which candidate fits best?

When several distributions are possible, compare them on the cumulative scale from Chapter 2.1. The data have their own cumulative curve, a staircase that rises by one step per value. A good distribution follows the staircase everywhere, and the largest gap is a measure of how badly it misses.

Code

CANDIDATE_FITS = {
    "normal": (stats.norm, {}),
    "exponential": (stats.expon, {"floc": 0}),
    "gamma": (stats.gamma, {"floc": 0}),
    "log-normal": (stats.lognorm, {"floc": 0}),
}


def aic(dist: stats.rv_continuous, fixed: dict) -> float:
    """AIC = 2 × free parameters − 2 × log-likelihood (lower is better)."""
    params = dist.fit(income, **fixed)
    free = len(params) - len(fixed)
    return 2 * free - 2 * dist.logpdf(income, *params).sum()


def plot_staircase(ax: Axes) -> None:
    """The data's cumulative staircase against two fitted distributions, with the largest gaps."""
    x = np.sort(income)
    steps = np.arange(1, len(x) + 1) / len(x)
    ax.step(x, steps, where="post", color=style.GREY, lw=1.5, label="data")
    for name, color in (("normal", style.CORAL), ("log-normal", style.BLUE)):
        dist, fixed = CANDIDATE_FITS[name]
        fitted = dist.cdf(x, *dist.fit(income, **fixed))
        ax.plot(x, fitted, color=color, lw=1.6, label=name)
        gaps = np.maximum(abs(steps - fitted), abs(steps - 1 / len(x) - fitted))
        i = int(np.argmax(gaps))
        ax.annotate(
            "",
            xy=(x[i], fitted[i]),
            xytext=(x[i], steps[i]),
            arrowprops=dict(arrowstyle="<->", color=color, lw=1.5),
        )
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8, loc="lower right")
    ax.set(title="Cumulative", xlabel="household income", ylabel="share below", xlim=(0, 3000))


def plot_aic(ax: Axes) -> None:
    """One bar per candidate distribution: its AIC."""
    names = list(CANDIDATE_FITS)
    scores = [aic(*CANDIDATE_FITS[name]) for name in names]
    best = int(np.argmin(scores))
    colors = [style.BLUE if i == best else style.GREY for i in range(len(names))]
    ax.barh(names, scores, color=colors, alpha=0.7)
    ax.set_xlim(min(scores) - 60, max(scores) + 20)
    ax.invert_yaxis()
    ax.set(title="AIC", xlabel="AIC (lower is better)")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_staircase(axes[0])
plot_aic(axes[1])
fig.tight_layout()
plt.show()

Figure 2.3.6: Left: the incomes’ own cumulative staircase (grey), with a fitted normal (red) and a fitted log-normal (blue); the arrows show the largest gap for each. Right: the AIC of four candidate distributions; lower is better.

The largest gap is the statistic of the Kolmogorov-Smirnov test,

$D = \max_x \big|\hat F_n(x) - F(x)\big|$

where $\hat F_n$ is the data’s staircase and $F$ the candidate’s cumulative curve. The test asks whether a gap that large could be chance. (When the candidate’s parameters were fitted to the same data, the gaps come out too small and the plain test is too lenient; the Lilliefors version corrects that for the normal.) A test built specially for normality, Shapiro-Wilk, comes in Chapter 3.8.

The yardstick for D can be built by simulation, which also shows why fitting matters. Draw many samples of 235 incomes from the fitted log-normal itself, so that the candidate is true by construction; refit the log-normal to each sample, as was done for the real data; and compute D each time:

Code

shape, loc, scale = stats.lognorm.fit(income, floc=0)
fitted = stats.lognorm(shape, loc, scale)
observed_d = stats.kstest(income, fitted.cdf).statistic


def refitted_gap(n: int, local: np.random.Generator) -> float:
    """D for a sample from the fitted log-normal, after refitting to that sample."""
    sample = fitted.rvs(size=n, random_state=local)
    refit = stats.lognorm(*stats.lognorm.fit(sample, floc=0))
    return float(stats.kstest(sample, refit.cdf).statistic)


local = np.random.default_rng(2)
null_d = np.array([refitted_gap(len(income), local) for _ in range(2000)])
fig, ax = plt.subplots(figsize=(8, 3.2))
style.null_panel(
    ax,
    observed_d,
    samples=null_d,
    dist=stats.kstwo(len(income)),
    tail="right",
    xlabel="largest gap D",
)
plt.show()

Figure 2.3.7: The Kolmogorov-Smirnov gap D for a log-normal fitted to the 235 incomes. Grey: D for 2,000 samples drawn from the fitted log-normal and refitted each time, the right yardstick. Blue: the textbook KS distribution, which assumes the candidate was fixed in advance. Because refitting pulls the curve towards each sample, the honest gaps are smaller, and the textbook curve sits too far right: it would call the fit good too easily. For the observed D (dashed), the grey yardstick gives p = 0.06; the textbook curve’s tail (red) would give about 0.43.

The AIC (Akaike’s information criterion) and the BIC (Bayesian information criterion) compare candidates that are not nested, by their log-likelihood with a penalty for every parameter k:

$\text{AIC} = 2k - 2\ell(\hat\theta), \qquad \text{BIC} = k \log n - 2\ell(\hat\theta)$

so a distribution cannot win just by being more flexible. The BIC punishes parameters harder once n is above 7, and prefers simpler models. Only differences matter; a difference of 10 or more is a clear win.

Two warnings. With a lot of data, every test finds some gap, because no real data are exactly normal or exactly anything. With little data, no test finds any gap. The question is never “is this distribution true?” but “is it close enough for what I want to do with it?” The pictures answer that better than any single number.

2.3.7 In one picture

Maximum likelihood picks the parameters that make the observed data most probable.
The curvature of the log-likelihood gives standard errors: more data, steeper hill.
The likelihood ratio test compares nested models: $2(\ell_1 - \ell_0)$ against a chi-square.
A kernel density estimate is a smooth histogram: a small bell on every point.
A QQ plot shows the tails: points on the diagonal mean the distribution fits; bends show where it does not.
AIC and BIC rank candidates, with a penalty for every parameter.
No distribution is exactly true. Choose one that is close where it matters for your question.

Next: several variables at once.