Non-parametric vs parametric bootstrap

Every estimate comes from a sample, and another sample would give another number. The central question of inference is how much that number would change. For the mean there is a ready-made formula; for the median, a quantile, a ratio or a correlation, almost never. The bootstrap answers the question by simulation, using only the data you have. This article builds the method from scratch in both its versions — non-parametric and parametric —, works every step by hand, derives the confidence intervals and compares the two versions in simulations where the answer is known. Every number comes from a reproducible script.

§1The problem: the sampling distribution

Consider observations \(X_1, \dots, X_n\), independent and all with the same distribution \(F\) (in short: i.i.d. from \(F\)). \(F\) is unknown — it is the “population”. We want to learn something about it: a parameter \(\theta = t(F)\), such as the mean, the median or the 90th percentile. To do so we compute an estimator, a function of the data:

Because the \(X_i\) are random, so is \(\hat\theta\). If we could repeat the data collection many times we would get many different \(\hat\theta\); the distribution of those values is the sampling distribution of \(\hat\theta\), which we will call \(G_F\). It depends on \(F\), on \(n\) and on the function \(s\). Almost everything we want to report about an estimate is a summary of \(G_F\):

• Standard error: \(\operatorname{se}(\hat\theta) = \sqrt{\operatorname{Var}_F(\hat\theta)}\) — the typical size of the error.
• Bias: \(\operatorname{bias}(\hat\theta) = \mathbb{E}_F[\hat\theta] - \theta\) — whether the estimator errs systematically to one side.
• Confidence interval: limits \(L, U\) computed from the data with \(\Pr_F(L \le \theta \le U) \approx 1 - \alpha\).

For the sample mean \(\bar X\), classical theory delivers all of it: \(\operatorname{se}(\bar X) = \sigma/\sqrt n\), zero bias and, by the central limit theorem, \(\bar X \approx \mathcal N(\mu, \sigma^2/n)\). But swap the mean for the median and the asymptotic formula becomes

which requires the density \(f\) at the median \(m\) — precisely what we do not know. For the coefficient of variation, a ratio of means, a Gini coefficient or the eigenvalue of a covariance matrix, the formulas are worse or do not exist.

If we knew \(F\), the solution would be simple: simulate. Draw many samples of size \(n\) from \(F\), compute \(\hat\theta\) on each and look at the histogram. That is Monte Carlo. There is only one obstacle: we do not know \(F\). The bootstrap idea (Efron, 1979) is so simple it feels like cheating: replace \(F\) by an estimate \(\hat F\) and simulate from that.

§2The bootstrap principle

Write the same experiment twice. In the real world, nature draws a sample \(\mathbf x = (x_1, \dots, x_n)\) from \(F\) and we compute \(\hat\theta = s(\mathbf x)\). In the bootstrap world, we draw samples \(\mathbf x^* = (x_1^*, \dots, x_n^*)\) from \(\hat F\) ourselves and compute \(\hat\theta^* = s(\mathbf x^*)\). The asterisk always marks “something from the bootstrap world”.

The principle is an analogy: \(\hat F\) is to \(F\) as \(\hat\theta^*\) is to \(\hat\theta\). If \(\hat F\) resembles \(F\), the variation of \(\hat\theta^*\) around \(\hat\theta\) mimics the variation of \(\hat\theta\) around \(\theta\). Formally, we approximate the distribution of \(\hat\theta - \theta\) under \(F\) by the distribution of \(\hat\theta^* - \hat\theta\) under \(\hat F\), conditional on the observed data:

The right-hand side is still a theoretical probability — the ideal bootstrap. It can rarely be computed exactly, so we approximate it by Monte Carlo: draw \(B\) samples from \(\hat F\) and use frequencies. This creates two layers of error that must not be confused:

1. Statistical error — \(\hat F \ne F\). Depends on \(n\) and on the quality of \(\hat F\). It does not shrink with \(B\).
2. Monte Carlo error — \(B\) is finite. Depends only on \(B\) and vanishes as \(B \to \infty\) (§8).

The whole difference between the two bootstraps lies in how \(\hat F\) is built.

