The Statistics of Evaluation
Real Signal, or Sampling Noise?

Chapter 1: The Eval Score IS a Random Variable

We watched the score jitter with our own eyes. Now we put a number on the jitter — and derive, from scratch, why quadrupling your test set only halves it. This one formula is the workhorse of the entire lesson.

Reframe the eval as a sequence of coin flips. On each test example the model is either right (1) or wrong (0). The accuracy is simply the average of those N zero-or-one outcomes. That single reframing is the whole trick: accuracy is a sum of random things, and sums of random things have completely understood behavior.

Before we go further, separate two quantities people constantly confuse. The estimand is p, the model's TRUE accuracy over the whole distribution of possible questions — a number you can never observe. The estimate is p-hat (written p̂), the accuracy you actually measured on your N questions — the number on your dashboard. The estimate wobbles around the estimand. Statistics is the discipline of saying how much.

Derive the variance of a single flip from scratch. Let X be 1 with probability p and 0 otherwise — a Bernoulli variable. Its mean is E[X] = p. Its variance is E[X²] − (E[X])². Because X is 0 or 1, X² = X, so E[X²] = p. Therefore the variance is p − p² = p(1−p).

Plug in a real number to feel it: at p = 0.85, the per-flip variance is 0.85 − 0.7225 = 0.1275. Notice this quantity peaks at p = 0.5 (a 50/50 coin is maximally uncertain) and shrinks toward 0 as p approaches 0 or 1 (a model that is right 99% of the time barely varies). Uncertainty is largest exactly when the model is most undecided.

Now average N independent flips. Two facts do all the work. First, variances of independent variables add: the variance of the sum of N flips is N × p(1−p). Second, dividing a quantity by a constant c divides its variance by c². We divide the sum by N to get the average, so we divide the variance by N². Net: Var(p̂) = N p(1−p) / N² = p(1−p)/N.

The standard error — the typical distance p-hat lands from p — is the square root of that variance:

Stare at the N under the square root, because it governs the economics of every eval you will ever run. Because N sits inside a square root, precision improves only as 1/√N. To halve your error bar you must quadruple your test set. To get 10× the precision you need 100× the examples — and 100× the labeling or inference cost. This single fact is why eval precision is genuinely expensive, and why Chapter 8's compute budgeting matters so much.

Why is the sampling distribution shaped like a bell? The Central Limit Theorem says that for a reasonable N, the average of many independent things is approximately normal — a bell curve centered at the true p with width SE. That is why the "±2 SE covers about 95%" rule works. Hold onto the word approximately: in Chapter 2 we will watch this approximation break when p is extreme or N is small.

Plug in real scales to build muscle memory. At 85% accuracy: SE is about ±1.6 points at N=500, ±1.1 at N=1000, ±0.5 at N=5000. Look at that first number and reread the last three papers you saw claim an "improvement." Many of those deltas are smaller than the ±1.6 that a 500-example benchmark simply cannot resolve.

One caveat that sets up several later chapters: this formula assumes the examples are independent draws. If your benchmark has ten questions all about the same document, those ten share a fate — get the document wrong and you likely miss all ten. That clustering inflates the true variance beyond p(1−p)/N. Anthropic's November 2024 error-bars note makes exactly this correction with clustered standard errors; the mental model is that your effective N is smaller than your nominal N.

python — from scratch, every intermediate
import numpy as np
outcomes = np.array([1]*425 + [0]*75)   # 425 right, 75 wrong
N = outcomes.size
phat = outcomes.mean()                    # 0.85
var1 = phat * (1 - phat)               # 0.1275  per-flip variance
var_mean = var1 / N                     # 0.000255
se = var_mean ** 0.5                   # 0.015969
print(phat, var1, round(se,6), round(1.96*se,4))
# Monte-Carlo check: empirical std of accuracies should match the formula
rng = np.random.default_rng(0)
sims = rng.binomial(N, phat, 100000) / N
print("empirical SE:", round(sims.std(),6))  # ~0.0160

python — the library one-liner
from scipy import stats
print(stats.sem(outcomes))          # standard error of the mean
# or the proportion shortcut: math.sqrt(p*(1-p)/N)

We can now quantify the jitter of a single score. The next step is to wrap it into the thing every eval report should show and almost none do: an honest interval.

Chapter 2: Confidence Intervals — Honest Brackets

We have a point score and a noise scale. Combine them and you get the single most under-used artifact in evaluation: the interval. A number without an interval is a claim without evidence attached; a number with one tells the reader how far to trust it.

