Foundations

Bias–Variance Tradeoff

How the expected error of a learned predictor or estimator splits into noise, squared bias and variance, how the split is derived, computed in closed form for linear smoothers and estimated by resampling, why the tradeoff is a tendency rather than a law, and how it changes under the 0–1 loss of classification and segmentation.

intermediate

The bias–variance tradeoff describes how the expected error of anything estimated from noisy data splits into a part that comes from systematic error and a part that comes from sensitivity to the particular data seen. A denoising filter that averages over a wide window removes most of the noise but smears edges: it has low variance and high bias. A narrow window keeps the edges but leaves the noise: low bias, high variance. A classifier, a regression network or a k-nearest-neighbor rule faces the same choice whenever its flexibility is set, by the window width, the number of neighbours, the depth of a tree or the strength of a penalty. The bias–variance decomposition makes the split exact for the squared error: expected error = irreducible noise + squared bias + variance.

The machine learning article states the decomposition at one input in a paragraph, and Mean Squared Error gives its estimator form in one equation. Overfitting and Underfitting covers the diagnosis from training and validation curves. This article owns the decomposition itself: its derivation for predictors and estimators and how the two views connect, closed forms for linear smoothers, how bias and variance are estimated in practice, why the “tradeoff” is a tendency and not a law, and what happens to it under the 0–1 loss, where Domingos’s (2000) unified decomposition shows that variance can lower the expected error. An example computes every term exactly, and by Monte Carlo, for Gaussian smoothing of a noisy synthetic image, and then for the segmentation obtained by thresholding it.

The words have other meanings on this site. The bias term of a neuron or linear model, dataset bias and inductive bias are different things, as are the explained variance of PCA and the variance of stochastic gradients.

Definition

Let a learning procedure turn a training set T\mathcal{T} into a predictor f^(x;T)\hat{f}(\mathbf{x}; \mathcal{T}), and let the target at an input x\mathbf{x} be y=f(x)+εy = f(\mathbf{x}) + \varepsilon with E[ε]=0\mathbb{E}[\varepsilon] = 0 and Var⁡(ε)=σ2\operatorname{Var}(\varepsilon) = \sigma^2. The training set is random, so the prediction at a fixed x\mathbf{x} is a random variable. Its average over training sets, fˉ(x)=ET[f^(x;T)]\bar{f}(\mathbf{x}) = \mathbb{E}_{\mathcal{T}}[\hat{f}(\mathbf{x}; \mathcal{T})], defines the two terms:

  • the bias at x\mathbf{x} is fˉ(x)−f(x)\bar{f}(\mathbf{x}) - f(\mathbf{x}), the error of the average prediction;
  • the variance at x\mathbf{x} is ET[(f^(x;T)−fˉ(x))2]\mathbb{E}_{\mathcal{T}}\big[(\hat{f}(\mathbf{x}; \mathcal{T}) - \bar{f}(\mathbf{x}))^2\big], the spread of the predictions around their average.

Goodfellow, Bengio and Courville (2016, §5.4.4) put it in words: bias measures the expected deviation from the true value, and variance the deviation from the expected estimator value that a particular sample of the data is likely to cause. The tradeoff is the observation that changes which lower one term often raise the other.

Intuition

Imagine fitting the same model many times, each time to a fresh training set drawn the same way: new photographs of the same scene, new noise on the same image, a new random sample of labelled patches. At one input, the predictions form a cloud. Its centre may sit away from the truth, which is bias; it has a width, which is variance. A rigid model gives a tight cloud in the wrong place; a flexible model gives a cloud centred near the truth but spread wide. The squared error, averaged over the clouds, is the squared distance from the centre to the truth plus the squared width, plus whatever noise sits in the target itself.

Smoothing an image makes this concrete. Each output pixel of a Gaussian filter is a weighted average of the noisy pixels around it. Wide weights average many independent noise values, so the output barely changes from one noisy copy to the next (low variance), but they also mix in pixels whose true values differ, so the average sits off the truth wherever the image is not flat (bias). On a flat patch the wide filter loses nothing; at an edge or on fine texture it loses a great deal. The best width therefore depends on the content, which the example below measures.

Formal Definition

Predictors. For squared error at a fixed input x0\mathbf{x}_0, with the expectation over training sets and over the noise in a new target yy independent of the training set,

E[(y−f^(x0))2]=σ2+(fˉ(x0)−f(x0))2+ET[(f^(x0)−fˉ(x0))2].\mathbb{E}\big[(y - \hat{f}(\mathbf{x}_0))^2\big] = \sigma^2 + \big(\bar{f}(\mathbf{x}_0) - f(\mathbf{x}_0)\big)^2 + \mathbb{E}_{\mathcal{T}}\big[(\hat{f}(\mathbf{x}_0) - \bar{f}(\mathbf{x}_0))^2\big].

The derivation (ours, the standard one) adds and subtracts f(x0)f(\mathbf{x}_0) and fˉ(x0)\bar{f}(\mathbf{x}_0): the cross term between ε\varepsilon and the prediction vanishes because the new noise is independent of the training set, and the cross term between the bias and f^−fˉ\hat{f} - \bar{f} vanishes because ET[f^−fˉ]=0\mathbb{E}_{\mathcal{T}}[\hat{f} - \bar{f}] = 0. This is the form the machine learning article quotes from Hastie, Tibshirani and Friedman (§7.3, eq. 7.9). Measured against the noise-free value f(x0)f(\mathbf{x}_0) instead of a noisy target, the σ2\sigma^2 term disappears and the error is squared bias plus variance. Which of the two is meant matters whenever numbers are compared: the example prints both, and they differ by exactly σ2\sigma^2.

Averaging over inputs gives the decomposition of the expected test error: Ex\mathbb{E}_{\mathbf{x}} of each term. Yang, Yu, You, Steinhardt and Ma (2020, §2.1) write it this way, defining Bias2=Ex,y[∥y−fˉ(x)∥22]\mathrm{Bias}^2 = \mathbb{E}_{\mathbf{x},\mathbf{y}}\big[\lVert \mathbf{y} - \bar{f}(\mathbf{x}) \rVert_2^2\big] (their eq. 1) and Variance=ExET[∥f(x,T)−fˉ(x)∥22]\mathrm{Variance} = \mathbb{E}_{\mathbf{x}}\mathbb{E}_{\mathcal{T}}\big[\lVert f(\mathbf{x}, \mathcal{T}) - \bar{f}(\mathbf{x}) \rVert_2^2\big] (eq. 2). Their bias term is measured against the noisy target y\mathbf{y}, so it includes the noise: in the three-term form above it corresponds to σ2\sigma^2 plus the squared bias, while the variance terms correspond directly.

Estimators. For an estimator θ^\hat{\theta} of a fixed quantity θ∗\theta^\ast, computed from a random data set D\mathcal{D}, Murphy (2022, eqs. 4.238–4.242) gives the same algebra without a noise term,

E[(θ^−θ∗)2]=Var⁡[θ^]+bias⁡2(θ^),bias⁡(θ^)=E[θ^]−θ∗,\mathbb{E}\big[(\hat{\theta} - \theta^\ast)^2\big] = \operatorname{Var}[\hat{\theta}] + \operatorname{bias}^2(\hat{\theta}), \qquad \operatorname{bias}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta^\ast,