§3Two ways to estimate \(F\)

Non-parametric: the empirical distribution function

Assuming nothing about the shape of \(F\), the natural estimator is the empirical distribution function (ECDF):

the proportion of observations less than or equal to \(t\). It is a staircase that rises \(1/n\) at each observed point — that is, a discrete distribution that puts probability \(1/n\) on each \(x_i\). It is a good estimate of \(F\) for precise reasons:

• For each fixed \(t\), \(nF_n(t) \sim \text{Binomial}(n, F(t))\); hence \(\mathbb E[F_n(t)] = F(t)\) (unbiased) and \(\operatorname{Var}[F_n(t)] = F(t)(1-F(t))/n \to 0\).
• Glivenko–Cantelli: \(\sup_t |F_n(t) - F(t)| \to 0\) almost surely — the staircase hugs the true curve everywhere at once.
• DKW inequality: \(\Pr\big(\sup_t |F_n(t) - F(t)| > \varepsilon\big) \le 2e^{-2n\varepsilon^2}\), a finite-\(n\) guarantee.
• \(F_n\) is the non-parametric maximum likelihood estimator of \(F\).

Drawing one observation from \(F_n\) means picking one of the \(x_i\) at random, each with probability \(1/n\). Drawing a sample of size \(n\) from \(F_n\) means repeating that \(n\) times, independently — therefore with replacement. That is the famous “resampling with replacement”: it is not a trick, it is simply what i.i.d. sampling from \(F_n\) means. (Without replacement, every resample would be a permutation of the same data and \(\hat\theta^*\) would always equal \(\hat\theta\).)

Parametric: a fitted model

If we have reasons to believe that \(F\) belongs to a family \(\{F_\eta\}\) indexed by a few parameters \(\eta\) (normal, exponential, Poisson, gamma…), we estimate \(\eta\) — usually by maximum likelihood — and use \(\hat F = F_{\hat\eta}\). Drawing from \(\hat F\) means generating random numbers from the fitted model: values that were never observed appear.

The example we will use throughout: waiting times modelled as exponential, density \(f(x) = \lambda e^{-\lambda x}\) for \(x \ge 0\), mean \(\mu = 1/\lambda\) and standard deviation also \(1/\lambda\). The log-likelihood of a sample is

(\(\ell''(\lambda) = -n/\lambda^2 < 0\), so it is a maximum.) The parametric bootstrap therefore draws from an exponential with mean \(\bar x\).

The data in the figure are the 15 times of the worked example (§7). The two estimates tell different stories: \(F_n\) is faithful to the data and to nothing else; the fitted exponential is smooth, extrapolates into the tail and assumes a shape. The whole comparison in this article boils down to that dilemma: flexibility without assumptions versus efficiency with assumptions.

§4Resampling by hand: a toy example

Take five observations, \(\mathbf x = (2.1,\; 0.4,\; 3.8,\; 1.2,\; 6.5)\), with mean \(2.8\) and median \(2.1\). For an NP resample, we draw five indices between 1 and 5 with replacement and take the corresponding values. For a P resample, we generate five values from an exponential with mean \(2.8\).

NP = non-parametric: indices 1–5 drawn with replacement. P = parametric: draws from an Exponential with mean 2.8 (the sample mean), rounded to one decimal. NP only repeats observed values; P produces new ones (11.1 was never observed).

Three observations that hold for any \(n\):

1. The number of distinct NP resamples is finite. A resample is determined by how many times each \(x_i\) appears, \((k_1, \dots, k_n)\) with \(\sum k_i = n\). By the “stars and bars” argument there are \(\binom{2n-1}{n}\) such combinations: 126 for \(n = 5\), 77,558,760 for \(n = 15\), about \(6.9 \times 10^{10}\) for \(n = 20\). They are not equally likely: the combination \((k_1, \dots, k_n)\) has multinomial probability \(\frac{n!}{k_1!\cdots k_n!}\,n^{-n}\). For \(n = 5\) one could enumerate them all and compute the ideal bootstrap exactly; in practice, we sample.
2. Each resample leaves observations out. The chance that a given observation is never drawn is \((1 - 1/n)^n\): 32.8% with \(n = 5\), 35.5% with \(n = 15\), and \(\to e^{-1} \approx 36.8\%\). On average a resample contains only about 63.2% of the distinct observations; the rest are repeats.
3. The NP never leaves the observed set; the P does. This will matter a lot for statistics that depend on the tail (§11).