Define a confidence interval (CI) by its GUARANTEE, not its formula. A 95% CI is a procedure: if you reran the whole experiment — drew a fresh test set — 100 times and built the interval each time, about 95 of those 100 intervals would contain the true p. The interval is the random thing; the truth is fixed. This is subtle and worth slowing down for.

The standard recipe, called the Wald or normal-approximation interval, is p̂ ± 1.96 × √(p̂(1−p̂)/N). Where does 1.96 come from? It is the point on the bell curve that leaves 2.5% of the area in each tail, so 95% sits in the middle. In your head, "2" is a fine shortcut for 1.96.

Now the reading that separates careful practitioners from casual ones. The RIGHT reading: "the procedure that produced this interval captures the truth 95% of the time." The WRONG reading: "there is a 95% probability the true accuracy is inside THIS interval." The distinction sounds pedantic — until you notice that a fixed interval either contains the fixed truth or it does not; there is no probability left once both are fixed. The 95% describes the method, across many hypothetical reruns.

Here is where the normal approximation quietly breaks: small N, or p near 0 or 1 — which is exactly the regime of the most interesting modern evals. A hard agentic benchmark where a model solves 17 of 20 tasks. A safety eval where the violation rate hovers near zero. In these corners the Wald recipe produces nonsense. Concretely, the Wald interval for 17/20 pokes above 100% accuracy — a bracket that includes impossible values.

The Wilson interval fixes this by asking a smarter question. Instead of approximating a bell curve around p-hat, it asks: "which true values of p would make my observation unsurprising?" — it inverts the test. The result has two beautiful properties. The center shrinks toward 0.5 (a built-in skepticism about extreme small-sample scores), and the width accounts for uncertainty in the variance itself. It stays inside [0, 1] by construction.

The rule of thumb for practice: Wald is fine when both N p̂ and N(1−p̂) exceed about 10; below that, or any time the interval nears 0 or 1, reach for Wilson. If you need guaranteed conservative coverage — a regulated safety claim — use Clopper-Pearson, which is wider but never under-covers. Wilson is the default in statsmodels and most eval libraries precisely because it behaves well across the whole range.

This is not academic hygiene — the field moved. OpenAI's SimpleQA-style and swe-bench-style reports and Anthropic's error-bars methodology (Nov 2024) publish CIs; when HuggingFace relaunched the Open LLM Leaderboard v2 in June 2024 it re-ran and normalized scores explicitly because small-N subscores were being over-read as rankings.

python — Wilson from scratch, every intermediate
import math
def wilson(k, n, z=1.96):
    phat = k / n
    zz = z * z
    center = (phat + zz / (2*n)) / (1 + zz/n)
    half   = z * math.sqrt(phat*(1-phat)/n + zz/(4*n*n)) / (1 + zz/n)
    return center - half, center + half
print(wilson(17, 20))     # (0.6396, 0.9476)
# and watch Wald misbehave:
phat, n = 0.85, 20
se = math.sqrt(phat*(1-phat)/n)
print(phat - 1.96*se, phat + 1.96*se)   # 0.6935 1.0065  (>1!)

python — the library one-liner
from statsmodels.stats.proportion import proportion_confint
print(proportion_confint(17, 20, method='wilson'))
# -> (0.6396..., 0.9476...)  same answer

Wilson handles proportions beautifully. But your dashboard also shows F1, BLEU, and judge win-rates — metrics with no clean closed-form SE. For those we need a tool that needs no formula at all: the bootstrap.

Chapter 3: The Bootstrap — Error Bars for Any Metric

The formula SE = √(p(1−p)/N) works for accuracy because accuracy is an average of independent 0/1 outcomes. But your dashboard also shows F1, BLEU, and judge win-rates. Where do their error bars come from? Chapter 1's derivation gives no answer — and that is the wall this chapter climbs over.

Why do those metrics resist the formula? F1 is a ratio of sums (2·precision·recall over their sum) — it does not decompose into a per-example average. BLEU is a geometric mean of clipped n-gram precisions multiplied by a brevity penalty, computed over the whole corpus at once. Neither is "the mean of N independent things," so no closed-form standard error exists in general.

The bootstrap idea, in one sentence: we cannot rerun reality to get fresh test sets, but we can fake fresh test sets by resampling our one test set with replacement. Each fake set is N draws from the empirical distribution — our best available stand-in for the true distribution — so the spread of the metric across fake sets mimics its spread across real reruns of the world.

