The bootstrap, permutation tests and A/B testing

Before this lesson: lesson 11 (p-values, power, the two-proportion z-test, multiple testing), lesson 04 (the variance of a sum, correlation), lesson 06 (standard errors and confidence intervals). Time: about 80 minutes, plus the notebook.

---

Why this matters #

• Error bars for any ML metric. Accuracy, AUC, F1, BLEU or an LLM benchmark score: the bootstrap over test examples gives an interval for each, and a paired bootstrap gives an interval for the difference of two models.
• A/B tests are how ML systems are launched. New ranking, recommendation, search and ads models, and new LLM features, go to real users through A/B tests. Data scientist and ML engineer interviews ask about sample size, peeking and CUPED.
• Bagging and random forests train each model on a bootstrap resample; the points left out of a resample (about 37%) give the out-of-bag error.
• Small randomised studies in science (a few patients, plants or field plots) use permutation tests, because they need no normality assumption and are exact.

1. Predict first #

2. The picture #

Permutation: shuffle the labels. Six students prepare for an exam. Three of them, chosen at random, use a new study app. Their score gains are 7, 9 and 12 points. The other three gain 3, 5 and 6 points. The statistic is the difference of the group means, "mean of the app group minus mean of the other group":

Under $H_0$ the app does nothing. Then each student would have the same gain with or without the app, and the labels "app" and "no app" are just a random split of six numbers into two groups of three. There are $\binom{6}{3} = 20$ such splits, all equally likely under $H_0$. Compute $T$ for each one:

Only the observed split reaches $T \ge 4.67$, so the one-sided p-value is $1/20 = 0.05$. Two splits reach $|T| \ge 4.67$, so the two-sided p-value is $2/20 = 0.10$. No formula and no normality assumption were used. Notice also the floor: with 3 against 3, no outcome can give a one-sided p-value below $1/20$.

Interactive figure: loading…

Bootstrap: resample the data. You want the spread of a statistic $\hat\theta$ (a mean, a median, an AUC) over repeated samples from the unknown distribution $F$. You cannot draw new samples from $F$. So draw them from the closest thing you have: the empirical distribution $\hat F_n$, which puts mass $1/n$ on each observation. A sample from $\hat F_n$ is $n$ draws with replacement from the data: some points appear twice or more, others not at all.

real world	bootstrap world
unknown $F$	known $\hat F_n$
sample $x_1, \dots, x_n$ from $F$	resample $x_1^*, \dots, x_n^*$ from $\hat F_n$ (with replacement)
estimate $\hat\theta$; true value $\theta$	replicate $\hat\theta^*$; "true" value $\hat\theta$
spread of $\hat\theta - \theta$: unknown	spread of $\hat\theta^* - \hat\theta$: simulate it

The bootstrap assumes that the right column behaves like the left column. It approximates the sampling distribution of $\hat\theta$, not the population.

Interactive figure: loading…

3. Definitions and main results #

Permutation tests #

Two samples $x_1, \dots, x_m$ and $y_1, \dots, y_n$ are pooled into $z = (z_1, \dots, z_N)$, $N = m + n$. A statistic $T$ compares the groups, for example the difference in means. The null hypothesis is that the pooled vector is exchangeable: its joint distribution does not change when you permute its entries. Two i.i.d. samples from the same distribution are exchangeable. So are units that a randomised experiment assigned at random to two treatments that have no effect.

Permutation p-values. For the exact test, compute $T$ on all $M = \binom{N}{m}$ splits and let $p$ be the fraction with $T \ge T_{\text{obs}}$ (the observed split is one of them). When $M$ is too large, draw $B$ random permutations, compute $T_1^*, \dots, T_B^*$, count $C = \#\{b : T_b^* \ge T_{\text{obs}}\}$, and use

Theorem (permutation tests are exact). Under exchangeability, both p-values satisfy $P(p \le \alpha) \le \alpha$ for every $\alpha$, every $N$, and (for the Monte Carlo version) every $B$.

Proof (Monte Carlo version; the exact version uses the same argument with all $M$ splits).

1. Condition on the multiset of pooled values. Under exchangeability, the observed labelling is a uniformly random labelling of these values. Each random permutation gives another uniformly random labelling, independent of the first.
2. So $T_{\text{obs}}, T_1^*, \dots, T_B^*$ are i.i.d. given the multiset, hence exchangeable: $T_{\text{obs}}$ is equally likely to take each of the $B + 1$ ranks.
3. Without ties, $P(C \le j) = (j + 1)/(B + 1)$ for $j = 0, \dots, B$, because $C \le j$ means that $T_{\text{obs}}$ is among the $j + 1$ largest. Ties can only make $C$ larger, so in general $P(C \le j) \le (j + 1)/(B + 1)$.
4. $p \le \alpha$ exactly when $C + 1 \le \alpha(B + 1)$. So $P(p \le \alpha) \le \lfloor \alpha(B + 1) \rfloor/(B + 1) \le \alpha$. $\blacksquare$

