Statistical Significance in Evaluation

A reported score is a measurement, and measurements are noisy. A classifier’s accuracy on a test set depends on which images happened to be drawn into that set; a network’s accuracy depends, in addition, on the random seed that set its initial weights, the order of its minibatches and its augmentations. Statistical significance in evaluation is the practice of measuring that noise and asking whether a difference between two models is larger than it: confidence intervals around a score, paired tests for two models scored on the same images, repeated training runs for the randomness of training, and corrections when many models or datasets are compared at once. In computer vision the differences that decide a comparison are often a few tenths of a point, and Bouthillier et al. (2021) found the variance of benchmark results to be of the same order as published improvements, so the question is a practical one. The roles of the training, validation and test sets are covered in Training, Validation, and Test Sets, and the binomial standard error of accuracy and McNemar’s test are introduced in Accuracy; this article covers where the variance comes from, how to put an interval around any metric, how to compare models, and how to report results.

Definition

An evaluation produces an estimate $\hat{\theta}$ of a quantity $\theta$, such as a model’s accuracy on the population of images it is meant for. A confidence interval at level $1 - \alpha$ is a range computed from the data by a procedure that, over repetitions of the experiment, contains $\theta$ with probability at least $1 - \alpha$ (typically 95%). A significance test of a null hypothesis $H_0$, such as “models A and B have equal accuracy”, returns a p-value: the probability, if $H_0$ were true, of a result at least as extreme as the one observed. A difference is called statistically significant at level $\alpha$ when $p < \alpha$; the rate of false detections under $H_0$ is the type I error, and the probability of detecting a true difference is the test’s power.

What counts as a repetition of the experiment has to be stated, because two different questions hide behind “is A better than B?”:

• Are these two trained models different? Here the models are fixed and the only randomness is the sample of test images. Intervals for accuracy, the bootstrap over test images, and McNemar’s test answer this question.
• Is training procedure A better than training procedure B? Here a repetition also retrains both models, with new seeds and ideally new data, so the randomness of training is part of the answer. Bouthillier et al. (2021) make the same distinction: when models are given rather than trained by the experimenter, the test data are the only source of variation left.

Intuition

Two archers each shoot ten arrows and one scores a point higher. Whether that says anything about their skill depends on how much a single archer’s score varies from one round of ten to the next; if it varies by five points, a one-point gap is noise. Evaluating a vision model is the same: the test set is one round, and a training run is one archer on one day. Confidence intervals ask how much the score would move on another round, seed spreads ask how much it would move on another day, and paired tests sharpen the comparison by having both archers shoot at the same targets, so that a target that is hard for both does not count against one of them.

Formal Definition

Two sources of variance. Let $\xi$ collect the random choices of training (initialization, data order, augmentation, dropout masks, nondeterministic arithmetic) and let $S$ be a test set of $n$ images drawn from the population. The test accuracy $\hat{a}(\xi, S)$ varies with both, and by the law of total variance

$$\operatorname{Var}[\hat{a}] = \mathbb{E}_{\xi}\big[\operatorname{Var}_{S}[\hat{a} \mid \xi]\big] + \operatorname{Var}_{\xi}\big[a(\xi)\big],$$

where $a(\xi)$ is the population accuracy of the network trained with $\xi$. The first term is test-set sampling, about $a(1-a)/n$ for independent test images; the second is training randomness, which no amount of test data removes.

Intervals for a proportion. With $k$ of $n$ test images correct and $\hat{a} = k/n$, the textbook (Wald) interval is $\hat{a} \pm z\sqrt{\hat{a}(1-\hat{a})/n}$, with $z = 1.96$ for 95%. Inverting the normal test with the standard error at the hypothesized value instead of the estimated one gives the Wilson score interval (Wilson, 1927; in the form of Brown, Cai & DasGupta, 2001, eq. 4):

$$\frac{k + z^2/2}{n + z^2} \;\pm\; \frac{z\sqrt{n}}{n + z^2}\sqrt{\hat{a}(1-\hat{a}) + \frac{z^2}{4n}} .$$

Bootstrap intervals. For a metric $M$ computed from $n$ test images, draw $B$ resamples of $n$ image indices with replacement, recompute $M^{*}_b$ on each, and take the $\alpha/2$ and $1 - \alpha/2$ quantiles of $M^{*}_1, \dots, M^{*}_B$ as the percentile interval (Efron & Tibshirani, 1993). To compare two models, compute $M_A^{*} - M_B^{*}$ on the same resample each time: this is the paired bootstrap.

