3.12 Bayesian statistics

Bayesian statistics treats an unknown quantity as uncertain and describes that uncertainty with a probability distribution. It starts from a prior, what we believed before the data, and Bayes’ rule turns it into a posterior, what we believe after.

Code: data and style used in this chapter

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from scipy.special import betaln
from sklearn import datasets

from statly import style

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

spector = (
    sm.datasets.spector.load_pandas().data
)  # 32 students: grade point average, grade improved
diabetes = datasets.load_diabetes()  # 442 patients: 10 measurements, disease progression
share = np.linspace(0, 1, 501)  # every possible value of a share


def plot_density(
    ax: Axes, dist: stats.rv_continuous, color: str, label: str, ls: str = "-"
) -> None:
    """A distribution over the share, as a curve."""
    ax.plot(share, dist.pdf(share), color=color, lw=2, ls=ls, label=label)

3.12.1 From prior to posterior

A new treatment helped 7 of the first 10 patients. What is its true success rate θ? The methods of Part 3 so far treat θ as a fixed unknown and ask how the data would vary. The Bayesian approach treats θ itself as uncertain and asks how plausible each value is, given the data. Bayes’ rule (Chapter 2.4) does the work:

$\underbrace{p(\theta \mid \text{data})}_{\text{posterior}} \;\propto\; \underbrace{p(\text{data} \mid \theta)}_{\text{likelihood}} \times \underbrace{p(\theta)}_{\text{prior}}$

The likelihood is the same as in Chapter 2.3: how probable the data are for each θ. The prior says how plausible each θ was before the data. Their product, rescaled to have area 1, is the posterior. For a share, a convenient prior is the beta distribution (Chapter 2.1): with a Beta(a, b) prior and k successes in n tries, the posterior is again a beta,

$\theta \mid \text{data} \sim \text{Beta}(a + k,\; b + n - k)$

so the prior acts like a successes and b failures seen before.

Code

PRIORS = [("flat prior, Beta(1, 1)", 1, 1), ("sceptical prior, Beta(10, 10)", 10, 10)]

fig, axes = plt.subplots(2, 2, figsize=(10, 5.2), sharex=True)
for row, (k, n) in enumerate(((7, 10), (70, 100))):
    likelihood = stats.binom.pmf(k, n, share)
    for col, (name, a, b) in enumerate(PRIORS):
        ax = axes[row, col]
        posterior = stats.beta(a + k, b + n - k)
        plot_density(ax, stats.beta(a, b), style.BLUE, "prior", ls="--")
        ax.plot(
            share,
            likelihood / likelihood.max() * posterior.pdf(share).max(),
            color=style.GREY,
            lw=1.5,
            label="likelihood",
        )
        plot_density(ax, posterior, style.CORAL, "posterior")
        low, high = posterior.interval(0.95)
        style.note(
            ax,
            f"{k} of {n}\nposterior mean {posterior.mean():.0%}\n95%: {low:.0%} to {high:.0%}",
            x=0.03,
            y=0.95,
            color=style.INK,
        )
        ax.set(title=name if row == 0 else "", yticks=[], ylabel="how plausible")
        ax.xaxis.set_major_formatter(style.PERCENT)
axes[0, 1].legend(fontsize=8, loc="upper right")
for ax in axes[1]:
    ax.set_xlabel("success rate θ")
fig.tight_layout()
plt.show()

Figure 3.12.1: Prior (dashed), likelihood (grey, rescaled) and posterior (solid) for a success rate. Top: 7 successes in 10 patients. With a flat prior (left), the posterior is the likelihood; with a sceptical prior centred on 50% (right, worth 20 earlier patients), the posterior is pulled towards 50%. Bottom: 70 successes in 100 patients. Now the data dominate, and both priors give nearly the same posterior.

Two lessons. With little data, the prior matters, so it must be stated and defended. With much data, reasonable priors agree, and the posterior is close to the likelihood.

3.12.2 Learning as the data arrive

Yesterday’s posterior is today’s prior. Updating patient by patient gives the same posterior as updating on all of them at once, and the curve narrows as evidence accumulates:

Code

outcomes = rng.uniform(size=500) < 0.65
STAGES = (0, 1, 5, 20, 100, 500)

fig, ax = plt.subplots(figsize=(8, 3.3))
for n, color in zip(STAGES, plt.cm.viridis(np.linspace(0, 0.9, len(STAGES)))):
    k = outcomes[:n].sum()
    density = stats.beta(1 + k, 1 + n - k).pdf(share)
    ax.plot(share, density, color=color, lw=2, label=f"after {n}")
