3.3 The bootstrap

We cannot draw new samples from the population, but we can draw them from the sample, which is our best picture of the population. How much an estimate varies across such resamples estimates how much it would vary across real samples.

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
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)
income, food = engel.income.to_numpy(), engel.foodexp.to_numpy()
sun = (
    sm.datasets.sunspots.load_pandas().data.SUNACTIVITY.to_numpy()
)  # sunspots per year, 1700 to 2008
RESAMPLES = 5_000

3.3.1 The idea

Every estimate in Chapter 3.2 came with a standard error: how much the estimate would change if we drew a new sample. The formulas for it exist only for simple estimates, and they assume things, such as a normal shape, that may not hold.

The bootstrap (Efron, 1979) replaces the formula with computing. The sample is our best picture of the population, so we treat it as one: draw a new “sample” of the same size from it, with replacement, and compute the estimate. Repeat thousands of times. The spread of these estimates stands in for the spread across real samples.

sample = np.round(rng.lognormal(0, 0.6, 12), 1)
CHIP_COLORS = plt.cm.tab20(np.linspace(0, 1, 12))


def plot_row(ax: Axes, y: float, picks: np.ndarray, label: str) -> None:
    """One row of chips: the values picked (repeats darker), with the row's mean."""
    counts = np.bincount(picks, minlength=12)
    for x, i in enumerate(np.sort(picks)):
        shade = min(1, 0.45 + 0.25 * (counts[i] - 1))
        ax.add_patch(plt.Rectangle((x, y), 0.9, 0.7, color=CHIP_COLORS[i], alpha=shade))
        ax.text(
            x + 0.45,
            y + 0.35,
            f"{sample[i]:.1f}",
            ha="center",
            va="center",
            fontsize=7.5,
            color="black",
        )
    ax.text(-0.3, y + 0.35, label, ha="right", va="center", fontsize=9)
    ax.text(12.3, y + 0.35, f"mean {sample[picks].mean():.2f}", va="center", fontsize=9)


def plot_chips(ax: Axes, rows: int = 3) -> None:
    """The sample and a few resamples as rows of coloured chips with their means."""
    plot_row(ax, rows + 0.6, np.arange(12), "sample")
    for r in range(rows):
        plot_row(ax, rows - 1 - r, rng.integers(0, 12, 12), f"resample {r + 1}")
    ax.set(xlim=(-2.2, 15), ylim=(-0.2, rows + 1.5))
    ax.set_axis_off()


fig, ax = plt.subplots(figsize=(8.5, 2.8))
plot_chips(ax)
plt.show()

Figure 3.3.1: A sample of 12 values (top row) and three bootstrap resamples. Each resample draws 12 times from the sample with replacement: some values appear twice or three times (darker), others not at all. Each resample has its own mean (right).

On average a resample contains about 63% of the distinct original values; the rest are repeats. That is exactly enough variation to mimic drawing new samples.

3.3.2 Standard errors without a formula

For the mean there is a formula, $s/\sqrt n$, so we can check the bootstrap against it. For the median there is no simple formula at all, and the bootstrap still works.

picks = rng.integers(0, len(income), (RESAMPLES, len(income)))
boot_means = income[picks].mean(axis=1)
boot_medians = np.median(income[picks], axis=1)


def plot_boot(ax: Axes, values: np.ndarray, title: str, formula_se: float | None) -> None:
    """A bootstrap distribution with its standard error, and the formula's curve if there is one."""
    width = style.hist(ax, values, bins=50)
    note = f"bootstrap SE {values.std(ddof=1):.1f}"
    if formula_se:
        grid = np.linspace(values.min(), values.max(), 200)
        ax.plot(
            grid,
            stats.norm.pdf(grid, values.mean(), formula_se) * len(values) * width,
            color=style.CORAL,
        )
        note += f"\nformula SE {formula_se:.1f}"
    style.note(ax, note, x=0.03, color=style.INK)
    ax.set(title=title, xlabel="household income", ylabel="resamples")


fig, axes = style.subplots(2, height=3.1, width=9.5)
plot_boot(axes[0], boot_means, "Mean", income.std(ddof=1) / np.sqrt(len(income)))
plot_boot(axes[1], boot_medians, "Median", None)
fig.tight_layout()
plt.show()

Figure 3.3.2: Bootstrap distributions for the 235 household incomes. Left: the mean, whose bootstrap standard error matches the formula s/√n, with the normal curve that formula implies (red). Right: the median, which has no simple formula; its resamples can only take values that are (or lie between two) observed incomes, so its histogram is spiky.

The same recipe gives standard errors for anything you can compute: a median, a ratio, a correlation (Chapter 3.2), a trimmed mean, the 90th percentile, the difference between two groups’ medians. A few thousand resamples are enough for a standard error; intervals need more, since they depend on the tails.

3.3.3 Bootstrap intervals

