Uncertainty Estimation

Uncertainty estimation gives a model’s prediction a statement of how far it may be wrong: a variance or interval around a regressed depth, flow vector or box edge, or a spread of class probabilities for a classified image. A softmax output already looks like such a statement, but its confidence need not match its accuracy, which is what expected calibration error measures, and it says nothing about whether the image resembles anything seen in training. Uncertainty estimation separates two reasons a prediction can be wrong, noise in the data and ignorance in the model, and offers ways to estimate each: likelihoods with a predicted variance, Bayesian neural networks, Monte Carlo dropout, ensembles and single-pass distance-aware models. The term “predictive uncertainty” is used for the combined quantity.

Purpose

A deployed vision system acts on its predictions. Ilg et al. (2018) argue that a driving system cannot be deployed without knowing how reliable the optical flow beneath it is, and their flow network marks the flow in a car’s shadow as unreliable while trusting the static background. Poggi et al. (2020) make the same case for depth from a single camera, where a wrong distance to a pedestrian can be dangerous. Uncertainty lets a system defer to a human or another sensor, filter detections, weigh measurements against each other, choose images to label, and notice that the input has drifted away from the training data.

Two different failures need two different signals. Kendall, Gal and Cipolla (2018, §3.1), citing Kendall and Gal (2017; not read), give the definitions:

• Epistemic uncertainty is uncertainty in the model, about what it does not know for lack of training data; more training data explains it away.
• Aleatoric uncertainty is uncertainty about what the data cannot explain, such as sensor noise or an object boundary that is genuinely ambiguous; it would be explained away only by observing all explanatory variables more precisely, not by more images of the same kind.

They split aleatoric uncertainty further. Heteroscedastic uncertainty depends on the input and is predicted as a model output; homoscedastic uncertainty is constant for all inputs of a task. Their paper uses the homoscedastic kind to weigh the losses of multi-task learning. Mukhoti et al. (2023) add a point that matters for out-of-distribution inputs: aleatoric uncertainty is only meaningful in distribution, because it describes the ambiguity between the classes that might be observed for that input.

When to Use It

Estimate uncertainty when its value changes a decision: a threshold for referral to a person, the weight of a measurement in a filter or a fusion, the choice of which images to annotate next, or a warning that the input is unlike the training data. A calibrated confidence for in-distribution classification alone can be served by temperature scaling, which expected calibration error covers; it changes no prediction and costs one scalar.

The two kinds point to different actions. High aleatoric uncertainty says that this input is hard and more data of the same kind will not help; high epistemic uncertainty says that the model has not seen enough inputs like this one, which more data, or a refusal to predict, addresses. Mukhoti et al. (2023) give active learning as the case where confusing them is costly: an acquisition rule should avoid inputs with high aleatoric but low epistemic uncertainty, and an out-of-distribution detector must not flag ambiguous in-distribution images.

Uncertainty estimates do not make a model robust. Izmailov et al. (2021) found Bayesian neural networks able to detect corrupted inputs but less accurate on them than ordinary training (see Evidence and Impact).

How It Works

Heteroscedastic regression: a predicted variance

For a regression output, the network predicts a mean $\hat{y}(x)$ and a variance $\hat\sigma^2(x)$ and is trained by the negative log-likelihood of a Gaussian, which up to constants is

$$\mathcal{L} = \sum_{i} \frac{\big(y_i - \hat{y}(x_i)\big)^2}{2\hat\sigma^2(x_i)} + \frac{1}{2}\log \hat\sigma^2(x_i),$$

the form Laves et al. (2020, Eq. 4) write (they drop the factors of one half) and Kendall et al. (2018, Eq. 5) write for a single noise parameter. Large residuals are discounted where the network predicts a large variance, and the log term stops it from predicting infinite variance everywhere; in practice the network outputs $\log\hat\sigma^2$ for numerical stability (Laves et al.; Poggi et al.). With a Laplace likelihood the squared error becomes an absolute error divided by a scale, the form Poggi et al. (2020, Eq. 12) use for depth. The derivation from maximum likelihood estimation is the same as for weighted least squares, except that the weights are now learned. Bounding box regression describes He et al.’s Gaussian box edges, an instance in detection.

The predicted variance captures aleatoric, heteroscedastic uncertainty. It is learned from training residuals, and Laves et al. (2020, §2.2) point out the catch: when the mean fits the training labels closely, the residuals understate the noise, so the maximum-likelihood variance is biased low and test-time uncertainty is systematically underestimated, the same bias that Bessel’s correction removes from a sample variance.

Bayesian prediction: a distribution over weights

Epistemic uncertainty is uncertainty about the weights $w$. Bayesian inference replaces the single fitted weight vector by a posterior $p(w \mid \mathcal{D})$ and predicts with the posterior predictive distribution (Goodfellow et al., 2016, §5.6, Eq. 5.68, written there for data without inputs)

$$p(y \mid x, \mathcal{D}) = \int p(y \mid x, w)\, p(w \mid \mathcal{D})\, dw .$$

For a model linear in features, $y = \phi(x)^\top w + \varepsilon$, with a Gaussian prior $w \sim \mathcal{N}(0, s_0^2 I)$ and Gaussian noise of known variance $\sigma^2(x)$, the posterior is Gaussian. The Deep Learning book (§5.6) derives it for unit noise; with per-example variances (our generalization, by the same algebra) it is

$$\Sigma_N = \Big(\Phi^\top W \Phi + s_0^{-2} I\Big)^{-1}, \qquad m_N = \Sigma_N \Phi^\top W \mathbf{y}, \qquad W = \operatorname{diag}\big(1/\sigma^2(x_i)\big),$$

and since a linear function of a Gaussian is Gaussian (see Gaussian distribution), the predictive distribution at a new image is

