Population F
The true density. In practice, out of reach.
The two worlds
On the left, what nature does and you never see: a population F generates a sample, and the sample generates an estimate. On the right, what the bootstrap does and you see in full: the observed sample becomes the population, and the whole thing repeats. The wager is that the right-hand world imitates the left-hand one.
Real world
F → X₁…Xₙ → θ̂One sample of size n
The only object you actually hold in your hands.
Sampling distribution of θ̂
Repeating the collection 4,000 times from F. It only exists because this is a simulation.
Bootstrap world
F̂ₙ → X*₁…X*ₙ → θ̂*Stand-in population F̂ₙ
Mass 1/n on each observed point. This is the population the bootstrap draws from.
One resample with replacement
Stack height = how many times that point was drawn.
Bootstrap distribution of θ̂*
Resampling the sample B times. This one you can compute with real data.
The bootstrap distribution is the result of two substitutions chained together. The first is the plug-in principle: the parameter is a functional of the distribution, θ = T(F), and you estimate it by applying the same functional to the empirical distribution, θ̂ = T(F̂ₙ). The second is the bootstrap substitution proper: the distribution of θ̂ − θ under F is approximated by the distribution of θ̂* − θ̂ under F̂ₙ.
What holds up the first step is Glivenko–Cantelli: sup|F̂ₙ(x) − F(x)| → 0 almost surely. The grey band in the F̂ₙ panel is the 95% Dvoretzky–Kiefer–Wolfowitz bound, of half-width √(ln(2/α)/2n): with probability 0.95 the entire true CDF fits inside it. Notice that it shrinks like n^(−1/2) — drag n and watch. But F̂ₙ being close to F only translates into closeness of the distributions of θ̂ when the functional T is smooth enough. When it is not — the maximum being the classic counterexample — the chain snaps at the second link, not the first.
The comparison that matters
The bootstrap does not try to imitate the distribution of θ̂, but that of the error θ̂ − θ. That is why the two curves below are shown centred: the true one at θ, the bootstrap one at θ̂. The more they overlap, the better any confidence interval built from the magenta one will work.
Estimate on the sample θ̂—
True parameter θ = T(F)—
Bootstrap standard error sê*—
True standard error—
Bootstrap bias θ̄* − θ̂—
True bias E[θ̂] − θ—
Skewness of the bootstrap distribution—
Kolmogorov distance between the two—
The two errors advanced
Every bootstrap number carries two errors of opposite natures, and conflating them is the source of most misuse of the method. Raising B kills one of them and leaves the other untouched.
Bootstrap standard error computed from the first b resamples, for b = 1…B (log scale). The solid line is the true target; the band is the expected Monte Carlo noise.
Monte Carlo error in sê* — vanishes with B—
Approximation error — does not vanish with B—
Ratio of the two—
The Monte Carlo error comes from having run B resamples instead of all nⁿ possible ones. For the standard deviation it is roughly sê*/√(2(B−1)) and goes to zero at rate B^(−1/2). It is cheap to kill: raise B.
The approximation error comes from F̂ₙ not being F. It is governed by n, not by B, and no amount of compute removes it. Running a million resamples on 12 observations still delivers an answer about 12 observations — with three decimal places of precision on the wrong quantity.
Rule of thumb: pick B large enough that the first error sits well below the second. Once the ratio above drops under 0.1, raising B has stopped buying you precision.
Anatomy of a resample advanced
Distinct points in this resample—
Distinct proportion—
Expected value 1 − (1 − 1/n)ⁿ—
Limit as n → ∞0.6321
Largest multiplicity observed—
θ̂* of this resample—
Each resample leaves out, on average, about 37% of the original observations. The probability that a given point escapes all n draws is (1 − 1/n)ⁿ → e^(−1) ≈ 0.3679, which is where the complementary fraction 0.632 comes from — the one that names the .632 estimator of prediction error in machine learning.
The multiplicity vector (M₁,…,Mₙ) is Multinomial(n; 1/n,…,1/n), and it helps to read a resample as a vector of random weights Mᵢ/n with mean 1/n applied to fixed data. That reading is the door into the Bayesian bootstrap, which swaps the multinomial weights for Dirichlet ones, and into the wild bootstrap, which swaps them for continuous weights applied to residuals. Both arrive in stage 5.
The confidence interval ladder
Six ways to turn one bootstrap distribution into an interval. They differ in which quantiles they reach for, and the differences only matter when the distribution is skewed. The bars sit under the histogram they were computed from; the tick marks on the histogram show the quantiles each percentile-type method actually pulled.
Normal throws away the shape of the bootstrap distribution and keeps only its spread. Basic keeps the shape but reflects it: it treats θ̂* − θ̂ as a pivot, so a bootstrap distribution with a long right tail produces an interval with a long left arm. Percentile does the opposite, inheriting the tail directly. The two disagree by exactly twice the distance from θ̂ to the bootstrap median, which is why they coincide when the distribution is symmetric and centred.
BC adds one number, ẑ₀ = Φ⁻¹(Ĝ(θ̂)), which measures how far off centre θ̂ sits in its own bootstrap distribution, and shifts the quantile levels to compensate. BCa adds a second, the acceleration â, estimated from the jackknife, which corrects for the standard error changing with θ. Bootstrap-t takes a different route entirely: it studentizes each replicate, t* = (θ̂* − θ̂)/sê*, and reads the quantiles of that, which requires a variance estimate inside every resample.
Both BCa and bootstrap-t are second-order accurate: their one-sided coverage error is O(n⁻¹) rather than O(n⁻¹⁄²). The coverage laboratory below measures that claim, and the panel after it explains where the extra order comes from.
Coverage laboratory
An interval is only as good as the fraction of the time it actually contains θ. Here the whole procedure is repeated on fresh datasets drawn from F, and the two tails are counted separately — because an interval can hit its nominal coverage while missing badly on one side and compensating on the other.
Nothing measured yet. Coverage has to be simulated, and the run takes a few seconds.
Read the chart as two separate failures. The left arm is P(θ < lower), the right arm is P(θ > upper), and each should equal α/2 — the dashed marks. A method whose arms are lopsided is systematically placing the interval on the wrong side of θ̂, and the total coverage hides it: 2% on one side and 8% on the other still adds up to a 90% interval that is wrong in a way that matters for every one-sided claim you make from it.
With a skewed population and the mean, the percentile and normal intervals will show exactly this lopsidedness, while BCa and bootstrap-t straighten it out. That asymmetry is the visible signature of the n⁻¹⁄² Edgeworth term, and it is what the next panel is about.
How fast the error dies advanced
The same measurement repeated across a grid of sample sizes, plotting the one-sided error — how far the left-tail miss is from α/2 — on log–log axes. On these axes a rate becomes a slope, and the two families separate.
Not run yet. This one is heavy — five sample sizes, each a full coverage study.
The two grey guides have slopes −1/2 and −1. First-order methods should run parallel to the shallow one, second-order methods parallel to the steep one. Expect noise: with a few hundred datasets the measured error has a Monte Carlo uncertainty of its own, and at large n the true error can fall below that noise floor.
Where the extra order comes from advanced
The central limit theorem says the studentized mean is approximately normal. The Edgeworth expansion says how approximately, and the leading correction term is a polynomial in x times φ(x), scaled by n⁻¹⁄² and by the population skewness. Plotting CDFs would show nothing — they look identical — so what is plotted is the residual, the simulated CDF minus Φ.
Not simulated yet. This panel is written for the mean; pick Mean as the functional to use it.
Both correction polynomials are even in x, and φ is even too, so the whole n⁻¹⁄² term takes the same value at x and −x. Subtract the two tails to get a two-sided probability and the term cancels exactly. That single fact explains why a plain normal interval has one-sided coverage error O(n⁻¹⁄²) but two-sided error O(n⁻¹), and it is why two-sided coverage is a poor diagnostic: it flatters every method equally.
Note also the sign flip between the two panels. Studentizing does not remove the skewness effect, it reverses and roughly doubles it — the s in the denominator is itself correlated with x̄ when the population is skewed. This is why you cannot fix a skewness problem by switching from σ to s.
The correction the bootstrap finds on its own advanced
Inverting the Edgeworth expansion gives the Cornish–Fisher expansion for the quantiles: tᵕ ≈ zᵕ − (γ/6√n)(2zᵕ² + 1). Plotted below is the correction itself, quantile − z, so that a normal approximation is the flat line at zero. The point of the panel is the magenta curve: the bootstrap reproduces the correction without ever being told what γ is.
Run the simulation in the panel above to populate this one.
Population skewness γ—
Skewness correction scale γ̂/(6√n)—
BCa bias correction ẑ₀—
BCa acceleration â—
For the mean in the smooth-function model, both BCa constants are asymptotically the same quantity: ẑ₀ ≈ â ≈ γ̂/(6√n). The three numbers above should sit in the same neighbourhood, and they drift apart only through Monte Carlo noise and terms of order n⁻¹. That coincidence is not a curiosity — it is the reason BCa is second-order accurate. The two corrections it applies are precisely the two things the Cornish–Fisher expansion says are missing from the naive percentile interval.
The delta method, for comparison advanced
Before resampling was cheap, standard errors came from a Taylor expansion: Var(g(θ̂)) ≈ g′(θ)²σ²/n. Where a closed form exists it is free and exact to first order, and the bootstrap should agree with it. Where it disagrees, one of them is telling you something.
Not run yet.
All three curves should decay like n⁻¹⁄², appearing as parallel straight lines on log–log axes. The delta method curve is computed from a single sample, so it carries that sample's noise; the bootstrap curve is computed from the same sample and should track it closely for smooth functionals. Systematic separation at small n is the Taylor remainder making itself felt — the delta method is an asymptotic statement, and it is exactly at small n that you wanted a standard error in the first place.
The delta method is available here for the mean, variance, standard deviation and coefficient of variation. For the median, quantiles and trimmed means there is a delta-type formula too, but it involves the population density at the quantile — a quantity that is itself hard to estimate, and whose estimation error swamps the gain. That difficulty is the practical reason the bootstrap took over.