§5The algorithms, step by step

Non-parametric bootstrap

1. Compute the estimate on the original data: \(\hat\theta = s(\mathbf x)\).
2. For \(b = 1, \dots, B\):
   1. draw \(n\) indices \(i_1, \dots, i_n\) uniformly from \(\{1, \dots, n\}\), with replacement;
   2. form \(\mathbf x^{*b} = (x_{i_1}, \dots, x_{i_n})\);
   3. compute \(\hat\theta^{*b} = s(\mathbf x^{*b})\) — exactly the same procedure used for \(\hat\theta\).
3. Summarise \(\hat\theta^{*1}, \dots, \hat\theta^{*B}\) (§6).

Parametric bootstrap

1. Choose a model \(\{F_\eta\}\) and check that it fits the data (ECDF against the model, QQ-plot, goodness-of-fit test).
2. Estimate \(\hat\eta\) (maximum likelihood) and compute \(\hat\theta = s(\mathbf x)\).
3. For \(b = 1, \dots, B\):
   1. generate \(x_1^{*b}, \dots, x_n^{*b}\) i.i.d. from \(F_{\hat\eta}\) — same size \(n\);
   2. compute \(\hat\theta^{*b} = s(\mathbf x^{*b})\), redoing inside the loop everything \(s\) does, including any parameter estimation.
4. Summarise \(\hat\theta^{*1}, \dots, \hat\theta^{*B}\) in the same way.

Note: the steps are identical except for the line that generates the data. In Python, with NumPy:

import numpy as np
rng = np.random.default_rng(42)                 # seed: reproducible results
x = np.array([0.6, 5.1, 1.7, 2.7, 3.9, 4.3, 5.1, 0.1, 3.6, 2.5, 6.8, 11.1, 5.5, 1.8, 0.7])
n, B = len(x), 10_000
s = np.median                                   # the statistic: swap at will

theta_hat = s(x)

# non-parametric: resample the data itself, with replacement
idx = rng.integers(0, n, size=(B, n))
theta_np = np.array([s(x[i]) for i in idx])

# parametric: generate from the fitted model (exponential, MLE = mean)
xs = rng.exponential(scale=x.mean(), size=(B, n))
theta_p = np.array([s(row) for row in xs])

for name, t in (("NP", theta_np), ("P", theta_p)):
    se   = t.std(ddof=1)
    bias = t.mean() - theta_hat
    lo, hi = np.quantile(t, [0.025, 0.975])
    print(f"{name}: se={se:.3f}  bias={bias:+.3f}  percentile CI=[{lo:.2f}, {hi:.2f}]")

And in R, where the boot package does both:

library(boot)
x <- c(0.6, 5.1, 1.7, 2.7, 3.9, 4.3, 5.1, 0.1, 3.6, 2.5, 6.8, 11.1, 5.5, 1.8, 0.7)

# NP: the statistic receives the data and the drawn indices
b_np <- boot(x, statistic = function(d, i) median(d[i]), R = 10000)

# P: sim = "parametric"; ran.gen generates one sample from the fitted model
b_p  <- boot(x, statistic = median, R = 10000, sim = "parametric",
             ran.gen = function(d, mle) rexp(length(d), rate = 1 / mle), mle = mean(x))

b_np; b_p                                       # bias and standard error
boot.ci(b_np, type = c("norm", "basic", "perc", "bca"))

§6What to get out of the \(B\) replicates

With the replicates \(\hat\theta^{*1}, \dots, \hat\theta^{*B}\) in hand, write \(\bar\theta^* = \frac1B\sum_b \hat\theta^{*b}\) and \(q^*_p\) for the \(p\)-quantile of the replicates.

Standard error