Make "with replacement" concrete. Each resample draws N examples from the original N, but after drawing an example we put it back, so some examples get picked twice or three times and others get picked zero times. On average about 36.8% of the original examples are left out of any given resample — that number is (1−1/N)^N, which converges to 1/e as N grows. Every fake universe is a slightly reweighted version of your data.

The recipe, mechanically:

Why does this work, honestly rather than by hand-waving? The fluctuation of the metric across bootstrap resamples mimics its fluctuation across true resamples of the world, because the empirical distribution converges to the true one as data accrues. This is called the plug-in principle — we plug the empirical distribution in where the (unknown) true distribution belongs. It is a principled approximation, not magic.

Be clear about what the bootstrap cannot do. It never invents examples you did not collect. If your 500-example set contains no adversarial phrasings, no bootstrap interval will reflect uncertainty about adversarial phrasings. The bootstrap quantifies sampling noise within your distribution; it is silent about distribution shift. It reveals how little you know from the data you have — not what you failed to collect.

Resample the RIGHT unit. For a QA benchmark, resample questions. For a dialogue benchmark, resample whole conversations, not individual turns — turns within a conversation share user, topic, and context, so they are correlated. For a multi-document eval, resample documents. Resampling below the independence unit pretends you have more independent data than you do, and gives intervals that are too narrow — the same clustering trap from Chapter 1.

You will meet this everywhere in the wild: HELM publishes bootstrapped CIs; machine-translation shared tasks have used paired bootstrap resampling since Koehn (2004); lm-eval-harness exposes a bootstrap standard error for aggregate metrics. And the paired version of the bootstrap returns as our best model-comparison tool in the very next chapter.

python — a 12-line bootstrap, any metric
import numpy as np, math
data = np.array([1,1,0,1,1,1,0,1,1,0])
print(data.mean())                       # 0.70
for idx in ([3,0,7,7,2,9,4,1,5,8],
            [6,6,1,3,0,2,9,5,4,8],
            [0,1,2,3,4,5,6,7,8,9]):
    print(data[np.array(idx)].mean(), end=' ')  # 0.8 0.6 0.7
print()
rng = np.random.default_rng(42)
boots = np.array([rng.choice(data, size=10, replace=True).mean()
                  for _ in range(10000)])
print(np.percentile(boots,2.5), np.percentile(boots,97.5),
      round(boots.std(),3), round(math.sqrt(0.7*0.3/10),3))
# -> 0.4 1.0 0.144 0.145

python — the library one-liner
from scipy.stats import bootstrap
res = bootstrap((data,), np.mean, n_resamples=10000,
                confidence_level=0.95, method='percentile')
print(res.confidence_interval)      # ~ (0.4, 1.0)

The bootstrap gives one model's error bar for any metric. But the question that actually decides launches is a comparison — and there is a magic trick hiding in the fact that both models took the same exam.

Chapter 4: Comparing Two Models — Paired Beats Unpaired

Here is a magic trick. The same 2-point gap, the same 1,000 examples — one test says p=0.22, the other says p=0.023. Nothing changed but the method. The difference is simply noticing that both models took the same exam. This chapter is that trick, fully explained.

First, hypothesis testing honestly, in one beat. The null hypothesis is "the two models are equally accurate." The p-value is the probability of seeing a gap at least this large IF the null were true. It is NOT the probability the null is true, and 0.05 is a social convention, not a law of physics.

The pairing insight. When models A and B answer the SAME questions, most questions come out identical — both right on the easy ones, both wrong on the hard ones. That shared component is pure noise for the comparison: it tells you the questions were easy or hard, not which model is better. Pairing subtracts it out. An unpaired test pays for that shared difficulty twice, once in each model's variance.

Build the 2×2 disagreement table. Label each question by how the two models did: both-correct (a), A-only-correct (b), B-only-correct (c), both-wrong (d). Here is the crucial fact — the entire comparison lives in b versus c, the discordant cells. The a and d cells, however enormous, contain zero information about which model is better.

McNemar's test from scratch. Under the null, each discordant example is a fair coin — equally likely to favor A or B — so b should be about (b+c)/2. The statistic (|b−c|−1)² / (b+c) measures how far the observed split is from fair, and it follows a chi-square (χ²) distribution with 1 degree of freedom. The "−1" is a continuity correction — a small adjustment for using a smooth curve to approximate discrete counts. When b+c is small, skip the approximation entirely and use the exact binomial tail.

Now run the magic trick with numbers. 1,000 examples: a=810, b=45, c=25, d=120. Marginals: A = 85.5%, B = 83.5%. The unpaired two-proportion z-test gives z=1.24, p=0.216 — "not significant." McNemar on the very same data: χ²=5.16, p=0.023 — significant. Same data, same 2-point gap. The paired test simply refuses to be distracted by the 930 agreements.

