Both bootstraps place a random probability distribution on the observed data and evaluate a statistic on it.
- Efron: resample with replacement; weights lie on a discrete lattice.
- Bayesian bootstrap (Rubin): draw Dirichlet weights on the simplex; weights are continuous.
Core claim:
A statistic is \(T(F)\). Both methods randomise the empirical distribution \(\hat F_n\) via weights \(w\). They differ only in the law of \(w\).
Derivations below are from first principles unless noted. Every closed form in prose is checked against simulation at build time.
Constraints
- Not investment advice.
- Return series are illustrative; monthly equity returns are not i.i.d. (volatility clusters). Both bootstraps assume i.i.d. sampling anyway.
- See Non-weightable estimators for dependent-data fixes.
1 Statistics as functionals
A statistic is a functional \(T(F)\) of the unknown data-generating distribution \(F\), not a property of one sample alone.
| you say | you mean |
|---|---|
| mean | \(T(F) = \int x \, dF(x)\) |
| median | \(T(F) = \inf\{x : F(x) \ge 1/2\}\) |
| mode | \(T(F) = \arg\max_x f(x)\) |
| maximum | \(T(F) = \sup\{x : F(x) < 1\}\) |
| minimum | \(T(F) = \inf\{x : F(x) > 0\}\) |
| tail risk | \(T(F) = F(-5\%)\) |
Empirical cdf with equal weights:
\[\hat F_n = \frac{1}{n}\sum_{i=1}^{n} \delta_{x_i}\]
Plug-in estimate: \(T(\hat F_n)\). Bootstrap question: how does \(T\) vary over plausible perturbations of \(\hat F_n\)?
1.1 Datasets
Four datasets recur below. Toy samples support exact enumeration; real series support reporting decisions.
| dataset | \(n\) | the functional \(T(F)\) | what the answer is for |
|---|---|---|---|
| toy sample \(x = (1,2,3,4,10)\) | 5 | mean, median, max, min, mode | exhaustive verification |
| three points \(x = (1,3,10)\) | 3 | mean, median | drawing the geometry |
| S&P 500 monthly returns | 59 | mean, median, \(F(-5\%)\) | market risk reporting |
| IBM monthly returns | 59 | the same, plus the maximum | is this stock different? |
| MNIST test set | 10,000 | accuracy, per-class recall | ship / don’t ship |
Toy sample \(x=(1,2,3,4,10)\): \(\bar x=4\), \(\tilde s^2=10\); \(n=5\) gives 126 resamples — all enumerated.
Three-point sample \(x=(1,3,10)\): simplex plots use \(n=3\).
Both toy sets are hand-picked for exact checks, not realism.
Returns and MNIST are introduced in Real return data and MNIST example.
2 Efron’s bootstrap
Resampling with replacement yields counts \(N \sim \text{Multinomial}(n; 1/n,\ldots,1/n)\).
\[N = (N_1, \dots, N_n) \sim \text{Multinomial}(n;\ 1/n, \dots, 1/n).\]
Now write down the empirical distribution of the resample. The value \(x_i\) appears \(N_i\) times out of \(n\), so it carries mass \(N_i / n\):
\[F^*_n = \sum_{i=1}^n \frac{N_i}{n}\,\delta_{x_i} = \sum_{i=1}^n w_i\,\delta_{x_i}, \qquad w = N/n.\]
The resample empirical cdf equals \(\sum_i (N_i/n)\,\delta_{x_i}\). Weights \(w=N/n\) lie on a lattice.
The moments follow from the multinomial. With \(p_i = 1/n\):
\[E[N_i] = n p_i = 1, \qquad \operatorname{Var}(N_i) = n p_i (1 - p_i) = 1 - \tfrac1n, \qquad \operatorname{Cov}(N_i, N_j) = -n p_i p_j = -\tfrac1n.\]
Divide by \(n\) to get the weights, remembering that variances pick up \(1/n^2\):
\[E[w_i] = \frac{1}{n}, \qquad \operatorname{Var}(w_i) = \frac{n-1}{n^3}, \qquad \operatorname{Cov}(w_i, w_j) = -\frac{1}{n^3}.\]
The covariance is negative, and it has to be. The weights sum to one, so if one goes up another must come down. That is not a quirk of the multinomial; it is a consequence of the constraint, and it will reappear for the Bayesian weights with the identical sign and a different scale.
How many distinct weight vectors are there? A weight vector is determined by the counts \((N_1, \dots, N_n)\), non-negative integers summing to \(n\). That is the classic stars-and-bars count: lay out \(n\) stars and \(n-1\) bars in a row, and each arrangement encodes one allocation. There are
\[\binom{2n-1}{n-1}\]
such arrangements.
At \(n = 5\) that is \(\binom{9}{4} = 126\) patterns, and enumerating them with their exact multinomial probabilities gives a total of 1.000000000000 — the enumeration is complete. At \(n = 59\), the size of the returns window used later, it is \(\binom{117}{58} \approx 1.218 \times 10^{34}\).
So Efron’s bootstrap draws from an astronomically large but strictly finite set of weight vectors. The weights are confined to a lattice inside some space. Naming that space is the next step, and it turns out to be the whole story.
3 The weight simplex
Weights satisfy \(w_i\ge 0\) and \(\sum_i w_i=1\). The feasible set is the standard simplex:
\[\Delta^{n-1} = \Big\{ w \in \mathbb{R}^n : w_i \ge 0,\ \sum_{i=1}^n w_i = 1 \Big\}.\]
The superscript is \(n-1\) because one equation removes one degree of freedom: fix any \(n-1\) of the weights and the last is determined.
Four features of this object carry the rest of the post.
The corners are degenerate distributions. At corner \(i\), \(w_i = 1\) and everything else is zero: all your mass on a single observation.
The edges and faces are deletions. A point lies on a face exactly when some \(w_i = 0\) — that observation has been removed from the distribution entirely. Remember this; it explains almost every failure discussed later.
The centre is the plug-in. \(w = (1/n, \dots, 1/n)\) is \(\hat F_n\) itself.
The volume is \(1/(n-1)!\), so the uniform density on the simplex is the constant \((n-1)!\). At \(n = 59\) the simplex is a \(58\)-dimensional polytope of volume \(4.254 \times 10^{-79}\), and Efron’s \(1.2 \times 10^{34}\) lattice points sit inside it.
The most important thing to internalise is that a point in this triangle is not a summary of a distribution. It is a distribution.
The level sets are straight because the mean \(\sum_i w_i x_i\) is a linear function of \(w\). Any functional that is linear in the weights will have flat level sets like these; the median’s will not be, and that difference is exactly what makes the median harder.
Now that the space is named, put a distribution on it.
4 Dirichlet posterior
Bayesian update on the simplex, with Dirichlet appearing from conjugacy.
Assumption
\(F\) is discrete and supported on the \(n\) observed values. The unknown parameter is then the probability vector \(p = (p_1, \dots, p_n)\) attached to those values, and the parameter space is exactly the simplex \(\Delta^{n-1}\).
This is a strong assumption, and both bootstraps make it. The Bayesian version is merely the one that forces you to write it down.
Step 1 — the likelihood. If the data are \(n\) draws from the distribution with masses \(p\), and observation \(i\) was seen \(n_i\) times, the likelihood is multinomial:
\[L(p) \propto \prod_{i=1}^n p_i^{\,n_i}.\]
With distinct data — no exact ties, which is the normal case for returns — every \(n_i = 1\), so
\[L(p) \propto \prod_{i=1}^n p_i .\]
Step 2 — a prior of matching shape. Try a prior of the form \(\pi(p) \propto \prod_i p_i^{\alpha_i - 1}\) on the simplex. This is the Dirichlet density, but nothing yet depends on knowing its name.
Step 3 — multiply. Posterior \(\propto\) prior \(\times\) likelihood:
\[\pi(p \mid \text{data}) \;\propto\; \prod_i p_i^{\alpha_i - 1} \cdot \prod_i p_i^{\,n_i} \;=\; \prod_i p_i^{\alpha_i + n_i - 1}.\]
The exponents simply increment. Prior and posterior have the same shape, with \(\alpha_i \mapsto \alpha_i + n_i\) — which is what conjugacy means, shown here rather than asserted.
Step 4 — take the limit. Set every \(\alpha_i = a\) and let \(a \to 0\). The prior \(\propto \prod_i p_i^{-1}\) is improper: it does not integrate to a finite number, so it is not a distribution at all, only a formula you can still multiply by. With \(n_i = 1\) the posterior is \(\propto \prod_i p_i^{a}\). Match that against the Dirichlet’s \(\prod_i p_i^{\alpha_i - 1}\) and the posterior parameters are \(\alpha_i = a + 1\), which go to 1 as \(a \to 0\). So the posterior is
\[p \mid \text{data} \;\sim\; \text{Dirichlet}(1, 1, \dots, 1),\]
whose density \(\propto \prod_i p_i^{0} = 1\) is constant: the uniform distribution on the simplex. That is the Bayesian bootstrap.
Posterior versus prior
- Prior (Haldane): \(\pi(p)\propto\prod_i p_i^{-1}\) — improper, edge-heavy.
- Posterior after one observation per category: \(\text{Dirichlet}(1,\ldots,1)\) — uniform on the simplex.
- “Uniform prior on the simplex” misstates the method.
Two sanity checks make the construction believable.
First, \(E[w_i \mid \text{data}] = 1/n\) by symmetry, so the posterior mean of \(F\) is exactly \(\hat F_n\). The method is centred on the plug-in estimate — it does not pull your answer anywhere.
Second, \(a \to \infty\) collapses the posterior onto the single point \((1/n, \dots, 1/n)\): infinite prior confidence that all masses are equal. So \(a\) is a smoothing knob, \(\infty\) is “I already know the answer”, and \(0\) is the let-the-data-speak end.
4.1 Why the Dirichlet is forced, not merely convenient
Here is the property that makes this family the only sensible choice. Suppose you merge two of the categories — you stop distinguishing observation 1 from observation 2 and ask only about the combined mass \(p_1 + p_2\). A coherent theory must give the same answer whether you bin the sample space coarsely from the start or merge afterwards.
The Dirichlet does. Merging coordinates keeps you inside the family with the parameters added:
\[(p_1 + p_2, p_3, \dots, p_n) \sim \text{Dirichlet}(\alpha_1 + \alpha_2, \alpha_3, \dots, \alpha_n).\]
Iterating this to a set \(S\) of size \(k\) with all \(\alpha_i = 1\) collapses everything to two categories and gives
\[\sum_{i \in S} w_i \;\sim\; \text{Beta}(k,\ n-k).\]
This one line is worth more than any other in the post. Accuracy, tail probabilities, per-class recall, any proportion at all — all of them are sums of subsets of the weights, so all of them have an exact Beta posterior with no simulation whatsoever. It is used repeatedly below.
Pushing the same consistency requirement to a continuous sample space gives Ferguson’s Dirichlet process, and the Bayesian bootstrap is its \(\alpha \to 0\) posterior.
5 Exponential weight construction
Normalized independent \(\text{Gamma}(1,1)=\text{Exp}(1)\) draws give \(\text{Dirichlet}(1,\ldots,1)\) weights.
Let \(G_1, \dots, G_n\) be independent with \(G_i \sim \text{Gamma}(\alpha_i, 1)\), so the joint density is
\[f(g) = \prod_{i=1}^n \frac{1}{\Gamma(\alpha_i)} g_i^{\alpha_i - 1} e^{-g_i}.\]
Change variables to the total and the direction: \(s = \sum_i g_i\) and \(w_i = g_i / s\). The inverse map is \(g_i = s\,w_i\). Because \(w\) has only \(n-1\) free coordinates, use \((w_1, \dots, w_{n-1}, s)\) as the new variables; the Jacobian of \(g_i = s w_i\) is
\[\left|\frac{\partial g}{\partial(w_1,\dots,w_{n-1},s)}\right| = s^{\,n-1}.\]
Substituting \(g_i = s w_i\) and \(\sum_i g_i = s\):
\[f(w, s) \;\propto\; \prod_i (s w_i)^{\alpha_i - 1} e^{-s} \cdot s^{n-1} \;=\; \underbrace{s^{\sum_i \alpha_i - 1} e^{-s}}_{\text{Gamma}(\sum \alpha_i,\,1) \text{ in } s} \;\times\; \underbrace{\prod_i w_i^{\alpha_i - 1}}_{\text{Dirichlet}(\alpha) \text{ in } w}.\]
The joint density factors into a function of \(s\) alone times a function of \(w\) alone, with no cross-term. That single factorisation delivers three facts at once: \(w \sim \text{Dirichlet}(\alpha)\), \(S \sim \text{Gamma}(\sum_i \alpha_i, 1)\), and \(w \perp S\).
Now set every \(\alpha_i = 1\). Since \(\text{Gamma}(1,1) = \text{Exp}(1)\), the recipe becomes: draw \(n\) exponentials, divide by their sum. The exponentials are not a trick; the normalisation costs nothing precisely because the direction carries no information about the length.
5.1 The same weights come from breaking a stick
Rubin’s original description is more vivid. Scatter \(n-1\) points uniformly on the interval \([0,1]\), sort them, and take the \(n\) gaps between consecutive points, counting the two ends. Those gaps — uniform spacings — are exactly \(\text{Dirichlet}(1,\dots,1)\).
The two constructions describe the same law: writing the sorted uniforms as running sums of scaled exponentials, the Rényi representation, turns one picture into the other. The spacings picture then hands you something the gamma picture hides. Because the gaps are laid out along the stick in order, their partial sums are the sorted uniforms themselves:
\[w_1 + \dots + w_j \;\stackrel{d}{=}\; U_{(j)},\]
the \(j\)-th order statistic of \(n-1\) uniforms. That fact does all the work for the median later.
Both constructions and the library agree, and normalised \(\text{Gamma}(4,1)\) variates give \(\text{Dirichlet}(4,\dots,4)\) as the factorisation predicts. Against the exact \(\text{Beta}(1, n-1)\) marginal at \(n = 6\), a Kolmogorov–Smirnov test gives \(p\) = 0.134 for the exponential construction and \(p\) = 0.735 for uniform spacings; the two samples are indistinguishable from each other (\(p\) = 0.520), and the \(\text{Gamma}(4,1)\) check against \(\text{Beta}(4, 4(n-1))\) gives \(p\) = 0.915.
5.2 Why the weights must sum to one
Three arguments, shallowest first.
1. It is the definition of the object. \(F_w\) is supposed to be a probability distribution. Total mass one is what makes it one. Without it, \(F_w(-5\%)\) is not a probability and there is nothing to interpret.
2. Otherwise the statistic stops behaving like a statistic. With \(\sum_i w_i = 1\), shifting every observation by \(c\) shifts the weighted mean by \(c\). Without the constraint, \(\sum_i w_i (x_i + c) = \sum_i w_i x_i + c\sum_i w_i\), so the answer scales with the total. And the total is random: \(\sum_i E_i \sim \text{Gamma}(n, 1)\) has standard deviation \(\sqrt n\) against a mean of \(n\), a relative fluctuation of \(1/\sqrt n\) — the same order as the standard error you are trying to measure. Skipping the normalisation roughly inflates your uncertainty by \(\sqrt 2\).
On the S&P window, dividing by \(n\) instead of by \(\sum_i E_i\) inflates the standard deviation of the bootstrapped mean by a factor of 1.12 — against the predicted \(\sqrt 2 \approx 1.41\).
3. The simplex is the parameter space. This is the real reason. The constraint is not something imposed on the posterior after the fact; it is the domain the posterior lives on. Sum-to-one and the negative correlation between weights are the same geometric fact seen from two sides.
5.3 Identical means, second moments off by a single factor
By symmetry \(E[w_i] = 1/n\) for \(\text{Dirichlet}(1,\dots,1)\). Merging gives \(w_i \sim \text{Beta}(1, n-1)\), whose variance is
\[\operatorname{Var}(w_i) = \frac{1 \cdot (n-1)}{(1 + n - 1)^2 (1 + n)} = \frac{n-1}{n^2 (n+1)},\]
and the constraint \(\operatorname{Var}(\sum_i w_i) = 0\) forces \(n \operatorname{Var}(w_i) + n(n-1)\operatorname{Cov}(w_i,w_j) = 0\), hence \(\operatorname{Cov}(w_i, w_j) = -1/[n^2(n+1)]\).
Put the two methods side by side.
Efron
\(E[w_i] = 1/n\)
\(\operatorname{Var}(w_i) = \dfrac{n-1}{n^3}\)
\(\operatorname{Cov}(w_i,w_j) = -\dfrac{1}{n^3}\)
\(P(w_i = 0) = (1-1/n)^n\)
Bayes
\(E[w_i] = 1/n\)
\(\operatorname{Var}(w_i) = \dfrac{n-1}{n^2(n+1)}\)
\(\operatorname{Cov}(w_i,w_j) = -\dfrac{1}{n^2(n+1)}\)
\(P(w_i = 0) = 0\)
Same mean, same sign of covariance, and every variance in the ratio
\[\frac{\operatorname{Var}_{\text{Bayes}}(w_i)}{\operatorname{Var}_{\text{Efron}}(w_i)} = \frac{n^3}{n^2(n+1)} = \frac{n}{n+1}.\]
Every one of those closed forms is checked against 400,000 simulated replicates in the cell above; the build fails if any disagrees.
6 Conjugate-pair duality
Efron and Bayes are opposite conditionings of one multinomial–Dirichlet pair.
\[\underbrace{\text{counts} \mid p = \hat F_n \;\sim\; \text{Multinomial}(n;\ 1/n, \dots, 1/n)}_{\textbf{Efron: fix the probabilities, randomise the counts}}\]
\[\underbrace{p \mid \text{counts} \;\sim\; \text{Dirichlet}(1, \dots, 1)}_{\textbf{Bayes: fix the counts, randomise the probabilities}}\]
Same pair of distributions, opposite arrows. Efron’s object is a sampling distribution: the randomness answers “what other datasets might I have drawn?” The Bayesian bootstrap’s is a posterior: the randomness answers “what might \(F\) be?” They produce nearly identical numbers and mean opposite things.
6.1 Does the resampling side secretly generate a posterior?
No — and the argument is short enough to give in full.
In this model every observed \(x_i\) contributes a factor \(p_i\) to the likelihood, so \(L(p) = 0\) whenever any \(p_i = 0\). Bayes’ theorem multiplies prior by likelihood. Therefore every posterior under every prior must assign probability zero to the faces of the simplex. You watched \(x_i\) happen; no amount of conditioning can leave you believing it was impossible.
But an Efron resample that omits an observation lands exactly on a face. And the share of resamples that omit at least one observation is \(1 - n!/n^n\).
At \(n=3\) it is 0.7778; at \(n=5\), 0.9616; at \(n=12\), 0.99995; at \(n=59\) it is 1 to within \(5\times10^{-25}\). Essentially every Efron replicate corresponds to a distribution that assigns probability zero to something you watched happen. At \(n = 59\) a typical resample omits about 21.5 of the 59 observations.
So Efron’s bootstrap is not a Bayesian procedure in disguise. It is a different kind of object that happens to give nearly the same numbers.
6.2 A second, independent obstruction
Take the prior \(\text{Dirichlet}(a, \dots, a)\), so the posterior is \(\text{Dirichlet}(a+1, \dots, a+1)\) and the symmetric-Dirichlet variance formula gives, for the mean functional,
\[\operatorname{Var}(\text{mean} \mid \text{data}) = \frac{\tilde s^2}{n a + n + 1}, \qquad \tilde s^2 = \frac1n \sum_i (x_i - \bar x)^2 .\]
This is largest at \(a = 0\). Every proper prior is tighter than the Bayesian bootstrap. To reach Efron’s \(\tilde s^2/n\) you would need \(na + n + 1 = n\), that is \(a = -1/n\) — a negative concentration, which is not a distribution.
| prior | \(a\) | denominator | sd |
|---|---|---|---|
| Haldane — the Bayesian bootstrap | 0 | 6 | 1.291 |
| Perks | 0.20 | 7 | 1.195 |
| Jeffreys | 0.50 | 8.5 | 1.085 |
| uniform | 1.00 | 11 | 0.953 |
| boundary of what a prior can do | |||
| Efron bootstrap | −1/n (invalid) | 5 | 1.414 |
| classical \(s^2/n\) | — | 4 | 1.581 |
6.3 What the Bayesian bootstrap’s prior actually believes
Lest this read as a free lunch. Under the stick-breaking construction of a Dirichlet process, the first weight is \(\text{Beta}(1, \alpha)\), which concentrates at 1 as \(\alpha \to 0\). So a prior draw from the Bayesian bootstrap’s prior is a point mass at a single random location: its opinion is that the world is deterministic and you have no idea where.
At \(\alpha = 0.01\) the largest atom averages 0.9932 and the effective number of atoms is 1.02; by \(\alpha = 20\) those are 0.1222 and 22.5. This single fact explains both the good behaviour and the bad: it is why the posterior mean is exactly \(\hat F_n\), and it is why the maximum degenerates.
6.4 Rubin’s actual point
Rubin’s 1981 paper is less an endorsement than a diagnosis. He built the Bayesian analogue of Efron’s method and then pointed at the prior it had forced on him: no unobserved value can ever occur. Efron’s bootstrap makes the identical assumption the moment it substitutes \(\hat F_n\) for \(F\) — it simply never has to say so out loud. The value of the Bayesian version is that it puts the assumption on the page where you can object to it.
If the two live on the same space with the same centre, what does the difference actually look like?
7 Lattice versus continuous weights
At \(n=3\), Efron has ten lattice points; Dirichlet\((1,1,1)\) has a density on the full triangle.
Look at one weight at a time. Efron’s marginal is a rescaled \(\text{Binomial}(n, 1/n)\) — a picket fence with a spike of height \((1-1/n)^n\) at zero. The Bayesian marginal is the continuous \(\text{Beta}(1, n-1)\).
7.1 One variance, three denominators
For the mean, do the algebra once. With weights having \(\operatorname{Var}(w_i) = v\) and \(\operatorname{Cov}(w_i,w_j) = c\),
\[\operatorname{Var}\Big(\sum_i w_i x_i\Big) = v \sum_i x_i^2 + c \sum_{i \ne j} x_i x_j = (v - c)\sum_i x_i^2 + c\Big(\sum_i x_i\Big)^2 .\]
For the symmetric Dirichlet with posterior concentration \(\alpha_0'\), \(v - c = \frac{1}{n(n\alpha_0'+1)}\) and \(c = -\frac{1}{n^2(n\alpha_0'+1)}\), so
\[\operatorname{Var}\Big(\sum_i w_i x_i\Big) = \frac{n\sum_i x_i^2 - (\sum_i x_i)^2}{n^2 (n\alpha_0' + 1)} = \frac{\tilde s^2}{n\alpha_0' + 1}.\]
Setting \(\alpha_0' = 1\) gives the Bayesian bootstrap; running the same computation with Efron’s \(v, c\) gives \(\tilde s^2 / n\). So all three familiar answers are the same quantity with a different fudge, and they are ordered:
\[\underbrace{\frac{\tilde s^2}{n+1}}_{\text{Bayes}} \;<\; \underbrace{\frac{\tilde s^2}{n}}_{\text{Efron}} \;<\; \underbrace{\frac{\tilde s^2}{n-1}}_{\text{textbook } s^2/n}\]
On the toy sample \(x = (1,2,3,4,10)\) with \(\tilde s^2 = 10\): variances 1.6667 < 2.0000 < 2.5000, standard deviations 1.2910 < 1.4142 < 1.5811, and a simulated check at \(B = 200\,000\) giving 1.2925 and 1.4147.
7.2 The asymptotic position, stated honestly
For smooth functionals both bootstraps converge to the same normal limit. That is Lo’s 1987 theorem for the Bayesian bootstrap, and Præstgaard and Wellner (1993) generalise it: the limit is the same for any weighting scheme that is exchangeable — meaning the joint law of the weights does not care which observation is which — and has the right first two moments. Efron’s multinomial weights, Dirichlet weights and Poisson weights are three members of one family.
This is the one result in the post I am not deriving — the proof is a functional central limit theorem and genuinely does not fit. What I can do is verify its consequence numerically, which the \(n = 437\) row of the sample-size table below does: by then the two standard errors agree to four significant figures.
Differences appear only at small \(n\), or for functionals the asymptotic theory does not cover. Which is the next question.
8 Mutual singularity
No \(\text{Dirichlet}(\alpha)\) with \(\alpha>0\) can equal Efron’s lattice law: continuous vs atomic measures are mutually singular.
8.1 But the first two moments match exactly
Set every \(\alpha_i = \alpha\) and equate the variances:
\[\frac{n-1}{n^2(n\alpha + 1)} = \frac{n-1}{n^3} \implies n\alpha + 1 = n \implies \alpha = 1 - \frac1n .\]
The covariance then matches automatically, because both families satisfy the same sum-to-one identity linking \(v\) and \(c\). So \(\text{Dirichlet}(1 - 1/n, \dots, 1 - 1/n)\) has the identical mean vector and covariance matrix to Efron’s weights, and for any linear functional the first two moments agree exactly on every dataset.
Note the value is below 1 — the edge-seeking regime. To imitate Efron you need a prior that likes deleting observations, which is conceptually exactly right. And \(\alpha \to 1\) as \(n \to \infty\).
8.2 Matching the variance does not match the shape
Two moments are not a distribution. Push past the variance and the agreement stops, which is easiest to see on the toy sample: run its bootstrapped mean through Efron’s weights and through Dirichlet weights at both \(\alpha = 0.8\) — the matching value at \(n=5\) — and \(\alpha = 1\), then compare not only the spread but the skewness and the excess kurtosis, which measure how lopsided and how heavy-tailed the bootstrap distribution is.
| weights | sd | skewness | excess kurtosis |
|---|---|---|---|
| Efron | 1.4128 | 0.5125 | -0.0300 |
| Dirichlet(0.8) | 1.4130 | 0.8486 | 0.4245 |
| Dirichlet(1) | 1.2930 | 0.8063 | 0.4172 |
The standard deviation matches Efron at \(\alpha \approx\) 0.80, exactly as the algebra predicted. The skewness matches at \(\alpha \approx\) 3.2. One parameter cannot do both jobs.
Worse, tuning to match the variance makes other things worse. The median law under \(\text{Dirichlet}(0.8)\) is [0.079, 0.247, 0.346, 0.248, 0.08], which is further from Efron’s [0.058, 0.26, 0.363, 0.26, 0.058] (total variation distance 0.0426) than plain \(\alpha = 1\) is (0.0207).
8.3 The construction that does contain both
Both methods are “draw i.i.d. non-negative weights from a unit-rate source, then fix up the total”.
| generator | fix-up | result | |
|---|---|---|---|
| Efron | \(W_i \sim \text{Poisson}(1)\) | condition on \(\sum_i W_i = n\) | exactly \(\text{Multinomial}(n; 1/n)/n\) |
| Bayes | \(W_i \sim \text{Exp}(1)\) | divide by \(\sum_i W_i\) | exactly \(\text{Dirichlet}(1,\dots,1)\) |
The Poisson claim is exact, and cheap to verify by brute force.
Enumerating every pattern for \(n = 3, 4, 5\) and comparing conditioned-Poisson probabilities against multinomial probabilities gives a maximum absolute discrepancy of 2.8e-17 — machine precision.
And now the punchline. \(\text{Exp}(1)\) and \(\text{Poisson}(1)\) are the interarrival times and the counts of the same unit-rate Poisson process. The chunky-versus-smooth distinction has been pushed all the way back to its source: Efron counts arrivals, the Bayesian bootstrap measures the gaps between them.
8.4 Two corners of one family
The Dirichlet-multinomial with resample size \(m\) and concentration \(\alpha\) contains Efron at \((m = n,\ \alpha \to \infty)\) and Bayes at \((m \to \infty,\ \alpha = 1)\). Both are corners; neither contains the other.
The naive “discretised Bayesian bootstrap” at \((m = n, \alpha = 1)\) is worth a warning: it is exactly uniform over all \(\binom{2n-1}{n-1}\) lattice patterns — verified above, all 35 patterns at \(n = 4\) carry probability 0.028571 \(= 1/35\) — and it has exactly twice the Bayesian bootstrap’s weight variance. It is not a good approximation to anything.
8.5 The reverse direction is cleaner
Run Efron’s resampling with a self-reinforcing urn. Draw an observation, then return it together with a copy — a Pólya urn — and let the resample size grow. The frequency vector converges to \(\text{Dirichlet}(1,\dots,1)\) exactly.
At \(n = 5\) with a resample of size 3,000, the marginal of a single frequency against the exact \(\text{Beta}(1, n-1)\) gives a Kolmogorov–Smirnov statistic \(D\) = 0.0170 with \(p\) = 0.197 — no evidence against the identity.
One line for the whole section: ordinary resampling forgets what it has drawn; the Bayesian bootstrap reinforces it.
9 Statistic choice
Bootstrap uncertainty is area on the simplex under Dirichlet\((1,\ldots,1)\) weights.
9.1 Mean — both work
Derived above: variances \(\tilde s^2/(n+1)\) against \(\tilde s^2/n\), in ratio \(n/(n+1)\). Smooth, linear, no surprises. This is the case the asymptotic theory covers and the case where the choice does not matter.
9.2 Median — an exact law, and a century-old interval falls out
First, a definitional point that is usually fudged. Define the median as \(F_w^{-1}(1/2) = \inf\{x : F_w(x) \ge 1/2\}\) for both methods. Without a shared definition the two are being scored on different statistics and the comparison is meaningless.
Now use the uniform-spacings construction. The partial sums of the weights are the order statistics of \(n-1\) uniforms, so the median sits at the rank
\[J = \min\{j : U_{(j)} \ge 1/2\} = 1 + \#\{i : U_i < 1/2\},\]
and each of the \(n-1\) uniforms falls below \(1/2\) independently with probability \(1/2\). Therefore
\[P(\text{median} = x_{(j)}) = \binom{n-1}{j-1} 2^{-(n-1)}.\]
A \(\text{Binomial}(n-1, \tfrac12)\) law over the ranks, completely free of the data values.
At \(n = 5\) the exact law is \((1,4,6,4,1)/16\) = [0.0625, 0.25, 0.375, 0.25, 0.0625], and \(200\,000\) simulated replicates give [0.0624, 0.2489, 0.3745, 0.2519, 0.0623].
Here is the payoff. At \(n = 59\), the central 95% of a \(\text{Binomial}(58, \tfrac12)\) runs from rank 23 to rank 37. So the Bayesian bootstrap’s 95% interval for the median is order statistics 23 and 37 — which is exactly the classical distribution-free sign-test interval for a median. A century-old nonparametric procedure drops out of the Dirichlet posterior with no approximation anywhere.
Quantile lumpiness
Continuous weights do not smooth quantiles. Both methods assume support on observed points only.
9.3 Max and min — mirror images, both useless
Under \(\text{Dirichlet}(1,\dots,1)\) every weight is strictly positive almost surely. So \(\sup \operatorname{supp}(F_w) = x_{(n)}\) with probability one: the posterior for the maximum is a point mass. No uncertainty at all, which is obviously wrong.
Efron is merely non-degenerate. The resample maximum is at most \(x_{(j)}\) exactly when no draw exceeded rank \(j\), and each of the \(n\) independent draws does so with probability \(j/n\):
\[P(\max{}^* \le x_{(j)}) = (j/n)^n \implies P(\max{}^* = x_{(j)}) = \Big(\frac{j}{n}\Big)^n - \Big(\frac{j-1}{n}\Big)^n .\]
The top atom is \(1 - (1-1/n)^n \to 1 - e^{-1} \approx 0.632\).
Simulation confirms \(P(\max^* \le x_{(j)}) = (j/n)^n\) to within 0.0004, and the top atom is 0.6723. Meanwhile 100% of Bayesian replicates return \(x_{(n)}\).
But notice that Efron’s distribution is also bounded above by \(x_{(n)}\). Neither bootstrap can say anything about an extremum, because both assume \(F\) gives zero probability to unseen values, and an extremum is entirely a statement about that region. The theory says the same thing: the bootstrap is consistent for functionals that respond smoothly to a small nudge in \(F\) — Hadamard-differentiable ones — and the maximum, which jumps the moment the largest point moves, is not one. For extremes the bootstrap is provably inconsistent.
There are three real fixes. Ask for a high quantile instead of the maximum. Or model the tail rather than resample it: keep only the observations above some high threshold and fit a generalised Pareto distribution to their excesses, the peaks-over-threshold approach, which can then say something about values you have not seen. Or stay Bayesian and use a full Dirichlet process prior with \(\alpha > 0\) and a continuous base measure \(F_0\) — a guess at where unseen values would fall — which gives them posterior predictive mass \(\alpha/(\alpha + n)\).
Demonstrating the third on IBM: with \(F_0\) a normal fitted to the returns (\(\mu\) = -0.222, \(\sigma\) = 5.027) and the observed maximum +15.584%, the probability that a fresh month beats the record goes from exactly \(0\) at \(\alpha = 0\), to \(1.386 \times 10^{-5}\) at \(\alpha = 1\), to \(6.499 \times 10^{-5}\) at \(\alpha = 5\). The degeneracy is a property of the prior, not of Bayesian inference.
9.4 Mode — here the estimator is broken, not the bootstrap
For a distribution supported on distinct observed values, the mode is \(\arg\max_i w_i\): the observation with the largest weight. By exchangeability, the Bayesian posterior over which observation that is must be uniform, \(1/n\) each — completely uninformative.
At \(n = 6\) the Bayesian law is [0.1664, 0.1666, 0.1665, 0.1669, 0.1673, 0.1663] against \(1/6\) = 0.1667.
Efron’s is worse in a more insidious way: it is dominated by ties. At \(n = 6\), 40.7% of resamples have a tied modal count, so the answer is decided by your tie-breaking rule rather than by the data. NumPy’s argmax breaks ties toward the first index, producing the spuriously non-uniform law [0.2519, 0.2048, 0.1741, 0.1485, 0.1223, 0.0983]. That is an artefact of the code, not a finding about the data; with random tie-breaking it is uniform by symmetry.
The lesson here is different from the maximum. There, the bootstrap was being asked an unanswerable question. Here, the estimator is broken: the mode of an empirical distribution over distinct points is meaningless. Fix the estimator — use a weighted kernel density estimate
\[\hat f_w(t) = \sum_i w_i K_h(t - x_i)\]
and take its argmax — and both bootstraps work.
On the IBM returns with a Silverman bandwidth of 2.36 percentage points, the KDE mode has a posterior standard deviation of 1.514 pp under the Bayesian bootstrap and 1.476 pp under Efron — both perfectly usable. Note that the Dirichlet weights slot into the weighted KDE with no modification at all, which is the first concrete instance of a theme that returns below.
9.5 The whole distribution — an exact Beta at every point
With \(k = \#\{x_i \le t\}\), the value \(F_w(t)\) is a sum of \(k\) of the weights. Aggregation gives it immediately:
\[F_w(t) \sim \text{Beta}(k,\ n-k) \quad \text{(Bayes)}, \qquad F^*(t) \sim \text{Bin}(n,\ k/n)/n \quad \text{(Efron)}.\]
Same mean \(k/n\); variances \(\frac{k(n-k)}{n^2(n+1)}\) and \(\frac{k(n-k)}{n^3}\), in ratio \(n/(n+1)\) once again.
The Beta is exactly the posterior you would get by treating “is the return below \(t\)?” as a coin flip and taking the same \(a \to 0\) limit in two categories rather than \(n\) — the Haldane prior. The construction agrees with itself, which is reassuring.
10 Real return data
Headline mean/median differences are tiny at \(n=59\); tail probabilities on sparse events differ materially.
10.1 Data provenance
The S&P 500 is a broad US equity index — roughly “the market” — and IBM is one large-cap stock, so the pair is the smallest interesting contrast: everything, against one thing.
The index figures come from the monthly S&P 500 series compiled by the economist Robert Shiller, who publishes it alongside his work on long-run market valuation; this post reads the datasets/s-and-p-500 repackaging of it. The IBM figures come from plotly’s five-year file of daily closes for S&P 500 constituents, assembled as a sample dataset for charting demos. Neither file was collected with a bootstrap in mind, which is rather the point: they are ordinary public data of the kind anyone actually has.
Both series are converted to monthly percentage changes, and both are cut to the same 59 months, March 2013 to January 2018, so that every comparison below is like-for-like rather than an artefact of one series covering a calmer period than the other.
One row of these datasets is one month. That matters: the resampling unit is a month, so everything here is a statement about month-to-month variation, and every interval below answers “if the next five years of months were drawn from the same process, how differently might this number have come out?”
10.2 Functionals reported
Four questions are asked of both series, and each is a functional of the return distribution \(F\):
| question | functional | why anyone asks |
|---|---|---|
| what does this earn on average? | \(T(F) = \int x\,dF\) | expected return, the input to any allocation |
| what does a typical month look like? | \(T(F) = F^{-1}(1/2)\) | a median is robust to one crash |
| how often is there a bad month? | \(T(F) = F(-5\%)\) | tail risk, drawdown budgeting |
| how good can a month get? | \(T(F) = \sup \operatorname{supp} F\) | best case — and, as it turns out, unanswerable |
Plus one two-sample question: is the mean of one series different from the mean of the other?
10.3 Decision stakes
Concretely, these are the numbers behind three decisions. How much of a portfolio to hold in equities rests on the mean and its uncertainty. How much loss to plan for rests on \(F(-5\%)\) — a risk limit, a capital buffer, a margin calculation. Whether a stock has beaten the market rests on the two-sample difference, and is the question behind every performance-attribution argument.
Getting any of those wrong is expensive in a specific direction: an interval that is too narrow understates how much of this is luck, and a lower endpoint of exactly zero on a tail probability says a \(-5\%\) month cannot happen, which is the one thing a capital buffer exists to survive.
That is also why a bootstrap rather than a formula. With 59 months there is no reliable normal approximation for a tail probability, the tail here rests on a single observed month, and the question is precisely what an interval does when the count underneath it is one. All three decisions are made from an interval rather than a point estimate, and the interval is exactly what the two methods disagree about — by nothing at all for the mean, and enormously for the tail.
Sample-size constraint
These series are here because they are small, public, and contain a genuine rare-event problem — not because five years of monthly data can settle an investment question. It cannot. The intervals below are wide, the tail estimate rests on a single month, and the i.i.d. assumption both methods make is false for returns. Read them as a worked example of how the uncertainty behaves, not as a finding about markets.
| series | n | mean % | sd % | median % | best | worst | months < -5% |
|---|---|---|---|---|---|---|---|
| S&P 500 | 59 | 1.068 | 2.244 | 1.268 | +6.171 (Mar 2016) | -6.596 (Jan 2016) | 1 |
| IBM | 59 | -0.222 | 5.027 | -0.471 | +15.584 (Mar 2016) | -12.165 (Oct 2014) | 10 |
Series definitions
The S&P figures come from the Shiller monthly series, in which each monthly level is the average of that month’s daily closes, not a month-end close. Averaging smooths, so this slightly understates monthly volatility. As a check on the download, the annualised volatility from Feb 1990 - Jun 2026 comes out at 12.3%.
The IBM figures come from unadjusted daily closes, so they are price returns: dividends are excluded, which understates total return by roughly 3–4%/yr over this window.
And, as stated at the top: these returns are not i.i.d., and both bootstraps assume they are.
| series | Bayes sd | Efron sd | closed Bayes | closed Efron | classical | ratio | Efron 95% | Bayes 95% |
|---|---|---|---|---|---|---|---|---|
| S&P 500 | 0.28732 | 0.28946 | 0.28723 | 0.28965 | 0.29214 | 0.99262 | [0.486, 1.625] | [0.481, 1.615] |
| IBM | 0.64373 | 0.64838 | 0.64340 | 0.64883 | 0.65440 | 0.99283 | [-1.490, 1.058] | [-1.479, 1.056] |
Lead with this: at \(n = 59\) the choice between the two methods changes your standard error in the third decimal place. The theoretical ratio is \(\sqrt{59/60}\) = 0.9916, and the simulated ratios are 0.9926 and 0.9928.
Monte Carlo precision
At \(B = 200\,000\), the relative Monte-Carlo standard error of a simulated standard deviation is about \(1/\sqrt{2B}\) = 0.00158, so a ratio of two simulated sds carries roughly ±0.0045. Several Efron-versus-Bayes comparisons in this post differ by less than that. Where they do, they are ties, and this post reports them as ties.
So spend the section on the two places where the choice does matter.
10.4 Where it matters, case 1 — a probability resting on one month
The question is \(F(-5\%)\): how often a month loses more than five per cent. Both series get the identical treatment, but they are in very different positions to answer it. IBM had ten such months in this window; the S&P had one. A number resting on a single observation is where the deletion property stops being a curiosity and starts changing the answer.
| series | k | plug-in | Bayes sd | Beta sd | Efron sd | Beta 95% | Efron 95% | P(Efron = 0) |
|---|---|---|---|---|---|---|---|---|
| S&P 500 | 1 | 0.01695 | 0.01671 | 0.01666 | 0.01683 | [0.00044, 0.0616] | [0.00000, 0.0508] | 0.36602 |
| IBM | 10 | 0.16949 | 0.04848 | 0.04844 | 0.04896 | [0.08590, 0.2742] | [0.08475, 0.2712] | 0.00003 |
IBM has \(k = 10\) such months and the two methods agree to within Monte-Carlo noise. The S&P has \(k = 1\) — a single month, January 2016 — and
\[P(\text{that month is omitted}) = \Big(1 - \tfrac{1}{59}\Big)^{59} = 0.365 .\]
So 36.5% of the time Efron’s answer is exactly zero: “a \(-5\%\) month is impossible.” Its 95% interval has a hard zero as its lower endpoint, which is a statement about the resampling scheme rather than about markets. The Bayesian bootstrap never says that. It returns the exact \(\text{Beta}(1, 58)\) posterior, with a lower endpoint of 0.04%.
This generalises, and it is the practical heart of the post: when a statistic rests on a handful of observations, Efron spends a large share of its replicates in a state where those observations do not exist. The same arithmetic reappears on MNIST’s rare classes below.
10.5 Where it matters, case 2 — one enormous month
The maximum was shown above to be unanswerable by either method. IBM’s window makes that concrete, because it contains one month far larger than any other — and a single observation like that decides an extremum on its own.
IBM’s best month, +15.58% in March 2016, is more than twice its second best, +6.72%. Efron reports a bimodal law with 63.4% of its mass on 15.58 and most of the rest near 6.72, and a 95% lower endpoint of +6.21% that is an artefact of the resampling scheme rather than a statement about markets. Bayes reports certainty. Both are useless, in the two different ways derived earlier.
10.6 A two-sample difference needs no new theory
Both failures above were about one series at a time. The question people actually argue over compares two, and it costs nothing extra: because the weight vectors for the two series are drawn independently, the difference of means is just the difference of two bootstrapped means, replicate by replicate.
| method | point | sd | 95% | P(Δ < 0) |
|---|---|---|---|---|
| Efron | 1.290 | 0.710 | [-0.109, 2.669] | 0.035 |
| Bayes | 1.290 | 0.704 | [-0.109, 2.660] | 0.035 |
The point estimate is +1.29 percentage points per month in the market’s favour, with a 95% interval that straddles zero. Nothing here establishes that IBM underperformed; the data are consistent with no difference at all.
The interpretive difference between the columns is the whole point of the duality section:
Efron licenses
“A one-sided test at the 3.5% level would reject.”
A statement about the procedure’s long-run behaviour over hypothetical repeat samples.
Bayes licenses
“There is a 3.5% posterior probability that IBM’s mean exceeded the market’s, given the model.”
A statement about the parameter.
Same numbers to two decimal places. Only one of them is a probability statement about the thing you wanted to know.
10.7 The gap closes as the window lengthens
Every number so far was computed at \(n = 59\). The ratio \(\sqrt{n/(n+1)}\) predicts that the disagreement widens on a shorter window and all but vanishes on a longer one, which is easy to check directly: run the same S&P series over twelve months, over the head-to-head 59, and over everything back to 1990.
| window | n | Bayes | Efron | classical | ratio |
|---|---|---|---|---|---|
| Jul 2025 - Jun 2026 | 12 | 0.6996 | 0.7281 | 0.7605 | 0.9608 |
| Mar 2013 - Jan 2018 | 59 | 0.2872 | 0.2896 | 0.2921 | 0.9916 |
| Feb 1990 - Jun 2026 | 437 | 0.1693 | 0.1695 | 0.1697 | 0.9989 |
11 Computational cost
Cost is \(O(B\times T_\text{estimator})\); weight generation is \(O(Bn)\).
Cheap to implement. The two methods differ by one line.
Code
n, B, rng = 59, 200_000, np.random.default_rng(20260726)
# Efron
W = rng.multinomial(n, np.full(n, 1 / n), size=B) / n
# Bayesian
E = rng.exponential(1.0, size=(B, n)); W = E / E.sum(axis=1, keepdims=True)Cheap to compute — but be precise about why.
- Weight generation is \(O(n)\) per replicate for both, so \(O(Bn)\) overall.
- Working memory is \(O(n)\) per replicate. Vectorising to a dense \((B, n)\) array makes it \(O(Bn)\) — 94 MB at \(B = 200\,000, n = 59\), and 16 GB for the MNIST test set. So chunk the replicate axis; every result in this post is computed in bounded memory that way.
- Evaluating the functional is \(O(n)\) for the mean and for any weighted cumulative-sum statistic.
- The sort-once trick, which falls out of the weights view. Order-statistic functionals — median, quantiles, max, min, the CDF — look like they need a re-sort per replicate, \(O(Bn\log n)\), if you think of Efron as producing a new dataset. Viewing it as weights lets you sort the data once, \(O(n\log n)\), and make every replicate a single \(O(n)\) cumulative sum. The same trick serves both methods. Viewing Efron as a weighted bootstrap is the optimisation.
Both are \(O(Bn)\). The differences are constants, and constants must be measured rather than assumed.
| operation | ms |
|---|---|
| Efron weights via rng.multinomial | 366.4 |
| Bayes weights via rng.exponential + normalise | 46.2 |
| Efron via index resampling (rng.integers) | 18.1 |
| Efron median, naive re-sort per replicate | 105.6 |
| Bayes median, sort once + cumulative sum | 81.5 |
The counter-intuitive result, which surprises most people: generating Dirichlet weights is about 7.9× faster than NumPy’s multinomial sampler at this \(n\). rng.multinomial is expensive when there are many categories. This is library-specific and version-specific; the table above was regenerated when this page was built, and you should regenerate it on yours.
Where the Bayesian bootstrap is asymptotically cheaper: sometimes it needs \(B = 0\). The CDF at any point is \(\text{Beta}(k, n-k)\); the median’s rank law is \(\text{Binomial}(n-1, \tfrac12)\); accuracy is \(\text{Beta}(k, n-k)\). Closed forms, exact, no simulation. Efron has no comparable shortcut and must simulate.
Degenerate-replicate waste. If your estimator fails on a resample that is missing a class, your effective \(B\) is lower than your nominal \(B\) — you paid for replicates you had to discard. The MNIST rare-class case below quantifies this.
Both methods are embarrassingly parallel. And for streaming settings there is the Poisson bootstrap: because its weights are i.i.d., it needs no knowledge of \(n\), uses \(O(1)\) memory and a single pass, which is why it is the version that runs at scale.
Not cheap, and this is the honest part: the total cost is \(B\) times the cost of the estimator. For a mean that is nothing. For a neural network it is everything.
12 MNIST example
Two resampling targets: test-set metrics (cheap) vs training-set refits (expensive).
12.1 Where MNIST came from, and why it is the right dataset here
MNIST is 70,000 greyscale images of handwritten digits, 28×28 pixels, labelled 0–9: 60,000 for training and a fixed, public test set of 10,000. The handwriting was not gathered for machine learning at all. It comes from forms filled in for the US National Institute of Standards and Technology — US Census Bureau staff and high-school students writing digits into boxes — which Yann LeCun, Corinna Cortes and Christopher Burges then size-normalised and centred into the images everyone now downloads.
It became the field’s standard first benchmark, and that is exactly what makes it useful here. Nobody is going to be misled about the modelling, so all the attention stays on the uncertainty.
The model that produced every number in this section is a 784–128–10 MLP — a plain feed-forward network, 784 pixels in, one hidden layer of 128 units, ten class scores out — trained for three epochs, on torch 2.13.0 on mps. Its test accuracy is 0.9687 — \(k = 9\,687\) correct out of \(n = 10\,000\) — and a single fit takes 4.2 seconds.
The crucial structural point is that there are two different populations you could resample here, and they answer different questions:
Resample the test set
One row = one test image. \(n = 10\,000\).
Answers: how much would this metric move if I had drawn a different test set from the same distribution?
Costs nothing — the model is already trained.
Resample the training set
One row = one training image. \(n = 60\,000\).
Answers: how much would this metric move if I had trained on a different sample?
Costs a full refit per replicate.
12.2 Every evaluation metric is a functional too
Nothing about the earlier machinery changes here. Every evaluation metric is a functional of the distribution of test examples, exactly as the mean and the tail probability were functionals of the distribution of months:
| metric | functional | effective denominator |
|---|---|---|
| overall accuracy | \(T(F) = F(\text{correct})\) | all 10,000 |
| per-class recall | \(F(\text{correct} \mid \text{true} = c)\) | ~1,000 |
| per-class precision | \(F(\text{correct} \mid \text{pred} = c)\) | ~1,000 |
| a confusion-matrix cell | \(F(\text{true}=c,\ \text{pred}=c')\) | can be single digits |
12.3 What a wrong number would ship
Three decisions rest on these numbers. Ship or don’t ship: is accuracy above the threshold in the requirements, including uncertainty? Model A versus model B: is a 0.3-point gap on the leaderboard real, or inside the noise of a 10,000-example test set? Where to spend the next labelling budget: which class is genuinely worst, as opposed to merely unlucky in this particular test set?
The third is where the choice of bootstrap starts to matter, because per-class and per-cell metrics have small denominators — and that is precisely the regime, established on the S&P tail above, where Efron’s resampling spends a large share of its replicates in a world where the relevant examples do not exist.
12.4 The cheap path — bootstrap the evaluation, not the training
Train once. Cache the per-example correctness vector. Now every evaluation metric is a functional of the empirical distribution of test examples, and bootstrapping is free.
Accuracy is \(\sum_i w_i \mathbb{1}\{\text{correct}_i\}\) — a sum of \(k\) of the weights. By the aggregation property derived earlier, the Bayesian posterior is
\[\text{accuracy} \mid \text{data} \;\sim\; \text{Beta}(k,\ n-k)\]
exactly, with no simulation at all.
Code
# The entire Bayesian bootstrap for accuracy, with B = 0.
from scipy import stats
posterior = stats.beta(k_te, n_te - k_te)
print(f"accuracy {k_te/n_te:.4f} sd {posterior.std():.6f} "
f"95% [{posterior.ppf(.025):.5f}, {posterior.ppf(.975):.5f}]")accuracy 0.9687 sd 0.001741 95% [0.96520, 0.97202]Monte Carlo agrees: 20,000 simulated replicates give a standard deviation of 0.001747 against the exact 0.001741. Efron gives \(\text{Bin}(n, k/n)/n\), indistinguishable at this \(n\).
Not every metric is that obliging. Per-class recall is a ratio of two weighted sums rather than a single one, and macro-F1 — the F1 score worked out separately for each digit and then averaged over the ten, so a rare class gets the same say as a common one — is a ratio of ratios. Those need simulation, and the interesting column is the last one.
| metric | point | Bayes sd | Efron sd | support |
|---|---|---|---|---|
| macro-F1 | 0.96844 | 0.00175 | 0.00172 | 10000 |
| recall, class 0 | 0.98776 | 0.00350 | 0.00343 | 980 |
| recall, class 1 | 0.98767 | 0.00327 | 0.00328 | 1135 |
| recall, class 2 | 0.96415 | 0.00580 | 0.00577 | 1032 |
| recall, class 3 | 0.95941 | 0.00626 | 0.00624 | 1010 |
| recall, class 4 | 0.97352 | 0.00519 | 0.00514 | 982 |
| recall, class 5 | 0.96300 | 0.00627 | 0.00632 | 892 |
| recall, class 6 | 0.97599 | 0.00492 | 0.00490 | 958 |
| recall, class 7 | 0.95817 | 0.00628 | 0.00625 | 1028 |
| recall, class 8 | 0.95380 | 0.00664 | 0.00669 | 974 |
| recall, class 9 | 0.96135 | 0.00611 | 0.00599 | 1009 |
12.5 Where the two diverge on real ML numbers: rare cells
Take any per-class or per-slice metric whose denominator is small. When an error count is \(k = 1\) out of \(n = 50\):
Efron’s bootstrap reports an error rate of exactly zero in 36.3% of replicates (theory: 36.4%) — the identical \(((n-k)/n)^n\) atom that appeared in the S&P tail. A 95% interval whose lower endpoint is a hard zero is a statement about the resampling scheme, not about the classifier. The Bayesian answer is the exact \(\text{Beta}(1, 49)\), with a 95% interval of [0.00052, 0.0725].
Recommendation, stated plainly: use the Bayesian bootstrap, or its exact Beta form, for any per-class or per-slice metric with a small denominator. The reason is one sentence — continuous weights cannot delete an observation.
12.6 The expensive path — bootstrap the training
This is the weighted likelihood bootstrap: draw \(w^{(b)} \sim \text{Dirichlet}(1,\dots,1)\) over the \(n\) training points and minimise
\[\sum_{i=1}^{n} w^{(b)}_i \,\ell(x_i, \theta).\]
It gives an approximate posterior for any estimator defined by minimising a loss — generalised linear models, boosted trees, neural nets — and it gets there without Markov chain Monte Carlo: no proposals to accept or reject, no chains to diagnose, only refits you already know how to run. In PyTorch it is a per-example loss multiplied by a weight tensor:
Code
criterion = nn.CrossEntropyLoss(reduction="none") # the whole trick
...
per_example = criterion(net(xb), yb)
loss = (per_example * wb).mean() # the weighted likelihood bootstrap
loss.backward()One line different from weighted training. The Efron equivalent is WeightedRandomSampler or integer duplication.
| what | fits | wall clock |
|---|---|---|
| single fit (the cheap path trains once) | 1 | 4.2 s |
| weighted likelihood bootstrap, B = 20 | 20 | 41.3 s |
| extrapolated to B = 200 | 200 | 6.9 min |
| a model taking 6 h per fit, B = 200 | 200 | 1200 h |
At \(B = 20\) the refits give a spread of test accuracy with standard deviation 0.00236, against the evaluation-only posterior standard deviation of 0.001741. The retraining interval is wider, and it is measuring something different — which brings us to the part that should not be softened.
12.7 The real drawbacks
Both paths measure one thing well and several things not at all, and it is worth being blunt about which is which. Every item below is a source of error the interval does not contain.
- The bootstrap prices sampling variability of the data. It does not price model misspecification, distribution shift, label noise, or the choice of architecture.
- It does not capture training randomness — initialisation, SGD ordering — unless you deliberately resample that too, in which case you are measuring a different quantity and should say which.
- Bootstrapping the test set tells you about test-set sampling variability, not about behaviour on a new deployment distribution. On a fixed public benchmark like MNIST everyone shares one test set, so this interval does not cover the thing people usually want when they ask how accurate a model is.
- Retraining-based intervals inherit every pathology of the training procedure, including non-convexity: two fits with different weights may land in different basins, and the resulting spread mixes parameter uncertainty with optimisation noise. The 0.00236 above contains both, inseparably.
- Both bootstraps still assume i.i.d. sampling. Grouped, temporal or hierarchical data need the resampling to happen at the level of the independent unit.
13 Non-weightable estimators
Weights must enter linearly as \(\sum_i w_i \ell(x_i,\theta)\). Failures fall into four classes.
13.1 (a) The weighting is ambiguous
Some statistics are not built from single observations at all, but from every pair of them: Gini, Kendall’s \(\tau\), Spearman, Mann–Whitney \(U\), the Hodges–Lehmann estimator, energy distance, maximum mean discrepancy. These are the U- and V-statistics, and they are where the weights view first runs into trouble. For a statistic \(\sum_{i<j} h(x_i, x_j)\) there is no single obvious weighted form. Using \(w_i w_j\) is defensible but is not the only choice, and different choices give different answers.
Concretely, on the S&P window the Hodges–Lehmann estimator (the median of the \(\binom{59}{2}\) pairwise averages) weighted by \(w_i w_j\) has posterior standard deviation 0.2537, while weighting by \(\min(w_i, w_j)\) gives 0.1378 — a 46% difference produced entirely by a modelling choice with no principled resolution. Efron sidesteps the question by producing an actual sample, which has an unambiguous set of pairs.
13.2 (b) A count is structurally required
\(k\)-nearest-neighbours: is \(k\) a count of points, or a mass of \(k/n\)? Decision trees with min_samples_leaf — a count, which is precisely why scikit-learn also offers min_weight_fraction_leaf, the weighted-world translation and the canonical example of the fix. Cross-validation fold construction. Anything with an integer threshold or a hard cardinality constraint. The number of distinct values. And the extrema and the mode from earlier, which are functionals of the support rather than of the masses.
13.3 (c) The objective does not decompose over observations
Time series and dependent data generally: an AR or GARCH likelihood does not factor into per-observation terms that reweighting respects, and reweighting individual months breaks exactly the dependence you care about. This applies directly to the returns used throughout this post — see the caveat at the top.
Also: ranking and contrastive losses defined over pairs or batches; a softmax normalised over a whole dataset; batch normalisation, whose statistics depend on batch composition rather than on any weighted sum; graph and network data; and clustered or matched designs, where the weights must attach to clusters rather than to rows. The fix in that last case is a Dirichlet-weighted block bootstrap — the same idea, one level up. That is also the fix for the non-i.i.d. returns problem.
13.4 (d) The software will not take weights
The mundane case, and by far the most common. There is a clean escape hatch, and it is derivable.
Draw \(w \sim \text{Dirichlet}(1,\dots,1)\), then resample \(m\) points i.i.d. with probabilities \(w\) instead of uniformly, and hand that ordinary-looking sample to the weight-blind code.
What does this cost? Use the law of total variance on the resample mean \(\bar X_m\), conditioning on \(w\). Given \(w\), the \(m\) draws are i.i.d. from \(F_w\), so \(E[\bar X_m \mid w] = \mu_w := \sum_i w_i x_i\) and \(\operatorname{Var}(\bar X_m \mid w) = \sigma^2_w/m\) where \(\sigma^2_w = \sum_i w_i x_i^2 - \mu_w^2\). Then
\[\operatorname{Var}(\bar X_m) = \underbrace{\operatorname{Var}(\mu_w)}_{=\ \tilde s^2/(n+1)} + \frac{1}{m}\,E[\sigma^2_w].\]
For \(\text{Dirichlet}(1,\dots,1)\), \(E[\sigma_w^2] = \frac{n}{n+1}\tilde s^2\), and substituting gives
\[\operatorname{Var}(\bar X_m) = \frac{\tilde s^2}{n+1} + \frac{n\,\tilde s^2}{m(n+1)} = \frac{\tilde s^2}{n+1}\Big(1 + \frac{n}{m}\Big).\]
The inflation factor is exactly \(1 + n/m\).
| m | m/n | MC variance | formula | inflation |
|---|---|---|---|---|
| 5.000 | 1.000 | 3.328 | 3.333 | 2.000 |
| 25.000 | 5.000 | 2.002 | 2.000 | 1.200 |
| 100.000 | 20.000 | 1.749 | 1.750 | 1.050 |
| 500.000 | 100.000 | 1.682 | 1.683 | 1.010 |
So \(m = n\) doubles the variance — that is exactly the \((m = n, \alpha = 1)\) corner warned about earlier — while \(m \approx 100n\) costs about 1%. This is the Pólya-urn bridge turned into a recipe, and it is the single most useful practical takeaway in the section.
14 Method selection
Reach for the Bayesian bootstrap when you want a posterior rather than a sampling distribution; when degenerate resamples would break the estimator — a missing class, a zero-variance feature, an empty stratum; when the statistic leans on few observations, which is the rare-cell case that recurs everywhere in machine learning; when you want to differentiate through the weights or drop them into a loss; and when you want an exact answer instead of a simulated one.
Reach for Efron when your estimator needs a sample rather than weights; when you want to make no Bayesian commitment at all; and when you want the surrounding machinery — BCa and bootstrap-\(t\) for better interval coverage, the block bootstrap for dependent data, the double bootstrap for correcting the correction — which mostly exists only for the Efron version.
And note that both fail identically on dependent data, extrema, heavy tails, and very small \(n\). Those failures matter far more than the choice between the two.
Bayesian bootstrap: posterior interpretation; no zero weights; exact Beta laws for proportions; weighted likelihood.
Efron: sampling distribution; BCa/block bootstrap tooling; estimators needing actual resamples.
15 Reference card
Everything derived above, on one page: each row a quantity you might actually need, with the Efron answer and the Bayesian answer side by side.
| quantity | Efron | Bayesian bootstrap |
|---|---|---|
| weight mean | \(1/n\) | \(1/n\) |
| weight variance | \((n-1)/n^3\) | \((n-1)/[n^2(n+1)]\) |
| weight covariance | \(-1/n^3\) | \(-1/[n^2(n+1)]\) |
| marginal law | \(\mathrm{Bin}(n,1/n)/n\) | \(\mathrm{Beta}(1,n-1)\) |
| \(P(w_i = 0)\) | \((1-1/n)^n \to e^{-1}\) | \(0\) |
| variance of the mean | \(\tilde s^2/n\) | \(\tilde s^2/(n+1)\) |
| median law | no closed form | \(P(x_{(j)}) = \binom{n-1}{j-1}2^{-(n-1)}\) |
| \(F(t)\), \(k=\#\{x_i \le t\}\) | \(\mathrm{Bin}(n,k/n)/n\) | \(\mathrm{Beta}(k,n-k)\) |
| a proportion / accuracy | \(\mathrm{Bin}(n,k/n)/n\) | \(\mathrm{Beta}(k,n-k)\), exact |
| maximum | \(P(\max \le x_{(j)}) = (j/n)^n\) | point mass at \(x_{(n)}\) |
| mode of \(\hat F\) | decided by tie-breaking | uniform, \(1/n\) each |
| support | a lattice of \(\binom{2n-1}{n-1}\) points | the whole simplex |
| generator | \(\mathrm{Poisson}(1)\) cond. \(\sum W_i = n\) | \(\mathrm{Exp}(1)\), normalised |
Each method is six lines.
Code
def efron_weights(n, B, rng):
return rng.multinomial(n, np.full(n, 1 / n), size=B) / n
def dirichlet_weights(n, B, rng):
E = rng.exponential(1.0, size=(B, n))
return E / E.sum(axis=1, keepdims=True)15.1 Reproducibility
The claim made at the top — that the page does not build if the algebra and the arithmetic disagree — is worth making checkable. What follows is the build’s own report on itself: a layout check over every figure, then the versions and seeds behind every number above.
fig-simplex ok
fig-point-is-dist ok
fig-dirichlet ok
fig-prior ok
fig-lattice ok
fig-marginal ok
fig-regions ok
fig-staircase ok
fig-returns ok
fig-mnist ok
------------------------------------------------
10 figures checked, 0 collisions| item | value |
|---|---|
| seed | 20260726 |
| replicates, headline | 200,000 |
| replicates, figures | 20,000 |
| python | 3.13.2 |
| numpy | 2.5.1 |
| scipy | 1.18.0 |
| pandas | 3.0.5 |
| matplotlib | 3.11.1 |
| torch | 2.13.0 |
| MNIST backend | torch 2.13.0 on mps |
| S&P 500 source | Shiller monthly series, datasets/s-and-p-500 |
| IBM source | plotly/datasets all_stocks_5yr.csv |
| data retrieved | 2026-07-26 |
Every figure above is reproducible from the single seed 20260726; each experiment derives its own independent generator from it by name, so adding an experiment never shifts an existing one. Every closed form stated in the prose is asserted against simulation in a code cell at build time, with tolerances set by the Monte-Carlo standard error — this page does not render if the algebra and the arithmetic disagree. The figure layout checker reported above walks every text artist in every figure and fails the build on any overlap.
Data provenance and its caveats are given where the two series are introduced: the S&P series is Shiller’s monthly averages of daily closes, and the IBM series is unadjusted closes and therefore price-only. Both were retrieved on 2026-07-26. The MNIST test set is fixed and public, so the intervals here describe test-set sampling variability and nothing else.
16 References
Efron, B. (1979). Bootstrap methods: another look at the jackknife. Annals of Statistics 7(1), 1–26.
Rubin, D. B. (1981). The Bayesian bootstrap. Annals of Statistics 9(1), 130–134.
Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics 1(2), 209–230.
Lo, A. Y. (1987). A large sample study of the Bayesian bootstrap. Annals of Statistics 15(1), 360–375.
Newton, M. A. and Raftery, A. E. (1994). Approximate Bayesian inference with the weighted likelihood bootstrap. JRSS B 56(1), 3–48.
Præstgaard, J. and Wellner, J. A. (1993). Exchangeably weighted bootstraps of the general empirical process. Annals of Probability 21(4), 2053–2086.
Blackwell, D. and MacQueen, J. B. (1973). Ferguson distributions via Pólya urn schemes. Annals of Statistics 1(2), 353–355.
Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica 4, 639–650.
Efron, B. (2012). Bayesian inference and the parametric bootstrap. Annals of Applied Statistics 6(4), 1971–1997.
The Poor Person’s Bayesian — survey-level view of bootstrap-as-cheap-Bayes.