It is simply the standard deviation of the replicates. As \(B \to \infty\) it converges to the ideal-bootstrap standard error, \(\operatorname{se}_{\hat F}(\hat\theta^*)\). A case you can check by hand: for the mean under NP, \(\operatorname{Var}_{F_n}(\bar x^*) = \hat\sigma^2/n\) with \(\hat\sigma^2 = \frac1n\sum(x_i - \bar x)^2\) (the variance of \(F_n\), divisor \(n\)), so the ideal bootstrap gives \(\sqrt{(n-1)/n}\; s/\sqrt n\) — practically the classical formula.

Bias

By principle (3), \(\mathbb E_F[\hat\theta] - \theta\) is estimated by \(\mathbb E_{\hat F}[\hat\theta^*] - \hat\theta\):

Be careful with \(\tilde\theta\): correcting bias increases variance, and the bias estimate itself is noisy. A rule of thumb (Efron and Tibshirani) is to worry only if \(|\widehat{\operatorname{bias}}|/\widehat{\operatorname{se}} > 0.25\). A large estimated bias is, first of all, a warning sign — we will see one in §7.

Confidence intervals

There are several constructions, all with the same nominal level \(1-\alpha\) and different qualities. We use \(z_p = \Phi^{-1}(p)\).

(a) Normal. Assumes \(\hat\theta\) is approximately normal and uses only the standard error: \(\hat\theta \pm z_{1-\alpha/2}\,\widehat{\operatorname{se}}_B\) (optionally centred at \(\tilde\theta\)).

(b) Basic (or “reverse percentile”). Applies principle (3) literally. If \(\hat\theta^* - \hat\theta\) mimics \(\hat\theta - \theta\), then the quantiles of \(\hat\theta - \theta\) are approximately \(a = q^*_{\alpha/2} - \hat\theta\) and \(b = q^*_{1-\alpha/2} - \hat\theta\), and

Note the flip: the upper quantile becomes the lower limit. If the bootstrap distribution has a long right tail, the basic interval is stretched to the left.

(c) Percentile. Uses the quantiles directly: \(\big[q^*_{\alpha/2},\; q^*_{1-\alpha/2}\big]\). It is the most used and most intuitive, but its justification is different: it is correct if there is a monotone transformation \(h\) (unknown — we never need to find it) such that \(h(\hat\theta) - h(\theta)\) is symmetric around zero and has the same distribution as \(h(\hat\theta^*) - h(\hat\theta)\). Being transformation-invariant, it is sensible for correlations, ratios and variances; it fails badly when the estimator is strongly biased (§11).

(d) Bootstrap-t (studentised). Instead of \(\hat\theta - \theta\), it mimics the t statistic, \(T = (\hat\theta - \theta)/\widehat{\operatorname{se}}\). In each resample it computes \(T^{*b} = (\hat\theta^{*b} - \hat\theta)/\widehat{\operatorname{se}}^{*b}\), where \(\widehat{\operatorname{se}}^{*b}\) is the standard error computed within the resample (by formula, like \(s^*/\sqrt n\) for the mean, or by a nested bootstrap). With \(t^*_p\) the quantiles of \(T^*\), by the same reasoning as (8):

It is Student's interval, except that the t table is replaced by a table built from the data itself — which captures skewness.

(e) BCa (bias-corrected and accelerated). A percentile interval with adjusted levels. Compute

where \(\hat\theta_{(i)}\) is the estimate without observation \(i\) (jackknife) and \(\bar\theta_{(\cdot)}\) their mean. \(\hat z_0\) measures median bias; \(\hat a\) measures how fast the standard error changes with the parameter (skewness). The interval is \(\big[q^*_{\alpha_1},\, q^*_{\alpha_2}\big]\) with

With \(\hat z_0 = \hat a = 0\) it reduces to the percentile.

For those who want the rigour. For smooth statistics, the coverage error of a one-sided limit is \(O(n^{-1/2})\) for the normal, basic and percentile intervals (“first-order accurate”) and \(O(n^{-1})\) for bootstrap-t and BCa (“second-order”) — Hall (1992). In equal-tailed two-sided intervals the \(n^{-1/2}\) terms of the two sides cancel and all of them have error \(O(n^{-1})\), but the constants differ a lot when the data are skewed. The simulation in §10 shows this in practice.

