Why Resampling Works
Imagine you want to know the sampling distribution of some statistic $\hat{\theta}$ — say, the median of a skewed distribution. The median has no simple closed-form standard error when the underlying distribution is unknown. The parametric route requires you to guess a distribution family and derive the variance analytically; the bootstrap avoids this entirely.
The Plug-In Principle
Every bootstrap argument starts with the plug-in principle: replace the unknown population distribution $F$ with the empirical distribution $\hat{F}_n$, which places mass $1/n$ on each observed data point $x_1, \ldots, x_n$.
A bootstrap sample $X^* = (X_1^*, \ldots, X_n^*)$ is drawn i.i.d. from $\hat{F}_n$ — equivalently, sampled with replacement from the original data. The bootstrap estimate of the statistic is $\hat{\theta}^* = \theta(X^*)$. Repeating this $B$ times produces a collection $\{\hat{\theta}^*_b\}_{b=1}^B$ whose distribution approximates the true sampling distribution of $\hat{\theta}$.
The key theoretical result (Efron, 1979; Bickel & Freedman, 1981) is that under mild regularity conditions,
In words: the bootstrap pivot distribution converges to the true pivot distribution as $n \to \infty$. This is a remarkable result — the data is being used both to estimate $\theta$ and to approximate the variation around that estimate.
A Minimal Implementation
The following code implements the bootstrap from scratch for any scalar statistic.
import numpy as np
def bootstrap(
data: np.ndarray,
stat_fn,
B: int = 10_000,
rng: np.random.Generator | None = None,
) -> np.ndarray:
"""Return B bootstrap replicates of stat_fn(data)."""
if rng is None:
rng = np.random.default_rng(42)
n = len(data)
replicates = np.empty(B)
for b in range(B):
sample = rng.choice(data, size=n, replace=True)
replicates[b] = stat_fn(sample)
return replicates
# Example: 95% CI for the median of a log-normal sample
rng = np.random.default_rng(0)
data = rng.lognormal(mean=1.5, sigma=0.8, size=80)
reps = bootstrap(data, np.median, B=10_000, rng=rng)
print(f"Sample median : {np.median(data):.3f}")
print(f"Bootstrap SE : {reps.std(ddof=1):.3f}")
print(f"95% percentile CI: ({np.percentile(reps, 2.5):.3f}, {np.percentile(reps, 97.5):.3f})")How Many Bootstrap Replicates?
For a 95% confidence interval, $B = 1{,}000$ is often sufficient for the interval endpoints to stabilize. For the BCa method (see the next section) or for tail quantiles beyond 99%, use $B \ge 10{,}000$. Doubling $B$ halves the Monte Carlo variance in the endpoint estimates, so there is a clear diminishing-returns curve. A practical default is $B = 9{,}999$ (odd, so the median is a data point).
The Bootstrap Standard Error
The bootstrap standard error is simply the standard deviation of the replicates:
This is a consistent estimator of the true standard error under the same regularity conditions that guarantee bootstrap validity. Unlike the delta method, it does not require you to compute derivatives of the statistic.
Pitfall: When the Bootstrap Fails
The bootstrap is not universally valid. It fails for:
- The sample maximum (and other extreme-order statistics) because <!--MATHBLOCK19--> has bounded support while the true extreme-value distribution may be heavy-tailed.
- Very small <!--MATHBLOCK20--> (rule of thumb: <!--MATHBLOCK21-->) because <!--MATHBLOCK22--> is a poor approximation of <!--MATHBLOCK23-->.
- Dependent data (time series, spatial data) where i.i.d. resampling destroys autocorrelation structure — use the block bootstrap instead.
- Nondifferentiable statistics such as quantiles of discrete distributions: the bootstrap still works in these cases but the coverage can be off by <!--MATHBLOCK24--> rather than the usual <!--MATHBLOCK25-->.