and draws the moral that a biased estimator can be the wiser choice “so long as it reduces our variance by more than the square of the bias” (§4.7.6.3). The predictor form is the estimator form applied to θ∗=f(x0)\theta^\ast = f(\mathbf{x}_0) at each input, with the noise of the test target added. Goodfellow et al. (§5.4.4) link the relation between bias and variance closely to capacity, underfitting and overfitting: increasing capacity tends to decrease bias and increase variance, which yields the U-shaped generalization error of their Fig. 5.6.

Random and fixed design. What “a new training set” means must be stated. In the random-design setting the inputs xi\mathbf{x}_i of the training set are drawn anew each time; in the fixed-design setting they are held fixed and only the targets are redrawn (Yang et al., §2, who use the random design for their experiments and attribute the fixed design to the theoretical work they compare with). The two settings give different decompositions. The image example below is fixed-design: the pixel positions never change, only the noise.

Linear smoothers. Many estimators are linear in the observations: y^=Sy\hat{\mathbf{y}} = S\mathbf{y} for a matrix SS that does not depend on y\mathbf{y}. Gaussian and box filters, ridge regression evaluated at the training inputs, and kernel regression all have this form. With y=f+ε\mathbf{y} = \mathbf{f} + \boldsymbol{\varepsilon} and Cov⁡(ε)=σ2I\operatorname{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 I, it follows directly (our derivation) that

bias⁡=(S−I) f,Var⁡(y^i)=σ2∑jSij2.\operatorname{bias} = (S - I)\,\mathbf{f}, \qquad \operatorname{Var}(\hat{y}_i) = \sigma^2 \sum_j S_{ij}^2 .

Both are exact, for every output, without simulation. The bias depends on the unknown signal; the variance depends only on the weights. For a filter, ∑jSij2\sum_j S_{ij}^2 is the sum of the squared weights that output pixel actually uses, so it changes near the image border, where the boundary rule folds or truncates the kernel. For a continuous 2-D Gaussian of standard deviation ss pixels the sum of squared weights is 1/(4πs2)1/(4\pi s^2) (our calculation, from ∫Gs2=1/(4πs2)\int G_s^2 = 1/(4\pi s^2)), so the noise variance falls as the square of the width; sampling the kernel on the pixel grid changes this at small ss, which the example prints. Murphy (§4.7.6.5, Fig. 4.25) illustrates the same tradeoff for ridge regression: strongly regularized fits (20 of 100 data sets are plotted) agree with one another but their average over the 100 misses the true curve, while lightly regularized fits scatter widely around an average close to it. The k-nearest-neighbor article gives the variance σ2/k\sigma^2/k of a kk-NN average, another smoother.

Properties

A tendency, not a law

The classical claim is monotone: as capacity grows, bias falls and variance rises, so test error is U-shaped (Goodfellow et al., Fig. 5.6). Yang et al. (§1) trace the claim to Geman, Bienenstock and Doursat (1992), who, as Yang et al. report, gave experimental evidence that the variance of neural networks increases with their width. The decomposition itself is an identity; the monotone curves are an empirical regularity, and measured curves can break them, as the studies below found.

  • Decision trees and k-NN. Domingos (2000, §4.1) measured zero-one bias and variance for C4.5 trees whose depth was limited, on 30 UCI data sets. The minimum of loss at an intermediate depth, the expected tradeoff, appeared in only 10 data sets; in 6 a decision stump was best and in 14 the unlimited tree. In 26 data sets variance grew with depth, but in several the bias also rose after the first levels. For k-NN (§4.3), as kk went from 1 to 21, an increase of bias and a decrease of variance with a minimum in between was rarely observed; in several data sets bias and variance moved in the same direction, and on average bias rose by 4.9% while variance fell by only 0.8%.
  • Wide neural networks. Yang et al. (2020) measured bias and variance of ResNet34 on CIFAR-10 as the width of its convolutional layers grows, training with the squared error between softmax outputs and one-hot labels on half of the 50,000 training images. The bias decreased monotonically, but the variance was unimodal: it rose and then fell (§3.1, Fig. 2). They report the same shape for ResNeXt29, a one-hidden-layer network on MNIST and Fashion-MNIST, and with the cross-entropy loss under a generalized decomposition for Bregman divergences that they attribute to Pfau (2013) (§3.2). Because the variance was small next to the bias, the risk mostly fell monotonically with width; adding label noise raised the variance and produced a double-descent risk curve (§3.3, Fig. 4). The regularization article describes double descent and the interpolation threshold. Murphy (§13.5.7) summarizes the same picture: models with fewer parameters than examples show the usual U-shaped curve that the tradeoff predicts, and test error falls again in the over-parameterized regime.
  • Depth and distribution shift. Yang et al. (§4) report that deeper ResNets had lower bias and higher variance, and that on the corrupted CIFAR-10-C test set the drop in accuracy came mainly from increased bias.

Ensembles reduce variance

Averaging MM predictors with the same bias leaves the bias unchanged. If their errors around the common mean have variance vv and pairwise correlation ρ\rho, the average has variance ρv+(1−ρ)v/M\rho v + (1 - \rho)v/M (our calculation): averaging removes the uncorrelated part and keeps the correlated one. Murphy (§18.2) states the consequence for ensembles: similar bias to the base models but lower variance. Bagging fits each base model to a bootstrap resample of the data, which, Murphy notes (§18.3), does not always improve performance: it relies on the base models being unstable, so that small changes of the data change the fit. Murphy (§18.2) calls decision trees high-variance estimators for this reason, and the random forest article describes how random feature choices lower the correlation between trees. Domingos (§4.2) found that boosting C4.5 lowered zero-one loss mostly through variance: averaged over the data sets, variance reduction contributed 2.5% to the reduction of loss and bias reduction 0.6%.

Estimating bias and variance

Bias needs the true function and variance needs many training sets, and in practice one has a single training set. Three routes are used.

  • A known generator. When the data are synthetic or simulated, as in the example below, new training sets can be drawn at will and the expectations replaced by averages over RR draws. For a linear smoother with known noise, the closed forms above make simulation unnecessary, and the Monte Carlo averages serve as a check.
  • Resampling one data set. Domingos (§4) split each data set into two thirds for training and one third for testing, drew 100 bootstrap replicates of the training part, fitted a model to each and computed the main prediction, bias and variance from the 100 predictions at every test example. Because the noise term is very difficult to estimate, he followed earlier authors in assuming it to be zero. Bootstrap replicates overlap, so they are not independent training sets; the estimates describe the learner’s behaviour on resamples of this data.
  • Disjoint splits. Yang et al. (§2.2) split the training set into NN disjoint parts, train one model on each, and estimate the variance at a test point by the unbiased sample variance of the NN predictions, 1N−1∑j∥f(x,Tj)−1N∑j′f(x,Tj′)∥2\frac{1}{N-1}\sum_j \lVert f(\mathbf{x}, \mathcal{T}_j) - \frac{1}{N}\sum_{j'} f(\mathbf{x}, \mathcal{T}_{j'}) \rVert^2 (an unnumbered display in §2.2). They average such estimates over several random splits, and obtain the bias by subtracting the variance from the risk. Splitting halves each training set, so the estimates describe a model trained on less data: their mainline experiment trained on 25,000 of CIFAR-10’s 50,000 images. The paper uses NN for the number of parts in §2.2 but writes “N=3N = 3 such random splits” in §3.1 for the number of repeated splits, which its Algorithm 1 calls kk.