§7Worked example: 15 waiting times

A service desk records the interval, in minutes, between the arrivals of 15 consecutive customers:

0.6 · 5.1 · 1.7 · 2.7 · 3.9 · 4.3 · 5.1 · 0.1 · 3.6 · 2.5 · 6.8 · 11.1 · 5.5 · 1.8 · 0.7

Mean \(\bar x = 3.70\), standard deviation \(s = 2.85\), median \(3.6\). These data were generated by me from an exponential with mean 4. That is on purpose: because we know \(F\), we can compute the true sampling distribution (by drawing 100,000 fresh samples) and grade each method. In real life you do not have that dashed line.

Inter-arrival times are the textbook case for the exponential model (a Poisson process), and the ECDF in §3 does not contradict the fit. We will estimate the mean and the median, with \(B = 10{,}000\).

The mean

Both bootstrap distributions have the right shape — slightly right-skewed, like the true one — but different widths. The NP is too narrow: \(\widehat{\operatorname{se}} = 0.71\) against a true 1.03. The reason is not the bootstrap, it is the sample: the observed standard deviation, 2.85, came out well below the true one, 4. The NP only knows what the sample shows. The P does much better (0.96) because it uses a piece of information from the model: in the exponential, \(\sigma = \mu\), so it estimates the spread by \(\bar x = 3.7\), which in this case is closer to 4 than \(s\) is. This is called “borrowing strength” from the model: with \(n = 15\), estimating \(\sigma\) directly from skewed data is very imprecise.

The median

Two lessons appear here.

The NP bootstrap of the median is discrete. With \(n = 15\) (odd), the median of a resample is always one of the observed values — in practice one of those in the middle of the sample. Across the 10,000 resamples only 12 different values appear. The histogram becomes a set of isolated bars; the standard error is still reasonable (0.93 against 1.05), but the quantiles, and therefore the percentile interval, jump from one data point to another.

The P is centred in the “wrong” place — and that is informative. In the parametric bootstrap world, the true median is that of \(\text{Exp}(\bar x)\), namely \(\bar x \ln 2 = 2.56\). The observed median is 3.6. Hence the estimated bias of \(-0.93\): huge, almost a full standard error. The bootstrap is saying that, if the model is true, the observed sample median is atypical. And indeed the true median is \(4\ln 2 = 2.77\). The practical lesson: if you are going to assume the exponential model, the natural estimator of the median is not the sample median, it is \(\hat m = \bar x \ln 2\) (by invariance of maximum likelihood). It is much more efficient: \(\operatorname{Var}(\bar x \ln 2) = (\ln 2)^2\mu^2/n \approx 0.48\,\mu^2/n\), whereas, by (2), the sample median has asymptotic variance \(1/(4nf(m)^2) = \mu^2/n\) (since \(f(m) = \lambda e^{-\ln 2} = 1/(2\mu)\)). The sample median would need about twice as much data for the same precision — provided the model is right.

Data: the 15 times above. Truth = 100,000 fresh samples from Exp(mean 4). For comparison: \(s/\sqrt n\) = 0.735; the exact NP value as \(B \to \infty\) is \(\sqrt{(n-1)/n}\; s/\sqrt n\) = 0.710; the P value is \(\bar x/\sqrt n\) = 0.955; the true one is \(4/\sqrt{15}\) = 1.033. Student t CI: [2.12, 5.28]. Exact CI under the exponential model (pivot \(\bar x/\mu \sim \text{Gamma}(n,1)/n\)): [2.36, 6.61].

In the table, both percentile CIs for the mean contain the true value 4. The P one, [2.05, 5.83], is 35% wider than the NP one, [2.41, 5.19] — again because the NP inherits this sample's low \(s\). Neither matches the exact interval of the exponential model, [2.36, 6.61], which is shifted to the right: the percentile is only a first-order approximation (§6). A single sample cannot decide which method is better: for that we need to repeat the experiment thousands of times (§10).

§8How many resamples? Choosing \(B\)