Paired tests on one test set. Let $n_{01}$ count the test images that only model A classifies correctly and $n_{10}$ those that only B classifies correctly. Under $H_0$ of equal accuracy, each of the $m = n_{01} + n_{10}$ disagreements is equally likely to favour either model, so $n_{01} \sim \mathrm{Binomial}(m, \tfrac12)$, and the exact two-sided McNemar p-value is

$$p = \min\Big(1,\; 2\sum_{j=0}^{\min(n_{01}, n_{10})} \binom{m}{j} 2^{-m}\Big).$$

McNemar (1947) gave formulas for the significance of a difference between correlated proportions, together with a chi-square equivalent; the accuracy article gives that $\chi^2$ form. A paired permutation test generalizes the idea to any metric: under $H_0$ the labels “A” and “B” are exchangeable on each image, so the null distribution of the difference is obtained by swapping them at random, image by image.

Multiple comparisons. When $m$ hypotheses are tested, the chance that at least one true $H_0$ is rejected, the family-wise error rate (FWER), grows with $m$. The Bonferroni correction tests each at $\alpha/m$. Holm’s (1979) step-down procedure orders the p-values $p_{(1)} \le \dots \le p_{(m)}$ and rejects $H_{(1)}, H_{(2)}, \dots$ while $p_{(i)} \le \alpha/(m - i + 1)$, stopping at the first failure; it controls the FWER and rejects at least as much as Bonferroni. Benjamini and Hochberg (1995) control instead the expected proportion of false rejections among all rejections, the false discovery rate (FDR), by rejecting $H_{(1)}, \dots, H_{(i)}$ for the largest $i$ with $p_{(i)} \le i q / m$, which they proved for independent test statistics.

Many models on many datasets. For $k$ models on $N$ datasets, the Friedman test ranks the models on each dataset (rank 1 for the best), averages the ranks $R_j$, and compares

$$\chi^2_F = \frac{12N}{k(k+1)}\Big[\sum_{j} R_j^2 - \frac{k(k+1)^2}{4}\Big]$$

with a $\chi^2$ distribution with $k - 1$ degrees of freedom, an approximation Demšar gives for roughly $N > 10$ and $k > 5$; he recommends the less conservative $F$ statistic of Iman and Davenport derived from $\chi^2_F$. If it rejects, the Nemenyi post-hoc test declares two models different when their average ranks differ by at least the critical difference $\mathrm{CD} = q_\alpha\sqrt{k(k+1)/(6N)}$, with $q_{0.05} = 2.569$ for four models (Demšar, 2006, Table 5(a)).

Properties

Where the variance comes from

Bouthillier et al. (2021) measured the sources separately in five case studies, two of them in vision: VGG11 on CIFAR-10 and an FCN with an ImageNet-pretrained ResNet-18 backbone for PASCAL VOC segmentation. They randomized one source at a time over 200 seeds, holding the others fixed. Resampling the data (bootstrap training sets, with the left-out examples as test set) stood out as the largest source; weight initialization generally contributed less than half of that variance, on a par with the order in which stochastic gradient descent visits the data; and the choice of hyperparameters by random search, a perturbed grid search or Bayesian optimization contributed on average about as much as initialization. For the classification tasks, the standard deviation they observed when resampling data matched a binomial model of the test set, so the data-sampling term is largely the test set’s limited size. In their PyTorch 1.2.0 setup, convolutional models gave reproducible results only with cuDNN’s deterministic mode on and its benchmarking off, and the segmentation pipeline could not be made exactly reproducible at all; they measured the residual “numerical noise” with every seed fixed. The PyTorch documentation still says, as of version 2.14, that completely reproducible results are not guaranteed across releases, platforms, or CPU and GPU, and it describes cuDNN’s benchmark mode as a source of nondeterminism across runs.

Training randomness on its own is large in vision benchmarks:

• CIFAR-10. Picard (2021), in an arXiv report rather than a peer-reviewed paper, trained a nine-layer ResNet from the DAWNBench results with 10,000 seeds in a short schedule: accuracy on what the report calls the validation set averaged 90.02% with a standard deviation of 0.23, but ranged from 89.01% to 90.83%. With a longer schedule over 500 seeds, the mean was 90.70% and the range 90.14–91.41%; the spread did not narrow once training had converged.
• ImageNet fine-tuning. With 50 seeds and pretrained backbones (only the classification layer initialized from scratch), Picard’s standard deviations were about 0.1 points and the gaps between the best and worst seed about half a point: 75.25–75.69% for a supervised ResNet-50, 74.98–75.35% for a self-supervised one and 76.63–77.09% for a self-supervised ViT. The author notes that a gap of this size is widely treated as an improvement in the field, and that his accuracies are below the state of the art.
• ImageNet from scratch. Wightman, Touvron and Jégou (2021) trained a ResNet-50 with their procedure A2 under 100 seeds; the ResNet article gives the spread on the ImageNet validation set (standard deviation 0.10), and the seed 0 weights they report elsewhere were 0.13 points above the 100-seed mean. On ImageNet-V2, a test set of 10,000 images instead of 50,000, the standard deviation was 0.23 (67.39–68.69%), and the seed that was best on the validation set was not the best on ImageNet-V2. They read this as measurement noise and recommend reporting on several datasets, with a clear distinction between validation and test.

The last point connects the two terms of the decomposition: a smaller test set makes the seed-to-seed spread look larger, because each network’s particular mistakes on a few thousand images add their own noise. The example below shows the same effect.

Intervals for a single score

The Wald interval is the one most often written, but its coverage is poor. Brown, Cai and DasGupta (2001) showed that its actual coverage oscillates erratically with $n$ and with the true proportion: at $n = 100$ a nominal 95% interval covers 95.2% of the time when the true proportion is 0.106 but 91.1% when it is 0.107, and the problem is worst near 0 and 1, which is where the error rates of strong classifiers live. They recommend the Wilson interval or an equal-tailed Jeffreys interval for $n \le 40$ and the Agresti–Coull interval, which is centred like Wilson’s, for larger $n$. For accuracy on a few hundred images of a near-perfect classifier, or for the error rate of one rare class, the choice matters; for 10,000 images at 80% accuracy it does not. The accuracy article’s example computes Wilson’s interval with SciPy.

These intervals assume independent test images. When test images come in correlated groups (frames of one video, crops of one photograph, scans of one patient), the effective number of images is smaller and the binomial interval is too narrow, as the split article demonstrates. The bootstrap then has to resample whole groups, not images, so that each resample keeps the dependence of the original test set (our reasoning).

The bootstrap for metrics that are not averages

Accuracy is an average over images, but many vision metrics are not. Average precision ranks all detections of a class across the test set before integrating the precision–recall curve, and the mean intersection over union of a segmentation benchmark sums intersections and unions over the whole test set before dividing. No simple formula gives their standard error; the bootstrap does, provided each resample recomputes the metric exactly as the benchmark does, from the resampled images’ predictions and ground truth, rather than averaging per-image scores (our reasoning). The PASCAL VOC organizers used this approach for the 2012 challenge (Everingham et al., 2015, §3.5): they resampled the test images to obtain percentile intervals for each method’s AP or segmentation score, percentile intervals of the paired difference to test two methods against each other, and intervals for each method’s rank. In detection on the bicycle class, for example, the leading entry scored 54.5 AP with a 95% interval of 51.7–57.3 and a rank range of 1–2 (their Table 10); their summary tables mark the leading methods that are not significantly different from each other. They point out that unpaired intervals can overlap even when one method beats the other on every resample, because resamples vary in difficulty; the paired interval of the difference does not have this problem.

The percentile interval is the simplest bootstrap interval. SciPy’s scipy.stats.bootstrap (version 1.16) notes that it is rarely used in practice and defaults to the bias-corrected and accelerated (BCa) interval, which adjusts the quantiles for bias and skewness in the bootstrap distribution; it also offers the basic (reverse percentile) interval, and draws 9,999 resamples by default. For the AUC there is also a nonparametric method for comparing classifiers scored on the same images, due to DeLong and colleagues, described in the ROC AUC article.

Paired comparisons

Two models scored on the same images make correlated errors: an image that is blurred, occluded or mislabelled typically defeats both. Comparing their accuracies as if they came from independent samples puts that shared difficulty into the variance of the difference and wastes the pairing. McNemar’s test, the paired bootstrap and the paired permutation test all work with per-image differences, so the images both models get right or both get wrong contribute nothing. For a 0–1 metric such as accuracy the permutation test is McNemar’s exact test in Monte Carlo form, since swapping the labels only matters on the disagreements.