With RR draws, the sample variance with denominator RR makes the empirical decomposition an exact identity, 1R∑k(uk−c)2=(uˉ−c)2+1R∑k(uk−uˉ)2\frac{1}{R}\sum_k (u_k - c)^2 = (\bar{u} - c)^2 + \frac{1}{R}\sum_k (u_k - \bar{u})^2, while the unbiased estimate of the variance uses R−1R - 1. The probability and statistics article shows the same choice of denominator as a bias–variance tradeoff of its own.

The 0–1 loss

For classification, and for segmentation scored pixel by pixel, the loss is 0 or 1 and the squared-error algebra no longer applies. The machine learning article reports Hastie and colleagues’ simulation (§7.3.1) in which a biased class-probability estimate that stays on the correct side of the decision boundary costs nothing. Murphy (§4.7.6.6) states that under the 0–1 loss bias and variance “combine multiplicatively”, citing an exercise in Hastie et al. and Friedman (1997): on the correct side of the boundary less variance lowers the error, but on the wrong side it pays to increase the variance. He concludes that the tradeoff “is not very useful for classification” and recommends focusing on the expected loss, estimated by cross-validation.

Domingos (2000) instead defines bias and variance for any loss LL, through the main prediction ym=arg⁡min⁡y′ED[L(y,y′)]y_m = \arg\min_{y'} \mathbb{E}_{\mathcal{D}}[L(y, y')], the prediction that differs least from all the predictions made on the training sets in D\mathcal{D}: the mean under squared loss, the median under absolute loss, the mode under 0–1 loss (Definition 1). With y∗y^\ast the optimal prediction and tt the true value:

  • bias B(x)=L(y∗,ym)B(\mathbf{x}) = L(y^\ast, y_m) (Definition 2);
  • variance V(x)=ED[L(ym,y)]V(\mathbf{x}) = \mathbb{E}_{\mathcal{D}}[L(y_m, y)] (Definition 3);
  • noise N(x)=Et[L(t,y∗)]N(\mathbf{x}) = \mathbb{E}_t[L(t, y^\ast)] (Definition 4).

His Theorem 1 shows that for two-class problems and a loss with L(y,y)=0L(y, y) = 0 the expected loss is

ED,t[L(t,y)]=c1N(x)+B(x)+c2V(x),\mathbb{E}_{\mathcal{D},t}[L(t, y)] = c_1 N(\mathbf{x}) + B(\mathbf{x}) + c_2 V(\mathbf{x}),

and for the symmetric 0–1 loss c1=2PD(y=y∗)−1c_1 = 2P_{\mathcal{D}}(y = y^\ast) - 1 and c2=1c_2 = 1 if B(x)=0B(\mathbf{x}) = 0, c2=−1c_2 = -1 if B(x)=1B(\mathbf{x}) = 1. Squared loss is the case c1=c2=1c_1 = c_2 = 1. On an example where the main prediction is already wrong, every training set that happens to give a different prediction gives the right one, so variance lowers the loss: “If a learner is biased on an example, increasing variance decreases loss” (§2). Averaged over examples, the loss is the c1c_1-weighted noise plus average bias plus the net variance Ex[c2V(x)]\mathbb{E}_{\mathbf{x}}[c_2 V(\mathbf{x})], the variance on unbiased examples minus that on biased ones (his eq. 6, and VuV_u, VbV_b in §4). Theorem 2 extends the decomposition to multiclass 0–1 loss, where only part of the variance on biased examples reduces the loss.

Domingos reads two consequences from this. Under 0–1 loss there is a much higher tolerance for variance than an additive decomposition would allow, which in his view helps explain why unstable learners such as decision trees do well in classification (§2). And Theorem 3 writes Schapire and colleagues’ margin of an ensemble as M(x)=±[2B(x)−1][2V(x)−1]M(\mathbf{x}) = \pm[2B(\mathbf{x}) - 1][2V(\mathbf{x}) - 1], so the margin and bias–variance explanations of boosting describe the same quantities (§3). As Domingos (§1) reports, earlier 0–1 decompositions by Kong and Dietterich, Breiman, Tibshirani, James and Hastie, Kohavi and Wolpert, and Friedman each had shortcomings, such as a variance that could be negative (Kong and Dietterich), terms undefined at a single example (Breiman) or a nonzero bias for the Bayes-optimal classifier (Kohavi and Wolpert). Murphy’s multiplicative statement and Domingos’s signed additive one describe the same effect with different definitions; neither is attributed to the other here.

Examples

Gaussian smoothing as an estimator

The example treats one noisy photograph of a known scene as the training set and a Gaussian filter as the learner. A 96×9696 \times 96 clean image CC with values in [0,1][0, 1] has three horizontal bands (rows and columns counted from 0 at the top left):

  • rows 0–31: a ramp in the column xx, C=0.1+0.8 x/95C = 0.1 + 0.8\,x/95;
  • rows 32–63: blocks of 0.2, 0.8 and 0.2 separated by vertical step edges at x=32x = 32 and x=64x = 64, with a 1-pixel-wide line of value 0.8 at x=12x = 12 and a 2-pixel-wide line at x=79x = 79–80;
  • rows 64–95: a flat block of 0.3 for x≤47x \le 47 and, for x≥48x \ge 48, a diagonal grating 0.5+0.3sin⁡(2π(x+r)/6+π/6)0.5 + 0.3\sin(2\pi(x + r)/6 + \pi/6) of period 6 pixels, which takes the values 0.2, 0.35, 0.65 and 0.8.

The observation is Y=C+σZY = C + \sigma Z with independent standard normal ZZ and σ=0.1\sigma = 0.1, without clipping or quantization. The estimate is scipy.ndimage.gaussian_filter(Y, s, mode="reflect", truncate=4.0) for widths s∈{0,0.5,1,1.5,2,3,4}s \in \{0, 0.5, 1, 1.5, 2, 3, 4\} pixels; these are SciPy’s defaults, written out. In SciPy 1.16.3’s source (scipy/ndimage/_filters.py), an axis whose width is at most 10−1510^{-15} is skipped, so s=0s = 0 returns the observation; the kernel radius is int(truncate * sigma + 0.5), while the docstring says round(truncate * sigma), which agree for these widths (radii 0, 2, 4, 6, 8, 12 and 16); the weights are normalized to sum 1; and reflect mirrors the image about its edge. Repeating the observation with new noise is repeating the training set, with the pixel positions fixed: a fixed-design setting.

Because the filter is separable, Y^=AYA⊤\hat{Y} = A Y A^\top with a 96×9696 \times 96 matrix AA whose columns are the filter’s own responses to unit impulses, so the closed forms above apply with the actual weights, borders included: bias ACA⊤−CA C A^\top - C, and variance σ2qrqx\sigma^2 q_r q_x at row rr and column xx, with qq the row sums of the squared entries of AA. The code checks this operator against SciPy’s 2-D filter on the clean image, on a noisy one and on impulses at six pixels including the four corners (C2). It then draws R=501R = 501 noisy observations and RR independent fresh ones, filters each at every width, and compares the Monte Carlo bias, variance and errors with the closed forms. Errors are reported against the clean image (squared bias plus variance) and against a fresh noisy observation (plus σ2\sigma^2). Six regions are scored: ramp (rows 0–15, columns 16–79, at least 16 pixels from other content and from the left and right borders), edges (the two pixels on each side of both step edges, rows 40–55), lines (the line pixels, rows 40–55), flat (rows 72–95, columns 0–39) and texture (rows 72–95, columns 56–95), and the whole image.