$$y \mid x, \mathcal{D} \sim \mathcal{N}\Big(\phi(x)^\top m_N,\; \underbrace{\sigma^2(x)}_{\text{aleatoric}} + \underbrace{\phi(x)^\top \Sigma_N\, \phi(x)}_{\text{epistemic}}\Big).$$

The two terms are the two kinds of uncertainty. The aleatoric term belongs to the noise model and does not change with the amount of data. The epistemic term shrinks as data accumulate in the directions $\phi(x)$ points along, and stays at the prior’s variance in directions no training image has explored. The Deep Learning book notes that the posterior mean is the weight decay (ridge) estimate with penalty $s_0^{-2}$; the covariance is what the point estimate throws away.

For a deep network the same predictive decomposition is used with samples. If $M$ weight samples (or ensemble members, or dropout masks) give means $\mu_i(x)$ and variances $\sigma_i^2(x)$, the predictive mean and variance are

$$\mu(x) = \frac{1}{M}\sum_i \mu_i(x), \qquad \sigma^2(x) = \underbrace{\frac{1}{M}\sum_i \big(\mu_i(x) - \mu(x)\big)^2}_{\text{epistemic}} + \underbrace{\frac{1}{M}\sum_i \sigma_i^2(x)}_{\text{aleatoric}},$$

the combination Poggi et al. (2020, Eqs. 16–17) and Laves et al. (2020, Eq. 11) take from Kendall and Gal (2017; not read) and Ilg et al. (2018). Gal and Ghahramani’s Monte Carlo dropout uses the same formula with a constant aleatoric term: the predictive variance is the sample variance of $T$ stochastic forward passes plus the inverse model precision $\tau^{-1}$, which they obtain from the weight decay, the probability that a unit is kept (the parameter of their Bernoulli masks), a prior length-scale and the number of training points (their Eq. 7).

For classification, Mukhoti et al. (2023, Eq. 1) write the analogous split of the predictive entropy:

$$\underbrace{\mathbb{H}[Y \mid x, \mathcal{D}]}_{\text{predictive}} = \underbrace{\mathbb{I}[Y; w \mid x, \mathcal{D}]}_{\text{epistemic}} + \underbrace{\mathbb{E}_{p(w \mid \mathcal{D})}\, \mathbb{H}[Y \mid x, w]}_{\text{aleatoric (in distribution)}} .$$

The mutual information is large when sampled models disagree about the class; the expected entropy is large when each model, individually, is unsure.

Variants

Bayesian neural networks

The posterior over a deep network’s weights has no closed form, and the methods differ in how they approximate it.

• Bayes by Backprop (Blundell et al., 2015) fits a diagonal Gaussian $q(w \mid \theta)$ by minimizing the variational free energy $\mathrm{KL}[q(w \mid \theta) \,\|\, P(w)] - \mathbb{E}_q[\log P(\mathcal{D} \mid w)]$ (their Eq. 1), a “complexity cost” plus a “likelihood cost”. Weights are sampled as $w = \mu + \log(1 + \exp\rho) \circ \epsilon$ with $\epsilon$ standard normal, so the gradient passes through the sample, and the prior is a scale mixture of two Gaussians. On MNIST with two hidden layers of 400 to 1,200 units, the scale-mixture version reached 1.32–1.36% test error against 1.33–1.51% for dropout (their Table 1), at about twice dropout’s cost per iteration. On a toy regression its confidence intervals widened away from the data, where a plain network fitted one extrapolation confidently.
• Laplace approximations fit a Gaussian at a trained network’s weights. Kristiadi, Hein and Hennig (2020) show for binary ReLU classifiers that a Gaussian over the weights bounds the confidence far from the training data away from one, and that a Gaussian over the last layer alone suffices. This supports the common practice of treating all but the last layer as a fixed feature map and putting a Bayesian linear model on top, the setting of the example below; with the probit approximation it leaves the decision boundary unchanged.
• Monte Carlo dropout (Gal and Ghahramani, 2016) casts a network trained with dropout as approximate variational inference in a deep Gaussian process and keeps dropout on at test time. Its mechanics and an example on digits belong to Dropout.
• Markov chain Monte Carlo. Izmailov et al. (2021) ran full-batch Hamiltonian Monte Carlo on modern architectures, parallelized over 512 TPUv3 devices for CIFAR-10 and IMDB, to see what the true posterior predicts; their findings are under Evidence and Impact.

Ensembles

An ensemble trains several networks and treats their spread as epistemic uncertainty. Sources describe deep ensembles (Lakshminarayanan et al., 2017; not read) differently. Van Amersfoort et al. (2020) describe them as networks trained from different random initializations with a different ordering of the data, while Ilg et al. (2018) and Poggi et al. (2020) build “bootstrapped ensembles” whose members train on different subsets of the data (67% of it per member in Ilg et al., a random 25% in Poggi et al.). Snapshot ensembles take members from one training run with a cyclic learning rate. Predictive ensembles combine members that each predict a variance, by the decomposition above.

An ensemble costs $M$ trainings and $M$ forward passes. Ilg et al. found that the best flow uncertainties came from bootstrapped ensembles of predictive networks and that a single network predicting several hypotheses, merged by a second network, matched them at interactive frame rates. Van Amersfoort et al. show a deep ensemble on the two-moons toy problem that is uncertain only along the decision boundary and confident far from the data, which they attribute to a lack of diversity among the members; the bootstrap rows of the example below show the linear version of the same failure.

Single-pass deterministic methods

Distance-aware models estimate epistemic uncertainty from where an input lands in feature space, in one forward pass. DUQ (van Amersfoort et al., 2020) replaces the softmax with radial basis functions around one centroid per class and measures uncertainty as the distance to the closest centroid; a gradient penalty keeps the features sensitive to the input, against what the authors call feature collapse, out-of-distribution inputs mapped onto in-distribution features. DDU (Mukhoti et al., 2023) keeps an ordinary softmax network made smooth and sensitive by residual connections and spectral normalization, then fits a Gaussian discriminant analysis to its features after training: the feature density is the epistemic signal and the softmax entropy the aleatoric one.