Dietterich (1998) compared five tests for deciding whether one learning algorithm beats another, measuring in experiments how often each detected a difference where none existed. In his results, the test for the difference of two proportions and the paired t-test over several random train/test splits had high type I error in some situations, and he advised never using them; the paired t-test over the folds of 10-fold cross-validation was somewhat elevated; McNemar’s test had low type I error; and his new 5×2cv test, a paired t-test built from five replications of twofold cross-validation, had acceptable type I error and slightly more power than McNemar’s. He recommended McNemar’s test for algorithms that can be run only once and 5×2cv for algorithms that can be run ten times. Why fold scores are not independent, and why no unbiased estimator of the variance of k-fold cross-validation exists, is the subject of cross-validation.

McNemar’s test treats the two trained models as fixed: it says whether these two networks differ on the population their test set was drawn from, not whether their training procedures would differ on another run. For deep networks trained once per paper, it is often the only test available; when training can be repeated, the runs themselves become the sample.

Significance is not effect size

With enough test images any difference, however small, becomes significant, and with too few a large one does not. A p-value therefore needs an effect size beside it: the difference itself with its confidence interval, which shows both how large the effect plausibly is and how precisely it is known. Bouthillier et al. (2021) frame the comparison of training procedures as the probability $P(A > B)$ that one run of A scores higher than one run of B, estimated over paired runs on random data splits, and call a result significant if the lower end of a bootstrap interval for $P(A > B)$ lies above 0.5 and meaningful if its upper end exceeds a threshold $\gamma$. They found $\gamma = 0.75$ to separate benchmark fluctuations from published improvements across their five case studies. In their simulations, comparing a single run of each procedure gave about 10% false detections and about 75% missed improvements, while comparing averages over 50 runs against a threshold set from typical published gains rarely made false detections but missed about 90% of real improvements; the probability criterion kept false detections near 5% and misses near 30%. They also recommend highlighting every method within the significance bounds of the best, not only the best.

Many comparisons

Every extra comparison is another chance for noise to look like an effect. Testing 20 variants against a baseline at $\alpha = 0.05$ yields one false “significant” result on average when none of them helps. Seed selection is the same problem in disguise: training many times and reporting the best run is a search over noise. Picard (2021) argues that trial-and-error development and the selection of papers by their scores amount to an implicit search over seeds across a field, and Wightman et al. (2021) avoided reporting the best intermediate checkpoint for the same reason, since picking the maximum over checkpoints would act like a seed selection. The Bonferroni and Holm corrections control the family-wise error; Demšar (2006) prefers Holm’s procedure, which is more powerful than the single-step Bonferroni–Dunn correction and makes no extra assumptions, and considers FDR methods less suitable for comparing algorithms because they require choosing an acceptable false discovery rate. Bouthillier et al. note that in competitions with hundreds of entrants the correction becomes so stringent that choosing a winner is partly arbitrary: a different test set might have reordered the top.

Many datasets

Claims of general superiority rest on several datasets, but scores on different datasets are not commensurable: one point on MNIST and one point on ImageNet are not the same amount of progress, and the scores’ variances differ. Demšar (2006) therefore recommends nonparametric tests whose sample is the set of datasets: the Wilcoxon signed-ranks test for two methods, which ranks the per-dataset differences by absolute size and compares the rank sums of positive and negative differences, and the Friedman test for several, followed by the Nemenyi test for all pairs or by Holm’s procedure against a control method. He introduced critical difference diagrams, which place the methods on an axis of average rank and join those whose ranks differ by less than the critical difference. When a new method is compared only with existing ones, comparisons against it as a control need fewer corrections and have more power than all pairwise comparisons. The weakness is sample size: Bouthillier et al., citing earlier work, note that machine learning papers typically include as few as three to five datasets, too few for these tests to have much power.

Reporting

The sources above agree on a short list. Picard (2021) asks authors to vary seeds, and if possible data splits, and to report the mean, standard deviation, minimum and maximum. Bouthillier et al. (2021) recommend randomizing as many sources of variation as possible, such as initialization, data order, augmentation and data splits, preferring multiple random splits to a single fixed test set, and accounting for variance in the decision. Wightman et al. (2021) report the final checkpoint rather than the best one and evaluate on more than one test set. Everingham et al. (2015) report bootstrap intervals and mark the methods that the leader does not significantly beat. To these, in our reading, belongs the test-set size and the unit of independence, so that readers can judge the binomial noise themselves; the number of runs, since the standard error of a mean over $k$ runs is the run-to-run standard deviation divided by $\sqrt{k}$ (with ImageNet’s 0.1-point seed deviation, five runs give a standard error of about 0.045 points); and whether a “test” score was measured on a benchmark’s validation split.