The paired bootstrap generalizes this beyond accuracy. Resample QUESTION INDICES (with replacement), recompute BOTH models' metrics on the same resampled indices, and record the delta. The percentile CI of those deltas — or the fraction of resamples where the delta flips sign — is your answer. It works for F1, BLEU, win rate, anything. The key is resampling the same indices for both systems, so the shared difficulty cancels.

Practical guidance. Same test set → always pair: McNemar for accuracy on 0/1 outcomes, paired bootstrap for everything else. Different test sets → you are stuck with unpaired, and you need far more data. And always report the CI on the DELTA, not just a p-value star, so readers see effect size and uncertainty in the same glance.

Field note: Anthropic's Adding Error Bars to Evals (Nov 2024) explicitly recommends paired analyses for model comparisons on shared eval sets, and MT/LLM tooling (sacrebleu, lm-eval-harness) implements paired bootstrap for exactly this reason. Judge-based win rates get their own paired treatment in GenAI Eval.

python — McNemar from scratch, then paired bootstrap
import math, numpy as np
from scipy import stats
b, c = 45, 25
chi2 = (abs(b - c) - 1)**2 / (b + c)      # 5.157
p_chi = stats.chi2.sf(chi2, 1)                # 0.0232
p_exact = 2 * stats.binom.sf(max(b,c)-1, b+c, 0.5)  # 0.0225
print(round(chi2,3), round(p_chi,4), round(p_exact,4))
# paired bootstrap on the delta (works for any metric):
rng = np.random.default_rng(0)
A = np.array([1]*855 + [0]*145)   # stand-in per-example correctness
B = np.array([1]*835 + [0]*165)
idx = np.arange(1000)
deltas = [ (A[s].mean() - B[s].mean())
           for s in (rng.choice(idx, 1000, replace=True) for _ in range(2000)) ]
print(np.percentile(deltas, [2.5, 97.5]))    # CI on the delta

python — the library one-liner
from statsmodels.stats.contingency_tables import mcnemar
tbl = [[810, 45], [25, 120]]
print(mcnemar(tbl, exact=True))     # pvalue ~ 0.0225
# or: scipy.stats.binomtest(45, 70, 0.5)

Pairing wrings the most out of the data you have. But no test can see a gap it is too small to resolve. How many examples do you actually need? That is power — Chapter 5.

Chapter 5: Power & Sample Size — Can Your Eval Even See the Gap?

Chapter 4 taught you to test properly. But a proper test on too few examples is a microscope with the lens cap on — you will see nothing and conclude "no difference." This chapter is about whether your eval could ever have seen the effect you are hunting.

Two ways a test can fail. A false alarm — declaring a gap that is not there — happens with probability α (alpha), which we set at 0.05 by convention. A miss — failing to detect a gap that IS there — happens with probability β (beta). The power of a test is 1−β: the chance you catch a real effect. Convention aims for 80% or more.

Build the intuition geometrically. Under the null, the observed gap is a bell curve centered at 0 with width SE_diff. Under the alternative (a real gap δ), it is a bell curve centered at δ. Power is an overlap question: how much of the alternative's curve lands beyond the significance threshold. Small δ or big SE → the curves overlap heavily → detection is a coin flip.

Derive the sample-size formula from that picture. To reject reliably, the true gap δ must exceed roughly (z_α/2 + z_β) standard errors. Squaring and solving for n gives:

Every symbol earns its place. 1.96 (that is z_α/2) buys the 5% false-alarm rate. 0.8416 (z_β) buys 80% detection. The variance sum is the noise. And δ² in the denominator is the killer.

That δ² denominator is the headline of the whole chapter: halving the gap you want to detect QUADRUPLES the required test set. Detecting 5-point gaps is cheap; detecting 1-point gaps costs about 19,000 examples; detecting 0.5-point gaps costs about 78,000. This single fact explains most eval frustration you will ever feel.

Run the flagship numbers (below). To detect 85% vs 86% with 80% power you need about 19,500 examples per model. Now hold that against reality: MMLU's smallest subjects have around 100 questions, GSM8K's test split is 1,319 problems, HumanEval is 164. Per-subject "wins" on such slices are astrology.

Flip the formula into a diagnostic — the minimum detectable effect (MDE): given your N, what is the smallest gap you can see with 80% power? For N=1,000 per model at ~85% accuracy, the MDE is about 4.4 points. Any claimed improvement below your MDE deserves the response: "your eval cannot support that claim."