ax.axvline(0.65, color=style.INK, ls="--", lw=1)
ax.xaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8, ncol=2)
ax.set(xlabel="success rate θ", ylabel="how plausible", yticks=[])
plt.show()

Figure 3.12.2: Posterior for a success rate whose true value is 65% (dashed line), starting from a flat prior and updated after 0, 1, 5, 20, 100 and 500 simulated patients. Each posterior becomes the prior for the next patients; the curve moves towards the truth and narrows.

3.12.3 Credible intervals and posterior probabilities

The posterior answers questions directly. A 95% credible interval is a range that holds θ with probability 95%, given the data and the prior. That is the statement people often wrongly make about a confidence interval (Chapter 3.2). The posterior also gives the probability of any claim, such as “the treatment works better than the old one, which helped 50%”:

Code

posterior = stats.beta(8, 4)
low, high = posterior.interval(0.95)
inside = (share >= low) & (share <= high)
above = share > 0.5

fig, axes = style.subplots(2, height=3, width=10)
for ax, mask, color, text in (
    (axes[0], inside, style.CORAL, f"95% credible interval\n{low:.0%} to {high:.0%}"),
    (axes[1], above, style.GREEN, f"P(θ > 50% | data) = {posterior.sf(0.5):.0%}"),
):
    plot_density(ax, posterior, style.INK, "")
    ax.fill_between(share[mask], posterior.pdf(share[mask]), color=color, alpha=0.35)
    style.note(ax, text, x=0.03, y=0.95, color=style.INK)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.set(xlabel="success rate θ", ylabel="how plausible", yticks=[])
axes[1].axvline(0.5, color=style.INK, ls="--", lw=1)
fig.tight_layout()
plt.show()

Figure 3.12.3: The flat-prior posterior after 7 of 10 successes. Left: the 95% credible interval (shaded) runs from 39% to 89%. Right: the posterior probability that the success rate is above 50%, the area to the right of the line, is 89%. The frequentist answers for the same data are a 95% confidence interval of 40% to 89% (Wilson) and a one-sided p-value of 0.17.

With a flat prior and plenty of data, credible and confidence intervals nearly coincide; the difference is mostly in what they mean.

3.12.4 Shrinkage: the prior pulls towards what is typical

For a mean of normal data with a normal prior, the posterior mean is a weighted average of the prior mean and the sample mean, each weighted by its precision (one over its variance):

$\text{posterior mean} = \frac{\mu_0 / \tau^2 + n \bar{x} / \sigma^2}{1 / \tau^2 + n / \sigma^2}$

where $\mu_0$ and $\tau$ are the prior’s mean and standard deviation and σ is the spread of single values. Small samples are pulled more towards the prior, large ones less. This is the same shrinkage as in mixed models (Chapter 4.7) and ridge regression (Chapter 6.4).

Code

local = np.random.default_rng(3)
operations = np.array([8, 12, 15, 20, 30, 40, 60, 80, 100, 150, 200, 300])
true_rates = local.normal(0.75, 0.06, 12)
successes = local.binomial(operations, true_rates)
raw = successes / operations
prior_mean, prior_sd = np.average(raw, weights=operations), 0.06
noise_var = raw.mean() * (1 - raw.mean()) / operations
shrunk = (prior_mean / prior_sd**2 + raw / noise_var) / (1 / prior_sd**2 + 1 / noise_var)

fig, ax = plt.subplots(figsize=(7, 3.6))
for r, s, n in zip(raw, shrunk, operations):
    ax.plot([0, 1], [r, s], color=style.BLUE, lw=0.5 + 3 * n / operations.max(), alpha=0.8)
    ax.text(-0.04, r, f"{n}", ha="right", va="center", fontsize=7.5)
ax.axhline(prior_mean, color=style.INK, ls="--", lw=1)
ax.set_xticks([0, 1], ["raw rate", "posterior mean"])
ax.yaxis.set_major_formatter(style.PERCENT)
ax.set(ylabel="success rate", xlim=(-0.2, 1.1))
plt.show()

Figure 3.12.4: Twelve hospitals’ success rates for an operation (simulated), each from a different number of operations. Raw rates (left end of each line) are pulled towards the overall average (dashed) by a normal prior fitted to all hospitals together; small hospitals (thin lines, sizes written on the left) move a lot, large ones (thick lines) much less. The raw rate of a small hospital says little about its true rate.

3.12.5 When there is no formula: sampling the posterior