\(B\) only controls the Monte Carlo error. For the standard error, a delta-method argument gives

where \(\hat\kappa\) is the excess kurtosis of the bootstrap distribution (\(\approx 0\) if it is nearly normal, giving \(1/\sqrt{2B}\)). With \(B = 200\), about 5%. Checking empirically, repeating the whole bootstrap 400 times for each \(B\):

NP bootstrap of the mean of the 15 times, repeated 400 times for each B. The target (B → ∞) is 0.710. Monte Carlo error falls like \(1/\sqrt B\); for quantiles it falls much more slowly than for the standard error, because a 2.5% quantile depends on the few resamples in the tail.

The “relative error” column follows \(1/\sqrt{2B}\) closely. But the last column shows the important point: extreme quantiles are much less stable than the standard deviation. Rules of thumb: \(B \approx 1{,}000\) is enough for a standard error; for intervals use at least 2,000, and preferably something like 9,999 (with \(B + 1\) a multiple of \(1/\alpha\), the quantiles land on exact replicates). Computing is cheap today: when in doubt, increase \(B\) until the result stops changing in the decimal place you will report. And fix the seed.

What \(B\) does not do: no value of \(B\) fixes the statistical error. With \(B = 10^6\), the NP bootstrap of the mean in §7 still says 0.71.

§9Try it yourself

The panel below runs the two algorithms of §5 on the 15 waiting times. Click once and follow a single resample; then click +1,000 and watch the histograms converge to the bootstrap distributions of §7.

§10Simulation: who gets it right, and when

To compare interval methods, the right metric is coverage: the fraction of times the 95% interval contains the true value when the whole experiment is repeated. The study design:

1. Fix a true distribution \(F\) and a size \(n \in \{10, 20, 50, 100, 250\}\).
2. Draw \(M = 2{,}000\) samples from \(F\). On each, run the bootstraps with \(B = 999\) and build seven 95% intervals for the mean: NP percentile, NP basic, NP bootstrap-t, NP BCa, P percentile, P basic and, as a reference, the classical Student t.
3. Estimated coverage = proportion of the \(M\) samples whose interval contains \(\mu\). Its Monte Carlo error is \(\sqrt{0.95 \times 0.05/2000} \approx 0.5\) pp.

The parametric model is always the exponential. Two scenarios: in the first, the data are exponential (model right); in the second, they come from a lognormal(0, 1) — also positive and skewed, plausible for times, but with coefficient of variation \(\sqrt{e - 1} \approx 1.31\) instead of the exponential's 1 (model wrong).

Each cell: coverage · mean interval width (in units of the true mean). An interval is only good if it covers close to 95% — a narrower interval with low coverage is not a virtue.

What the numbers say:

1. With the model right and small \(n\), the parametric wins. With \(n = 10\), the P percentile covers 91.3%, against 86.5% for the NP percentile and 89.8% for Student t. The NP underestimates variability because small samples of skewed data rarely show the tail; the P recovers it through the model. From \(n \approx 50\) on the differences vanish: both converge to 95%.
2. The interval method matters as much as the NP × P choice. The NP bootstrap-t is the best in both scenarios (94.4% with \(n = 10\) when the model is right, 90.0% when it is wrong), exactly as second-order theory predicts. BCa improves on the percentile, but only a little with \(n\) this small. The basic interval is the worst: it flips the skewness (8) to the wrong side.
3. With the model wrong, the parametric converges to the wrong answer. We can compute where to. Under the exponential model, the bootstrap standard error of the mean is \(\bar x/\sqrt n \to \mu/\sqrt n\), but the true one is \(\sigma/\sqrt n = 1.31\,\mu/\sqrt n\). The interval is 1.31 times narrower than it should be, and its coverage tends to

no matter how much data you have. The simulation confirms it: 88.5% with \(n = 50\), 87.8% with 100, 85.0% with 250 — already at the limit. Meanwhile the NP methods climb with \(n\) into the 92–93% range (the percentile goes from 80.0% to 92.5%). This is the heart of the comparison: the NP error shrinks with \(n\); the error of a wrong model does not. More data only make the P more confident in the wrong answer.