The second part thresholds each smoothed image at 0.5, a segmentation, and scores it against the clean labels [C>0.5][C > 0.5] with the 0–1 loss. Because the labels are clean, Domingos’s noise term is zero and y∗=ty^\ast = t. Over the RR observations the code computes the main prediction (the mode; RR is odd, so it is never tied), BB, VV and the expected loss LL at every pixel, and checks Theorem 1, L=B+c2VL = B + c_2 V, as an identity between integer counts (C9). Because each smoothed pixel is exactly Gaussian with mean mm and variance vv, the probability that it lands above the threshold is Φ((m−0.5)/v)\Phi((m - 0.5)/\sqrt{v}), which gives the exact expected loss and the exact bias (the main prediction is the more probable class); the Monte Carlo loss is checked against it (C10).

The plan fixed these checks and hypotheses before any code ran:

  • C1–C4: the data have their planned properties, including a pooled noise standard deviation within 6σ/2N6\sigma/\sqrt{2N} of σ\sigma over the N=R⋅9,216N = R \cdot 9{,}216 draws (C1); the operator equals SciPy’s filter (C2); the bias is zero on the ramp (C3), since a symmetric kernel that sums to one maps a linear function to itself when it sees nothing else; and the sample decomposition with denominator RR holds to rounding (C4); float tolerance 10−1210^{-12}.
  • H1, Monte Carlo against the closed forms with tolerances computed from RR: (a) the bias error divided by its exact standard error v/R\sqrt{v/R} is at most 6 at every pixel; (b) the variance estimate, rescaled by R/(R−1)R/(R-1), deviates from the exact variance by at most 6 of its standard deviations 2/(R−1)\sqrt{2/(R-1)} at every pixel; (c) and (d) the region means of the error against the clean image and against fresh observations are within 4 standard errors, computed from the RR per-draw region means.
  • H2: (a) the variance falls strictly from each width to the next at every pixel, and so does its image mean; (b) the mean squared bias rises strictly in the image, edges, lines and texture regions.
  • H3: the MSE-optimal width (exact) is larger in the flat region than at the edges and than in the texture; the same comparison on the Monte Carlo means is reported too.
  • H4: under 0–1 loss, every pixel of the 1-pixel line is biased for s≥1s \ge 1 and unbiased for s≤0.5s \le 0.5, and every pixel of the 2-pixel line is biased for s≥1.5s \ge 1.5 and unbiased for s≤1s \le 1.
  • C9, C10: Domingos’s identity holds in integers at every pixel; the Monte Carlo expected loss is within 6L(1−L)/R+1/R6\sqrt{L(1-L)/R} + 1/R of the exact one at every pixel.

H3 and H4 are illustrations, not tests. In the plan’s words, the noise level, contrasts and grating period “were chosen at planning, with the hand arithmetic in §4, so that a difference shows”, and H4 was “derived by hand arithmetic at planning (§4), so presented as an illustration”. The code uses NumPy 2.3.5 and SciPy 1.16.3 and runs in about a second.

import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS",
           "VECLIB_MAXIMUM_THREADS"):
    os.environ[_v] = "1"                      # one BLAS thread, before NumPy
import numpy as np
from scipy import ndimage
from scipy.special import ndtr

SEED = 0                                      # first run (pilot: seed 1)
H = W = 96
SIGMA, R, T = 0.1, 501, 0.5                   # noise std, draws, threshold
WIDTHS = [0, 0.5, 1, 1.5, 2, 3, 4]            # blur widths s in pixels

# Clean image: ramp band, step-edge band with two thin lines, flat + texture band.
r, x = np.mgrid[0:H, 0:W]
C = np.empty((H, W))
C[:32] = 0.1 + 0.8 * x[:32] / 95                              # ramp in x
C[32:64] = np.where((x[32:64] >= 32) & (x[32:64] <= 63), 0.8, 0.2)
C[32:64, 12] = 0.8                                            # 1-px line
C[32:64, 79:81] = 0.8                                         # 2-px line
C[64:, :48] = 0.3                                             # flat
C[64:, 48:] = 0.5 + 0.3 * np.sin(2 * np.pi * (x[64:, 48:] + r[64:, 48:]) / 6
                                 + np.pi / 6)                 # texture

def box(r0, r1, cols):
    m = np.zeros((H, W), bool)
    m[r0:r1 + 1, cols] = True
    return m

REG = {"ramp": box(0, 15, slice(16, 80)),
       "edges": box(40, 55, [30, 31, 32, 33, 62, 63, 64, 65]),
       "lines": box(40, 55, [12, 79, 80]),
       "flat": box(72, 95, slice(0, 40)),
       "texture": box(72, 95, slice(56, 96)),
       "image": np.ones((H, W), bool)}
LINE1, LINE2 = box(40, 55, [12]), box(40, 55, [79, 80])

def gf(img, s):
    return ndimage.gaussian_filter(img, sigma=s, mode="reflect", truncate=4.0)

def op(s):                                    # 1-D operator from impulse responses
    if s == 0:
        return np.eye(H)
    return ndimage.gaussian_filter1d(np.eye(H), s, axis=0, mode="reflect",
                                     truncate=4.0)

def hf(ok):
    return "held" if ok else "failed"

rng = np.random.default_rng(SEED)
Z = rng.standard_normal((R, H, W))            # noise of the R observations
Zf = rng.standard_normal((R, H, W))           # noise of R fresh observations
Y = C + SIGMA * Z
Yf = C + SIGMA * Zf
ystar = C > T

print(f"seed {SEED}, R = {R}, sigma = {SIGMA}, image {H}x{W}, threshold {T}")
# ---- C1 data checks
sizes = {k: int(m.sum()) for k, m in REG.items()}
tex = np.unique(np.round(C[64:, 48:], 9))
N = R * H * W
tol_sig = 6 * SIGMA / np.sqrt(2 * N)
sig_y, sig_f = (SIGMA * Z).std(), (SIGMA * Zf).std()
c1 = [sizes == {"ramp": 1024, "edges": 128, "lines": 48, "flat": 960,
                "texture": 960, "image": 9216},
      np.abs(C[REG["flat"]] - 0.3).max() == 0,
      np.abs(np.diff(C[:32], 2, axis=1)).max() <= 1e-12,
      tex.size == 4 and not np.any(tex == 0.5),
      np.abs(C - T).min() > 1e-3,
      abs(sig_y - SIGMA) <= tol_sig, abs(sig_f - SIGMA) <= tol_sig]
print("C1 data: sizes", sizes)
print(f"  flat max|C-0.3| {np.abs(C[REG['flat']] - 0.3).max():.1e}, ramp max|2nd diff| "
      f"{np.abs(np.diff(C[:32], 2, axis=1)).max():.1e}, texture values {tex.tolist()}")
print(f"  min|C-0.5| {np.abs(C - T).min():.6f}, foreground px {int(ystar.sum())}, "
      f"sigma-hat {sig_y:.6f} / fresh {sig_f:.6f} (tol {tol_sig:.3e})  C1 {hf(all(c1))}")