The simplest interval takes the middle 95% of the bootstrap estimates: the percentile interval. It works well when the bootstrap distribution is roughly symmetric and centred on the estimate. When it is skewed, or shifted, two refinements help. The basic interval reflects the percentiles around the estimate. The BCa interval (bias-corrected and accelerated) shifts and stretches the percentiles to correct for bias and skew, and is usually the most accurate.

How accurate is accurate? Simulate: draw many small samples from a skewed population whose true mean is known, compute each kind of interval, and count how often it contains the truth:

TRUE_MEAN = np.exp(0.8**2 / 2)
SIMULATIONS, N = 300, 20


def interval_hits() -> dict[str, float]:
    """For each method, the share of intervals that catch the true mean."""
    hits = {"t-interval": 0, "percentile": 0, "basic": 0, "BCa": 0}
    for _ in range(SIMULATIONS):
        x = rng.lognormal(0, 0.8, N)
        half = stats.t.ppf(0.975, N - 1) * x.std(ddof=1) / np.sqrt(N)
        hits["t-interval"] += x.mean() - half <= TRUE_MEAN <= x.mean() + half
        for method in ("percentile", "basic", "BCa"):
            ci = stats.bootstrap(
                (x,),
                np.mean,
                n_resamples=1_000,
                method=method.lower() if method != "BCa" else "BCa",
                random_state=rng,
            ).confidence_interval
            hits[method] += ci.low <= TRUE_MEAN <= ci.high
    return {name: count / SIMULATIONS for name, count in hits.items()}


coverage = interval_hits()


def plot_hits(ax: Axes) -> None:
    """One bar per method: how often its 95% interval caught the truth."""
    names, values = list(coverage), list(coverage.values())
    ax.barh(names, values, color=[style.GREY, style.BLUE, style.BLUE, style.GREEN], alpha=0.7)
    for i, value in enumerate(values):
        ax.text(value + 0.005, i, f"{value:.0%}", va="center", fontsize=9)
    ax.axvline(0.95, color=style.INK, ls="--", lw=1)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.invert_yaxis()
    ax.set(xlabel="intervals that contained the true mean", xlim=(0.75, 1))


fig, ax = plt.subplots(figsize=(8, 2.6))
plot_hits(ax)
plt.show()

Figure 3.3.3: 300 samples of 20 values each from a skewed (log-normal) population with a known mean. Share of the 95% intervals that contained the true mean, for the t-interval of Chapter 3.2 and three bootstrap intervals. With small, skewed samples none reaches 95% (dashed); for a plain mean, the bootstrap intervals do no better than the formula.

The bootstrap is not magic: from 20 values of a skewed variable, the sample simply does not show the long tail well, and every method is too optimistic. For a plain mean, where the formula already exists, the bootstrap brings nothing extra. Its value is for the estimates where no formula exists, and with larger samples, where all the methods come close to 95% and BCa is usually the most accurate.

3.3.4 Bootstrapping a regression

To bootstrap a regression slope, resample whole rows, keeping each household’s income and food spending together, refit the line, and keep the slope. This case resampling makes no assumption about the noise. The textbook standard error of a slope (Chapter 4.1) assumes the noise is equally large everywhere, and the food spending fans out (Chapter 4.3):

X = np.c_[np.ones(len(income)), income]
slopes = np.array([np.linalg.lstsq(X[p], food[p], rcond=None)[0][1] for p in picks])
ols = sm.OLS(food, X).fit()


def plot_slopes(ax: Axes) -> None:
    """Bootstrap slopes with the formula's normal curve."""
    width = style.hist(ax, slopes, bins=50)
    grid = np.linspace(slopes.min(), slopes.max(), 200)
    ax.plot(
        grid,
        stats.norm.pdf(grid, ols.params[1], ols.bse[1]) * len(slopes) * width,
        color=style.CORAL,
    )
    style.note(
        ax,
        f"bootstrap SE {slopes.std(ddof=1):.4f}\nformula SE {ols.bse[1]:.4f}",
        x=0.03,
        color=style.INK,
    )
    ax.set(xlabel="slope: food spending per unit of income", ylabel="resamples")


fig, ax = plt.subplots(figsize=(8, 3.1))
plot_slopes(ax)
plt.show()

Figure 3.3.4: Slopes of food spending on income in 5,000 case resamples (bars), against the normal curve of the textbook standard error (red). The bootstrap distribution is wider: the formula, which assumes equal noise everywhere, is too optimistic for these fanning data.

Another way, residual resampling, keeps the incomes fixed and reshuffles the residuals onto the fitted line. It assumes the noise is the same everywhere, so it would repeat the formula’s mistake here; case resampling is the safer default.

3.3.5 Bootstrap tests

The bootstrap can also test. Is the average income 900? The bootstrap distribution of the mean is centred on the observed mean, 982. To see what chance produces if the true mean were 900, shift the data so their mean is exactly 900, resample from the shifted data, and see how often the resampled means land as far from 900 as 982 does.

NULL_MEAN = 900
shifted = income - income.mean() + NULL_MEAN
null_means = shifted[picks].mean(axis=1)
distance = abs(income.mean() - NULL_MEAN)
p_boot = (abs(null_means - NULL_MEAN) >= distance).mean()