The "+1" counts the observed labelling as one of the permutations, which it is. Without it, $C/B$ can be 0, a p-value that claims the impossible, and the test can reject too often. The smallest possible p-value is $1/(B + 1)$. The Monte Carlo error of $p$ is about $\sqrt{p(1 - p)/B}$, so use $B = 999$ for a quick look, $B = 9{,}999$ or more for a final answer, and many more for a tiny threshold (Bonferroni over many tests).

Practical points.

• Two-sided p-values. Count $|T_b^*| \ge |T_{\text{obs}}|$, or take twice the smaller one-sided p-value (SciPy does this). They agree when the permutation distribution is symmetric, as in section 2.
• Any statistic works: a difference in medians, a t-statistic, a correlation. To test the independence of two variables, permute one column and keep the other fixed.
• Paired data: under $H_0$ the sign of each within-pair difference is random, so flip signs instead of shuffling labels.
• What $H_0$ really says: the two groups have the same distribution, not only the same mean. If the variances differ under $H_0$, the labels are not exchangeable, and the permutation test of means is no longer exact. Use Welch (lesson 11), or permute a studentised statistic.
• Dependent data (time series, users with many sessions) are not exchangeable unit by unit. Permute whole blocks or whole clusters.

The bootstrap #

Let $\hat\theta = s(x_1, \dots, x_n)$ be a statistic. The bootstrap distribution is the distribution of $\hat\theta^* = s(x_1^*, \dots, x_n^*)$ when the $x_i^*$ are drawn i.i.d. from $\hat F_n$. You approximate it with $B$ resamples $\hat\theta_1^*, \dots, \hat\theta_B^*$.

• Bootstrap standard error: the standard deviation of $\hat\theta_1^*, \dots, \hat\theta_B^*$.
• Percentile interval at level $1 - \alpha$: the $\alpha/2$ and $1 - \alpha/2$ quantiles of the $\hat\theta_b^*$.
• How many resamples: about 1,000 for a standard error, and 10,000 or more for interval endpoints, so that two analysts with the same data get the same answer.

Who is in a resample? A given point is missed by one draw with probability $1 - 1/n$, so it appears in the resample with probability $1 - (1 - 1/n)^n$. This is 0.651 for $n = 10$, 0.634 for $n = 100$, and it tends to $1 - 1/e = 0.632$. A resample contains about 63% of the distinct points; the rest of its entries are repeats.

The mean, exactly. For $\hat\theta = \bar x$, the bootstrap variance is $\hat\sigma^2/n$ with the plug-in variance $\hat\sigma^2 = \frac1n \sum_i (x_i - \bar x)^2$. So the bootstrap standard error is $\sqrt{(n-1)/n}$ times the usual $S/\sqrt n$: slightly too small. This "narrowness bias" matters only for small $n$.

Why the bootstrap works, and how well

For means (with finite variance), smooth functions of means, and quantiles of a distribution with a continuous positive density, the distribution of $\sqrt n(\hat\theta^* - \hat\theta)$ given the data converges to the limit distribution of $\sqrt n(\hat\theta - \theta)$. So the bootstrap gets the standard error and the shape right to first order. The percentile interval would be exact if some monotone transformation $g$ made $g(\hat\theta) - g(\theta)$ symmetric about 0, with the same distribution in both worlds. Otherwise each of its two tails has a coverage error of order $1/\sqrt n$; for skewed statistics most misses fall on one side. The bootstrap-t and BCa (bias-corrected and accelerated) intervals reduce this error to order $1/n$. BCa is the default of scipy.stats.bootstrap.

Checking coverage by simulation. Pick a known $F$ and a true $\theta$. Repeat many times: draw a sample of size $n$ from $F$, compute the interval, record whether it contains $\theta$. The share of hits estimates the coverage; with $R$ repetitions its standard error is about $\sqrt{0.95 \cdot 0.05/R}$. Section 5 does this: the 95% percentile interval covers the mean of Exp(1) data with $n = 30$ in about 93% of data sets, the median with $n = 60$ in about 95%, and the maximum of uniform data in 0%.

When the bootstrap fails, and what to use instead.