4. With \(n = 10\) and very skewed data, nobody is good. With the wrong model, the P percentile (86.7%) still beats the NP percentile (80.0%) — with 10 observations the sample almost never contains the heavy tail of the lognormal, and a slightly wrong model beats no model. The bootstrap does not create information the sample does not have.

§11Where the non-parametric bootstrap breaks: the maximum

Let \(X_i \sim \text{Uniform}(0, \theta)\) and \(\hat\theta = X_{(n)} = \max_i X_i\), the maximum likelihood estimator of \(\theta\). Its distribution is exact: \(\Pr(X_{(n)} \le t) = (t/\theta)^n\) for \(0 \le t \le \theta\), with bias \(-\theta/(n+1)\). Since \(\hat\theta < \theta\) always, the sign of the error is known.

In the NP bootstrap, a resample contains the observed maximum whenever it is drawn at least once:

The true distribution of \(n(\theta - \hat\theta)/\theta\) converges to an Exponential(1), which is continuous. The bootstrap distribution of \(n(\hat\theta - \hat\theta^*)/\hat\theta\) has an atom of mass 63.2% at zero, for every \(n\). The NP bootstrap is inconsistent here: more data do not help. The reason is structural: the maximum depends on the right tail of \(F\), and \(F_n\) has no information beyond the largest observed point.

The parametric bootstrap solves it elegantly. Under \(\text{U}(0, \hat\theta)\), the ratio \(\hat\theta^*/\hat\theta\) has a \(\text{Beta}(n, 1)\) distribution — exactly the same as \(\hat\theta/\theta\) in the real world. That is a pivot: a quantity whose distribution does not depend on the parameter. Using it as in (8), but with a ratio instead of a difference,

For the sample in the figure (\(\hat\theta = 0.965\), \(n = 30\)): [0.966, 1.092], which contains \(\theta = 1\). The coverage of each construction over 2,000 repetitions:

2,000 samples from U(0, 1) with n = 30, B = 999.

Two general lessons. First: the percentile interval has zero coverage — in both bootstraps — because every replicate lies below \(\hat\theta\), which lies below \(\theta\). Strong bias destroys the percentile interval. Second: the NP bootstrap requires the statistic to be a “smooth” function of \(F\). It fails (or needs corrections) for maxima and minima, extreme quantiles, parameters on the boundary of the parameter space, means of distributions with infinite variance and statistics that are not differentiable at critical points. A common fix is the m-out-of-n bootstrap, which resamples \(m \ll n\) observations.

§12Why it works: the theory on one page

Consistency. The bootstrap is consistent if the conditional distribution of \(\sqrt n(\hat\theta^* - \hat\theta)\) and the distribution of \(\sqrt n(\hat\theta - \theta)\) get close as \(n \to \infty\) (for example, the Kolmogorov distance between them goes to zero in probability). For the mean this holds whenever \(\operatorname{Var}(X) < \infty\) (Bickel and Freedman, 1981): by the CLT applied in the bootstrap world, \(\sqrt n(\bar x^* - \bar x) \Rightarrow \mathcal N(0, \hat\sigma^2)\), and \(\hat\sigma^2 \to \sigma^2\). For statistics that are smooth (Hadamard-differentiable) functionals of \(F\), the functional delta method extends the result. The maximum of §11 is not smooth, which is why it fails.

Why the bootstrap can beat the normal approximation. The Edgeworth expansion of the standardised mean is

where \(\gamma\) is the skewness of \(F\). The normal approximation ignores the second term and errs by \(O(n^{-1/2})\). The bootstrap reproduces the same expansion with \(\hat\gamma\) in place of \(\gamma\); since \(\hat\gamma - \gamma = O_p(n^{-1/2})\), the difference between the two expansions is \(O_p(n^{-1})\). This holds in full only for pivotal (studentised) statistics, like \(T\) in the bootstrap-t: for \(\bar X - \mu\) without dividing by the standard error, the gap between \(\hat\sigma\) and \(\sigma\) reintroduces an error of order \(n^{-1/2}\). That is why the bootstrap-t was the best in the simulation.