What underpowered comparisons do to a whole FIELD is worth pausing on. Only lucky-large deltas cross the significance bar, so published and leaderboard gaps are systematically inflated — a phenomenon called the winner's curse — and follow-up reproductions "mysteriously" shrink. This needs no fraud, just small N and selective attention. It sets up Chapter 6 directly.

Practical levers before you buy more data: pairing (Ch4) shrinks the effective noise and can cut required N dramatically; stratifying by difficulty helps; and for generative models, more samples per prompt reduces decode variance (Ch8). Online-experiment power with real users, and CUPED-style variance reduction, live in Experiment Design — same math, different battlefield.

python — required_n from scratch, with an empirical check
import math, numpy as np
from scipy import stats
def required_n(p1, p2, alpha=0.05, power=0.8):
    za = stats.norm.ppf(1 - alpha/2)          # 1.96
    zb = stats.norm.ppf(power)                # 0.8416
    varsum = p1*(1-p1) + p2*(1-p2)
    return (za + zb)**2 * varsum / (p2 - p1)**2
n = required_n(0.85, 0.86)
print(math.ceil(n))                          # 19458
# empirical: simulate at that N, fraction of z-tests that reject ~ 0.80
rng = np.random.default_rng(0); N = math.ceil(n); hits = 0
for _ in range(10000):
    a = rng.binomial(N, 0.86)/N; bb = rng.binomial(N, 0.85)/N
    se = math.sqrt(a*(1-a)/N + bb*(1-bb)/N)
    if abs(a-bb)/se > 1.96: hits += 1
print(hits/10000)                            # ~0.80

python — the library one-liner
from statsmodels.stats.power import NormalIndPower
from statsmodels.stats.proportion import proportion_effectsize
es = proportion_effectsize(0.86, 0.85)
print(NormalIndPower().solve_power(es, power=0.8, alpha=0.05))

Power protects a single comparison. But real teams do not run one test — they sweep twenty configs and report the best. That is a different, sneakier way to fool yourself, and it is Chapter 6.

Chapter 6: Multiple Comparisons & Benchmark Shopping

You tried 20 prompt variants and one beat baseline at p=0.03. Champagne? Before you pop it, compute the probability that pure noise hands you at least one such "win." The answer is going to be uncomfortable.

This is the garden of forking paths, quantified. Each test at α=0.05 has a 5% false-alarm rate. But run 20 independent tests where NOTHING works, and P(at least one false win) = 1 − 0.95²⁰ = 64%. The dashboard almost guarantees a champagne moment. The family of tests, not each individual test, is what needs the 5% guarantee.

Name the everyday versions of this sin, because you commit them without noticing: sweeping 20 hyperparameter configs and reporting the best; evaluating on 15 benchmark subsets and highlighting the 2 that improved; rerunning a flaky eval "to be sure" and keeping the good run; trying 5 judge prompts and picking the favorable one. All are multiple comparisons wearing different clothes.

Bonferroni, the blunt instrument. To keep the FAMILY false-alarm rate at 5% across m tests, require each individual test to clear α/m. For 20 configs: 0.05/20 = 0.0025. It is simple, assumption-free, and it works — but it is brutally conservative when m is large, because it guards against even a single false positive anywhere in the family.

Benjamini-Hochberg (BH), the practical instrument, controls the false discovery rate (FDR) instead: rather than "probably zero false wins," it promises "at most ~5% of my declared wins are false." Recipe: sort the m p-values ascending, find the LARGEST k such that p_(k) ≤ (k/m)·α, and reject hypotheses 1..k. Early p-values face a near-Bonferroni bar; later ones face progressively gentler bars — a rising staircase rather than a flat wall.

Seed variance is the same disease in disguise. Train the identical config with 5 random seeds and scores land at [84.1, 83.5, 84.8, 83.9, 84.4] — mean 84.14, std 0.49. Reporting your best seed (84.8) against a rival's mean is benchmark shopping against yourself. The honest unit of comparison is mean ± std over seeds. "Quantifying Variance in Evaluation Benchmarks" (Madaan et al., June 2024) measured exactly this and found seed noise swamps many claimed gaps.

The field-scale consequence, why leaderboard tops are inflated: hundreds of teams submit to the same benchmark, and the top slot is selected as the MAX over many noisy scores, so it is biased upward even if everyone is honest — the winner's curse from Chapter 5, now at ecosystem scale. The Leaderboard Illusion study (April 2025) documented the industrial version: labs privately testing many Chatbot Arena variants and surfacing only the best.