• Extremes. For the maximum $M = \max_i x_i$ of $\text{Uniform}(0, \theta)$ data, every resample has $\hat\theta^* \le M < \theta$, so the percentile interval never contains $\theta$. Also $P(\hat\theta^* = M) = 1 - (1 - 1/n)^n \to 0.632$, while the true distribution of $M$ is continuous. Instead: the exact pivot $M/\theta$, whose CDF is $u^n$ on $[0, 1]$; it gives the interval $[M,\ M/\alpha^{1/n}]$ with coverage exactly $1 - \alpha$. Or extreme-value theory, or an $m$-out-of-$n$ bootstrap.
• Heavy tails. For the mean of data with infinite variance (Cauchy, Pareto with tail index below 2), the bootstrap is inconsistent. Even with finite but huge variance (revenue per user), one extreme value decides each resample. Instead: a robust statistic (median, trimmed mean), capping (winsorising) the metric, or a parametric tail model.
• Dependent data. Resampling single time steps of an autocorrelated series destroys the correlation, so the standard error is far too small. For the mean of an AR(1) series $x_t = 0.7\,x_{t-1} + e_t$ of length 500, the true SE is 0.150 and the i.i.d. bootstrap gives 0.062 (notebook, exercise 3). Instead: a moving block bootstrap: resample blocks of $\ell$ consecutive values, with $\ell$ much longer than the correlation time (blocks of 25 give 0.135 here). Or a cluster bootstrap that resamples whole units (users, stations).
• Very small $n$. With $n = 5$ the empirical distribution has at most 5 points, and the percentile interval undercovers. Instead: a t interval, a bootstrap-t interval, an exact or parametric method, or more data.

Permutation vs bootstrap. The permutation test is exact for finite samples, because it uses a symmetry that $H_0$ implies. The bootstrap is approximate, because $\hat F_n$ is close to $F$ only for large $n$. Use a permutation test to test "no difference" in a randomised comparison; use the bootstrap for standard errors and intervals.

Comparing two ML models on one test set #

Two classifiers are scored on the same $n$ test examples. The two accuracies are strongly correlated (both models find the same examples easy), so treat them as paired data:

1. For each example $i$, compute the per-example score difference $d_i$ (for accuracy: $d_i \in \{-1, 0, 1\}$).
2. Draw $B$ bootstrap resamples of the example indices. Use the same resample for both models.
3. The percentile interval of the resampled mean difference is the interval for the gain. Or run a permutation test that swaps the two models' outputs on each example (a sign flip of $d_i$).

Resample the independent unit (documents, users, patients), not tokens or pixels. The test set is only one source of noise: retrain with several seeds too, and report the spread over seeds. A gain smaller than the seed-to-seed variation is not a result.

A/B tests: sample size #

Users are randomised to control (arm 1) and treatment (arm 2), $n$ users per arm. The metric has standard deviation $\sigma$ per user. The minimum detectable effect (MDE) $\delta$ is the smallest true difference you want to detect, with a two-sided test at level $\alpha$ and power $1 - \beta_{\text{II}}$.

Theorem (users per arm, equal variances). With the normal approximation,

Proof.

1. The difference in means $D$ has standard error $\operatorname{SE} = \sigma\sqrt{2/n}$. The test rejects when $|D| > z_{1-\alpha/2}\,\operatorname{SE}$.
2. If the true difference is $\delta > 0$, then $D \approx \mathcal N(\delta, \operatorname{SE}^2)$. Ignore the tiny chance of rejecting in the wrong direction: $\text{power} \approx \Phi\big((\delta - z_{1-\alpha/2}\,\operatorname{SE})/\operatorname{SE}\big)$.
3. Set the power to $1 - \beta_{\text{II}}$: $\delta/\operatorname{SE} - z_{1-\alpha/2} = z_{1-\beta_{\text{II}}}$, so $\delta = (z_{1-\alpha/2} + z_{1-\beta_{\text{II}}})\,\sigma\sqrt{2/n}$.
4. Solve for $n$. $\blacksquare$

Rule of thumb. For $\alpha = 0.05$ and power 0.8, $2(1.960 + 0.842)^2 = 15.7$, so

(Kohavi et al., 2009). For power 0.9 the constant is 21.

Two proportions. For conversion rates $p_1$ (baseline) and $p_2 = p_1 + \delta$, the variance of one user is $p(1 - p)$, and it is not the same under $H_0$ and under $H_1$. The test uses the pooled variance under $H_0$, with $\bar p = (p_1 + p_2)/2$; the power uses the variances under $H_1$. The same four steps give

with $s_0 = \sqrt{2\bar p(1 - \bar p)}$ (from $H_0$) and $s_1 = \sqrt{p_1(1 - p_1) + p_2(1 - p_2)}$ (from $H_1$).

Interactive figure: loading…

What the formula says.

• $n \propto 1/\delta^2$: half the MDE needs about 4 times the users.
• $n \propto \sigma^2$: a metric with less noise (conversion instead of revenue) or less variance (CUPED, below) needs fewer users.
• An MDE is absolute (4.0% → 4.5%, $\delta = 0.005$) or relative (+12.5%). Keep $\delta$ and $\sigma$ in the same units.
• An unequal split costs power: with a share $f$ of users in treatment, the test needs about $1/(4f(1 - f))$ times as many users as a 50/50 test (2.8 times for 90/10).

A/B tests: the main pitfalls #