Examples

The script generates 16×16 grayscale images of four shapes (disc, square, cross, triangle) on a zero background, with floating-point intensities: each shape is shifted by up to two pixels, scaled by a contrast drawn from 0.5 to 1.0, and covered by Gaussian noise of standard deviation 0.8. Part 1 compares two fixed template matchers that differ only in blur: both correlate the image with every shifted template of each class and pick the best class, A without blurring and B after blurring image and templates with a Gaussian of $\sigma = 0.5$ pixels. Their population accuracies are measured on 200,000 images; the evaluation itself uses one test set of 1,000. Part 1 computes Wald and Wilson intervals and their exact coverage, the paired standard error and bootstrap interval of the difference, McNemar’s exact test and a paired permutation test, cross-checking the Wilson interval, the exact test and the bootstrap against SciPy. Part 2 trains a small network, one hidden ReLU layer on blurred pixels with plain stochastic gradient descent, with 32 or 64 hidden units, 20 times each with seeds 0–19 that set the initialization and minibatch order. Training and test images stay fixed; each network is also scored on 100,000 fresh images to approximate its true accuracy. The script runs in under ten seconds with NumPy 2.3.5 and SciPy 1.16.3, and two runs give identical output.

import numpy as np
from scipy import stats
from scipy.ndimage import gaussian_filter

K, S = 4, 16                                    # classes, image size (pixels)
yy, xx = np.mgrid[:S, :S] - (S - 1) / 2
SHAPES = np.stack([(xx**2 + yy**2 < 20).astype(float),                 # disc
                   ((abs(xx) < 4) & (abs(yy) < 4)).astype(float),      # square
                   ((abs(xx) < 1.5) | (abs(yy) < 1.5)).astype(float),  # cross
                   ((yy > -4) & (yy < 5) & (abs(xx) < (yy + 4) / 2)).astype(float)])  # triangle
SHIFTS = [(a, b) for a in range(-2, 3) for b in range(-2, 3)]

def make_images(n, rng):
    """16x16 grayscale shapes, intensities around 0: random shift of up to 2 pixels,
    contrast 0.5-1.0, additive Gaussian noise with standard deviation 0.8."""
    y = rng.integers(0, K, n)
    X = np.stack([np.roll(SHAPES[c], tuple(rng.integers(-2, 3, size=2)), axis=(0, 1)) for c in y])
    return rng.uniform(0.5, 1.0, (n, 1, 1)) * X + rng.normal(scale=0.8, size=(n, S, S)), y

def template_matcher(sigma):
    """Blur image and templates by sigma, pick the class whose shifted template correlates best."""
    T = np.stack([gaussian_filter(np.roll(s, d, axis=(0, 1)), sigma).ravel()
                  for s in SHAPES for d in SHIFTS])
    T -= T.mean(1, keepdims=True)
    T /= np.linalg.norm(T, axis=1, keepdims=True)
    def predict(X):
        F = gaussian_filter(X, (0, sigma, sigma)).reshape(len(X), -1)
        F -= F.mean(1, keepdims=True)
        return (F @ T.T).reshape(len(X), K, len(SHIFTS)).max(2).argmax(1)
    return predict

def wald(k, n, z=1.96):
    p = k / n
    h = z * np.sqrt(p * (1 - p) / n)
    return p - h, p + h

def wilson(k, n, z=1.96):                       # Wilson (1927) score interval
    p = k / n
    centre = (k + z**2 / 2) / (n + z**2)
    h = z * np.sqrt(n) / (n + z**2) * np.sqrt(p * (1 - p) + z**2 / (4 * n))
    return centre - h, centre + h

