14 - Bootstrap and Resampling Methods

You have fifteen response times of a service and you want to say how uncertain their median is. For the mean there is a formula, the standard error s/√n; for the median, for a ratio between two means or for a concentration index the formula either does not exist or requires assumptions about the distribution that you cannot justify. The bootstrap, proposed by Efron in 1979, sidesteps the problem with a concrete idea: if the sample is the best description you have of the population, draw new samples from it and see how much the statistic varies. This article, the fourteenth in the Probability and Statistics series, builds the bootstrap by hand, compares three kinds of interval by measuring their coverage, shows a case in which the method fails, and arrives at the permutation test and the block bootstrap for dependent data. It assumes the sampling distribution of an estimator, seen in article 6.

The algorithm: resampling with replacement

Given a sample of n values, a bootstrap sample is obtained by drawing n values from the sample with replacement: some values appear several times, others never. On each one you recompute the statistic of interest (here the median), and repeating B times you obtain the bootstrap distribution, which approximates the sampling distribution of the statistic. Its standard deviation is the bootstrap standard error. The listing generates fifteen response times from a lognormal distribution (asymmetric, as latencies often are) with a fixed seed, and uses B = 9999. From that distribution it derives three 95% confidence intervals. The percentile takes the 2.5% and 97.5% quantiles as they are. The basic reflects them around the estimate: 2t − q_97.5 and 2t − q_2.5. The BCa (bias-corrected and accelerated, Efron 1987) corrects the quantiles for bias and for the skewness of the error, and here it is computed by scipy.stats.bootstrap.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2026)
# 15 response times (ms) of a service: skewed distribution, with long tails
x = np.round(rng.lognormal(mean=5.0, sigma=0.6, size=15)).astype(int)
print("sorted data:", np.sort(x))
n, B = len(x), 9999
med = np.median(x)

# Resample n values WITH replacement, B times, and recompute the median on each sample
idx = rng.integers(0, n, size=(B, n))
boot = np.median(x[idx], axis=1)
print(f"observed median = {med:.1f}   bootstrap standard error = {boot.std(ddof=1):.2f}")
print("quantiles of the bootstrap distribution (2.5%, 50%, 97.5%):", np.round(np.quantile(boot, [.025, .5, .975]), 1))

q_lo, q_hi = np.quantile(boot, [.025, .975])
perc = (q_lo, q_hi)                               # percentile: the quantiles as they are
basic = (2 * med - q_hi, 2 * med - q_lo)          # basic: reflects the quantiles around the estimate
bca = stats.bootstrap((x,), np.median, confidence_level=0.95, n_resamples=B, method="BCa", random_state=1).confidence_interval
print(f"percentile 95% [{perc[0]:.1f}, {perc[1]:.1f}]   basic [{basic[0]:.1f}, {basic[1]:.1f}]   BCa [{bca.low:.1f}, {bca.high:.1f}]")
print(f"how many bootstrap medians are distinct values: {len(np.unique(boot))} out of {B}")

# How much noise does a small B add? Repeat the standard error estimate 200 times for each B
for b in (50, 200, 1000, 5000):
    estimates = [np.median(x[rng.integers(0, n, size=(b, n))], axis=1).std(ddof=1) for _ in range(200)]
    print(f"B = {b:5d}: estimated standard error {np.mean(estimates):6.2f}  (variability across repetitions: sd {np.std(estimates):.2f})")

# Output:
# sorted data: [ 48  92 123 125 126 130 141 143 163 171 178 202 218 229 343]
# observed median = 143.0   bootstrap standard error = 18.51
# quantiles of the bootstrap distribution (2.5%, 50%, 97.5%): [125. 143. 178.]
# percentile 95% [125.0, 178.0]   basic [108.0, 161.0]   BCa [126.0, 202.0]
# how many bootstrap medians are distinct values: 13 out of 9999
# B =    50: estimated standard error  18.40  (variability across repetitions: sd 1.94)
# B =   200: estimated standard error  18.42  (variability across repetitions: sd 0.99)
# B =  1000: estimated standard error  18.46  (variability across repetitions: sd 0.41)
# B =  5000: estimated standard error  18.46  (variability across repetitions: sd 0.20)