Defenses that actually work: preregister the ONE primary metric and test before running the sweep (demote everything else to exploratory); hold out a confirmation set the sweep never touches and validate the single winner on it once; apply BH to any screen you do report; and always report seed mean±std. The release-gate mechanics that operationalize this live in Regression Testing & Quality Gates.

python — Benjamini-Hochberg from scratch + seed report
import numpy as np
print(round(1 - 0.95**20, 4))            # 0.6415
def benjamini_hochberg(ps, alpha=0.05):
    order = np.argsort(ps)
    m = len(ps)
    passing = [i+1 for i, j in enumerate(order)
               if ps[j] <= (i+1)/m * alpha]
    kmax = max(passing) if passing else 0
    return set(order[:kmax])              # indices of rejections
ps = [0.001,0.004,0.012,0.021,0.041,0.09,0.13,0.22,0.45,0.78]
print(len(benjamini_hochberg(ps)))         # 3
seeds = np.array([84.1,83.5,84.8,83.9,84.4])
print(round(seeds.mean(),2), round(seeds.std(ddof=1),3))  # 84.14 0.493

python — the library one-liner
from statsmodels.stats.multitest import multipletests
reject, p_adj, _, _ = multipletests(ps, alpha=0.05, method='fdr_bh')
print(reject.sum())                       # 3

Every trap so far — jitter, weak tests, underpowered slices, selection bias — converges on one artifact everyone trusts and almost nobody error-bars: the leaderboard. Let's watch one lie in real time.

Chapter 7: Showcase — The Leaderboard Lie

Everything from Chapters 1–6 collides here, in the single artifact your whole industry stares at and almost never error-bars: the leaderboard. This is the showcase. YOU control the truth — the models' real skills — and you get to watch the rankings lie.

Set the scene: 8 models. In the default preset, all eight have IDENTICAL true accuracy 85%. Every "Run eval" draws a fresh N-example test set per model and ranks them. The podium reshuffles every single run. Gold medals are being awarded by a random number generator, and the medal ceremony looks exactly as authoritative either way.

Why are ranks noisier than scores? A rank is an order statistic — a max/ordering over several noisy quantities. Each pairwise ordering flips whenever the score gap between two models is small relative to SE_diff, and with 8 models there are 28 pairs, so SOME pair flips almost every run. The rank-1 slot is claimed by whoever got luckiest — the winner's curse from Chapter 5, made visible.

The CI-overlap lens makes it honest. Turn on confidence intervals and at N=500 all eight bars overlap heavily, saying "these are statistically indistinguishable" — yet the RANK column keeps confidently printing 1 through 8. The number that looks most authoritative (the rank) carries the least statistical support on the page.