Recalibration

Uncertainty estimates can be calibrated after training on held-out data. For classification, temperature scaling (Guo et al., 2017) is described in expected calibration error. For regression, Kuleshov, Fenner and Ermon (2018) call a forecaster calibrated when the fraction of targets below the predicted $p$-quantile tends to $p$ for every $p$ (their Eq. 3), and recalibrate any model by fitting a monotone map, such as isotonic regression, from predicted to observed cumulative probability. Laves et al. (2020) instead scale the predicted standard deviation by one scalar $s$ fitted by likelihood on a validation set (“σ scaling”, their Eqs. 7–10), applied to the total of MC-dropout epistemic and predicted aleatoric variance. On their medical regression tasks a learned auxiliary recalibration network could overfit a small calibration set, and the single scalar gave the better test likelihood in those cases (their Table 1).

How to Apply It

Bayesian linear regression with an exact posterior

The example uses a model whose posterior is exact, so that every uncertainty it reports can be checked: Bayesian linear regression on fixed features of synthetic images. This is the last-layer setting that Kristiadi et al. (2020) motivate, with fixed hand-made features in place of a trained network.

Each $24 \times 24$ image holds one disk, rendered with $4 \times 4$ supersampling so that each pixel stores the fraction of its area inside the disk. The task is to estimate the disk’s area in pixels from 37 features: a bias and the sums of the 36 blocks of $4 \times 4$ pixels. The ground-truth area is the rendered area, which equals the sum of the block sums, so a linear model with all block weights 1 is exactly right. Labels carry noise with standard deviation half the radius, $\sigma(x) = 0.5\sqrt{\text{area}/\pi}$, as if annotators traced the boundary and erred in proportion to its length: aleatoric and heteroscedastic. Training images have their centre in $[10, 14]^2$ and radius 2 to 4, so they light only the 16 inner blocks; the 20 border blocks are zero in every training image. Two shifted test sets probe the model off its training distribution: shifted disks of the same sizes in a random corner, which light the border blocks, and large disks of radius 5 to 6 at the centre, which light only inner blocks but with areas of 78 to 113 against at most 50 in training. The prior is $w \sim \mathcal{N}(0, s_0^2 I)$ with $s_0 = 1$.

Every result is averaged over 20 independent training sets, each used at sizes $N$ = 10, 40, 160, 640 and 2,560 (nested), and measured on 2,000 test images per set. The script:

• checks the data and its components: border blocks unlit in training, labels with standardized noise of mean 0 and standard deviation 1, exact recovery of the area from noise-free labels, a fitted noise model against the true one, and an identity used by the perturbed ensemble;
• A. with the known noise model, reports mean aleatoric and epistemic variance on in-distribution test images as $N$ grows, and the coverage of nominal 95% intervals built from the total, the aleatoric part alone and the epistemic part alone;
• B. compares the known noise model with a fitted heteroscedastic one ($\log\sigma^2$ linear in the log of the image’s total intensity, fitted by maximum likelihood, which is convex in its two parameters, alternating with the posterior mean) and a fitted homoscedastic one, by the ratio of mean fitted to true variance and by coverage overall and in the smallest and largest third of the objects;
• C. at $N = 640$, reports both terms, the RMSE of the posterior mean and the coverage on each test set, and the ROC AUC for telling in-distribution images from each shifted set by aleatoric, epistemic and total variance. The prior scale fixes the epistemic variance in unexplored directions and no source fixes it, so this part runs $s_0$ = 0.3, 1 and 3;
• D. compares two 100-member ensembles with the exact posterior at $N$ = 40 and 640. Bootstrap members refit to images resampled with replacement. Ordinary least squares is undefined here, because the border blocks are zero in every training image, so each member is the same weighted ridge (MAP) fit that the posterior mean is, with the prior as its penalty. Perturbed members fit the same MAP objective to labels with fresh noise of variance $\sigma^2(x_i)$ added, pulled towards a weight vector drawn from the prior instead of towards zero. By our derivation each perturbed member is an exact posterior sample: its weights are $\Sigma_N\big(\Phi^\top W(\mathbf{y} + \boldsymbol\varepsilon) + w_0/s_0^2\big)$, whose mean is $m_N$ and whose covariance is $\Sigma_N(\Phi^\top W \Phi + s_0^{-2} I)\Sigma_N = \Sigma_N$, the identity the script checks. The table gives each ensemble’s variance as a ratio to the posterior’s epistemic variance, and the ROC AUC for detecting shifted images.

The script takes from a few seconds to about half a minute on a shared machine, with NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2. Passing three integers on the command line replaces the seeds; the planning runs used other seeds to set the expected values before the run shown.

import sys
import numpy as np
from scipy.stats import norm
from sklearn.metrics import roc_auc_score

# Seeds: test sets, training replicates, ensembles (the planning run passes other seeds).
SEED_TEST, SEED_TRAIN, SEED_ENS = (int(a) for a in (sys.argv[1:4] if len(sys.argv) > 3 else (0, 100, 200)))
S, SUB, G = 24, 4, 6                 # 24x24 images, 4x4 supersampling, 6x6 grid of 4x4 blocks
coords = (np.arange(S * SUB) + 0.5) / SUB
ring = np.zeros((G, G), bool); ring[[0, -1], :] = True; ring[:, [0, -1]] = True
ring = ring.ravel()                   # the 20 border blocks

def render(cy, cx, r):
    """Anti-aliased disks: each pixel holds the fraction of its area inside the disk."""
    dy = coords[None, :] - cy[:, None]; dx = coords[None, :] - cx[:, None]
    inside = dy[:, :, None] ** 2 + dx[:, None, :] ** 2 <= r[:, None, None] ** 2
    return inside.reshape(len(r), S, SUB, S, SUB).mean(axis=(2, 4))