# ---- Part 1: two template matchers (no blur, blur 0.5) on one test set of 1,000 images
A, B = template_matcher(0.0), template_matcher(0.5)
rng = np.random.default_rng(0)
Xpop, ypop = make_images(200_000, rng)          # stands in for the population
pa, pb = np.mean(A(Xpop) == ypop), np.mean(B(Xpop) == ypop)
X, y = make_images(1000, rng)                   # the test set we actually have
ca, cb = A(X) == y, B(X) == y                   # per-image correctness
n, ka = len(y), int(ca.sum())
print(f"population accuracy  A {pa:.4f}  B {pb:.4f}  difference {pa - pb:+.4f}")
print(f"test accuracy (n={n})  A {ca.mean():.3f}  B {cb.mean():.3f}  difference {ca.mean() - cb.mean():+.3f}")
lo, hi = wald(ka, n); print(f"A: Wald 95% [{lo:.3f}, {hi:.3f}]", end="  ")
lo, hi = wilson(ka, n); print(f"Wilson 95% [{lo:.3f}, {hi:.3f}]", end="  ")
ci = stats.binomtest(ka, n).proportion_ci(0.95, method="wilson")
print(f"scipy Wilson [{ci.low:.3f}, {ci.high:.3f}]")

print("exact coverage of nominal 95% intervals (Wald / Wilson) for a true accuracy p:")
for p in (pa, 0.98):
    row = []
    for m in (20, 50, 100, 1000):
        k = np.arange(m + 1)
        pmf = stats.binom.pmf(k, m, p)    # probability of each possible number correct
        row += [f"n={m}: " + " / ".join(f"{pmf[(lo <= p) & (p <= hi)].sum():.3f}"
                                         for lo, hi in (wald(k, m), wilson(k, m)))]
    print(f"  p={p:.4f}  " + "   ".join(row))

# Standard error of the difference: ignoring the pairing vs using it.
d = ca.astype(float) - cb                       # per-image difference: -1, 0 or +1
se_unpaired = np.sqrt(ca.mean() * (1 - ca.mean()) / n + cb.mean() * (1 - cb.mean()) / n)
se_paired = d.std(ddof=1) / np.sqrt(n)
print(f"SE of A - B: treating the accuracies as independent {se_unpaired:.4f}, paired {se_paired:.4f}")

# Paired bootstrap: resample test images, recompute both accuracies on the same resample.
boot = rng.integers(0, n, (10_000, n))
diffs = d[boot].mean(1)
lo, hi = np.percentile(diffs, [2.5, 97.5])
out = [f"percentile [{lo:+.3f}, {hi:+.3f}]"]
for method in ("percentile", "BCa"):
    r = stats.bootstrap((ca.astype(float), cb.astype(float)), lambda a, b: np.mean(a - b),
                        paired=True, vectorized=False, n_resamples=10_000, method=method,
                        rng=np.random.default_rng(1))
    out.append(f"scipy {method} [{r.confidence_interval.low:+.3f}, {r.confidence_interval.high:+.3f}]")
print("paired bootstrap 95% CI of A - B: " + "  ".join(out))

# McNemar's exact test: only the images where the two classifiers disagree carry information.
n01, n10 = int(np.sum(ca & ~cb)), int(np.sum(~ca & cb))
m = n01 + n10
p_exact = min(1.0, 2 * stats.binom.cdf(min(n01, n10), m, 0.5))   # two-sided exact binomial
p_scipy = stats.binomtest(n01, m, 0.5).pvalue
print(f"McNemar: A right/B wrong {n01}, A wrong/B right {n10}; exact p = {p_exact:.4f} (scipy {p_scipy:.4f})")

# Paired permutation test: under H0 the labels A and B can be swapped on each image.
flips = rng.choice([-1.0, 1.0], size=(20_000, n))
null = np.abs((flips * d).mean(1))
p_perm = (1 + np.sum(null >= abs(d.mean()) - 1e-12)) / (1 + len(null))
print(f"paired permutation test (20,000 sign flips): p = {p_perm:.4f}")