The parametric case. If the model is true and \(\hat\eta\) is consistent, \(F_{\hat\eta} \to F_\eta\), and the bootstrap distribution converges to the right one under regularity conditions (continuity of \(G_{F_\eta}\) in \(\eta\)). If the model is false, \(\hat\eta\) converges to the pseudo-true parameter \(\eta^\dagger\), the one minimising the Kullback–Leibler divergence to \(F\), and the bootstrap converges to \(G_{F_{\eta^\dagger}}\) — the sampling distribution of a world that does not exist. That is exactly what (13) quantified.

§13Side by side

Between the two extremes there are intermediate options: the residual bootstrap in regression (the model for the mean is parametric, the errors are resampled without assuming a distribution), the wild bootstrap (for heteroscedasticity) and the smoothed bootstrap (resampling from a kernel density estimate instead of \(F_n\), useful precisely for quantiles and medians).

§14Common mistakes

• Thinking the bootstrap creates data. It estimates variability from what the sample shows. If the sample is small or unrepresentative, the bootstrap inherits that (§7, §10).
• Resampling the wrong unit. The NP assumes i.i.d. observations. With clustered data (students in schools, repeated measures on the same patient), resample whole clusters. With time series, use a block bootstrap or a model for the dependence. Resampling dependent observations underestimates the standard error.
• Not redoing the whole procedure inside the loop. If the analysis includes outlier removal, variable selection, hyperparameter tuning or parameter estimation, all of it goes inside \(s(\cdot)\). Otherwise the variability of those steps disappears from the standard error.
• Using the percentile interval with a strongly biased estimator (§11). Look at the estimated bias first: if \(|\widehat{\operatorname{bias}}|/\widehat{\operatorname{se}}\) is large, prefer BCa, bootstrap-t or a pivot.
• Running a parametric bootstrap without checking the model. The result is only as good as the model. Compare the ECDF with the fitted model, make a QQ-plot and, if possible, run the NP alongside: if the two disagree a lot, investigate before reporting.
• Small \(B\) for intervals. \(B = 100\) gives a reasonable standard error but unstable interval limits (§8).
• Confusing the standard deviation of the data with that of the replicates. The standard error is the standard deviation of \(\hat\theta^*\), not of the \(x_i\).
• Statistics without finite variance. The mean of very heavy-tailed data (e.g. Cauchy, Pareto with index \(\le 2\)): the bootstrap of the mean does not converge to anything useful.

§15Summary

1. The goal of the bootstrap is to approximate the sampling distribution of an estimator — and from it get the standard error, bias and intervals — when there is no formula.
2. The principle: replace the population \(F\) by an estimate \(\hat F\) and simulate. There is a statistical error (\(\hat F \ne F\)) and a Monte Carlo error (finite \(B\)); only the second shrinks with \(B\).
3. Non-parametric: \(\hat F = F_n\), resample the data with replacement. Almost no assumptions; needs a reasonable \(n\) and smooth statistics.
4. Parametric: \(\hat F = F_{\hat\eta}\), generate from the fitted model. More accurate with small \(n\) and for tail statistics if the model is right; converges to the wrong answer if the model is wrong.
5. The interval type matters: percentile and basic are first-order; bootstrap-t and BCa second-order. With biased estimators the percentile can fail completely.
6. When in doubt, run both. If they agree, good. If they disagree, the disagreement is the most valuable information in the analysis.

§16Further reading

• B. Efron (1979). Bootstrap methods: another look at the jackknife. Annals of Statistics, 7(1). The original paper.
• B. Efron and R. Tibshirani (1993). An Introduction to the Bootstrap. Chapman & Hall. The textbook reference.
• A. C. Davison and D. V. Hinkley (1997). Bootstrap Methods and their Application. Cambridge. Practical and thorough; the basis of R's boot package.
• P. Hall (1992). The Bootstrap and Edgeworth Expansion. Springer. The second-order theory.
• P. J. Bickel and D. A. Freedman (1981). Some asymptotic theory for the bootstrap. Annals of Statistics, 9(6).
• L. Wasserman (2004). All of Statistics, ch. 8. A good summary in a few pages.