The observed median is 143 ms and the bootstrap standard error is 18.51. The three intervals are different: percentile [125, 178], basic [108, 161], BCa [126, 202]. This is not a mistake: they are different constructions that coincide only when the bootstrap distribution is symmetric and centred on the estimate, and here it is not. There is also a warning sign: the 9999 bootstrap medians take only 13 distinct values, because the median of 15 numbers can only be one of the 15 observed values; with samples this small, the bootstrap distribution of the median is «stepped», and the quantiles fall on those steps.

In the last part the listing repeats the standard error estimate 200 times for each B: with B = 50 the variability across repetitions is 1.94, with B = 1000 it drops to 0.41 and with B = 5000 to 0.20, while the average estimate (18.4-18.5) hardly changes. The number of resamples does not make the estimate better, only less noisy: beyond a few thousand you pay time to reduce an error that is already small compared with the one due to the fifteen observations alone.

Which interval? Measuring coverage

A 95% interval should contain the true value in 95% of samples. With real data you never know, but in simulation you do: the program below uses a lognormal population with true median 1 and true mean e^0.5 = 1.6487, draws 1000 samples of 30 observations, builds the intervals for the median and for the mean with each method (B = 1000) and counts how many times they contain the true value. It also includes the classic t interval for the mean as a reference.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
n, B, R = 30, 1000, 1000                       # sample size, resamples, repetitions
true = {"median": 1.0, "mean": np.exp(0.5)}    # lognormal(0, 1) population: known true values
statistics = {"median": np.median, "mean": np.mean}

coverage = {}
def count(key, lo, hi, truth):
    coverage[key] = coverage.get(key, 0) + (lo <= truth <= hi)

for _ in range(R):
    x = rng.lognormal(0.0, 1.0, n)
    for name, f in statistics.items():
        t0 = f(x)
        boot = f(x[rng.integers(0, n, size=(B, n))], axis=1)   # B resamples with replacement
        lo, hi = np.quantile(boot, [.025, .975])
        count((name, "percentile"), lo, hi, true[name])
        count((name, "basic"), 2 * t0 - hi, 2 * t0 - lo, true[name])
        ci = stats.bootstrap((x,), f, n_resamples=B, method="BCa", random_state=rng).confidence_interval
        count((name, "BCa"), ci.low, ci.high, true[name])
    lo, hi = stats.t.interval(0.95, n - 1, loc=x.mean(), scale=stats.sem(x))   # classic t interval
    count(("mean", "t"), lo, hi, true["mean"])

print(f"coverage over {R} samples (n={n}, B={B}); nominal 0.95, error ~0.007")
for (s, m), c in coverage.items():
    print(f"  {s:8s} {m:10s} {c / R:.3f}")

# Output:
# coverage over 1000 samples (n=30, B=1000); nominal 0.95, error ~0.007
#   median   percentile 0.952
#   median   basic      0.814
#   median   BCa        0.954
#   mean     percentile 0.871
#   mean     basic      0.845
#   mean     BCa        0.890
#   mean     t          0.876

For the median, percentile (0.952) and BCa (0.954) are close to the nominal 0.95 (the simulation error is about 0.007), while the basic is at 0.814: the method that reflects the quantiles can be badly wrong when the distribution is asymmetric. For the mean of lognormal data with n = 30 no method reaches 95%: BCa 0.890, t 0.876, percentile 0.871, basic 0.845. This is the limit of all methods on the mean of very skewed data with few observations: the sampling distribution of the mean is not yet close to the normal (article 5), and the bootstrap, which estimates the distribution from the data, inherits the same slowness. In this experiment BCa is the best or tied for best in every case, but it does not save a sample that is too small. The bootstrap does not create information that the data do not contain.

Where the bootstrap fails: the maximum