# ---- closed forms, Monte Carlo and checks per width
EX, MC = {}, {}
c2_ok = c3_ok = c4_ok = c5_ok = c6_ok = True
pix = [(0, 0), (0, 95), (95, 0), (95, 95), (5, 50), (48, 48)]
print("\nper width: C2 operator residuals, C3 ramp |bias|, C4 identity, "
      "C5/C6 max |z|, interior variance factor vs 1/(4 pi s^2)")
for s in WIDTHS:
    A = op(s)
    m = A @ C @ A.T                           # exact mean
    b = m - C                                 # exact bias
    q = (A ** 2).sum(1)
    v = SIGMA ** 2 * np.outer(q, q)           # exact variance, borders included
    d_clean = np.abs(m - gf(C, s)).max()
    d_noisy = np.abs(A @ Y[0] @ A.T - gf(Y[0], s)).max()
    d_imp = 0.0
    for (pr, px) in pix:
        delta = np.zeros((H, W)); delta[pr, px] = 1.0
        d_imp = max(d_imp, np.abs(gf(delta, s) - np.outer(A[:, pr], A[:, px])).max())
    Yh = np.stack([gf(Y[k], s) for k in range(R)])
    M = Yh.mean(0)
    bh = M - C
    vh = ((Yh - M) ** 2).mean(0)              # R denominator
    e = (Yh - C) ** 2                         # error against the clean image
    ef = (Yf - Yh) ** 2                       # error against a fresh observation
    d_id = np.abs(e.mean(0) - (bh ** 2 + vh)).max()
    z5 = np.abs((bh - b) / np.sqrt(v / R)).max()
    z6 = np.abs((vh * R / (R - 1) / v - 1) / np.sqrt(2 / (R - 1))).max()
    ramp_b = np.abs(b[REG["ramp"]]).max()
    c2_ok &= max(d_clean, d_noisy, d_imp) <= 1e-12
    c3_ok &= ramp_b <= 1e-12
    c4_ok &= d_id <= 1e-12
    c5_ok &= z5 <= 6
    c6_ok &= z6 <= 6
    info = "" if s == 0 else f"  {q[48] ** 2:.4e} vs {1 / (4 * np.pi * s ** 2):.4e}"
    print(f"s={s:<3} C2 {d_clean:.1e} {d_noisy:.1e} {d_imp:.1e}  C3 {ramp_b:.1e}  "
          f"C4 {d_id:.1e}  C5 {z5:.2f}  C6 {z6:.2f}{info}")
    EX[s] = dict(m=m, b=b, v=v)
    MC[s] = {}
    for k, msk in REG.items():
        pe, pf = e[:, msk].mean(1), ef[:, msk].mean(1)   # per-draw region means
        MC[s][k] = dict(e=pe.mean(), se=pe.std(ddof=1) / np.sqrt(R),
                        f=pf.mean(), sf=pf.std(ddof=1) / np.sqrt(R))
    # 0-1 loss: Domingos's decomposition over the R draws, and the exact values
    y = Yh > T
    n1 = y.sum(0)
    ym = n1 > R / 2                           # main prediction (mode; R odd)
    B = ym != ystar
    err = (y != ystar).sum(0)
    dis = (y != ym[None]).sum(0)
    c2s = np.where(B, -1, 1)
    p1 = ndtr((m - T) / np.sqrt(v))
    L = np.where(ystar, 1 - p1, p1)
    Bx = (m > T) != ystar
    tol10 = 6 * np.sqrt(L * (1 - L) / R) + 1 / R
    MC[s]["01"] = dict(
        L=err.mean() / R, B=B.mean(), net=(c2s * dis).mean() / R,
        Vu=((~B) * dis).mean() / R, Vb=(B * dis).mean() / R,
        nB=int(B.sum()), nBx=int(Bx.sum()), ndiff=int((B != Bx).sum()),
        nlow=int((B & (dis > 0)).sum()),
        c9=int((err != R * B + c2s * dis).sum()),
        c10=(np.abs(err / R - L) / tol10).max(),
        l1=int(Bx[LINE1].sum()), l2=int(Bx[LINE2].sum()))
    del Yh, e, ef, y

# ---- squared-loss tables
print("\nregion   s    bias2(ex)   var(ex)     mse(ex)     mse(MC)     |z|C7  "
      "fresh(ex)   fresh(MC)   |z|C8")
c7_ok = c8_ok = True
for k in REG:
    for s in WIDTHS:
        msk = REG[k]
        b2 = (EX[s]["b"][msk] ** 2).mean()
        vv = EX[s]["v"][msk].mean()
        d = MC[s][k]
        z7 = abs(d["e"] - (b2 + vv)) / d["se"]
        z8 = abs(d["f"] - (SIGMA ** 2 + b2 + vv)) / d["sf"]
        c7_ok &= z7 <= 4
        c8_ok &= z8 <= 4
        EX[s][k] = dict(b2=b2, v=vv, mse=b2 + vv)
        print(f"{k:8s} {s:<4} {b2:.4e}  {vv:.4e}  {b2 + vv:.4e}  {d['e']:.4e}  "
              f"{z7:5.2f}  {SIGMA ** 2 + b2 + vv:.4e}  {d['f']:.4e}  {z8:5.2f}")

# ---- H2, H3
viol = sum(int((EX[t]["v"] >= EX[s]["v"]).sum()) for s, t in zip(WIDTHS, WIDTHS[1:]))
vmean = [EX[s]["v"].mean() for s in WIDTHS]
h2a = viol == 0 and all(a > c for a, c in zip(vmean, vmean[1:]))
h2b = {}
for k in ("image", "edges", "lines", "texture"):
    seq = [EX[s][k]["b2"] for s in WIDTHS]
    h2b[k] = all(c > a for a, c in zip(seq, seq[1:]))
opt_ex = {k: WIDTHS[int(np.argmin([EX[s][k]["mse"] for s in WIDTHS]))] for k in REG}
opt_mc = {k: WIDTHS[int(np.argmin([MC[s][k]["e"] for s in WIDTHS]))] for k in REG}
print(f"\nH2a variance strictly decreasing: pixel-step violations {viol}, "
      f"image means strictly decreasing {all(a > c for a, c in zip(vmean, vmean[1:]))}"
      f"  -> {hf(h2a)}")
print("H2b bias2 strictly increasing: " +
      ", ".join(f"{k} {hf(ok)}" for k, ok in h2b.items()))
print("MSE-optimal s (exact):", opt_ex)
print("MSE-optimal s (MC):   ", opt_mc)
for lab, o in (("exact", opt_ex), ("MC", opt_mc)):
    print(f"H3 ({lab}): flat > edges {hf(o['flat'] > o['edges'])}, "
          f"flat > texture {hf(o['flat'] > o['texture'])}")

# ---- 0-1 loss table and H4
print("\n0-1 loss (MC over R draws; exact from Phi)")
print("s    L        B        net var   Vu       Vb       biased(MC) biased(ex) "
      "differ var-lowers C9fail C10   1px(16) 2px(32)")
h4 = True
c9_ok = c10_ok = True
for s in WIDTHS:
    d = MC[s]["01"]
    c9_ok &= d["c9"] == 0
    c10_ok &= d["c10"] <= 1
    e1 = 16 if s >= 1 else 0
    e2 = 32 if s >= 1.5 else 0
    h4 &= d["l1"] == e1 and d["l2"] == e2
    print(f"{s:<4} {d['L']:.5f}  {d['B']:.5f}  {d['net']:.5f}  {d['Vu']:.5f}  "
          f"{d['Vb']:.5f}  {d['nB']:10d} {d['nBx']:10d} {d['ndiff']:6d} "
          f"{d['nlow']:10d} {d['c9']:6d} {d['c10']:.2f}  {d['l1']:7d} {d['l2']:7d}")