def make(n, rng, kind):
    """kind: 'id' (centred, r 2-4), 'shifted' (in a corner, r 2-4), 'large' (centred, r 5-6)."""
    if kind == "large":
        r = rng.uniform(5, 6, n); cy, cx = rng.uniform(11, 13, (2, n))
    else:
        r = rng.uniform(2, 4, n)
        if kind == "id":
            cy, cx = rng.uniform(10, 14, (2, n))
        else:                          # a random corner, 0.2-1.5 px from both near edges
            u = rng.uniform(0.2, 1.5, (2, n)); far = rng.random((2, n)) < 0.5
            cy, cx = np.where(far, S - r - u, r + u)
    img = render(cy, cx, r)
    blocks = img.reshape(n, G, 4, G, 4).sum(axis=(2, 4)).reshape(n, G * G)
    area = img.sum(axis=(1, 2))                       # ground-truth area (px^2)
    sigma = 0.5 * np.sqrt(area / np.pi)               # label noise sd: half the radius
    y = area + sigma * rng.standard_normal(n)
    Phi = np.hstack([np.ones((n, 1)), blocks])        # 37 features: bias + 36 block sums
    return Phi, y, area, sigma

def posterior(Phi, y, var, s0):
    """Exact Gaussian posterior for y ~ N(Phi w, var), prior w ~ N(0, s0^2 I)."""
    A = Phi.T @ (Phi / var[:, None]) + np.eye(Phi.shape[1]) / s0 ** 2
    Sig = np.linalg.inv(A); Sig = (Sig + Sig.T) / 2
    return Sig @ (Phi.T @ (y / var)), Sig

def map_fit(Phi, y, var, s0, w0=None):
    """Weighted ridge (MAP) fit, optionally pulled towards w0 instead of 0."""
    A = Phi.T @ (Phi / var[:, None]) + np.eye(Phi.shape[1]) / s0 ** 2
    b = Phi.T @ (y / var) + (0 if w0 is None else w0 / s0 ** 2)
    return np.linalg.solve(A, b)

def fit_noise(Phi, y, psi, s0, iters=5):
    """Fit log sigma^2 = psi @ v by maximum likelihood (convex in v), alternating with the mean."""
    m = map_fit(Phi, y, np.ones(len(y)), s0)            # start: unweighted fit
    v = np.zeros(psi.shape[1]); v[0] = np.log(np.mean((y - Phi @ m) ** 2))
    for _ in range(iters):
        m, _ = posterior(Phi, y, np.exp(psi @ v), s0)
        r2 = (y - Phi @ m) ** 2
        nll = lambda v: 0.5 * np.sum(psi @ v + r2 * np.exp(-psi @ v))
        for _ in range(50):                            # damped Newton (backtracking line search)
            e = r2 * np.exp(-psi @ v)
            g = 0.5 * psi.T @ (1 - e); H = 0.5 * psi.T @ (psi * e[:, None])
            step, t = np.linalg.solve(H + 1e-9 * np.eye(len(v)), g), 1.0
            while nll(v - t * step) > nll(v) - 1e-4 * t * (g @ step) and t > 1e-10:
                t /= 2
            v = v - t * step
    return v

def cover(err, var, z=norm.ppf(0.975)):
    return np.mean(np.abs(err) <= z * np.sqrt(var))

psi_of = lambda Phi: np.column_stack([np.ones(len(Phi)), np.log(Phi[:, 1:].sum(1))])

# ---------------- data and checks ----------------
rng = np.random.default_rng(SEED_TEST)
test = {k: make(2000, rng, k) for k in ("id", "shifted", "large")}
R, NS, S0 = 20, (10, 40, 160, 640, 2560), 1.0
train = [make(max(NS), np.random.default_rng(SEED_TRAIN + i), "id") for i in range(R)]

print("data checks")
for k, (Phi, y, area, sig) in test.items():
    ink = Phi[:, 1:][:, ring].sum(1)
    print(f"  {k:8s} area {area.min():5.1f}-{area.max():5.1f}  sigma {sig.min():.2f}-{sig.max():.2f}  "
          f"border-block ink min {ink.min():.2f} max {ink.max():.2f}")
Phi, y, area, sig = np.vstack([t[0] for t in train]), *(np.concatenate([t[j] for t in train]) for j in (1, 2, 3))
print(f"  training pool: {len(y)} images, border-block ink max {Phi[:, 1:][:, ring].sum(1).max():.2f}, "
      f"area - sum of blocks max |diff| {np.abs(area - Phi[:, 1:].sum(1)).max():.1e}")
z = (y - area) / sig
print(f"  standardized label noise: mean {z.mean():+.3f}, sd {z.std():.3f}; "
      f"noise floor RMS {np.sqrt(np.mean(sig ** 2)):.2f} vs sd of y {y.std():.2f}")
Pt, yt, at, st = test["id"]
m, Sig = posterior(train[0][0], train[0][1], train[0][3] ** 2, S0)
print(f"  exact posterior, N = {max(NS)}: test RMSE {np.sqrt(np.mean((yt - Pt @ m) ** 2)):.3f} "
      f"vs noise floor {np.sqrt(np.mean(st ** 2)):.3f}")