The ordinary bootstrap does not work for statistics that depend on the extremes of the sample. The classic case is the maximum of a uniform sample on (0, θ), used to estimate θ. In a resample of n values, the original maximum is present with probability 1 − (1 − 1/n)ⁿ, which for n = 50 is 0.6358 and tends to 1 − 1/e ≈ 0.632 for large n: in 63% of resamples the bootstrap maximum is exactly the observed one. In the real world, on the other hand, the maximum of a second independent sample coincides with that of the first with probability zero. The bootstrap distribution has a big atom where the true distribution is continuous, and it does not approximate it.

import numpy as np

rng = np.random.default_rng(11)

# 1) Where the bootstrap fails: the maximum of a uniform sample on (0, theta)
n = 50
x = rng.uniform(0, 10, n)
maximum = x.max()
boot_max = x[rng.integers(0, n, size=(100_000, n))].max(axis=1)
exact = 1 - (1 - 1 / n) ** n                   # P(the sample maximum is in the resample)
print(f"P(bootstrap maximum = observed maximum): exact {exact:.4f}   simulated {np.mean(boot_max == maximum):.4f}")
# in the real world, the maximum of another sample is almost never equal to the observed one
real = rng.uniform(0, 10, size=(100_000, n)).max(axis=1)
print(f"in the real world P(maximum = observed value) = {np.mean(real == maximum):.4f}")

# Output:
# P(bootstrap maximum = observed maximum): exact 0.6358   simulated 0.6345
# in the real world P(maximum = observed value) = 0.0000

The two probabilities computed by the listing, the exact and the simulated, are 0.6358 and 0.6345, and in the real world the probability of coincidence is 0.0000 over 100,000 samples. The practical consequence is that bootstrap intervals for a maximum, a minimum or a very extreme quantile do not have the stated coverage, whatever B is. For these statistics other approaches are needed, such as the exact distribution of the maximum or the bootstrap with subsamples of size m < n.

The jackknife, the older brother

The jackknife (Quenouille, then Tukey) recomputes the statistic leaving out one observation at a time, obtaining n partial values. From those values bias and standard error are estimated with closed formulas, with no random numbers. The listing verifies two exact properties. The first: applied to the variance with divisor n (biased), the jackknife correction returns exactly the variance with divisor n − 1: 1.749077 against 1.749077. The second: the jackknife standard error of the mean coincides with s/√n, 0.341475 in both cases. For smooth statistics like these the jackknife is exact and fast.

import numpy as np

rng = np.random.default_rng(12)

# 2) Jackknife: leave out one observation at a time
y = rng.lognormal(0, 1, 15)
m = len(y)
def jackknife(stat, y):
    partial = np.array([stat(np.delete(y, i)) for i in range(len(y))])
    t0 = stat(y)
    bias = (len(y) - 1) * (partial.mean() - t0)
    se = np.sqrt((len(y) - 1) / len(y) * np.sum((partial - partial.mean()) ** 2))
    return t0, bias, se

t0, bias, se = jackknife(lambda a: a.var(ddof=0), y)
print(f"variance with 1/n: {t0:.6f}  jackknife bias {bias:.6f}  corrected {t0 - bias:.6f}  variance with 1/(n-1): {y.var(ddof=1):.6f}")
print(f"standard error of the mean: jackknife {jackknife(np.mean, y)[2]:.6f}   formula s/sqrt(n) {y.std(ddof=1) / np.sqrt(m):.6f}")

# 3) Standard error of the median (n = 15): the truth and two estimators, averaged over 300 samples
true_se = np.median(rng.lognormal(0, 1, size=(200_000, m)), axis=1).std()
jk, bs = [], []
for _ in range(300):
    z = rng.lognormal(0, 1, m)
    jk.append(jackknife(np.median, z)[2])
    bs.append(np.median(z[rng.integers(0, m, size=(500, m))], axis=1).std(ddof=1))
print(f"standard error of the median, true {true_se:.3f}")
print(f"  jackknife: mean {np.mean(jk):.3f}, sd across samples {np.std(jk):.3f}, maximum {np.max(jk):.3f}")
print(f"  bootstrap: mean {np.mean(bs):.3f}, sd across samples {np.std(bs):.3f}, maximum {np.max(bs):.3f}")