print(f"\nChecks: C1 {hf(all(c1))}, C2 {hf(c2_ok)}, C3 {hf(c3_ok)}, C4 {hf(c4_ok)}, "
      f"C9 {hf(c9_ok)}, C10 {hf(c10_ok)}")
print(f"H1a (C5) {hf(c5_ok)}, H1b (C6) {hf(c6_ok)}, H1c (C7) {hf(c7_ok)}, "
      f"H1d (C8) {hf(c8_ok)}, H4 {hf(h4)}")

Output:

seed 0, R = 501, sigma = 0.1, image 96x96, threshold 0.5
C1 data: sizes {'ramp': 1024, 'edges': 128, 'lines': 48, 'flat': 960, 'texture': 960, 'image': 9216}
  flat max|C-0.3| 0.0e+00, ramp max|2nd diff| 2.2e-16, texture values [0.2, 0.35, 0.65, 0.8]
  min|C-0.5| 0.004211, foreground px 3424, sigma-hat 0.099985 / fresh 0.099914 (tol 1.974e-04)  C1 held

per width: C2 operator residuals, C3 ramp |bias|, C4 identity, C5/C6 max |z|, interior variance factor vs 1/(4 pi s^2)
s=0   C2 0.0e+00 0.0e+00 0.0e+00  C3 0.0e+00  C4 4.3e-17  C5 3.77  C6 3.89
s=0.5 C2 3.3e-16 3.3e-16 1.4e-17  C3 2.2e-16  C4 2.7e-16  C5 3.67  C6 3.98  4.1134e-01 vs 3.1831e-01
s=1   C2 3.3e-16 3.3e-16 1.4e-17  C3 2.2e-16  C4 5.6e-16  C5 3.69  C6 3.96  7.9595e-02 vs 7.9577e-02
s=1.5 C2 4.4e-16 3.3e-16 6.9e-18  C3 2.2e-16  C4 5.8e-16  C5 3.49  C6 3.40  3.5369e-02 vs 3.5368e-02
s=2   C2 3.3e-16 3.3e-16 2.8e-17  C3 3.3e-16  C4 6.9e-16  C5 3.78  C6 3.56  1.9896e-02 vs 1.9894e-02
s=3   C2 4.4e-16 4.4e-16 6.9e-18  C3 3.3e-16  C4 6.4e-16  C5 3.58  C6 3.35  8.8429e-03 vs 8.8419e-03
s=4   C2 5.6e-16 4.4e-16 6.9e-18  C3 5.6e-16  C4 6.1e-16  C5 3.10  C6 3.31  4.9743e-03 vs 4.9736e-03

region   s    bias2(ex)   var(ex)     mse(ex)     mse(MC)     |z|C7  fresh(ex)   fresh(MC)   |z|C8
ramp     0    0.0000e+00  1.0000e-02  1.0000e-02  1.0009e-02   0.45  2.0000e-02  1.9956e-02   1.13
ramp     0.5  4.9232e-33  4.1806e-03  4.1806e-03  4.1848e-03   0.43  1.4181e-02  1.4171e-02   0.34
ramp     1    6.0704e-33  8.4002e-04  8.4002e-04  8.4040e-04   0.09  1.0840e-02  1.0870e-02   1.39
ramp     1.5  7.6992e-33  3.8308e-04  3.8308e-04  3.8341e-04   0.11  1.0383e-02  1.0416e-02   1.59
ramp     2    1.4340e-32  2.2100e-04  2.2100e-04  2.2135e-04   0.16  1.0221e-02  1.0251e-02   1.51
ramp     3    7.9904e-33  1.0312e-04  1.0312e-04  1.0331e-04   0.12  1.0103e-02  1.0130e-02   1.35
ramp     4    4.8927e-32  6.0763e-05  6.0763e-05  6.0850e-05   0.08  1.0061e-02  1.0086e-02   1.27
edges    0    0.0000e+00  1.0000e-02  1.0000e-02  9.9596e-03   0.68  2.0000e-02  1.9879e-02   1.06
edges    0.5  2.0499e-03  4.1134e-03  6.1632e-03  6.1423e-03   0.56  1.6163e-02  1.6060e-02   1.13
edges    1    1.6874e-02  7.9595e-04  1.7670e-02  1.7709e-02   1.14  2.7670e-02  2.7658e-02   0.10
edges    1.5  2.8518e-02  3.5369e-04  2.8872e-02  2.8938e-02   2.07  3.8872e-02  3.8923e-02   0.33
edges    2    3.7888e-02  1.9896e-04  3.8087e-02  3.8142e-02   2.08  4.8087e-02  4.8147e-02   0.35
edges    3    5.0883e-02  8.8429e-05  5.0972e-02  5.0994e-02   1.21  6.0972e-02  6.1016e-02   0.23
edges    4    5.9083e-02  4.9743e-05  5.9133e-02  5.9142e-02   0.74  6.9133e-02  6.9176e-02   0.21
lines    0    0.0000e+00  1.0000e-02  1.0000e-02  9.9346e-03   0.74  2.0000e-02  1.9969e-02   0.16
lines    0.5  8.2129e-03  4.1134e-03  1.2326e-02  1.2368e-02   0.39  2.2326e-02  2.2340e-02   0.06
lines    1    7.4298e-02  7.9595e-04  7.5094e-02  7.5403e-02   1.43  8.5094e-02  8.5264e-02   0.38
lines    1.5  1.2982e-01  3.5369e-04  1.3017e-01  1.3057e-01   1.69  1.4017e-01  1.4040e-01   0.41
lines    2    1.7050e-01  1.9896e-04  1.7070e-01  1.7108e-01   1.67  1.8070e-01  1.8090e-01   0.33
lines    3    2.2200e-01  8.8429e-05  2.2209e-01  2.2240e-01   1.49  2.3209e-01  2.3221e-01   0.18
lines    4    2.5142e-01  4.9744e-05  2.5147e-01  2.5170e-01   1.25  2.6147e-01  2.6151e-01   0.06
flat     0    0.0000e+00  1.0000e-02  1.0000e-02  9.9987e-03   0.06  2.0000e-02  1.9975e-02   0.63
flat     0.5  1.1549e-32  4.1853e-03  4.1853e-03  4.1865e-03   0.12  1.4185e-02  1.4157e-02   1.01
flat     1    0.0000e+00  8.4361e-04  8.4361e-04  8.4757e-04   1.02  1.0844e-02  1.0808e-02   1.62
flat     1.5  1.2840e-34  3.8569e-04  3.8569e-04  3.8804e-04   0.89  1.0386e-02  1.0350e-02   1.71
flat     2    8.1210e-34  2.2312e-04  2.2312e-04  2.2448e-04   0.68  1.0223e-02  1.0188e-02   1.72
flat     3    1.6156e-08  1.0475e-04  1.0477e-04  1.0540e-04   0.45  1.0105e-02  1.0069e-02   1.75
flat     4    1.1077e-06  6.2149e-05  6.3257e-05  6.3375e-05   0.11  1.0063e-02  1.0027e-02   1.79
texture  0    0.0000e+00  1.0000e-02  1.0000e-02  1.0032e-02   1.61  2.0000e-02  1.9981e-02   0.49
texture  0.5  1.8225e-03  4.1853e-03  6.0078e-03  6.0178e-03   0.73  1.6008e-02  1.5997e-02   0.33
texture  1    1.9533e-02  8.4361e-04  2.0376e-02  2.0375e-02   0.07  3.0376e-02  3.0401e-02   0.53
texture  1.5  3.7073e-02  3.8569e-04  3.7459e-02  3.7460e-02   0.11  4.7459e-02  4.7505e-02   0.77
texture  2    4.3574e-02  2.2312e-04  4.3797e-02  4.3798e-02   0.34  5.3797e-02  5.3850e-02   0.83
texture  3    4.4977e-02  1.0475e-04  4.5081e-02  4.5082e-02   0.11  5.5081e-02  5.5136e-02   0.85
texture  4    4.4999e-02  6.2149e-05  4.5061e-02  4.5061e-02   0.15  5.5061e-02  5.5115e-02   0.85
image    0    0.0000e+00  1.0000e-02  1.0000e-02  9.9969e-03   0.46  2.0000e-02  1.9973e-02   2.14
image    0.5  5.8383e-04  4.1583e-03  4.7421e-03  4.7413e-03   0.24  1.4742e-02  1.4726e-02   1.64
image    1    5.4143e-03  8.2560e-04  6.2399e-03  6.2426e-03   1.06  1.6240e-02  1.6241e-02   0.07
image    1.5  9.6009e-03  3.7355e-04  9.9744e-03  9.9780e-03   1.95  1.9974e-02  1.9980e-02   0.46
image    2    1.1697e-02  2.1392e-04  1.1910e-02  1.1913e-02   1.78  2.1910e-02  2.1916e-02   0.46
image    3    1.3586e-02  9.8496e-05  1.3685e-02  1.3686e-02   0.90  2.3685e-02  2.3690e-02   0.40
image    4    1.5010e-02  5.7361e-05  1.5067e-02  1.5068e-02   0.33  2.5067e-02  2.5073e-02   0.40