# ---- Part 2: training randomness. A tiny MLP on blurred pixels, 32 vs 64 hidden units.
def train_mlp(X, y, seed, hidden, epochs=10, lr=0.05, batch=20):
    """One-hidden-layer ReLU network, softmax output, plain SGD.
    The seed sets the weight initialization and the order of the minibatches."""
    r = np.random.default_rng(seed)
    nn, dd = X.shape
    W1, b1 = r.normal(scale=np.sqrt(2 / dd), size=(dd, hidden)), np.zeros(hidden)
    W2, b2 = r.normal(scale=np.sqrt(1 / hidden), size=(hidden, K)), np.zeros(K)
    Y = np.eye(K)[y]
    for _ in range(epochs):
        for idx in np.array_split(r.permutation(nn), nn // batch):
            h = np.maximum(X[idx] @ W1 + b1, 0)
            z = h @ W2 + b2
            p = np.exp(z - z.max(1, keepdims=True)); p /= p.sum(1, keepdims=True)
            g = (p - Y[idx]) / len(idx)          # gradient of the mean cross-entropy wrt z
            gh = (g @ W2.T) * (h > 0)
            W2 -= lr * h.T @ g; b2 -= lr * g.sum(0)
            W1 -= lr * X[idx].T @ gh; b1 -= lr * gh.sum(0)
    return lambda Z: (np.maximum(Z @ W1 + b1, 0) @ W2 + b2).argmax(1)

blur = lambda X: gaussian_filter(X, (0, 1, 1)).reshape(len(X), -1)
rng = np.random.default_rng(1)
Xtr, ytr = make_images(2000, rng)               # one fixed training set
Xte, yte = make_images(1000, rng)               # one fixed test set
Xp, yp = make_images(100_000, rng)              # fresh images: each network's true accuracy
Ftr, Fte, Fp = blur(Xtr), blur(Xte), blur(Xp)
res = {}
for hidden in (32, 64):
    runs = []
    for seed in range(20):
        f = train_mlp(Ftr, ytr, seed, hidden)
        runs.append((np.mean(f(Fte) == yte), np.mean(f(Fp) == yp)))
    res[hidden] = np.array(runs)
print("\n20 seeds per model; same training and test images")
print("hidden  test acc: mean  sd     min-max       | fresh acc: mean  sd")
for h, r in res.items():
    print(f"{h:6d}  {r[:, 0].mean():.3f}  {r[:, 0].std(ddof=1):.4f}  {r[:, 0].min():.3f}-{r[:, 0].max():.3f}"
          f"  |  {r[:, 1].mean():.4f}  {r[:, 1].std(ddof=1):.4f}")
a = res[32][:, 0].mean()
print(f"binomial SE of one test accuracy (n=1000): {np.sqrt(a * (1 - a) / 1000):.4f}")
# Every single-run comparison: each of the 20 runs of one model against each of the other's.
pair = res[64][:, 0][:, None] - res[32][:, 0][None, :]
print(f"single-run comparisons on the test set (400 pairs): 64 higher {np.mean(pair > 0):.2f}, "
      f"tied {np.mean(pair == 0):.2f}, lower {np.mean(pair < 0):.2f}")
allr = np.vstack([res[32], res[64]])
print(f"correlation of test and fresh accuracy over the 40 networks: {np.corrcoef(allr.T)[0, 1]:.2f}")
for j, name in ((0, "test set"), (1, "fresh images")):
    t = stats.ttest_ind(res[64][:, j], res[32][:, j], equal_var=False)
    print(f"Welch t-test, 20 runs vs 20 runs, {name}: difference of means "
          f"{res[64][:, j].mean() - res[32][:, j].mean():+.4f}, p = {t.pvalue:.4f}")

Output:

population accuracy  A 0.8588  B 0.8426  difference +0.0163
test accuracy (n=1000)  A 0.845  B 0.829  difference +0.016
A: Wald 95% [0.823, 0.867]  Wilson 95% [0.821, 0.866]  scipy Wilson [0.821, 0.866]
exact coverage of nominal 95% intervals (Wald / Wilson) for a true accuracy p:
  p=0.8588  n=20: 0.948 / 0.948   n=50: 0.928 / 0.936   n=100: 0.943 / 0.940   n=1000: 0.953 / 0.949
  p=0.9800  n=20: 0.332 / 0.940   n=50: 0.635 / 0.922   n=100: 0.866 / 0.949   n=1000: 0.950 / 0.947
SE of A - B: treating the accuracies as independent 0.0165, paired 0.0085
paired bootstrap 95% CI of A - B: percentile [+0.000, +0.033]  scipy percentile [+0.000, +0.033]  scipy BCa [-0.001, +0.033]
McNemar: A right/B wrong 44, A wrong/B right 28; exact p = 0.0764 (scipy 0.0764)
paired permutation test (20,000 sign flips): p = 0.0779

20 seeds per model; same training and test images
hidden  test acc: mean  sd     min-max       | fresh acc: mean  sd
    32  0.639  0.0095  0.620-0.665  |  0.6385  0.0040
    64  0.643  0.0071  0.634-0.658  |  0.6438  0.0053
binomial SE of one test accuracy (n=1000): 0.0152
single-run comparisons on the test set (400 pairs): 64 higher 0.64, tied 0.06, lower 0.30
correlation of test and fresh accuracy over the 40 networks: 0.40
Welch t-test, 20 runs vs 20 runs, test set: difference of means +0.0042, p = 0.1277
Welch t-test, 20 runs vs 20 runs, fresh images: difference of means +0.0053, p = 0.0010

The unblurred matcher really is better, by 1.6 points on 200,000 images, and the test set happens to reproduce that gap almost exactly. Yet on 1,000 images the evidence is borderline: the paired bootstrap interval of the difference just touches zero and McNemar’s exact test gives $p \approx 0.08$, so a significance test at the 5% level would not detect a real improvement. This is a matter of power, not of the improvement being absent. The two classifiers disagree on only 72 images, and those disagreements are all the information the test set holds about their difference. Using the pairing halves the standard error of the difference compared with treating the two accuracies as independent; without it, the same data would look even less conclusive. The permutation test reproduces McNemar’s p-value up to Monte Carlo error, as it must for a 0–1 metric, and the three bootstrap intervals agree to within 0.001.

At 86% accuracy the Wald and Wilson intervals are nearly identical and both cover close to 95% of the time. At 98% the Wald interval fails on small test sets: on 20 images it covers the true accuracy only a third of the time, because a test set with no errors yields a zero-width interval at 100%, while Wilson’s interval stays between 92% and 95% at every size shown.

In Part 2 the network with 64 hidden units is better on fresh images by about half a point, a difference the 100,000-image evaluation detects clearly. On the fixed 1,000-image test set, however, single runs disagree about the order in 30% of the 400 pairings of runs, the 20-run averages are not significantly different, and the best run of the smaller network (66.5%) beats every run of the larger one. Two effects combine. Seeds change each network’s true accuracy only by about 0.4–0.5 points, but they also change which of the 1,000 test images it gets wrong, so the test-set accuracy varies more from seed to seed than the true accuracy does (0.95 against 0.40 points for the smaller network), the effect Wightman et al. saw between ImageNet’s 50,000 validation images and ImageNet-V2’s 10,000. Hence the weak correlation of 0.40 between a network’s test score and its fresh score: picking the best seed by its test score mostly picks test-set luck. And the Welch test on test-set scores accounts only for seed variance on this one test set; the binomial standard error of 1.5 points, larger than either seed spread, describes how much both models’ scores would move together on a different test set of the same size, which no number of seeds reveals.

Common Misconceptions

• “Overlapping confidence intervals mean no significant difference.” Separate intervals for two models scored on the same images ignore the correlation between their errors. In the example the paired standard error of the difference is half the unpaired one, and Everingham et al. note that unpaired intervals can overlap even when one method wins on every bootstrap resample. Compare models with an interval or test for the difference.
• “Not significant means no difference.” It means the test set or the number of runs could not detect it. The template matchers above differ by 1.6 points in truth and $p \approx 0.08$ on 1,000 images.
• “Significant means important.” A large enough test set makes a 0.05-point difference significant. Whether it matters is a separate judgment about effect size, which Bouthillier et al. formalize as the probability that one procedure outperforms the other in a run.
• “One seed is enough if the test set is large.” A large test set shrinks the first term of the variance decomposition but leaves the second. Picard’s 10,000 CIFAR-10 seeds span 1.8 points on the same evaluation images.
• “Significance on a benchmark means the model generalizes.” A test only quantifies noise relative to the population the test set was drawn from. If the community has reused that test set for years, or if deployment images differ from it, the interval says nothing about either problem; benchmark overfitting covers the first, and the second is a question of distribution shift.

Where It Is Used

Every benchmark comparison in vision faces these questions, from image classification on ImageNet, where the half-point gap between the best and worst of 50 fine-tuning seeds matches gaps the field treats as improvements (Picard, 2021), to detection and segmentation, where the bootstrap over test images gives intervals for average precision and mean IoU as PASCAL VOC did. Classifiers scored on the same medical images can be compared through their AUCs with DeLong’s method (ROC AUC). Small datasets are evaluated by cross-validation, whose fold scores need the tests designed for correlated estimates, such as Dietterich’s 5×2cv. Repeated training with different seeds belongs in the training workflow whenever two training recipes, augmentations or architectures are compared, as in the ResNet training-procedure studies.