Play the three sliders as three lessons. (1) N drives everything — churn melts as N grows (Chapter 1's law). (2) More models = more churn at fixed gaps (Chapter 6's multiplicity). (3) The true-gap slider shows the detectability threshold — gaps below about 2 SE_diff churn, gaps above lock in (Chapter 5's power).

Now the real world, 2024–2025. The Leaderboard Illusion study (April 2025) showed Chatbot Arena's dynamics rewarded labs that privately tested many variants and surfaced the best — Chapter 6's max-selection, industrialized. The Llama 4 Maverick incident (April 2025) — the Arena-submitted variant differed from the released model — showed rank-optimizing pressure in action. And GSM1k (May 2024) showed several models dropping many points on freshly rewritten test items: some leaderboard positions were measuring memorization, a bias no error bar catches.

What an honest leaderboard looks like: a score with a CI per model; ranks only where separations are significant (cluster ties into bands); paired tests on shared test sets for adjacent entries; multi-run averages for stochastic decoding. Some of this exists — HELM ships CIs; the Open LLM Leaderboard aggregates carefully — but rank-1 headlines still ignore all of it.

Carry three questions into every leaderboard you ever read: how big is N, how big are the gaps in SE units, and how many models or configs were tried? You now own the math behind all three.

python — the leaderboard simulator + flip probability
import math, numpy as np
from scipy import stats
for N in (500, 5000):
    sed = math.sqrt(0.86*0.14/N + 0.85*0.15/N)
    print(N, round(sed,6), round(stats.norm.cdf(-0.01/sed),3))
# 500 0.022267 0.327   |   5000 0.007041 0.078
rng = np.random.default_rng(0)
a = rng.binomial(500, 0.86, 100000)/500
b = rng.binomial(500, 0.85, 100000)/500
print(round(np.mean(b > a) + 0.5*np.mean(b == a), 3))  # ~0.328

The leaderboard treated "noise" as a single blob. But rerun the same model on the same test set with a different seed or temperature and the score still moves. Time to split the noise into its parts.

Chapter 8: Variance Decomposition — Where Does the Noise Come From?

So far "noise" meant test-set sampling. But rerun the SAME model on the SAME test set with a different training seed or a different temperature and the score moves anyway. The single blob called "noise" is actually a stack of independent sources, and the compute budget should attack the biggest one — not the most convenient one.

Enumerate the sources an LLM eval actually has. (1) Training-seed variance — retrain the same config, get a different model (weight init, data order, dropout). (2) Test-set sampling variance — Chapter 1's √(p(1−p)/N). (3) Decode-sampling variance — temperature > 0 makes the SAME model give different answers per run. Plus honorable mentions: prompt-template sensitivity, quantization/hardware nondeterminism, and judge variance for LLM-judged metrics.

The key mathematical fact that makes this tractable: independent noise sources add in VARIANCE, not in standard deviation. Total variance = var_seed + var_test + var_decode; total std = √ of the sum. The consequence is huge: the total is dominated by the largest term, and shrinking a small term is nearly invisible. The square-root arithmetic makes "fixing the wrong noise" provably futile.

How to MEASURE each term with a designed grid: train S seeds; for each seed run the eval R times with sampling on. Variance across seeds (of per-seed means) estimates seed noise; variance across repeats within a seed estimates decode noise; the analytic p(1−p)/N estimates test-set noise. This is a one-way random-effects ANOVA (analysis of variance) wearing ML clothes — Madaan et al. (June 2024) ran exactly this design across popular benchmarks and found seed variance rivals or exceeds many published deltas.

Averaging attacks each term differently — this is the budgeting lever. Averaging over s seeds divides var_seed by s (test/decode untouched). Enlarging the test set k-fold divides var_test by k. Averaging r decode runs (or r samples per prompt) divides var_decode by r. And each remedy has a different unit price: a seed costs a full training run; test examples cost labeling; decode repeats cost only inference.

The temptation of a deterministic eval and its price: temperature=0 sets var_decode to zero and makes evals reproducible — but greedy decoding measures a slightly different quantity than the sampled behavior users actually experience, and single-sample pass@1 at t=0 can misrank models whose sampled distributions differ. The sensible middle ground is sampling multiple generations and averaging (the pass@k machinery itself belongs to GenAI Eval).

Reporting discipline that follows directly: state WHICH sources your error bar includes. "acc 84.1 ± 0.5 (test-set CI, single seed, t=0)" and "acc 84.1 ± 1.2 (3 seeds × 5 decode runs, all sources)" are different claims — most papers publish the first and implicitly defend the second. Regression gates need the ALL-sources number or they will fire on seed noise (Regression Testing owns that machinery).

python — variance decomposition, four scenarios
import math
vs, vt, vd = 0.36, 1.00, 0.16          # seed, test, decode variances
print(round(math.sqrt(vs + vt + vd), 3))       # 1.233 — all sources
print(round(math.sqrt(vs/3 + vt + vd), 3))     # 1.131 — average 3 seeds
print(round(math.sqrt(vs + vt/4 + vd), 3))     # 0.877 — 4x test set
print(round(math.sqrt(vs + vt), 3))            # 1.166 — kill decode entirely

python — recover the components from a designed S×R grid
import numpy as np
rng = np.random.default_rng(0)
S, R = 8, 6                                # seeds, decode repeats per seed
seed_eff = rng.normal(0, 0.6, S)             # true seed std 0.6
grid = seed_eff[:, None] + rng.normal(0, 0.4, (S, R))
seed_means = grid.mean(axis=1)
print("seed var est:", round(seed_means.var(ddof=1), 3))     # ~ 0.36 + small
print("decode var est:", round(grid.var(axis=1, ddof=1).mean(), 3))  # ~ 0.16

You now own the full toolkit. The next chapter drills it into interview reflexes; the last files everything into one cheat sheet and points to where each thread continues.

Chapter 9: The Interview Arsenal

Eval-statistics questions are staff-interview gold because they expose whether a candidate reasons about uncertainty or just reads dashboards. This chapter drills the reflexes and then hands you ten worked Q&A and three debates to argue both sides of.

The meta-skill interviewers probe: given any eval claim, instinctively locate N, the noise scale, the pairing structure, and the number of comparisons tried. Rehearse the 10-second triage — "delta vs 2×SE_diff, paired or not, how many shots were taken?"

Reflex 1 — the back-of-envelope SE. Memorize SE = √(p(1−p)/N) and its worst case √(0.25/N) = 0.5/√N. Instant anchors: N=100 → ±5pt; N=1,000 → ±1.6pt; N=10,000 → ±0.5pt. Producing these without a calculator is the single highest-value move in this topic.

Reflex 2 — the pairing question. "Same test set?" If yes, the disagreement cells b and c are everything (McNemar); paired bootstrap handles any metric. Practice re-deriving (|b−c|−1)²/(b+c) from "discordant items are fair coins under the null" — interviewers love the derivation more than the formula.

Reflex 3 — power sizing on a whiteboard. n ≈ 8 × variance-sum / δ² for 80% power at α=0.05 (since (1.96+0.84)² ≈ 7.85 ≈ 8). One-point gap near 85% → n ≈ 8×0.25/0.0001 ≈ 20k. Say the number BEFORE they ask.

Reflex 4 — multiplicity instinct. Hear "we tried k variants" and immediately compute 1−0.95^k out loud; know Bonferroni (α/m) as the conservative bar and BH's staircase as the practical one; know the confirmation-set escape hatch.

Reflex 5 — variance accounting. Any single-run claim gets the three-source checklist (seed / test-set / decode); components add in variance; attack the biggest. Anchor: seed noise alone rivals many published gaps (Madaan et al., 2024).

How to structure answers at staff level: state the framework first ("this is a paired comparison with multiplicity"), do ONE crisp calculation, name the correction, then land on a DECISION recommendation. Interviewers grade the landing, not the formula recall.

python — the eval_triage helper you'd keep in your snippets
import math
from scipy import stats
def eval_triage(k1, n1, k2, n2, discordant=None):
    p1, p2 = k1/n1, k2/n2
    se_diff = math.sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)
    z = (p1 - p2) / se_diff
    print(f"unpaired: gap={p1-p2:.3f} SE={se_diff:.4f} z={z:.2f} p={2*stats.norm.sf(abs(z)):.3f}")
    if discordant:
        b, c = discordant
        p_ex = 2 * stats.binom.sf(max(b,c)-1, b+c, 0.5)
        print(f"paired (McNemar exact): {b} vs {c}, p={p_ex:.3f}")