The beta and normal cases have formulas. Most models do not: the posterior is known only up to its scaling constant. Markov chain Monte Carlo (MCMC) gets around this by drawing samples from the posterior instead of writing it down. The simplest version, the Metropolis algorithm, walks through parameter space:

From the current value, propose a nearby value at random.
If the posterior is higher there, move. If it is lower, move with probability equal to the ratio of the two posterior values; otherwise stay.
Repeat thousands of times. After a warm-up, the visited values are draws from the posterior.

Only the ratio of posterior values is needed, so the unknown constant cancels.

Code

gpa, improved = spector.GPA.to_numpy(), spector.GRADE.to_numpy()


def log_posterior(params: np.ndarray) -> float:
    """Log likelihood of the logistic model plus log of the normal priors (up to a constant)."""
    eta = params[0] + params[1] * gpa
    log_likelihood = np.sum(improved * eta - np.log1p(np.exp(eta)))
    return log_likelihood - np.sum(params**2) / (2 * 10**2)


def metropolis(steps: int = 20_000, jump: float = 0.8) -> np.ndarray:
    """A random-walk Metropolis sampler started at zero."""
    local = np.random.default_rng(0)
    current, current_lp = np.zeros(2), log_posterior(np.zeros(2))
    chain = np.empty((steps, 2))
    for t in range(steps):
        proposal = current + local.normal(0, jump, 2) * [3, 1]
        proposal_lp = log_posterior(proposal)
        if np.log(local.uniform()) < proposal_lp - current_lp:
            current, current_lp = proposal, proposal_lp
        chain[t] = current
    return chain


chain = metropolis()
WARM_UP = 2000
kept = chain[WARM_UP:]
mle = sm.Logit(improved, sm.add_constant(gpa)).fit(disp=0).params
low, high = np.percentile(kept[:, 1], [2.5, 97.5])

fig, axes = style.subplots(3, height=3, width=11)
axes[0].plot(np.arange(WARM_UP), chain[:WARM_UP, 1], color=style.GREY, lw=0.5)
axes[0].plot(np.arange(WARM_UP, len(chain)), kept[:, 1], color=style.BLUE, lw=0.3)
axes[0].set(xlabel="step", ylabel="slope", title="The chain")
axes[1].scatter(kept[::10, 0], kept[::10, 1], s=3, color=style.BLUE, alpha=0.4)
axes[1].set(xlabel="intercept", ylabel="slope", title="Joint posterior")
axes[2].hist(kept[:, 1], bins=50, color=style.BLUE, alpha=0.6)
axes[2].axvline(mle[1], color=style.INK, ls="--", lw=1.2)
for edge in (low, high):
    axes[2].axvline(edge, color=style.CORAL, lw=1.5)
style.note(axes[2], f"95%: {low:.1f} to {high:.1f}", x=0.55, y=0.95, color=style.CORAL)
axes[2].set(
    xlabel="slope (log odds per GPA point)", ylabel="draws", title="Posterior of the slope"
)
fig.tight_layout()
plt.show()

Figure 3.12.5: A Bayesian logistic regression of whether 32 students’ grades improved on their grade point average, with normal priors (standard deviation 10) on intercept and slope, sampled by the Metropolis algorithm. Left: the path of the slope over 20,000 steps; the first 2,000 (grey) are discarded as warm-up. Middle: the draws of intercept and slope; they are strongly correlated. Right: the posterior of the slope, with its 95% credible interval; the dashed line is the maximum likelihood estimate (Chapter 4.4).

Checking an MCMC run is part of the job: the chain should wander freely over the whole posterior (no long trends, no getting stuck), and several chains from different starts should agree. Modern samplers, such as Hamiltonian Monte Carlo in Stan or PyMC, use gradients to take long, well-aimed steps and need far fewer draws.

3.12.6 Priors on coefficients are penalties

In a regression with normal noise, a normal prior with mean 0 and standard deviation τ on every coefficient gives a posterior whose most likely value is exactly the ridge estimate of Chapter 6.4, with penalty $\lambda = \sigma^2 / \tau^2$ . A Laplace (double exponential) prior gives the lasso. A tight prior is a strong penalty; a vague prior is almost none.

Code

X = (diabetes.data - diabetes.data.mean(axis=0)) / diabetes.data.std(axis=0)
y = diabetes.target - diabetes.target.mean()
sigma2 = np.var(y - X @ np.linalg.lstsq(X, y, rcond=None)[0])
taus = np.geomspace(0.3, 300, 60)
paths = np.array([np.linalg.solve(X.T @ X + sigma2 / t**2 * np.eye(10), X.T @ y) for t in taus])