Pc, _, ac, sc = train[0]                               # noise-free labels: posterior recovers area
mc, _ = posterior(Pc, ac, np.full(len(ac), 1e-6), 1e3)
print(f"  noise-free labels: max |prediction - area| on test {np.abs(Pt @ mc - at).max():.1e}")
vt = fit_noise(Pc, ac + sc * np.random.default_rng(SEED_ENS).standard_normal(len(ac)), psi_of(Pc), S0)
q = np.exp(psi_of(Pt) @ vt) / st ** 2
print(f"  noise-model fit, N = {max(NS)}: fitted / true variance on test images {q.min():.3f}-{q.max():.3f}")
A = Pc.T @ (Pc / (sc ** 2)[:, None]) + np.eye(37) / S0 ** 2
_, Sg = posterior(Pc, ac, sc ** 2, S0)
print(f"  perturbed-fit covariance identity: max |Sig A Sig - Sig| {np.abs(Sg @ A @ Sg - Sg).max():.1e}")

# ---------------- A: more data, epistemic falls; aleatoric stays ----------------
print("\nA. in-distribution test images, known noise model, mean over 20 training sets")
print("      N  aleatoric  epistemic  ratio   coverage of 95% intervals: total  aleatoric-only  epistemic-only")
prev = None
for n in NS:
    al, ep, ct, ca, ce = [], [], [], [], []
    for P, yy, _, s in train:
        m, Sig = posterior(P[:n], yy[:n], s[:n] ** 2, S0)
        e = yt - Pt @ m; va = st ** 2; ve = np.einsum("ij,jk,ik->i", Pt, Sig, Pt)
        al.append(va.mean()); ep.append(ve.mean())
        ct.append(cover(e, va + ve)); ca.append(cover(e, va)); ce.append(cover(e, ve))
    ep_m = np.mean(ep)
    print(f"  {n:5d}  {np.mean(al):9.3f}  {ep_m:9.3f}  {'' if prev is None else f'{prev / ep_m:5.2f}':>5s}"
          f"   {np.mean(ct):38.3f}  {np.mean(ca):14.3f}  {np.mean(ce):14.3f}")
    prev = ep_m

# ---------------- B: the noise model, known or fitted ----------------
small, big = at <= np.quantile(at, 1 / 3), at >= np.quantile(at, 2 / 3)
print("\nB. noise model (in-distribution test images, mean over 20 training sets)")
print("      N  noise model        mean var / true  coverage: all  smallest third  largest third")
for n in (40, 640):
    for name in ("known", "heteroscedastic", "homoscedastic"):
        rat, c_all, c_s, c_b = [], [], [], []
        for P, yy, _, s in train:
            P, yy, s = P[:n], yy[:n], s[:n]
            if name == "known":
                var_tr, var_te = s ** 2, st ** 2
            else:
                psi = psi_of(P) if name == "heteroscedastic" else np.ones((n, 1))
                v = fit_noise(P, yy, psi, S0)
                var_tr = np.exp(psi @ v)
                var_te = np.exp((psi_of(Pt) if name == "heteroscedastic" else np.ones((len(Pt), 1))) @ v)
            m, Sig = posterior(P, yy, var_tr, S0)
            tot = var_te + np.einsum("ij,jk,ik->i", Pt, Sig, Pt); e = yt - Pt @ m
            rat.append(var_te.mean() / (st ** 2).mean())
            c_all.append(cover(e, tot)); c_s.append(cover(e[small], tot[small])); c_b.append(cover(e[big], tot[big]))
        print(f"  {n:5d}  {name:17s}  {np.mean(rat):15.3f}  {np.mean(c_all):13.3f}  {np.mean(c_s):14.3f}  {np.mean(c_b):13.3f}")

# ---------------- C: shifted images ----------------
print("\nC. N = 640, known noise model, mean over 20 training sets; AUROC: in-distribution vs that set")
print("  prior sd  test set  aleatoric  epistemic     RMSE  coverage  AUROC: aleatoric  epistemic  total")
for s0 in (0.3, 1.0, 3.0):
    rows = {k: [] for k in test}
    for P, yy, _, s in train:
        m, Sig = posterior(P[:640], yy[:640], s[:640] ** 2, s0)
        out = {}
        for k, (Pk, yk, ak, sk) in test.items():
            ve = np.einsum("ij,jk,ik->i", Pk, Sig, Pk); e = yk - Pk @ m
            out[k] = (sk ** 2, ve, e)
        lab = np.r_[np.zeros(2000), np.ones(2000)]
        va0, ve0, _ = out["id"]
        for k in test:
            va, ve, e = out[k]
            au = [np.nan] * 3 if k == "id" else [roc_auc_score(lab, np.r_[va0, va]),
                                                  roc_auc_score(lab, np.r_[ve0, ve]),
                                                  roc_auc_score(lab, np.r_[va0 + ve0, va + ve])]
            rows[k].append([va.mean(), ve.mean(), np.sqrt(np.mean(e ** 2)), cover(e, va + ve), *au])
    for k in test:
        a = np.mean(rows[k], axis=0)
        auc = "" if k == "id" else f"{a[4]:16.3f}  {a[5]:9.3f}  {a[6]:5.3f}"
        print(f"  {s0:8.1f}  {k:8s}  {a[0]:9.2f}  {a[1]:9.2f}  {a[2]:7.2f}  {a[3]:8.3f}  {auc}")