eval_triage(180, 200, 174, 200, discordant=(14, 4))
# unpaired: gap=0.030 SE=0.0319 z=0.94 p=0.347
# paired (McNemar exact): 14 vs 4, p=0.031

Ten interview questions with model answers

Expand each for a staff-grade answer and a follow-up the interviewer will actually ask next.

Three debates — argue both sides

The reflexes are yours. The last chapter files every tool into a one-screen cheat sheet and shows where each thread continues across the family.

Chapter 10: Connections & The Cheat Sheet

You started with "84.3 vs 83.9 — real or noise?" You can now answer it five different ways. Let's file the tools, draw the decision tree, and point to where each thread continues.

Retell the lesson as one story: a score is a random variable (Ch1) → honest brackets via Wilson (Ch2) → brackets for any metric via bootstrap (Ch3) → comparisons via pairing (Ch4) → test sets sized via power (Ch5) → sweeps disciplined via multiplicity corrections (Ch6) → leaderboards read with rank-skepticism (Ch7) → noise budgeted via decomposition (Ch8). One toolkit, one thread.

The one-screen cheat sheet

The decision tree in prose. One model → CI (Wilson if small/extreme, bootstrap if the metric is exotic). Two models, same test set → McNemar / paired bootstrap. Two models, different test sets → unpaired z, and check power first. Many configs → correct for multiplicity or confirm on held-out data. Any single-run claim → variance decomposition.

Honest limitations of everything taught here. All of it quantifies SAMPLING uncertainty within a fixed distribution — none of it detects distribution shift, test-set contamination (GSM1k showed models dropping points on freshly rewritten items, May 2024), metric gaming, or a metric that measures the wrong thing entirely. Statistics tells you the measurement is precise; it cannot tell you it is the right measurement.

Where each thread continues

Further reading, all load-bearing above: Miller, Adding Error Bars to Evals (arXiv 2411.00640, Nov 2024); Madaan et al., Quantifying Variance in Evaluation Benchmarks (June 2024); Singh et al., The Leaderboard Illusion (April 2025); Zhang et al., GSM1k, A Careful Examination of LLM Performance on Grade School Arithmetic (May 2024); Koehn, Statistical Significance Tests for Machine Translation Evaluation (2004).

"An approximate answer to the right question is worth a good deal more than an exact answer to the wrong question." — John Tukey

An error bar is exactly that: approximate, honest, and worth more than any exact-looking point score.