• Peeking (optional stopping). If you check the p-value repeatedly and stop at the first $p < 0.05$, the type I error grows with the number of looks. With equally spaced looks at a fixed planned sample it is about 0.08 with 2 looks, 0.14 with 5, 0.19 with 10 and 0.25 with 20 (section 5). The z-statistics at different looks share most of their data, so they are strongly but not perfectly correlated; each new look adds a new chance to cross 1.96. With unlimited looks the error tends to 1. Prevent: fix the sample size in advance and analyse once, or use a sequential design built for many looks (group-sequential boundaries such as Pocock or O'Brien–Fleming, or always-valid p-values).
• Sample-ratio mismatch (SRM). The observed split differs from the design (say 50/50) more than chance allows. Test the assignment counts with a χ² goodness-of-fit test and a strict threshold such as $p < 0.001$. Example: 50,600 vs 49,400 users gives $X^2 = 2 \cdot 600^2/50{,}000 = 14.4$ and $p = 0.00015$. An SRM means a bug in assignment, logging or filtering (for example, the treatment crashes for some browsers, and those users are never logged). The remaining users of the two arms are then not comparable. Do not trust the result until you know the cause.
• Novelty and primacy effects. Users react to the change itself: they click a new feature because it is new (novelty), or they are slower with a new layout at first (primacy). Detect: plot the effect by day, and compare new with returning users. Prevent: run whole weeks (this also handles day-of-week effects), long enough for the effect to settle.
• Many metrics or many variants. With 20 independent metrics at $\alpha = 0.05$, the chance of at least one false positive is $1 - 0.95^{20} = 0.64$. Prevent: choose one primary metric (the overall evaluation criterion, OEC) before the test, correct the others for multiplicity (Holm or BH, lesson 11), and treat surprising secondary results as ideas for a new test.
• Interference. The analysis assumes that one user's outcome does not depend on the other users' assignments (SUTVA, the stable unit treatment value assumption). Social networks, two-sided marketplaces and shared ad budgets break it. Prevent: randomise clusters (friend groups, regions), use switchback designs (switch the whole system between versions over time), or randomise markets.
• Unit of randomisation vs unit of analysis. You randomise users but compute the metric per page view. The views of one user are correlated, so the naive variance, which treats them as independent, is too small, and so are the p-values. With $m$ views per user and intra-user correlation $\rho$, the variance is larger by the factor $1 + (m - 1)\rho$. Prevent: analyse at the unit of randomisation (one value per user), or use the delta method or a user-level bootstrap for ratio metrics such as clicks per view.

A healthy experimentation system also runs A/A tests: about 5% of them should be significant at $\alpha = 0.05$, and their splits should match the design.

Interactive figure: loading…

CUPED: variance reduction with pre-experiment data #

Let $Y$ be the outcome of a user during the experiment, and $X$ a covariate measured before the experiment started, often the same metric in the weeks before. CUPED (controlled-experiment using pre-experiment data; Deng et al., 2013) replaces $Y$ by

Theorem (CUPED). Let $\rho = \operatorname{Corr}(X, Y)$, and let $\tau$ be the true treatment effect on the mean of $Y$. (a) $\operatorname{Var}(Y_{\text{cv}})$ is smallest at $\theta^* = \operatorname{Cov}(X, Y)/\operatorname{Var}(X)$, where it equals $\operatorname{Var}(Y)(1 - \rho^2)$. (b) For every fixed $\theta$, the estimate $\hat\tau_{\text{cv}} = \bar Y_{\text{cv},T} - \bar Y_{\text{cv},C}$ is unbiased for $\tau$.

Proof.

1. $\operatorname{Var}(Y_{\text{cv}}) = \operatorname{Var}(Y) - 2\theta\operatorname{Cov}(X, Y) + \theta^2\operatorname{Var}(X)$, a parabola in $\theta$ (lesson 04: the variance of a difference).
2. Its minimum is at $\theta^* = \operatorname{Cov}(X, Y)/\operatorname{Var}(X)$, with value $\operatorname{Var}(Y) - \operatorname{Cov}(X, Y)^2/\operatorname{Var}(X) = \operatorname{Var}(Y)(1 - \rho^2)$. This proves (a).
3. In the difference of the two arms, $\mathbb{E}[X]$ cancels: $\hat\tau_{\text{cv}} = (\bar Y_T - \bar Y_C) - \theta(\bar X_T - \bar X_C)$. So you never need to know $\mathbb{E}[X]$.
4. $X$ is measured before the random assignment, so it has the same distribution in both arms: $\mathbb{E}[\bar X_T] = \mathbb{E}[\bar X_C]$.
5. Hence $\mathbb{E}[\hat\tau_{\text{cv}}] = \mathbb{E}[\bar Y_T - \bar Y_C] = \tau$. This proves (b). $\blacksquare$

The variance of the estimate also falls by the factor $1 - \rho^2$, so the users needed fall by the same factor. With $\rho = 0.7$ the factor is 0.51: about half the users, or half the duration. Deng et al. report variance reductions of about 50% for some metrics at Bing with 1–2 weeks of pre-experiment data.

Interactive figure: loading…

Details that matter.

• Use one $\theta$ for both arms, estimated from the pooled data. Estimating $\theta$ adds a bias of order $1/n$ only.
• Never use a covariate that the treatment can change. Then $\mathbb{E}[\bar X_T] \ne \mathbb{E}[\bar X_C]$, and step 4 fails: the correction removes part of the real effect, or even flips its sign (section 5).
• Users without pre-period data: add an indicator "has pre-period data" as a second covariate.
• The same idea elsewhere: CUPED is a control variate (lesson 14), with $\mathbb{E}[X]$ cancelled by randomisation instead of being known. It is also close to regression adjustment (ANCOVA): regress $Y$ on the treatment indicator and $X$ (lesson 13).

Answers to section 1. (1) About 63: each point appears with probability $1 - (1 - 1/100)^{100} = 0.634$. (2) About 25%, not 5% (the peeking simulation in section 5 gives 0.247). (3) About 4 times as many, because $n \propto 1/\delta^2$: halving $\delta$ multiplies $n$ by 4 (the exact formula gives 3.9 at this baseline).

4. Worked example #

A. A bootstrap by hand. Take $n = 3$ observations $x = (1, 2, 6)$, with mean 3. There are $3^3 = 27$ equally likely ordered resamples. Their means and counts:

resample mean	1	4/3	5/3	2	8/3	3	10/3	13/3	14/3	6
count (of 27)	1	3	3	1	3	6	3	3	3	1

• Step 1, the mean: the bootstrap distribution has mean 3, the sample mean.
• Step 2, the variance: the plug-in variance is $\hat\sigma^2 = \frac{(1-3)^2 + (2-3)^2 + (6-3)^2}{3} = \frac{14}{3}$, so the bootstrap variance is $\hat\sigma^2/n = 14/9$, and the bootstrap SE is $\sqrt{14/9} = 1.247$.
• Step 3, compare: the usual SE is $S/\sqrt n = \sqrt{7/3} = 1.528$ (with $S^2 = 14/2 = 7$). The ratio is $\sqrt{2/3} = 0.816$: the narrowness bias.
• Step 4, who is in a resample: each point appears with probability $1 - (2/3)^3 = 19/27 = 0.704$.
• Step 5, the interval: the extremes 1 and 6 each have probability $1/27 = 0.037 > 0.025$. So the 2.5% and 97.5% quantiles of the exact bootstrap distribution are 1 and 6, and the 95% percentile interval is the whole range of the data. With $n = 3$ the bootstrap has almost nothing to work with.

B. An A/B test from plan to decision. A new ranking model should raise the conversion rate from the baseline $p_1 = 4.0\%$ to at least $p_2 = 4.5\%$ (absolute MDE $\delta = 0.005$), with a two-sided test at $\alpha = 0.05$ and power 0.8 ($z_{0.975} = 1.960$, $z_{0.8} = 0.842$).

Plan. $\bar p = 0.0425$, so $s_0 = \sqrt{2\bar p(1 - \bar p)} = \sqrt{0.08139} = 0.2853$, and $s_1 = \sqrt{p_1(1-p_1) + p_2(1-p_2)} = \sqrt{0.0384 + 0.04298} = 0.2853$. So

The rule of thumb gives $16 \cdot 0.04 \cdot 0.96/0.005^2 = 24{,}576$. With 4,000 new users per arm per day, the test needs 6.4 days. Plan 7 days: a full week, and some time for novelty effects to fade. At 28,000 users per arm the power is 0.83. Fix this plan before the start.

Check the split. After 7 days, control has 28,112 users and treatment 27,888. Under 50/50, $E = 28{,}000$ each: $X^2 = 2 \cdot 112^2/28{,}000 = 0.896$ and $p = 0.34$. No sample-ratio mismatch.

Analyse once. Control converts 1,125 users ($\hat p_1 = 4.002\%$), treatment 1,255 ($\hat p_2 = 4.500\%$).

• Pooled rate: $\hat p = 2{,}380/56{,}000 = 0.0425$.
• Standard error under $H_0$: $\operatorname{SE}_0 = \sqrt{0.0425 \cdot 0.9575 \cdot (1/28{,}112 + 1/27{,}888)} = 0.001705$.
• Statistic: $z = 0.004983/0.001705 = 2.92$, so $p = 0.0035$.
• 95% interval (the unpooled SE is also $0.001705$): $0.498 \pm 1.96 \cdot 0.1705 = 0.498 \pm 0.334$ percentage points, that is $[0.16, 0.83]$ pp. Relative to the baseline, this is a lift of $+12.5\%$ (about $+4\%$ to $+21\%$).
• Decision: launch the change, if the guardrail metrics (latency, errors, complaints) are fine.

With CUPED. Suppose the pre-period conversion of each user has correlation $\rho = 0.5$ with the outcome. The variance factor is $1 - 0.25 = 0.75$. The plan would then need $0.75 \cdot 25{,}551 \approx 19{,}164$ users per arm (4.8 days). Or keep 7 days and detect an MDE $\sqrt{0.75} = 0.87$ times as large.

5. In code #

SciPy and simulation check the claims of sections 2–4. The CUPED simulation and the block bootstrap are in the notebook.

import numpy as np
from scipy import stats

# exact permutation test: app (7, 9, 12) vs no app (3, 5, 6)
a, b = np.array([7.0, 9.0, 12.0]), np.array([3.0, 5.0, 6.0])
def diff(x, y, axis):
    return np.mean(x, axis=axis) - np.mean(y, axis=axis)
for alt in ("greater", "two-sided"):
    r = stats.permutation_test((a, b), diff, permutation_type="independent",
                               alternative=alt, n_resamples=np.inf)   # inf = all splits
    print(alt, round(r.statistic, 3), round(r.pvalue, 3))
# greater 4.667 0.05
# two-sided 4.667 0.1

# coverage of 95% percentile bootstrap intervals, 1000 simulated data sets each
rng = np.random.default_rng(0)
cases = {"mean, Exp(1), n=30":   (lambda: rng.exponential(size=30), np.mean, 1.0),
         "median, Exp(1), n=60": (lambda: rng.exponential(size=60), np.median, np.log(2)),
         "max, U(0,1), n=50":    (lambda: rng.uniform(size=50), np.max, 1.0)}
for name, (draw, stat, truth) in cases.items():
    hits = 0
    for _ in range(1000):
        ci = stats.bootstrap((draw(),), stat, n_resamples=999, method="percentile",
                             rng=rng).confidence_interval
        hits += ci.low <= truth <= ci.high
    print(name, hits / 1000)
# mean, Exp(1), n=30 0.928
# median, Exp(1), n=60 0.95
# max, U(0,1), n=50 0.0

# users per arm for 4.0% -> 4.5%, and an SRM check
z = stats.norm.ppf([0.975, 0.8])
p1, p2 = 0.040, 0.045
pb = (p1 + p2) / 2
n = (z[0] * np.sqrt(2 * pb * (1 - pb)) + z[1] * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2 / (p2 - p1) ** 2
print(np.ceil(n), 16 * p1 * (1 - p1) / (p2 - p1) ** 2)   # 25551.0 24576.0
print(stats.chisquare([50_600, 49_400]))                  # statistic 14.4, pvalue 0.00015

With 1,000 data sets, each coverage estimate has a standard error of about 0.007. The mean of skewed data at $n = 30$ undercovers a little (0.928 instead of 0.95), the median is fine, and the maximum is never covered.

Peeking in an A/A test, and CUPED with a good and a bad covariate:

import numpy as np

# A/A test (H0 true) with k equally spaced looks; stop at the first |z| > 1.96
rng = np.random.default_rng(4)
for looks in (1, 2, 5, 10, 20, 100):
    S = np.cumsum(rng.normal(size=(100_000, looks)), axis=1)   # data arrive in equal batches
    Z = S / np.sqrt(np.arange(1, looks + 1))                    # z-statistic at each look
    print(looks, np.mean(np.any(np.abs(Z) > 1.96, axis=1)))
# 1 0.04946   2 0.08159   5 0.14084   10 0.19272   20 0.24672   100 0.36978

# CUPED: 2000 simulated experiments, n users per arm, true effect tau
rng = np.random.default_rng(5)
reps, n, tau = 2000, 5000, 0.1
treat = np.arange(2 * n) >= n                     # users are i.i.d., so a fixed split is random
X = rng.gamma(2.0, 1.0, size=(reps, 2 * n))       # pre-experiment activity
Y = 1.0 + X + rng.normal(0.0, 1.44, size=X.shape) + tau * treat     # corr(X, Y) = 0.70
Z = X + 0.3 * treat + rng.normal(0.0, 0.2, size=X.shape)            # in-experiment: moved by the treatment

def estimate(Y, cov=None):
    if cov is not None:                           # CUPED with one theta from the pooled data
        c, y = cov - cov.mean(1, keepdims=True), Y - Y.mean(1, keepdims=True)
        theta = (c * y).sum(1) / (c * c).sum(1)
        Y = Y - theta[:, None] * c
    return Y[:, treat].mean(1) - Y[:, ~treat].mean(1)

for name, d in (("plain", estimate(Y)), ("CUPED with X", estimate(Y, X)), ("CUPED with Z", estimate(Y, Z))):
    print(f"{name:13s} mean = {d.mean():+.4f}   SD = {d.std():.4f}")
# plain         mean = +0.1016   SD = 0.0406
# CUPED with X  mean = +0.1007   SD = 0.0292
# CUPED with Z  mean = -0.1914   SD = 0.0291

With the pre-period covariate, the estimate stays centred on $\tau = 0.1$, and the standard deviation falls by the factor $0.0292/0.0406 = 0.72 \approx \sqrt{1 - 0.70^2}$. The in-experiment covariate is just as precise, but it reports a negative effect: the correction removed the real effect.

6. Where it shows up #

• Evaluating ML models. A bootstrap interval over test examples for accuracy, AUC or a benchmark score; a paired bootstrap or a sign-flip permutation test for the difference of two models. Leaderboards of many models on one test set also need a multiple-testing view.
• Online experiments. New ranking, recommendation and search models are launched through A/B tests with a primary metric, guardrail metrics (latency, errors) and SRM checks. Rankers are also compared by interleaving the result lists of two models for the same user, which needs far fewer users.
• Bagging and random forests. Each tree is trained on a bootstrap resample; the 36.8% of points left out give the out-of-bag error.
• Weather and climate. Differences in forecast skill between two systems are paired by date and autocorrelated in time, so use a block bootstrap with blocks of a week or more.
• Science. Randomised field trials and lab experiments use permutation tests for small studies, and covariate adjustment (ANCOVA, the CUPED idea) with baseline measurements.

7. Common mistakes #

• $C/B$ instead of $(1 + C)/(1 + B)$. It can report $p = 0$ and can reject too often.
• Too few permutations for a small threshold. With $B = 999$ the smallest p-value is 0.001, so a Bonferroni threshold of $10^{-5}$ can never be reached.
• Permuting when $H_0$ does not give exchangeability: unequal variances, paired data (flip signs instead), time series (permute blocks).
• Resampling the wrong unit: page views instead of users, tokens instead of documents, time steps instead of blocks.
• The percentile bootstrap for a maximum, a minimum or an extreme quantile, or for the mean of heavy-tailed data.
• Thinking the bootstrap creates information. It estimates the sampling distribution from the one sample you have; with $n = 3$ it cannot do much.
• Peeking: stopping at the first $p < 0.05$, extending a test "a few more days" because $p = 0.06$, or stopping before a full week.
• Ignoring an SRM because the imbalance looks small: 1.2% of 100,000 users is far outside chance.
• CUPED with a covariate measured after the start, or with a different $\theta$ in each arm.
• Mixing units in the sample size: a relative MDE with an absolute $\sigma$, or users per arm vs users in total.

8. Check yourself #

1. Group A has 4 values and group B has 4 values, all 8 distinct, and A holds the 4 largest. What are the exact one-sided and two-sided permutation p-values for the difference in means?
   Answer

   There are $\binom{8}{4} = 70$ equally likely splits. Only the observed one reaches $T \ge T_{\text{obs}}$, so the one-sided $p = 1/70 = 0.014$. The mirror split (B holds the 4 largest) gives $-T_{\text{obs}}$, so the two-sided p-value is $2/70 = 0.029$.
2. A Monte Carlo permutation test with $B = 999$ finds 23 permutations with $T^* \ge T_{\text{obs}}$. Give the p-value and its Monte Carlo standard error. Why not $23/999$?
   Answer

   $p = 24/1000 = 0.024$, with Monte Carlo SE $\approx \sqrt{0.024 \cdot 0.976/1000} = 0.005$. The observed labelling is itself one of the equally likely permutations. Counting it gives a p-value with $P(p \le \alpha) \le \alpha$ for every $B$, and never $p = 0$.
3. Why does the percentile bootstrap fail for $\theta$ in $\text{Uniform}(0, \theta)$, and which interval should you use?
   Answer

   Every resample maximum is at most the sample maximum $M$, and $M < \theta$, so the interval never covers $\theta$. The bootstrap distribution also puts mass $1 - (1 - 1/n)^n \approx 0.632$ on $M$ itself, while the true distribution of $M$ is continuous. Use the pivot: $P(M/\theta \le u) = u^n$, so $P(\alpha^{1/n} \le M/\theta \le 1) = 1 - \alpha$, and the interval $[M,\ M/\alpha^{1/n}]$ has coverage exactly $1 - \alpha$.
4. Baseline conversion 20%, MDE 2 percentage points (to 22%), two-sided $\alpha = 0.05$, power 0.8. How many users per arm? Compare with the rule of thumb.
   Answer

   $\bar p = 0.21$, $\sqrt{2 \cdot 0.21 \cdot 0.79} = 0.5760$ and $\sqrt{0.16 + 0.1716} = 0.5758$. So $n = (1.960 \cdot 0.5760 + 0.842 \cdot 0.5758)^2/0.02^2 = 1.6136^2/0.0004 \approx 6{,}510$ per arm. Rule of thumb: $16 \cdot 0.16/0.0004 = 6{,}400$.
5. A pre-period covariate has correlation 0.8 with the outcome. By what factors do CUPED's variance, standard error and required sample size change? Why is the estimate still unbiased?
   Answer

   The variance and the sample size both scale by $1 - 0.64 = 0.36$ (2.8 times fewer users); the standard error by $\sqrt{0.36} = 0.6$. It is unbiased because the covariate was measured before randomisation, so $\mathbb{E}[\bar X_T] = \mathbb{E}[\bar X_C]$ and the correction term has mean zero.
6. Users are randomised, but the metric is the click rate per page view; each user has about 5 views, with intra-user correlation 0.3. What goes wrong with the naive analysis?
   Answer

   It treats the views as independent. The true variance is larger by $1 + (5 - 1) \cdot 0.3 = 2.2$, so the naive standard error is too small by the factor $\sqrt{2.2} = 1.48$, and the p-values are too small. Analyse per user, or use the delta method or a user-level bootstrap.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Why is a permutation test exact under $H_0$?

   Show answer

2. Give the Monte Carlo permutation p-value with $B$ shuffles and $C$ shuffles at least as extreme as observed.

   Show answer

3. State the bootstrap principle in one sentence.

   Show answer

4. What share of the distinct data points does one bootstrap resample contain, for large $n$?

   Show answer

5. How is the 95% percentile bootstrap interval computed?

   Show answer

6. Why does the percentile bootstrap fail for the maximum of $\text{Uniform}(0, \theta)$ data?

   Show answer

7. What is the rule of thumb for the users per arm of an A/B test ($\alpha = 0.05$, power 0.8)?

   Show answer

8. What is peeking, and what does it do to an A/B test?

   Show answer

9. By what factor does CUPED reduce the variance, and with which $\theta$?

   Show answer

10. Why must the CUPED covariate be measured before the experiment?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

See also the interview track, A/B testing.

• "Design an A/B test for a new ranking model." — The primary metric and guardrails, the unit of randomisation, the MDE and the sample size, the duration in whole weeks, SRM and A/A checks, one pre-registered analysis, no peeking.
• "The test is significant after 3 days. The PM wants to launch." — Peeking inflates the type I error; the planned power assumed the full duration; novelty effects fade. Wait for the planned end, or use a sequential method that was planned for early stopping.
• "Explain the bootstrap. When does it fail?" — Resample with replacement from the data to mimic sampling from $F$. It fails for extremes, heavy tails, dependent data and very small $n$.
• "Permutation test or t-test?" — The permutation test is exact under exchangeability and needs no normality; the t-test is cheaper and gives an interval directly. With large $n$ they agree.
• "How do you reduce variance in an experiment?" — CUPED or regression adjustment, stratification, analysing only the triggered users, a less noisy metric (conversion instead of revenue, capped revenue).
• "Model B beats model A by 0.4% accuracy on 10,000 test examples. Is it real?" — Paired bootstrap of the per-example differences for an interval, or McNemar's test; also check the spread over training seeds.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Write permutation_test (difference in means, $p = (1 + C)/(1 + B)$) and check it on the exact six-student example.
2. Write bootstrap_percentile_ci, compare it with scipy.stats.bootstrap(method="percentile"), and measure its coverage for the mean, the median and the maximum.
3. Write a moving block bootstrap for the mean of an AR(1) series and compare it with the i.i.d. bootstrap and the true standard error.
4. Write ab_sample_size and two_proportion_power, and check the plan of section 4B by simulating whole A/B tests.
5. Write cuped and check on simulated A/B data that the variance falls by about $1 - \rho^2$ and that the estimate stays unbiased.

11. Watch and read #

Bootstrapping Main Ideas!!!StatQuest with Josh Starmer · 9 min

Resampling with replacement, step by step, and how the histogram of resampled means gives an interval.

Power Analysis, Clearly Explained!!!StatQuest with Josh Starmer · 17 min

How the effect size, the noise and the power decide the sample size, as in the A/B sample-size figure.

Statistical Learning: 5.4 The BootstrapStanford Online · 12 min

The bootstrap world vs the real world, from the authors of An Introduction to Statistical Learning.

What Exactly is CUPED?Eppo by Datadog · 2 min

A two-minute summary of CUPED: use each user's pre-experiment data to reduce the variance.

• Tim Hesterberg, What Teachers Should Know about the Bootstrap (2015): the permutation test, the bootstrap idea, how many resamples, and the accuracy of bootstrap intervals.
• Ron Kohavi, Roger Longbotham, Dan Sommerfield and Randal Henne, Controlled experiments on the web: survey and practical guide (2009): sample size, the $16\sigma^2/\Delta^2$ rule, A/A tests, and limitations.
• Alex Deng, Ya Xu, Ron Kohavi and Toby Walker, Improving the sensitivity of online controlled experiments by utilizing pre-experiment data (CUPED, 2013).
• Evan Miller, "How not to run an A/B test": why peeking breaks the significance level.