# ---------------- D: ensembles against the exact posterior ----------------
print("\nD. 100-member ensembles, known noise model, mean over 20 training sets")
print("      N  ensemble     variance / posterior epistemic: id  shifted  large   AUROC id vs shifted")
for n in (40, 640):
    res = {"bootstrap": [], "perturbed": [], "posterior": []}
    for i, (P, yy, _, s) in enumerate(train):
        P, yy, s = P[:n], yy[:n], s[:n]
        m, Sig = posterior(P, yy, s ** 2, S0)
        g = np.random.default_rng(SEED_ENS + 1000 * i + n)
        W = {"bootstrap": [], "perturbed": []}
        for b in range(100):
            idx = g.integers(0, n, n)                                     # resample images
            W["bootstrap"].append(map_fit(P[idx], yy[idx], s[idx] ** 2, S0))
            W["perturbed"].append(map_fit(P, yy + s * g.standard_normal(n), s ** 2, S0,
                                          w0=S0 * g.standard_normal(P.shape[1])))
        post = {k: np.einsum("ij,jk,ik->i", test[k][0], Sig, test[k][0]) for k in test}
        for name in W:
            Wb = np.array(W[name])
            var = {k: (test[k][0] @ Wb.T).var(axis=1, ddof=1) for k in test}
            lab = np.r_[np.zeros(2000), np.ones(2000)]
            res[name].append([var[k].mean() / post[k].mean() for k in test] +
                             [roc_auc_score(lab, np.r_[var["id"], var["shifted"]])])
        res["posterior"].append([1.0, 1.0, 1.0, roc_auc_score(lab, np.r_[post["id"], post["shifted"]])])
    for name, r in res.items():
        a = np.mean(r, axis=0)
        print(f"  {n:5d}  {name:10s}  {a[0]:34.3f}  {a[1]:7.3f}  {a[2]:5.3f}  {a[3]:20.3f}")

Output:

data checks
  id       area  12.4- 50.2  sigma 0.99-2.00  border-block ink min 0.00 max 0.00
  shifted  area  12.4- 50.6  sigma 0.99-2.01  border-block ink min 11.69 max 35.19
  large    area  78.4-113.2  sigma 2.50-3.00  border-block ink min 0.00 max 0.00
  training pool: 51200 images, border-block ink max 0.00, area - sum of blocks max |diff| 0.0e+00
  standardized label noise: mean -0.005, sd 1.002; noise floor RMS 1.53 vs sd of y 11.05
  exact posterior, N = 2560: test RMSE 1.576 vs noise floor 1.527
  noise-free labels: max |prediction - area| on test 3.3e-12
  noise-model fit, N = 2560: fitted / true variance on test images 0.906-1.069
  perturbed-fit covariance identity: max |Sig A Sig - Sig| 7.3e-16

A. in-distribution test images, known noise model, mean over 20 training sets
      N  aleatoric  epistemic  ratio   coverage of 95% intervals: total  aleatoric-only  epistemic-only
     10      2.331      5.390                                           0.953           0.757           0.827
     40      2.331      1.187   4.54                                    0.951           0.908           0.618
    160      2.331      0.236   5.04                                    0.945           0.936           0.366
    640      2.331      0.055   4.27                                    0.944           0.941           0.192
   2560      2.331      0.014   3.90                                    0.944           0.944           0.096

B. noise model (in-distribution test images, mean over 20 training sets)
      N  noise model        mean var / true  coverage: all  smallest third  largest third
     40  known                        1.000          0.951           0.955          0.946
     40  heteroscedastic              0.635          0.867           0.883          0.844
     40  homoscedastic                0.612          0.881           0.950          0.811
    640  known                        1.000          0.944           0.953          0.934
    640  heteroscedastic              0.961          0.938           0.949          0.927
    640  homoscedastic                0.959          0.937           0.992          0.877

C. N = 640, known noise model, mean over 20 training sets; AUROC: in-distribution vs that set
  prior sd  test set  aleatoric  epistemic     RMSE  coverage  AUROC: aleatoric  epistemic  total
       0.3  id             2.33       0.05     1.59     0.944  
       0.3  shifted        2.35      23.35    30.85     0.000             0.506      1.000  1.000
       0.3  large          7.57       0.81     3.33     0.915             1.000      0.997  1.000
       1.0  id             2.33       0.06     1.58     0.944  
       1.0  shifted        2.35     253.99    30.21     0.719             0.506      1.000  1.000
       1.0  large          7.57       1.10     2.91     0.953             1.000      0.996  1.000
       3.0  id             2.33       0.06     1.59     0.944  
       3.0  shifted        2.35    2073.62    28.72     1.000             0.506      1.000  1.000
       3.0  large          7.57       2.49     2.94     0.960             1.000      0.995  1.000

D. 100-member ensembles, known noise model, mean over 20 training sets
      N  ensemble     variance / posterior epistemic: id  shifted  large   AUROC id vs shifted
     40  bootstrap                                0.794    0.001  0.879                 0.331
     40  perturbed                                1.014    1.019  1.036                 1.000
     40  posterior                                1.000    1.000  1.000                 1.000
    640  bootstrap                                0.935    0.014  0.738                 0.927
    640  perturbed                                0.995    1.012  0.991                 1.000
    640  posterior                                1.000    1.000  1.000                 1.000

The checks. Training images never light a border block, shifted images always do, and the labels have the noise the generator intends. The posterior at $N = 2{,}560$ predicts in-distribution areas with an RMSE of 1.576 against a noise floor of 1.527, and with noise-free labels it recovers the area to rounding error, so the features determine the target and the remaining error is the label noise. The fitted noise model, with plenty of data, comes within about 10% of the true variance on every test image.

More data, less epistemic uncertainty (A). The aleatoric column does not move with $N$; it is the noise in the labels. The epistemic column falls by a factor of about four each time $N$ grows fourfold, the $1/N$ rate expected once the data outweigh the prior, and by $N = 160$ it is a tenth of the aleatoric term. Intervals from the total variance cover close to 95% of the test labels at every size. Leaving out the epistemic part costs most when data are few, 76% coverage at $N = 10$, and nothing at $N = 2{,}560$; intervals from the epistemic part alone, which is all that a spread of models measures without a noise model, cover fewer and fewer labels as the model improves. In this example neither term is a substitute for the other.