# Output:
# variance with 1/n: 1.632471  jackknife bias -0.116605  corrected 1.749077  variance with 1/(n-1): 1.749077
# standard error of the mean: jackknife 0.341475   formula s/sqrt(n) 0.341475
# standard error of the median, true 0.344
#   jackknife: mean 0.326, sd across samples 0.265, maximum 1.947
#   bootstrap: mean 0.399, sd across samples 0.191, maximum 1.137

For the median, which is not a smooth function of the data, the picture changes. The true standard error (estimated with 200,000 samples) is 0.344. Over 300 samples of 15 observations, the jackknife has mean 0.326 but a standard deviation across samples of 0.265 and a maximum of 1.947: it is very unstable. The bootstrap has mean 0.399 (slightly high) but is more stable, with deviation 0.191. You should therefore not rely on the jackknife for the median, while the bootstrap is more regular; neither is precise with 15 data points. The very title of Efron's article (1979) presents the bootstrap as a new look at the jackknife.

Permutation test: the test when the null hypothesis is exchange

The bootstrap estimates a sampling distribution; the permutation test directly builds the distribution of a statistic under the null hypothesis. If two packing methods are equivalent, the labels «A» and «B» attached to the 12 times are arbitrary, and every way of distributing them among the data is equally likely. You recompute the difference between the means for each redistribution and see how extreme the observed difference is. With 6 and 6 parcels the redistributions are 924 and all are enumerated; with 20 parcels per group there would be about 138 billion and they are sampled at random, adding 1 to numerator and denominator so that a Monte Carlo p-value is never zero.

import itertools
import numpy as np

# Permutation test: two packing methods, time (minutes) for 6 parcels each
a = np.array([12.1, 14.3, 13.8, 15.2, 12.9, 14.7])
b = np.array([15.9, 14.8, 16.4, 13.9, 17.2, 15.5])
everyone = np.concatenate([a, b])
obs = a.mean() - b.mean()

# Exact: all C(12,6) = 924 ways of choosing which 6 values to call "A"
diff = np.array([everyone[list(i)].mean() - np.delete(everyone, list(i)).mean()
                 for i in itertools.combinations(range(12), 6)])
p_exact = np.mean(np.abs(diff) >= abs(obs) - 1e-12)
print(f"observed difference {obs:.3f}   exact p = {int(p_exact * len(diff))}/{len(diff)} = {p_exact:.4f}")

# Monte Carlo: 20000 random permutations (the +1 avoids a p-value of exactly zero)
rng = np.random.default_rng(3)
mc = np.array([(lambda p: p[:6].mean() - p[6:].mean())(rng.permutation(everyone)) for _ in range(20_000)])
p_mc = (1 + np.sum(np.abs(mc) >= abs(obs) - 1e-12)) / (1 + len(mc))
print(f"Monte Carlo p = {p_mc:.4f} +- {np.sqrt(p_mc * (1 - p_mc) / len(mc)):.4f}")

# Output:
# observed difference -1.783   exact p = 26/924 = 0.0281
# Monte Carlo p = 0.0269 +- 0.0011

The observed difference is −1.783 minutes. Only 26 of the 924 permutations produce a difference at least this large in absolute value: exact p 26/924 = 0.0281. The Monte Carlo with 20,000 permutations gives 0.0269 ± 0.0011, consistent. The test is «exact» in the sense that it uses no approximations, but it answers a precise question: if the groups were exchangeable, how surprising would this result be? It does not allow you to infer anything about the population beyond that comparison. The difference between permutation and bootstrap for tests comes down to this: permutation resamples without replacement and imposes the null hypothesis; the bootstrap resamples with replacement and, for a test, requires first «shifting» the data so that the null hypothesis is true.

When the data are not independent: the block bootstrap

The ordinary bootstrap assumes independent observations: shuffling the data of a time series destroys the order and therefore the dependence. The block bootstrap (Künsch, 1989) resamples consecutive blocks of length L, which preserve the dependence within them. The listing uses an autoregressive AR(1) series with coefficient ρ = 0.7, in which each value depends on the previous one: its true standard error of the mean is known by simulation (20,000 series) and also by the asymptotic formula, √(σ_x²/n · (1 + ρ)/(1 − ρ)). The standard errors estimated by the bootstrap are then compared for five block lengths, L = 1 being the ordinary bootstrap.