def plot_test(ax: Axes) -> None:
    """Resampled means under the null hypothesis, with the observed mean."""
    style.hist(ax, null_means, bins=50, color=style.GREY)
    ax.axvline(income.mean(), color=style.CORAL, ls="--", lw=1.6)
    ax.axvline(NULL_MEAN - distance, color=style.CORAL, ls=":", lw=1.2)
    style.note(
        ax, f"observed mean {income.mean():.0f}\np = {p_boot:.4f}", x=0.75, color=style.CORAL
    )
    ax.set(xlabel="mean income in a resample, if the true mean were 900", ylabel="resamples")


fig, ax = plt.subplots(figsize=(8, 3.1))
plot_test(ax)
plt.show()

Figure 3.3.5: A bootstrap test of ‘the mean income is 900’. The data are shifted to have mean 900 and resampled; the bars show the means that chance gives if the null hypothesis is true. The observed mean (dashed) is far out: hardly any resamples get that far from 900 in either direction.

Permutation tests (Chapter 3.1) and bootstrap tests are siblings. A permutation test makes the null hypothesis true by shuffling labels; a bootstrap test makes it true by shifting the data. Permutation tests are exact when groups were randomised; bootstrap tests handle questions where there is nothing to shuffle, such as a single mean.

3.3.6 When the bootstrap fails

The bootstrap assumes the sample represents the population, and that the units are independent. It breaks when either fails:

uniform = rng.uniform(0, 1, 50)
boot_max = uniform[rng.integers(0, 50, (RESAMPLES, 50))].max(axis=1)
BLOCK = 11


def block_resample(series: np.ndarray, block: int) -> np.ndarray:
    """A moving-block bootstrap resample: random blocks of consecutive values, glued together."""
    starts = rng.integers(0, len(series) - block, len(series) // block + 1)
    return np.concatenate([series[s : s + block] for s in starts])[: len(series)]


def lag_one(series: np.ndarray) -> float:
    """Correlation between each value and the next."""
    return float(np.corrcoef(series[:-1], series[1:])[0, 1])


def plot_max(ax: Axes) -> None:
    """Bootstrap maxima pile up on the sample maximum."""
    style.hist(ax, boot_max, bins=40, color=style.GOLD)
    ax.axvline(uniform.max(), color=style.CORAL, ls="--", lw=1.4)
    share = (boot_max == uniform.max()).mean()
    style.note(ax, f"{share:.0%} equal the\nsample maximum", x=0.03, color=style.INK)
    ax.set(title="The maximum", xlabel="largest value in a resample", ylabel="resamples")


def plot_memory(ax: Axes) -> None:
    """Lag-one autocorrelation in single-year and block resamples, against the real series."""
    single = [lag_one(sun[rng.integers(0, len(sun), len(sun))]) for _ in range(1_000)]
    blocks = [lag_one(block_resample(sun, BLOCK)) for _ in range(1_000)]
    edges = np.linspace(-0.3, 1, 60)
    ax.hist(single, bins=edges, color=style.GREY, alpha=0.6, label="single years")
    ax.hist(blocks, bins=edges, color=style.BLUE, alpha=0.6, label=f"blocks of {BLOCK} years")
    ax.axvline(lag_one(sun), color=style.CORAL, ls="--", lw=1.4, label="real series")
    ax.legend(fontsize=8, loc="upper center")
    ax.set(
        title="Memory in time",
        xlabel="correlation of each year with the next",
        ylabel="resamples",
    )


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

Figure 3.3.6: Left: bootstrapping the largest value of 50 draws from a uniform population. A resample can never exceed the sample’s own maximum, so 63% of resamples return exactly that value: the bootstrap cannot see beyond the edge of the data. Right: sunspots have strong year-to-year memory. Resampling single years (grey) destroys it; resampling blocks of 11 years (blue) keeps the memory within each block.

• Extremes: the bootstrap cannot produce values beyond the sample, so it fails for maxima, minima and very high percentiles.
• Very small samples: 10 values give a poor picture of a population, and so do their resamples.
• Dependent data: for time series, resample whole blocks of consecutive values (the block bootstrap); for clustered data, such as pupils in schools, resample whole clusters.

Two relatives are worth knowing. The parametric bootstrap fits a distribution first (say a log-normal) and draws new samples from the fitted distribution instead of from the data; it is better with small samples if the distribution is right. The jackknife leaves out one case at a time instead of resampling; it is older, cheaper and less flexible, but shows directly how much each case moves the estimate.

3.3.7 In one picture

• Resample the data with replacement, recompute the estimate, repeat thousands of times.
• The spread of the resampled estimates is the bootstrap standard error; the middle 95% is a percentile interval; BCa corrects for bias and skew.
• For regressions, resample cases to stay honest when the noise is uneven.
• A bootstrap test shifts the data to make the null hypothesis true.
• It fails for extremes, tiny samples and dependent data; use blocks or clusters for the last.

Next: comparing two groups with the classical tests.