H2a variance strictly decreasing: pixel-step violations 0, image means strictly decreasing True  -> held
H2b bias2 strictly increasing: image held, edges held, lines held, texture held
MSE-optimal s (exact): {'ramp': 4, 'edges': 0.5, 'lines': 0, 'flat': 4, 'texture': 0.5, 'image': 0.5}
MSE-optimal s (MC):    {'ramp': 4, 'edges': 0.5, 'lines': 0, 'flat': 4, 'texture': 0.5, 'image': 0.5}
H3 (exact): flat > edges held, flat > texture held
H3 (MC): flat > edges held, flat > texture held

0-1 loss (MC over R draws; exact from Phi)
s    L        B        net var   Vu       Vb       biased(MC) biased(ex) differ var-lowers C9fail C10   1px(16) 2px(32)
0    0.04463  0.00195  0.04268  0.04363  0.00095          18          0     18         18      0 0.95        0       0
0.5  0.02533  0.00141  0.02392  0.02457  0.00065          13          4      9         13      0 0.86        0       0
1    0.02167  0.00825  0.01343  0.01429  0.00087          76         75      1         71      0 0.73       16       0
1.5  0.05378  0.02355  0.03023  0.03414  0.00391         217        214      7        148      0 0.69       16      32
2    0.09118  0.04015  0.05104  0.05755  0.00651         370        354     40        203      0 0.69       16      32
3    0.10424  0.09852  0.00572  0.03403  0.02831         908        728    456        659      0 0.86       16      32
4    0.10644  0.10493  0.00151  0.02795  0.02644         967        935    426        640      0 0.68       16      32

Checks: C1 held, C2 held, C3 held, C4 held, C9 held, C10 held
H1a (C5) held, H1b (C6) held, H1c (C7) held, H1d (C8) held, H4 held

Checks. C1, C2, C3, C4, C9 and C10 held in this run. Check C10 failed in the pilot, which ran the same script with seed 1 (see Process), so its outcome depends on the seed. The operator built from impulse responses reproduces SciPy’s 2-D filter to 6×10−166 \times 10^{-16}, corners included, so the closed-form variances use exactly the weights the filter applies. On the ramp the mean squared bias is below 5×10−325 \times 10^{-32} at every width: the smoother is exact on linear intensity, which is why smoothing does no harm to shading that varies linearly. Domingos’s identity holds at every pixel and width with no failures, and the largest C10 ratio is 0.95, at s=0s = 0, below its bar of 1 but not far below it.

Monte Carlo against closed form (H1). All four parts held. The largest bias statistic is 3.78 and the largest variance statistic 3.98, against a bar of 6, and no region mean is more than 2.14 standard errors from its exact value. The fresh-observation column exceeds the clean-image column by σ2=0.01\sigma^2 = 0.01 in every exact row, the noise term of the predictor form; the Monte Carlo columns agree with it within the H1 tolerances. The closed form also shows where the flat region’s bias starts: it is at rounding level up to s=2s = 2 and appears at s=3s = 3 (1.6156×10−81.6156 \times 10^{-8}), the first width whose kernel radius, 12 pixels, reaches past the region’s 9-pixel margin to other content.

Variance and bias against width (H2). Both held. The variance falls at every pixel at every step, and its image mean falls from 10−210^{-2} to 5.7361×10−55.7361 \times 10^{-5}, a factor of about 174. The interior variance factor matches the continuous formula 1/(4πs2)1/(4\pi s^2) to within 0.03% from s=1s = 1 on, but at s=0.5s = 0.5 it is 0.41134 against 0.31831: a kernel of half a pixel, sampled at integer offsets, is far from the continuous Gaussian, and the continuous formula underestimates the remaining noise variance by about 23%. Bias rises at every step in all four tested regions, but at different rates. On the grating it saturates at 4.4999×10−24.4999 \times 10^{-2}, which is the mean square of the grating’s deviations from 0.5 over its six phases, (4×0.152+2×0.32)/6=0.045(4 \times 0.15^2 + 2 \times 0.3^2)/6 = 0.045: a width of 3 or 4 pixels has averaged the texture away completely.

Where smoothing pays (H3). It held, on the exact and on the Monte Carlo means alike: the optimal width is 4 in the flat region and 0.5 at the edges and in the texture. The flat region gains a factor of about 158 in mean squared error (from 1.0000×10−21.0000 \times 10^{-2} to 6.3257×10−56.3257 \times 10^{-5}), the edges and texture less than a factor of two, and for the thin lines no smoothing at all is best: at s=0.5s = 0.5 the squared bias on the lines, 8.2129×10−38.2129 \times 10^{-3}, already exceeds the variance saved, 10−2−4.1134×10−310^{-2} - 4.1134 \times 10^{-3}. The image as a whole is best at s=0.5s = 0.5. One global width is thus a compromise between regions that want opposite settings, which is one motivation for filters that adapt their averaging to the content, such as the bilateral filter and non-local means. The image filtering article measures the same compromise at one width on a step edge, as noise standard deviation against edge width.