import numpy as np

# AR(1) time series: x_t = rho * x_{t-1} + noise. The observations are NOT independent.
rho, n = 0.7, 200
rng = np.random.default_rng(21)

def series(r, k, burn=200):
    e = r.standard_normal((k, n + burn))
    x = np.zeros((k, n + burn))
    for t in range(1, n + burn):
        x[:, t] = rho * x[:, t - 1] + e[:, t]
    return x[:, burn:]

# Truth: standard error of the mean over 20,000 independent series, and the asymptotic theoretical value
means = series(rng, 20_000).mean(axis=1)
sigma2_x = 1 / (1 - rho**2)
theoretical = np.sqrt(sigma2_x / n * (1 + rho) / (1 - rho))
print(f"standard error of the mean: simulated {means.std():.4f}   asymptotic theoretical {theoretical:.4f}   if the data were independent {np.sqrt(sigma2_x / n):.4f}")

def se_bootstrap(x, L, B=500):
    k = int(np.ceil(n / L))                      # number of consecutive blocks of length L
    starts = rng.integers(0, n - L + 1, size=(B, k))
    samples = np.concatenate([x[starts + j][:, :, None] for j in range(L)], axis=2).reshape(B, -1)[:, :n]
    return samples.mean(axis=1).std(ddof=1)

all_series = series(rng, 200)
for L in (1, 5, 10, 20, 40):                      # L = 1 is the ordinary bootstrap
    estimates = [se_bootstrap(x, L) for x in all_series]
    print(f"block L = {L:2d}: estimated standard error (mean over 200 series) {np.mean(estimates):.4f}  sd {np.std(estimates):.4f}")

# Output:
# standard error of the mean: simulated 0.2319   asymptotic theoretical 0.2357   if the data were independent 0.0990
# block L =  1: estimated standard error (mean over 200 series) 0.0980  sd 0.0091
# block L =  5: estimated standard error (mean over 200 series) 0.1689  sd 0.0227
# block L = 10: estimated standard error (mean over 200 series) 0.1932  sd 0.0334
# block L = 20: estimated standard error (mean over 200 series) 0.2018  sd 0.0445
# block L = 40: estimated standard error (mean over 200 series) 0.1886  sd 0.0591

The true standard error is 0.2319, and the asymptotic formula gives 0.2357. If the data were independent it would be 0.0990: ignoring the dependence underestimates the uncertainty by a factor of 2.3, and the ordinary bootstrap (L = 1) indeed gives 0.0980. With growing blocks the estimate rises to 0.1689 (L = 5), 0.1932 (L = 10) and 0.2018 (L = 20), then falls back to 0.1886 with L = 40, when the blocks are few and the variability across series increases (sd 0.0591). No length reaches the true value: the moving-block bootstrap tends to underestimate with n = 200. The point is not precision, but direction: the choice of L is a compromise between preserving the dependence (long blocks) and having enough distinct blocks (short blocks).

A worked exercise and a test to do on your own

What is the probability that the maximum of a sample of n = 10 appears in a bootstrap resample? Each draw picks the maximum with probability 1/10, so it misses it with probability 9/10; in 10 independent draws it always misses it with probability (9/10)¹⁰ = 0.3487, and includes it with probability 1 − 0.3487 = 0.6513. Checked with the calculation in Python: 1 - 0.9**10 returns 0.6513215599.

What you can do now

You can build a bootstrap by hand, choose between percentile, basic and BCa intervals knowing that the basic can have very low coverage on asymmetric data, and recognise the cases in which it does not work, such as the maximum of a sample. You can apply the jackknife, the permutation test and the block bootstrap, with the limitations measured above. The bootstrap is the frequentist alternative to the credible intervals of article 11; for tests without distributional assumptions see also the non-parametric tests. The last article in the series, exploratory data analysis, is the step that comes before any of these calculations.