The noise model (B). At $N = 640$ the fitted heteroscedastic model’s variance is within 4% of the true one on average and its intervals cover 0.938 of the labels overall and 0.949 in the smallest third. The planned check that its coverage in the largest third of objects stays within 0.95 ± 0.02 failed: it printed 0.927. The known noise model reached 0.934 on the same images. The homoscedastic model has about the right average variance but spends it in the wrong places: intervals too wide for small objects (0.992) and too narrow for large ones (0.877). At $N = 40$, with 37 weights, both fitted models underestimate the noise by more than a third and cover under 89% of the labels; in our reading, the posterior mean fits the 40 training labels closely, so their residuals understate the noise, the bias Laves et al. (2020) describe. The known model at $N = 40$ is still calibrated, because its epistemic term honestly reports the small sample.

Shifted images (C). On the corner disks the aleatoric variance is the same as in distribution, and its ROC AUC of 0.506 means it cannot tell them apart: the noise model describes how labels scatter, not whether the image is familiar. The epistemic variance separates the shifted set perfectly under every prior. The posterior mean is wrong by about 30 on these images, because it puts the prior’s zero weight on the border blocks that carry part of the disk. Whether the intervals cover that error depends entirely on the prior scale: none of the labels at $s_0 = 0.3$, 72% at $s_0 = 1$ and all of them at $s_0 = 3$. In the directions no training image explored, the epistemic variance is whatever the prior says, a modelling choice and not something the data can estimate; in distribution, coverage is 0.944 whatever the prior. The large disks show a different kind of shift. They light only blocks the model has learned, the fitted weights extrapolate correctly, and epistemic variance grows by more than an order of magnitude yet stays small beside the aleatoric term, so coverage stays near 95% at $s_0$ = 1 and 3 and falls to 0.915 under the tightest prior. Their aleatoric variance does separate them perfectly, but only because the noise model gives bigger objects more noise, and no large disk is as small as the biggest training disk.

Ensembles against the posterior (D). The perturbed ensemble reproduces the posterior’s epistemic variance on every test set to within a few percent, as the derivation says it must; in our reading the remaining differences are the Monte Carlo error of 100 members. The bootstrap ensemble comes within about 7% in distribution at $N = 640$, but on shifted images it reports 1% of the posterior variance or less. Every resample has zero in the border blocks, so every member fits the same zero weight there, and the members disagree only about the inner blocks. Ranked by this variance, shifted images are still mostly separable at $N = 640$ (ROC AUC 0.927), because their inner blocks differ too; at $N = 40$ they look more certain than in-distribution images (0.331). In this linear model, resampling the data cannot produce disagreement in directions the data never cover; drawing the prior into each member can.

All settings, sizes, seeds and printed values were fixed before the first run and not changed after it. One component, the Newton fit of the noise model, diverged in a planning run on other seeds and was given a line search and a new starting point before the pilot, and the run shown is the first run.

Evaluation

Calibration. A regressor’s intervals should contain the target as often as they claim; Kuleshov et al. (2018) define this through predicted quantiles, and coverage of a 95% interval, as in the example, is one point of that curve. Calibration alone is not enough: Kuleshov et al. note that a forecast ignoring the input, the marginal distribution of the target, is calibrated but useless, so intervals should also be sharp, as narrow as possible. The example’s homoscedastic rows show a further limit: intervals can be roughly calibrated on average while badly calibrated for small and large objects. Laves et al. (2020) bin predictions by predicted variance and compare each bin’s mean squared error with its mean variance, an uncertainty calibration error built like expected calibration error. For classifiers, expected calibration error, the Brier score and the negative log-likelihood apply; that article covers them.

Sparsification. For dense predictions such as flow and depth, Ilg et al. (2018) sort pixels by predicted uncertainty, remove them in order of decreasing uncertainty and plot the error of the rest. The oracle removes pixels in order of true error; the difference between the two curves is the sparsification error, and its area is the AUSE (lower is better). For their best flow network, removing the 20% most uncertain pixels halved the average endpoint error on Sintel train clean. Poggi et al. (2020) add the area under the random gain (AURG), the improvement over removing pixels at random, which can be negative.

Out-of-distribution detection. An uncertainty score is used as a detector of inputs from another dataset, scored by the ROC AUC and the false-positive rate at 95% true-positive rate; the ROC AUC article describes Hendrycks and Gimpel’s maximum-softmax baseline. DDU’s comparisons report the ROC AUC.

Under shift. Calibration on the test split of the training distribution says little about shifted images. Ovadia et al. (2019; not read), as the expected calibration error article reports them, found accuracy and calibration degrading with shift for every method they tested on corrupted CIFAR-10 and ImageNet, with better calibration on the in-distribution test set not usually carrying over; dataset bias describes the kinds of shift.

Best Practices

• Model both kinds when both matter. A predicted variance captures aleatoric noise; a posterior, ensemble or feature-space density captures epistemic uncertainty; add the two for a predictive interval. In the example, intervals from either part alone undercover.
• Calibrate on held-out data and check by group. Fit temperatures or variance scales on a validation set (Laves et al., 2020; Guo et al., 2017), and report calibration for subgroups such as object size, not only on average.
• Choose the prior deliberately. Far from the data, a Bayesian model’s uncertainty is its prior’s; the example’s shifted coverage moved from 0% to 100% with the prior scale.
• Test on the shift you fear. Report calibration and detection under the corruptions, sites or sensors expected in deployment, as Izmailov et al. do with corrupted CIFAR-10, following Ovadia et al.’s set-up.
• Count the cost. Ensembles and Monte Carlo methods multiply training, memory or inference; DDU’s Pascal VOC table gives 275 ms for one softmax pass, 876 ms for a three-member ensemble and 1,577 ms for five dropout passes (Mukhoti et al., 2023, Table 3).

Pitfalls