0–1 loss (H4 and the counts). H4 held: the 1-pixel line is biased at all 16 pixels from s=1s = 1 on and the 2-pixel line at all 32 pixels from s=1.5s = 1.5 on, and neither is biased at smaller widths. The image-mean expected loss is lowest at s=1s = 1 (0.02167), not at the squared-error optimum of 0.5, and the two decompositions behave differently. The squared bias grows smoothly with the width, while the 0–1 bias jumps where a structure’s smoothed value crosses the threshold: the exact count of biased pixels goes 0, 4, 75, 214, 354, 728 and 935. On biased pixels the variance enters with a minus sign, and the net variance at s=3s = 3 is only 0.00572, the difference between 0.03403 on unbiased pixels and 0.02831 on biased ones. The column “var-lowers” counts biased pixels whose Monte Carlo variance is positive, where by the identity the expected loss is below 1: there, noise that sometimes pushes a pixel back to the correct side of the threshold (a blurred line pixel back above it, for instance) makes the segmentation better on average, Domingos’s point in pixels.

The Monte Carlo and exact bias counts differ on 18, 9, 1, 7, 40, 456 and 426 pixels. In our reading this follows from the construction. A pixel whose smoothed mean lies very close to the threshold has a probability of either class near one half, and the mode of 501 draws can then fall on either side. The construction puts such pixels in at least two places: 64 ramp pixels (columns 47 and 48 of rows 0–31) whose clean value is 0.0042 from the threshold, which matters at small widths, and the grating, whose mean is driven to 0.5 at widths of 3 and 4. Domingos’s quantities computed from a finite set of training sets describe that set; where the learner is nearly undecided, the main prediction itself is uncertain.

Process. The plan, with every check and hypothesis above, was written before any code ran (16:50 on 12 November 2026), and its checksum was recorded before the pilot. Between the two, a reachable-value script ran once (16:50) and printed only check residuals and tail probabilities. From its output, two statements in the plan were corrected before the checksum: the texture has four distinct values, where the draft said six, and the tolerance on the pooled noise standard deviation is 1.974×10−41.974 \times 10^{-4}, where the draft said of order 10−510^{-5}. A planning diagnostic without image data (16:51) computed the tail probabilities of C10’s tolerance. The pilot ran with seed 1, and this first run uses seed 0, as planned; it is the only difference between the two scripts. In the pilot, check C10 failed: at some widths the Monte Carlo 0–1 loss differed from the exact one by more than its tolerance at some pixel. Every other check and hypothesis held in the pilot. The failure was reported to the editor, no print was added to locate it and nothing was changed, and the editor approved running the planned script with seed 0, in which C10 held. No change followed the first run. The example departs from the batch plan’s description in ten ways, all fixed before the pilot: (1) the two thin lines and the lines region were added, because blurring moves steps, flats, ramps and a symmetric grating toward the threshold but not across it, which would leave the 0–1 part with almost nothing to count; (2) R=501R = 501, odd so that the mode is never tied; (3) the exact 0–1 quantities with check C10 were added; (4) hypothesis H4 was added, from hand arithmetic (3 and 4 were proposed in the plan and approved by the editor after the pilot); (5) H2 was made precise (every pixel and the image mean for variance; four regions for bias, leaving out flat and ramp, whose bias is zero or at rounding level at small widths), and the optimal width is decided on the exact means with the Monte Carlo argmin also reported; (6) H1 was split into the four statistics above; (7) the checks C1, C2, C4 and C9 were added; (8) two printed values that decide nothing were added, the continuous variance factor and the count of pixels where the Monte Carlo and exact bias differ; (9) the image layout, regions, σ\sigma, threshold and seeds were fixed; (10) the fixed-design setting is stated.

Common Misconceptions

  • “The tradeoff is a theorem.” The decomposition is an identity for squared error. That bias falls and variance rises with capacity is a tendency: Domingos found bias and variance moving together for k-NN on several data sets, and Yang et al. found unimodal variance in wide networks.
  • “Unbiased estimators are best.” Under squared error a little bias can buy a large reduction in variance; ridge regression and shrinkage estimators (Murphy, §4.7.6.4–4.7.6.5) are built on this, and in the example smoothing lowered the error on the edges, the texture and the whole image only by accepting bias (on the ramp, and in the flat region up to s=2s = 2, it cost no bias at all).
  • “The noise term belongs to the model.” σ2\sigma^2 is in the target, not the predictor, and no learner removes it. Against noise-free ground truth it does not appear; against noisy test targets it does. Comparing an error of one kind with a decomposition of the other is off by σ2\sigma^2, as the two columns of the example show.
  • “Variance always adds to the error.” Under 0–1 loss variance subtracts on biased examples (Domingos, Theorem 1), as the thin lines show.
  • “The generalization gap is the variance.” The gap compares training and test error of one fitted model; the variance compares fits to different training sets. Both grow with capacity in the classical picture, but they are different quantities, and Overfitting and Underfitting works with the gap.
  • “Bias here is dataset bias.” Statistical bias is the error of the average prediction under a fixed data distribution. Dataset bias is a mismatch between the training distribution and the world, which shows up as error that more data from the same source does not remove.

Where It Is Used

Related

  • Machine Learning

    How a program learns a task such as recognizing objects from examples instead of hand-written rules, why doing well on the training images is not the goal, and how features, models, data and evaluation fit together in computer vision.

  • k-Nearest Neighbors

    The classifier and regressor that predict from the k most similar stored examples, with no training step; how the vote, the number of neighbours, the distance and the feature scaling decide its accuracy, the Cover–Hart bound on its error, why it struggles in high dimensions, and where nearest-neighbour methods are used in computer vision.

  • Linear Regression

    The model that predicts a real-valued target as a weighted sum of input features, fitted by ordinary least squares, ridge, lasso or elastic net; what each estimator does to the coefficients, how scikit-learn scales its penalties, and how lasso selects pixels from image patches.

  • Random Forest

    How a random forest combines many decision trees, each grown on a bootstrap sample with a random subset of features tried at every split, how Breiman's strength and correlation explain its error, how out-of-bag votes estimate error without a test set, and how scikit-learn's averaged probabilities and feature importances differ from Breiman's original procedure.

  • Regularization

    Any change to how a network is trained that aims to lower its error on new images rather than on the training images, from weight penalties and data augmentation to early stopping and the implicit bias of SGD.

  • Gaussian Filter

    A linear smoothing filter that replaces each pixel by an average of its neighbors weighted by a Gaussian, whose standard deviation sets the scale of the detail removed.

References

  1. Domingos, P. (2000). A Unified Bias-Variance Decomposition and its Applications. Proceedings of the Seventeenth International Conference on Machine Learning (ICML), 231–238. Morgan Kaufmann.
  2. Yang, Z., Yu, Y., You, C., Steinhardt, J. & Ma, Y. (2020). Rethinking Bias-Variance Trade-off for Generalization of Neural Networks. Proceedings of the 37th International Conference on Machine Learning (ICML), PMLR 119, 10767–10777.
  3. Goodfellow, I., Bengio, Y. & Courville, A. (2016). Deep Learning. MIT Press. §5.2 and §5.4.
  4. Murphy, K. P. (2022). Probabilistic Machine Learning: An Introduction, §4.7.6, §13.5.7 and §18.2–§18.3. MIT Press. Online version of 18 April 2025.
  5. Geman, S., Bienenstock, E. & Doursat, R. (1992). Neural Networks and the Bias/Variance Dilemma. Neural Computation, 4(1), 1–58. (not read)
  6. SciPy developers. scipy/ndimage/_filters.py (gaussian_filter, gaussian_filter1d), source code of SciPy 1.16.3.