fig, ax = plt.subplots(figsize=(8, 3.4))
for k, name in enumerate(diabetes.feature_names):
    ax.plot(taus, paths[:, k], lw=1.6)
    ax.text(taus[-1] * 1.1, paths[-1, k], name, fontsize=7.5, va="center")
ax.set(
    xscale="log",
    xlim=(taus[0], taus[-1] * 3),
    xlabel="prior standard deviation τ (log scale)",
    ylabel="posterior mean",
)
plt.show()

Figure 3.12.6: Posterior means of the ten diabetes coefficients (standardised predictors) as the prior standard deviation τ on each coefficient grows from very tight (left: all coefficients held near 0) to very vague (right: the least squares estimates). This is the ridge path of Chapter 6.4, read as a family of priors.

3.12.7 Bayes factors

To compare two hypotheses, Bayesians compare how well each predicted the data. The Bayes factor is the ratio of the probabilities of the data under each, averaging over the prior within each hypothesis:

$\text{BF}_{01} = \frac{p(\text{data} \mid H_0)}{p(\text{data} \mid H_1)}$

A Bayes factor of 10 means the data were ten times as probable under $H_0$ . Unlike a p-value, it can support the null hypothesis, not only fail to reject it. It also behaves differently as data pile up:

Code

tosses = np.unique(np.geomspace(100, 20_000, 60).astype(int))
heads = np.round(0.52 * tosses).astype(int)
p_values = np.array([stats.binomtest(h, n, 0.5).pvalue for h, n in zip(heads, tosses)])
log_bf01 = (
    heads * np.log(0.5) + (tosses - heads) * np.log(0.5) - betaln(heads + 1, tosses - heads + 1)
)

fig, axes = style.subplots(2, height=3, width=10)
axes[0].plot(tosses, p_values, color=style.CORAL, lw=2)
axes[0].axhline(0.05, color=style.INK, ls="--", lw=1)
axes[0].set(
    xscale="log",
    yscale="log",
    xlabel="tosses (log scale)",
    ylabel="p-value (log scale)",
    title="p-value",
)
axes[1].plot(tosses, np.exp(log_bf01), color=style.BLUE, lw=2)
axes[1].axhline(1, color=style.INK, ls="--", lw=1)
axes[1].set(
    xscale="log",
    yscale="log",
    xlabel="tosses (log scale)",
    ylabel="Bayes factor for H0 (log scale)",
    title="Bayes factor",
)
fig.tight_layout()
plt.show()

Figure 3.12.7: A coin that lands heads 52% of the time in every sample, for sample sizes from 100 to 20,000. H0: the coin is fair; H1: any rate is possible (flat prior). Left: the two-sided p-value falls below 0.05 near 2,500 tosses and keeps falling. Right: the Bayes factor for H0 favours the fair coin, most strongly (about 11 to 1) near 1,000 tosses; at 2,500 tosses, where the p-value is already below 0.05, it still favours H0 about 5 to 1, and it only turns to favour H1 beyond about 5,000. A small but real difference is ‘significant’ long before it is convincing evidence against a fair coin under a flat prior; how to weigh such cases is a choice, not a fact.

Bayes factors depend strongly on the prior under $H_1$ : a very vague prior spreads $H_1$ ’s predictions thinly and so favours $H_0$ . Many Bayesians therefore prefer to report the posterior, and the probability of effects large enough to matter, rather than a test.

3.12.8 Choosing priors

Weakly informative priors rule out absurd values (a success rate of exactly 0, an effect of a thousand standard deviations) but let the data speak otherwise. They are the usual default.
Informative priors bring in real knowledge, from earlier studies or the physics of a problem; they should be stated and justified like any other assumption.
A sensitivity analysis refits with a few different priors. If the conclusion changes, the data alone do not settle the question, which is itself worth knowing.
Hierarchical priors learn the prior from the data, as in Section 3.12.4: hospitals share a common distribution whose centre and spread are estimated together with each hospital’s rate.

3.12.9 In one picture

Posterior ∝ likelihood × prior; the posterior is the full answer.
Beta prior + binomial data → Beta(a + k, b + n − k); the prior counts as earlier observations.
A credible interval contains θ with the stated probability; posteriors give the probability of any claim directly.
Posteriors are precision-weighted averages: shrinkage towards the prior, strongest for small samples.
MCMC draws samples from posteriors without a formula; check the chains.
Normal and Laplace priors on coefficients are ridge and lasso penalties.
Bayes factors compare hypotheses by how well they predicted the data, and depend on the priors.

Next: choosing among all these tests, and how to report them.