• Reading softmax confidence as epistemic uncertainty. Kristiadi et al. (2020), citing Hein et al. (2019), note that ReLU classifiers give arbitrarily high confidence as an input is scaled away from the data; Mukhoti et al. (2023) show softmax entropy taking arbitrary values on out-of-distribution data.
• Reading aleatoric uncertainty as novelty. The example’s corner disks have in-distribution aleatoric variance. A learned noise model answers for inputs it has never seen by extrapolation.
• Trusting fitted variances from an overfitted model. Variances fitted on training residuals underestimate the noise (Laves et al., 2020; the example’s $N = 40$ rows).
• Assuming an ensemble explores. Members that agree outside the data report certainty there: van Amersfoort et al.’s two-moons ensemble and the example’s bootstrap ensemble.
• Feature collapse. Distance and density methods fail when the feature extractor maps unfamiliar inputs onto familiar features (van Amersfoort et al., 2020); DDU needs its spectral normalization and residual connections, and in Mukhoti et al.’s Figure 1 networks without such a regularized feature space overlapped in-distribution and out-of-distribution densities.
• Assuming a method transfers between tasks. Poggi et al. (2020) found that dropout sampling made self-supervised monocular depth worse in every setting and ranked errors worse than at random under stereo supervision (below).

Evidence and Impact

What the posterior buys. Wenzel et al. (2020) sampled posteriors with stochastic-gradient MCMC for a ResNet-20 on CIFAR-10 and a CNN-LSTM on IMDB reviews and found the posterior predictive at temperature $T = 1$ systematically worse than point estimates from SGD, and that “cold” posteriors with $T < 1$, which overcount the data, improved it. Izmailov et al. (2021), with full-batch HMC on a ResNet-20 with filter response normalization (and Wenzel et al.‘s CNN-LSTM on IMDB), found the opposite: the Bayesian model average beat SGD and deep ensembles in accuracy and log-likelihood on CIFAR-10, CIFAR-100 and IMDB at $T = 1$, with little evidence of a cold-posterior effect once data augmentation was turned off, which they identify as its main cause. They also found deep ensembles’ predictive distributions about as close to HMC’s as standard SGLD’s and closer than mean-field variational inference’s. Under the 16 corruptions of CIFAR-10-C at five intensities, however, HMC was the least robust method they tested; deep ensembles and SGLD were consistently more robust, and at high intensities even a single SGD model beat the HMC ensemble, although the Bayesian model could detect the shift.

Single-pass against ensembles. With a Wide-ResNet-28-10 trained on CIFAR-10, Mukhoti et al. (2023, Table 1, averages over 25 runs) report ROC AUCs for detecting SVHN of 94.44 with softmax entropy, 97.86 with DDU’s feature density, 97.73 with a five-member ensemble’s predictive entropy and 97.18 with its mutual information. For ImageNet against ImageNet-O with a ResNet-50 (Table 2), the values were 51.42 for softmax entropy, 71.29 for DDU and 62.43 for a three-member ensemble’s mutual information. Van Amersfoort et al. (2020) report that DUQ improves on or matches deep ensembles for Fashion-MNIST against MNIST and CIFAR-10 against SVHN.

Optical flow. Ilg et al. (2018) compared dropout, bootstrapped and snapshot ensembles, each trained either on the endpoint error (“empirical”) or on a likelihood with a predicted variance (“predictive”), for FlowNetC on Sintel train clean with eight samples. Predictive versions generally outperformed empirical ones (their Table 2), the best uncertainties came from bootstrapped ensembles of predictive networks, and a single multi-hypothesis network with a merging network matched them at interactive frame rates.

Monocular depth. Poggi et al. (2020) evaluated the same families for self-supervised monocular depth estimation on KITTI (Eigen split, improved ground truth, Monodepth2). With supervision from monocular sequences, where pose and scale are unknown, the bootstrapped and snapshot ensembles alone did no better than chance at ranking errors, and dropout sampling ranked errors better but produced much worse depth; a variance learned by likelihood did better, and their self-teaching variant, which learns the variance against a teacher network’s depth, did better still. With stereo supervision, where the pose is known, the bootstrapped and snapshot ensembles worked well and the likelihood variance beat the self-teaching one, while dropout’s uncertainty was worse than random (negative AURG). Dropout sampling gave worse depth than the baseline under every kind of supervision (their Tables 1–3). Combining an ensemble with a predicted variance often improved on either alone.

Detection. Gaussian YOLOv3 (Choi et al., 2019) predicts a mean and variance for each YOLOv3 box coordinate, trains them with the loss attenuation of Kendall and Gal, and multiplies the detection score by one minus the average localization uncertainty. It raised mAP by 3.09 on its KITTI validation split and by 3.5 on the BDD test set, cut false positives by 41.40% and 40.62% on the KITTI and BDD validation sets, and ran at over 42 frames per second at 512 × 512. Bounding box regression covers He et al.’s KL loss and variance voting for two-stage detectors.

Segmentation. On Pascal VOC 2012 with DeepLab-v3+ and a ResNet-101, DDU matched the mIoU of the plain softmax model (78.53) and outperformed five-pass dropout and a three-member ensemble on all three of the uncertainty metrics it uses, which Mukhoti et al. take from earlier work by Mukhoti and Gal; in their examples the softmax entropy was high along object edges, an aleatoric signal, while low feature density marked a whole wrongly predicted region. Medical image segmentation lists estimating where a label map is unreliable, and relating it to inter-rater variation, among its open problems.

Medical imaging. Laves et al. (2020) combined MC dropout with a predicted variance for four medical regression tasks, among them tumour cellularity in breast histology (BreastPathQ), bone age from the RSNA pediatric bone age data and surgical instrument tracking in endoscopy (EndoVis). The raw predictive uncertainty was underestimated; The authors report that σ scaling recalibrated it in all experiments and kept accuracy; it lowered the negative log-likelihood of every model in their Table 2, although the table shows the uncertainty calibration error rising slightly in one row (ResNet-101 on EndoVis, 0.03 to 0.04). The auxiliary recalibration network did about as well in many rows but overfitted the calibration set in some (their Table 1 and Figure 3).