Evaluating generative models: bits per dimension, FID, KID, precision and recall, CRPS and calibration

Before this lesson: lesson 01 (the model families and maximum likelihood as a KL divergence) and lesson 04 (exact log-likelihoods). You also need the multivariate Gaussian, the symmetric matrix square root and the KL divergence. Helpful: probabilistic forecasting and Gaussian processes and uncertainty (scoring rules). Time: about 110 minutes, plus the notebook.

---

Why this matters #

You cannot improve what you measure badly. Image generators are ranked by FID; likelihood models by bits per dimension; weather and climate ensembles (GenCast, the ECMWF ensemble) by the CRPS, the spread–skill ratio and rank histograms. Papers compare numbers computed with different sample sizes, feature networks and estimators, and draw wrong conclusions.

For a scientific generative model, all three questions apply. A model of weather fields must give realistic samples (question 2), and its ensembles must be calibrated forecasts (question 3). After few-step distillation (lesson 12), a model can become sharper but lose diversity. Only the right metric shows it.

1. Predict first #

2. The picture #

Three kinds of evaluation #

• Likelihood treats the model as a density $p_\theta$ and averages $-\log p_\theta(x)$ over test data. This is a cross-entropy: $\mathbb E_{p_{\text{data}}}[-\log p_\theta] = H(p_{\text{data}}) + \mathrm{KL}(p_{\text{data}} \,\|\, p_\theta)$ (lesson 01). It needs a model with a density or a bound on it (flows, autoregressive models, VAEs, diffusion ELBOs).
• Sample-based distances compare a set of generated samples with a set of real data in a feature space. They need no density, so they also work for GANs.
• Forecast scores compare a predictive distribution $F$ (often given by $M$ samples, an ensemble) with the single outcome $y$ that happened, and average over many cases.

The simplest case: 1-D Gaussians #

Data $\mathcal N(0, 1)$, model $\mathcal N(\mu, \sigma^2)$.

• The expected negative log-likelihood is $\tfrac12\log(2\pi\sigma^2) + (1 + \mu^2)/(2\sigma^2)$. It is smallest at $\mu = 0$, $\sigma = 1$.
• The Fréchet distance is $\mu^2 + (\sigma - 1)^2$. It is the squared 2-Wasserstein distance between the two Gaussians. FID computes it between Gaussians fitted to features of images.
• For one outcome $y$, draw the forecast CDF $F(x)$ and the step function $\mathbb 1\{x \ge y\}$. The CRPS is the integral of the squared gap between the two curves. It is small when $F$ is sharp and centered near $y$.

Likelihood and samples measure different things #

A model can have a good likelihood and bad samples, or good samples and a bad likelihood. The figure shows three ways this happens in 1-D.

Interactive figure: loading…

3. Definitions and main results #

3.1 Likelihood and bits per dimension #

The problem. A continuous density $p$ on $\mathbb R^D$ cannot be compared directly with a probability mass function. It can put arbitrarily high density on the integer points, so the "log-likelihood" of discrete data can go to $+\infty$ and means nothing.

Uniform dequantization. Train on $x + u$ with $u \sim \mathcal U[0, 1)^D$. The density $p$ then defines a discrete model: the mass of the unit cube at $x$,

Theorem 1 (the dequantization bound) #

$P$ sums to at most 1 (to exactly 1 if $p$ lives on $[0, 256)^D$), and for every $x$

Proof.

1. The cubes $x + [0, 1)^D$, $x \in \{0, \dots, 255\}^D$, do not overlap and fill $[0, 256)^D$. So $\sum_x P(x) = \int_{[0,256)^D} p \le 1$.
2. $u$ has density 1 on the unit cube, so $P(x) = \mathbb E_u[p(x + u)]$.
3. $\log$ is concave. Jensen's inequality gives $\mathbb E_u \log p(x + u) \le \log \mathbb E_u\, p(x + u) = \log P(x)$. $\blacksquare$

So the dequantized negative log-likelihood $-\mathbb E_u \log p(x + u)$ is an upper bound on $-\log P(x)$, the code length (in nats) of a real discrete model. Maximizing the continuous likelihood on dequantized data maximizes a lower bound on the discrete log-likelihood. Variational dequantization (Flow++) learns $q(u \mid x)$ and tightens the bound to $\mathbb E_q[\log p(x + u) - \log q(u \mid x)] \le \log P(x)$.

Rescaling and bits per dimension. Models usually see $y = (x + u)/256 \in [0, 1)^D$. By the change of variables, $p_y(y) = 256^D\, p(x + u)$, so $\log p(x + u) = \log p_y(y) - D\log 256$. Divide the nats per image by $D \ln 2$ to get bits per dimension (bpd):

averaged over $u$ and over test images.

• The uniform density on $[0, 1]^D$ has $\log p_y = 0$, so $\mathrm{bpd} = 8$: the trivial code for 8-bit pixels.
• On CIFAR-10 test data: Glow 3.35, PixelCNN++ 2.92, the variational diffusion model (VDM) bound 2.65 bits/dim.
• For data scaled to $[-1, 1]$, the bin width is $2/256$, so the constant is $D \log 128$.
• A model that only gives a lower bound on $\log p_y$ (a VAE ELBO, a diffusion ELBO) gives an upper bound on the bpd. It is still valid to report, as "$\le$".

3.2 The Fréchet distance and FID #

Definition. Fit Gaussians $\mathcal N(\mu_r, \Sigma_r)$ and $\mathcal N(\mu_g, \Sigma_g)$ to the features of the real and the generated images (sample means; sample covariances with $1/(N-1)$). Then

This is the squared 2-Wasserstein distance between the two Gaussians. FID is this distance on Inception-v3 features (Heusel et al. 2017). Lower is better. The standard protocol is "FID-50k": 50,000 generated images against the statistics of the training set.

• 1-D: $\mathrm{FD} = (\mu_1 - \mu_2)^2 + (\sigma_1 - \sigma_2)^2$.
• Commuting covariances (same eigenvectors, eigenvalues $\lambda_j$ and $\kappa_j$): $\mathrm{FD} = \|\mu_1 - \mu_2\|^2 + \sum_j (\sqrt{\lambda_j} - \sqrt{\kappa_j})^2$.

Why the symmetric form. Older code computes $\mathrm{tr}\,(\Sigma_r\Sigma_g)^{1/2}$ with a general matrix square root. The product $\Sigma_r\Sigma_g$ is not symmetric. If $\Sigma_r$ is invertible, $\Sigma_r\Sigma_g = \Sigma_r^{1/2}\big(\Sigma_r^{1/2}\Sigma_g\Sigma_r^{1/2}\big)\Sigma_r^{-1/2}$ is similar to a symmetric positive semidefinite (PSD) matrix. So both have the same eigenvalues, all $\ge 0$, and the traces of their square roots agree. The symmetric form can use an eigendecomposition (eigh), which is stable and real (linear algebra lesson 16). A general sqrtm of a non-symmetric, nearly singular $2048 \times 2048$ matrix can return complex parts and lose accuracy.

Interactive figure: loading…

Theorem 2 (the small-sample bias of FID) #

Draw two independent sets of $N$ samples from the same distribution $\mathcal N(\mu, I_d)$ and compute the plug-in estimate $\widehat{\mathrm{FD}}_N$. To second order in the sampling errors,

Why this is true (intuition). The estimate is a squared distance between two noisy summaries. The noise in the two sets does not cancel; it adds. A sum of squares of noise has a positive mean.

Proof

1. Mean term. $\hat\mu_1 - \hat\mu_2 \sim \mathcal N(0, \tfrac2N I_d)$, so $\mathbb E\|\hat\mu_1 - \hat\mu_2\|^2 = 2d/N$ (exact).
2. Covariance term near $I$. Write $\hat\Sigma_k = I + E_k$. For commuting matrices the trace term is $\sum_j(\sqrt{1 + e_j} - \sqrt{1 + f_j})^2 \approx \tfrac14\sum_j(e_j - f_j)^2$, because $\sqrt{1 + e} \approx 1 + e/2$. In general, to second order, it is $\tfrac14\|E_1 - E_2\|_F^2$.
3. For the sample covariance of $\mathcal N(0, I_d)$: $\mathrm{Var}(\hat\Sigma_{ii}) = 2/(N-1)$ and $\mathrm{Var}(\hat\Sigma_{ij}) = 1/(N-1)$ for $i \ne j$. So $\mathbb E\|E_k\|_F^2 = (2d + d(d-1))/(N-1) = d(d+1)/(N-1)$.
4. $E_1, E_2$ are independent with mean 0, so $\mathbb E\|E_1 - E_2\|_F^2 = 2d(d+1)/(N-1)$. Multiply by $\tfrac14$ and add step 1. $\blacksquare$

For a general covariance $\Sigma$ with eigenvalues $\lambda_i$, the same argument (with the second-order Bures metric $\tfrac12\sum_{ij} E_{ij}^2/(\lambda_i + \lambda_j)$ in the eigenbasis) gives $2\,\mathrm{tr}\,\Sigma/N + \big(\mathrm{tr}\,\Sigma + \sum_{i \ne j}\lambda_i\lambda_j/(\lambda_i + \lambda_j)\big)/(N-1)$. The widget prints this value; a simulation with $N = 50$ gives $0.1204$ against the predicted $0.1205$.

Consequences.

• The bias decays like $1/N$ and grows like $d^2$ for white features. Real Inception features are not white, so the formula with $d = 2048$ does not give the real bias. But the real bias is still large: between the CIFAR-10 train and test sets, Bińkowski et al. (2018) measure FID ≈ 8.1 at $N = 10{,}000$, still decreasing.
• Compare FID only at the same $N$, feature network and preprocessing. No unbiased estimator of FID exists, and the bias depends on the distributions, so at small $N$ it can even change the ranking of two models. Chong and Forsyth (2020) extrapolate linearly in $1/N$ ($\mathrm{FID}_\infty$).
• FID assumes Gaussian features, but Inception features are ReLU outputs, many of them exactly zero.

3.3 KID: an unbiased kernel distance #

Maximum mean discrepancy (MMD). For a kernel $k$ and independent $X, X' \sim P$, $Y, Y' \sim Q$:

It is the squared distance between the mean feature vectors of $P$ and $Q$ in the feature space of the kernel. KID (Bińkowski et al. 2018) uses the same Inception features and the cubic polynomial kernel $k(a, b) = (a^\top b/d + 1)^3$. Its feature space contains all monomials up to degree 3, so KID compares means, covariances and third moments. (This kernel is not characteristic: $\mathrm{MMD} = 0$ only says that the moments up to order 3 agree.)

The unbiased estimator. With samples $x_1, \dots, x_m \sim P$ and $y_1, \dots, y_n \sim Q$:

Theorem 3 (KID is unbiased) #

$\mathbb E\,\widehat{\mathrm{MMD}}{}^2_u = \mathrm{MMD}^2(P, Q)$ for every $m, n \ge 2$.

Proof

1. For $i \ne j$, $x_i$ and $x_j$ are independent copies, so $\mathbb E\,k(x_i, x_j) = \mathbb E\,k(X, X')$. The first sum has $m(m-1)$ such terms. The same holds for the second sum.
2. Every cross pair $(x_i, y_j)$ is independent, so each of the $mn$ terms has mean $\mathbb E\,k(X, Y)$.
3. Linearity of expectation gives the result. The diagonal terms $k(x_i, x_i)$ are left out because $\mathbb E\,k(X, X) \ne \mathbb E\,k(X, X')$ in general. $\blacksquare$

Consequences. When $P = Q$ the expectation is 0, so single estimates scatter around 0 and are often negative. A negative KID is normal, not a bug. Values at different $N$ are comparable. The cost is $O(N^2 d)$, so in practice KID is averaged over many random subsets (for example of 1,000 samples each).

3.4 Precision and recall #

A model can fail in two ways: low fidelity (samples outside the data manifold) and low diversity (modes of the data are missing). One number such as FID cannot show which.

• Precision: the share of generated samples that lie in the support of the real distribution.
• Recall: the share of real samples that lie in the support of the model distribution.

The k-NN estimator (Kynkäänniemi et al. 2019). In a feature space, approximate the real support by the union of balls around the real features, each with radius equal to the distance to its $k$-th nearest other real feature. Do the same for the generated features. Precision is the share of generated points inside the real union; recall is the share of real points inside the generated union. The paper uses VGG-16 features, $k = 3$ and 50,000 samples per set.

Interactive figure: loading…

Weakness. An outlier has a far $k$-th neighbor, so it gets a huge ball. Outliers among the generated samples inflate recall, and outliers among the real samples inflate precision (section 4 shows a 1-D case). "Density and coverage" (Naeem et al. 2020) are more robust variants.

3.5 Why likelihood and sample quality can disagree #

• Good likelihood, bad samples (Theis, van den Oord and Bethge 2016). Mix a good model $p$ with a bad one $q$: $0.01\,p + 0.99\,q$. Then $\log(0.01\,p(x) + 0.99\,q(x)) \ge \log p(x) - \log 100$. The log-likelihood drops by at most $\log 100 = 4.61$ nats per image, which is $0.0022$ bits/dim for $D = 3072$. But 99% of the samples come from the bad model. (The first figure shows this.)
• Good samples, bad likelihood. A model that outputs near-copies of training images has perfect samples but a very poor test likelihood. GANs have no likelihood at all.
• The likelihood is dominated by low-level detail. Most bits encode local texture and noise, not semantics. Flows trained on CIFAR-10 give higher likelihood to SVHN images than to CIFAR-10 test images (Nalisnick et al. 2019).
• FID mixes both failure types. Mode dropping (low recall) and poor fidelity (low precision) can give the same FID. Report precision and recall next to it.

3.6 Proper scoring rules and the CRPS #

Setting. A forecast is a distribution $F$; then the outcome $y$ happens. A scoring rule $S(F, y)$ is a number, here negatively oriented. Let $S(F, G) = \mathbb E_{Y \sim G}\, S(F, Y)$ be the expected score when the truth is $G$. $S$ is proper if $S(G, G) \le S(F, G)$ for all $F$, and strictly proper if equality holds only for $F = G$. In words: reporting what you believe is the best strategy.

• The log score $-\log f(y)$ is strictly proper, but it needs a density $f$.
• The squared error of the forecast mean is proper but not strictly proper: it ignores the spread.
• The continuous ranked probability score (CRPS) is strictly proper on distributions with a finite mean:

It has the units of $y$, and it becomes the absolute error $|y - \mu|$ for a point forecast at $\mu$.

Theorem 4 (the kernel form of the CRPS) #

If $X, X'$ are independent with CDF $F$ and $\mathbb E|X| < \infty$, then

Proof.

1. For real $a, b$: $\int (\mathbb 1\{a \le x\} - \mathbb 1\{b \le x\})^2\, dx = |a - b|$, because the integrand is 1 exactly for $x$ between $a$ and $b$, and 0 elsewhere.
2. Take $(a, b) = (X, y)$, then the expectation. The integrand is non-negative, so we may swap $\mathbb E$ and $\int$ (Tonelli). With $\mathbb 1^2 = \mathbb 1$ and $\mathbb E\,\mathbb 1\{X \le x\} = F(x)$: $\mathbb E|X - y| = \int \big[F(x) - 2F(x)\mathbb 1\{y \le x\} + \mathbb 1\{y \le x\}\big]\, dx$.
3. Take $(a, b) = (X, X')$. By independence $\mathbb E[\mathbb 1\{X \le x\}\mathbb 1\{X' \le x\}] = F(x)^2$, so $\mathbb E|X - X'| = \int \big[2F(x) - 2F(x)^2\big]\, dx$.
4. Both integrals are finite because $\mathbb E|X| < \infty$, so we may subtract them: $\mathbb E|X - y| - \tfrac12\mathbb E|X - X'| = \int \big[F^2 - 2F\,\mathbb 1\{y \le x\} + \mathbb 1\{y \le x\}\big]\, dx = \int (F - \mathbb 1\{y \le x\})^2\, dx$. $\blacksquare$

Propriety in two lines. For $Y \sim G$ and fixed $x$, $\mathbb E(F(x) - \mathbb 1\{Y \le x\})^2 = (F(x) - G(x))^2 + G(x)(1 - G(x))$ (bias plus variance of a Bernoulli). Integrate: $S(F, G) = \int (F - G)^2\, dx + \int G(1 - G)\, dx$. This is smallest exactly when $F = G$.

The Gaussian closed form. For $F = \mathcal N(\mu, \sigma^2)$ and $z = (y - \mu)/\sigma$, let $\varphi, \Phi$ be the standard normal density and CDF. $X - y \sim \mathcal N(-\sigma z, \sigma^2)$, whose absolute value has mean $\sigma[z(2\Phi(z) - 1) + 2\varphi(z)]$ (the folded normal). $X - X' \sim \mathcal N(0, 2\sigma^2)$, so $\mathbb E|X - X'| = 2\sigma/\sqrt\pi$. With Theorem 4:

Ensemble estimators. For $M$ members $x_1, \dots, x_M$:

The standard estimator is exactly the CRPS of the empirical CDF of the members (Theorem 4 applied to the empirical distribution). The fair one is from Ferro (2014).

Interactive figure: loading…

Theorem 5 (the standard estimator is biased; the fair one is not) #

If the members are i.i.d. draws from $F$, independent of $y$, then

Proof

1. $\mathbb E\,\tfrac1M\sum_i |x_i - y| = \mathbb E|X - y|$.
2. In $\sum_{i,j}|x_i - x_j|$, the $M$ diagonal terms are 0, and each of the $M(M-1)$ other terms has mean $\mathbb E|X - X'|$.
3. So the standard pair term has mean $\tfrac{M-1}{2M}\mathbb E|X - X'| = \tfrac12\mathbb E|X - X'| - \tfrac{1}{2M}\mathbb E|X - X'|$, and the fair pair term has mean $\tfrac12\mathbb E|X - X'|$. Use Theorem 4. $\blacksquare$

For a Gaussian $F$ the bias is $\sigma/(M\sqrt\pi)$: $0.113\,\sigma$ for $M = 5$.

• The standard estimator punishes small ensembles, so a 10-member and a 50-member ensemble are compared unfairly. As a training loss, its expected spread term is $\tfrac{M-1}{M}\cdot\tfrac12\mathbb E|X - X'|$ instead of $\tfrac12\mathbb E|X - X'|$. It rewards spread too little, so its optimum is under-dispersed.
• The fair version has a degenerate case: with $M = 2$ it is 0 whenever $y$ lies between the two members. Some weather models (for example ECMWF's AIFS-CRPS) train with an "almost fair" mix of the two.

3.7 Multivariate fields: energy score and variogram score #

For a field $y \in \mathbb R^D$ (values at $D$ sites), the average per-site CRPS sees only the marginal distributions. An ensemble with correct marginals but wrong spatial correlation gets the same score.

Energy score (Gneiting and Raftery 2007, section 4.3): $\mathrm{ES}(F, y) = \mathbb E\|X - y\| - \tfrac12\mathbb E\|X - X'\|$ with the Euclidean norm. It is strictly proper on distributions with a finite mean, and it equals the CRPS for $D = 1$. Its ensemble estimators have the same form as the CRPS estimators, with norms.

Weakness. The energy score detects errors in the mean well, errors in the variance less well, and errors in the correlation structure only weakly, more so in high dimension (Pinson and Tastu 2013; Scheuerer and Hamill 2015). The reason: $\|X - y\|$ adds up the errors at all $D$ sites, and a change in correlation changes this sum very little compared with its sampling noise.

Variogram score of order $p$ (Scheuerer and Hamill 2015):

It compares the expected differences between pairs of sites (the variogram) with the observed differences. $p = 0.5$ is the usual choice, and the weights $w_{ij}$ can decay with distance. It is proper but not strictly proper: it sees only differences, so it is blind to a shift of the whole field. Use it together with the CRPS or the energy score.

3.8 Calibration: rank histograms and the spread–skill ratio #

Calibration means that forecasts and outcomes are statistically consistent. For an ensemble: if the outcome behaves like one more member (it is exchangeable with the members), then its rank $r$ = (number of members below $y$) is uniform on $\{0, \dots, M\}$. Reason: $y$ is equally likely to be in any of the $M + 1$ positions of the sorted list of all $M + 1$ values.

Rank histogram (Talagrand diagram): count the ranks over many cases, in $M + 1$ bins.

• Flat: consistent with calibration.
• U shape: the outcome often falls outside the ensemble: under-dispersion (too narrow), or conditional biases.
• Dome: over-dispersion (too wide).
• Slope: a bias. Many low ranks mean that the members are too high.

A flat histogram is necessary, not sufficient (Hamill 2001). A climatological ensemble is calibrated but not sharp, and opposite biases in different regions can cancel. The goal is maximal sharpness subject to calibration (Gneiting, Balabdaoui and Raftery 2007), so look at proper scores too.

Theorem 6 (spread and skill of a perfect ensemble) #

Given the information at forecast time, let $y, x_1, \dots, x_M$ be i.i.d. with mean $\mu$ and variance $\sigma^2$. Let $\bar x$ be the ensemble mean and $s^2$ the ensemble variance with $1/(M-1)$. Then

Proof

$\bar x - y = (\bar x - \mu) - (y - \mu)$. The two parts are independent with mean 0 and variances $\sigma^2/M$ and $\sigma^2$, so the variance of the difference is $\sigma^2(1 + 1/M)$. The variance estimator with $1/(M-1)$ is unbiased. $\blacksquare$

So spread / (RMSE of the ensemble mean) tends to $\sqrt{M/(M+1)}$ for a perfect ensemble: $0.913$ for $M = 5$, $0.953$ for $M = 10$. The corrected spread–skill ratio (Fortin et al. 2014) is

over $N$ cases. It is about 1 for a calibrated ensemble, below 1 when under-dispersed and above 1 when over-dispersed. Average the variances and the squared errors over cases before taking roots. If the members have the right mean but a spread $c\sigma$ instead of $\sigma$, then $\mathbb E\,s^2 = c^2\sigma^2$ and $\mathbb E(\bar x - y)^2 = \sigma^2 + c^2\sigma^2/M$, so $\mathrm{SSR} = c\sqrt{(M+1)/(M + c^2)}$.

Interactive figure: loading…

Answers to section 1. (1) Positive. A Fréchet distance is never negative, and the estimate is not always 0, so its mean is positive: about $2d/N + d(d+1)/(2(N-1))$ for white features (Theorem 2). KID is unbiased (Theorem 3): its expectation is 0, and single values are often negative. (2) Above, by $\mathbb E|X - X'|/(2M)$, which is $0.113\,\sigma$ for a Gaussian with $M = 5$ (Theorem 5). The uncorrected spread–skill ratio is about $\sqrt{5/6} = 0.913$ (Theorem 6).

The family table of lesson 01, with how to evaluate

family	likelihood	how to evaluate
autoregressive models, normalizing flows	exact	bits/dim on dequantized test data; plus sample metrics
VAEs, diffusion models	lower bound (ELBO)	bits/dim as an upper bound ("≤"); FID/KID, precision and recall
GANs, few-step distilled models	none	FID/KID, precision and recall; nearest-neighbor checks for memorization
energy-based models	unnormalized	sample metrics; likelihood only with an estimate of $Z_\theta$
forecast ensembles (any family)	—	fair CRPS, energy and variogram scores, rank histograms, spread–skill ratio

4. Worked example #

(a) Bits per dimension. 8-bit grey images with $D = 784$, rescaled to $[0, 1]$. A flow gives $\log p_y(y) = 3800$ nats (averaged over $u$). A positive log-density is normal on $[0, 1]^D$.

1. $D \ln 256 = 784 \cdot 5.5452 = 4347.4$.
2. $-\log P(x) \le -3800 + 4347.4 = 547.4$ nats.
3. $\mathrm{bpd} = 547.4/(784 \ln 2) = 547.4/543.4 = 1.007$.

(b) Fréchet distance. $\mathcal N(0, I)$ against $\mathcal N\big((1, 1), \begin{pmatrix}2&1\\1&2\end{pmatrix}\big)$.

1. Mean term: $\|(1, 1)\|^2 = 2$.
2. $\Sigma_1 = I$ commutes with $\Sigma_2$, whose eigenvalues are 3 and 1. Trace term: $(1 - \sqrt3)^2 + (1 - 1)^2 = 0.5359$.
3. $\mathrm{FD} = 2.5359$.

(c) FID bias in 1-D. Two sets of $N = 50$ from the same $\mathcal N(0, 1)$. Mean term $2/N = 0.04$; standard-deviation term $\approx d(d+1)/(2(N-1)) = 1/49 = 0.0204$. So $\mathbb E\,\widehat{\mathrm{FD}} \approx 0.060$ (a simulation with 200,000 repeats gives $0.0602$).

(d) KID by hand. $d = 1$, so $k(a, b) = (ab + 1)^3$. Real $\{0, 1\}$, generated $\{1, 2\}$.

1. Within real ($i \ne j$): $k(0, 1) = 1$, mean 1.
2. Within generated: $k(1, 2) = 27$, mean 27.
3. Cross: $k(0,1) = 1$, $k(0,2) = 1$, $k(1,1) = 8$, $k(1,2) = 27$; mean $37/4 = 9.25$.
4. $\widehat{\mathrm{MMD}}{}^2_u = 1 + 27 - 2 \cdot 9.25 = 9.5$.

(e) Precision and recall in 1-D, with $k = 1$. Real $\{0, 1, 2, 3\}$: every radius is 1, so the real union is $[-1, 4]$. Generated $\{0.5, 1.5, 10\}$: radii 1, 1 and 8.5 (the nearest neighbor of 10 is 1.5), so the generated union is $[-0.5, 18.5]$.

1. Precision: 0.5 and 1.5 are inside $[-1, 4]$, 10 is not: $2/3$.
2. Recall: 0, 1, 2, 3 are all inside $[-0.5, 18.5]$: $1$.
3. Without the outlier 10, the generated union is $[-0.5, 2.5]$ and the real point 3 is outside: recall $3/4$. The outlier raised recall from 0.75 to 1.

(f) CRPS. Forecast $\mathcal N(1, 2^2)$, outcome $y = 0$, so $z = -0.5$, $\Phi(-0.5) = 0.3085$, $\varphi(0.5) = 0.3521$: $\mathrm{CRPS} = 2\,[-0.5\,(2 \cdot 0.3085 - 1) + 2 \cdot 0.3521 - 0.5642] = 2\,[0.1915 + 0.7041 - 0.5642] = 0.6628$.

Ensemble $\{-1, 1, 2, 4\}$, $y = 0$:

1. $\tfrac1M\sum_i |x_i - y| = (1 + 1 + 2 + 4)/4 = 2$.
2. Distances of the 6 unordered pairs: 2, 3, 5, 1, 3, 2, sum 16; over ordered pairs, 32.
3. Standard: $2 - 32/(2 \cdot 16) = 1.000$. Fair: $2 - 32/(2 \cdot 4 \cdot 3) = 0.667$.
4. The difference $0.333$ is the mean pairwise distance $32/12$ divided by $2M$. For Gaussian members with $\sigma = 2$ its expectation is $\sigma/(M\sqrt\pi) = 0.282$.

(g) Calibration. A perfect 10-member ensemble has an uncorrected ratio of $0.953$ and a corrected ratio of 1. With half the right spread ($c = 0.5$): $\mathrm{SSR} = 0.5\sqrt{11/10.25} = 0.518$. Its rank histogram is U-shaped. For $M = 4$ the bin probabilities are $(0.314, 0.130, 0.113, 0.130, 0.314)$ instead of $0.2$ each (code below).

5. In code #

NumPy and SciPy check the theory and the worked example: the dequantization bound, the Fréchet distance in its symmetric form, the bias of Theorem 2, the CRPS and its estimators, and exact rank histograms. The notebook has these cells and the exercises.

import numpy as np
from scipy.integrate import quad
from scipy.linalg import sqrtm
from scipy.stats import truncnorm

# Theorem 1 for a 4-level variable: a density on [0, 4) (a normal N(1.7, 0.6^2) truncated to [0, 4)).
p = truncnorm((0 - 1.7) / 0.6, (4 - 1.7) / 0.6, loc=1.7, scale=0.6)
for x in range(4):
    P = quad(p.pdf, x, x + 1)[0]                        # induced discrete model P(x)
    bound = quad(lambda u: p.logpdf(x + u), 0, 1)[0]    # E_u log p(x + u), u ~ U[0, 1)
    print(x, f"{bound:.4f} <= {np.log(P):.4f}")
# 0 -2.5215 <= -2.1232
# 1 -0.5770 <= -0.5601
# 2 -1.4104 <= -1.2238
# 3 -5.0215 <= -4.1929

def frechet(m1, S1, m2, S2):
    """Squared Frechet (2-Wasserstein) distance between two Gaussians, symmetric form with eigh."""
    w, Q = np.linalg.eigh(S1)
    R = Q @ np.diag(np.sqrt(np.clip(w, 0, None))) @ Q.T          # S1^{1/2}
    ev = np.clip(np.linalg.eigvalsh(R @ S2 @ R), 0, None)          # eigenvalues of S1^{1/2} S2 S1^{1/2}
    return np.sum((m1 - m2) ** 2) + np.trace(S1) + np.trace(S2) - 2 * np.sum(np.sqrt(ev))

S2 = np.array([[2.0, 1.0], [1.0, 2.0]])
print(f"{frechet(np.zeros(2), np.eye(2), np.ones(2), S2):.4f}")   # 2.5359 (worked example b)

S1 = np.array([[4.0, 0.0], [0.0, 1.0]])                           # does not commute with S2
general = np.trace(sqrtm(S1 @ S2)).real                           # what older FID code does
sym = np.trace(sqrtm(sqrtm(S1) @ S2 @ sqrtm(S1))).real
print(np.allclose(S1 @ S2, (S1 @ S2).T), f"{general:.6f} {sym:.6f}")   # False 4.114390 4.114390

# Theorem 2: two samples of size N from the SAME distribution N(0, I_d).
rng = np.random.default_rng(0)
for d, N in [(1, 50), (16, 200), (16, 1000)]:
    fd = [frechet(a.mean(0), np.atleast_2d(np.cov(a, rowvar=False)), b.mean(0), np.atleast_2d(np.cov(b, rowvar=False)))
          for a, b in (rng.standard_normal((2, N, d)) for _ in range(2000))]
    print(d, N, f"{np.mean(fd):.3f} {2 * d / N + d * (d + 1) / (2 * (N - 1)):.3f}")
# 1 50 0.059 0.060
# 16 200 0.844 0.843
# 16 1000 0.168 0.168

The CRPS, its two ensemble estimators, and rank histograms computed exactly (without simulation):

import numpy as np
from scipy.integrate import quad
from scipy.stats import binom, norm

mu, s, y = 1.0, 2.0, 0.0
F = lambda x: norm.cdf(x, mu, s)
crps_int = quad(lambda x: F(x) ** 2, -np.inf, y)[0] + quad(lambda x: (F(x) - 1) ** 2, y, np.inf)[0]
z = (y - mu) / s
crps_cf = s * (z * (2 * norm.cdf(z) - 1) + 2 * norm.pdf(z) - 1 / np.sqrt(np.pi))
print(f"{crps_int:.4f} {crps_cf:.4f}")                      # 0.6628 0.6628

def crps_ens(x, y, fair=False):
    M = len(x)
    pair = np.abs(x[:, None] - x[None, :]).sum()
    return np.mean(np.abs(x - y)) - pair / (2 * M * (M - 1) if fair else 2 * M * M)

print(crps_ens(np.array([-1.0, 1.0, 2.0, 4.0]), 0.0), round(crps_ens(np.array([-1.0, 1.0, 2.0, 4.0]), 0.0, fair=True), 4))  # 1.0 0.6667

# Theorem 5 by simulation: M = 5 members from F = N(1, 2^2), outcome y = 0.
rng = np.random.default_rng(1)
X = rng.normal(mu, s, size=(400_000, 5))
a, pair = np.abs(X - y).mean(1), np.abs(X[:, :, None] - X[:, None, :]).sum((1, 2))
print(f"{np.mean(a - pair / 50):.3f} {np.mean(a - pair / 40):.3f} {crps_cf + s / (5 * np.sqrt(np.pi)):.3f}")
# 0.889 0.663 0.888: the standard mean is biased up by s/(M sqrt(pi)); the fair mean is the exact 0.663

# Rank histograms without simulation: P(rank = r) for M = 4 members ~ N(b, c^2), outcome ~ N(0, 1).
for name, b, c in [("calibrated", 0.0, 1.0), ("too narrow", 0.0, 0.5),
                   ("too wide", 0.0, 2.0), ("biased high", 0.5, 1.0)]:
    probs = [quad(lambda u: binom.pmf(r, 4, norm.cdf(u, b, c)) * norm.pdf(u), -12, 12)[0] for r in range(5)]
    print(f"{name:12s}", np.round(probs, 3))
# calibrated   [0.2 0.2 0.2 0.2 0.2]
# too narrow   [0.314 0.13  0.113 0.13  0.314]
# too wide     [0.113 0.24  0.294 0.24  0.113]
# biased high  [0.335 0.235 0.183 0.143 0.104]

6. Where it shows up #

• Image and video generation. FID-50k with Inception-v3 is the standard table entry, with precision and recall (or density and coverage) next to it. Resizing and JPEG compression change FID (clean-fid, Parmar et al. 2022). Stein et al. (2023) propose DINOv2 features instead of Inception ($\mathrm{FD}_{\text{DINOv2}}$, reported for example by EDM2).
• Likelihood models. Flows, autoregressive models and variational diffusion models report bits/dim; a diffusion ELBO gives an upper bound (lesson 08).
• Weather and climate. WeatherBench 2 uses the fair CRPS and the spread–skill ratio. GenCast, a diffusion ensemble, is compared with the ECMWF ensemble by the CRPS, the ensemble-mean RMSE, spread–skill ratios and rank histograms.
• Spatial fields in science. Per-site CRPS of the quantities of interest against observations; energy and variogram scores for the whole field; correlation against distance; rank histograms per station; power spectra of generated fields.
• After distillation (lesson 12). A sharper few-step student can lose diversity: check recall, the spread–skill ratio and rank histograms, not only FID.
• Practice. Report bootstrap confidence intervals, check memorization with nearest training neighbors, and keep $N$, features and preprocessing fixed within one table.

7. Common mistakes #

• Forgetting $D\log 256$, or using it for data on $[-1, 1]$ (there it is $D\log 128$). Reporting nats as bits. A continuous likelihood on data that was not dequantized means nothing.
• Comparing FID across sample sizes, feature networks or preprocessing pipelines.
• sqrtm of $\Sigma_1\Sigma_2$ with the imaginary part silently dropped. Use the symmetric PSD form.
• Clipping negative KID values to 0. That makes the estimator biased again.
• Sign conventions. Gneiting and Raftery (2007) define the CRPS and the energy score with the opposite sign (positively oriented). Most code and this lesson use "lower is better".
• The standard ensemble CRPS for ensembles of different sizes, or as a training loss.
• Per-site CRPS alone for spatial fields. It cannot see correlation errors.
• Reading an uncorrected SSR of 0.91 for $M = 5$ as under-dispersion, or a flat rank histogram as proof of a good forecast.

8. Check yourself #

1. A flow is trained on 8-bit images scaled to $[-1, 1]^D$ and reports $\log p(y)$. Write the bits/dim.
   Answer

   The bin width is $2/256 = 1/128$, so $p(y) = 128^D p(x + u)$ and $\mathrm{bpd} = (-\log p(y) + D\log 128)/(D\ln 2)$.
2. Compute KID with $k(a, b) = (ab + 1)^3$ between the set $\{0, 1\}$ and itself. Why is a negative value possible?
   Answer

   Within terms: $k(0,1) = 1$, mean 1, twice. Cross term: $(k(0,0) + k(0,1) + k(1,0) + k(1,1))/4 = (1 + 1 + 1 + 8)/4 = 2.75$. Result $1 + 1 - 5.5 = -3.5$. The within terms leave out the diagonal pairs, but the cross term contains them. The estimator is unbiased, so for equal distributions its mean is 0 and single values scatter on both sides of 0.
3. Two independent 1-D samples of size $N = 200$ come from the same Gaussian. What FID do you expect? And for white 16-dimensional features?
   Answer

   1-D: $2/N + 1/(N-1) = 0.010 + 0.005 = 0.015$. For $d = 16$: $2 \cdot 16/200 + 16 \cdot 17/(2 \cdot 199) = 0.16 + 0.683 = 0.843$. Both are pure bias: the true distance is 0.
4. Two members $\{-1, 1\}$, outcome $y = 0$. Compute the standard and the fair ensemble CRPS.
   Answer

   $\tfrac1M\sum|x_i - y| = 1$; the ordered pair sum is $4$. Standard: $1 - 4/(2 \cdot 4) = 0.5$. Fair: $1 - 4/(2 \cdot 2 \cdot 1) = 0$. This is the degenerate case of the fair score: it is 0 whenever $y$ lies between two members.
5. A 10-member ensemble has the right mean but half the right spread. Which corrected and uncorrected spread–skill ratios do you expect, and what does the rank histogram look like?
   Answer

   Corrected: $0.5\sqrt{11/10.25} = 0.518$. Uncorrected: $0.5/\sqrt{1 + 0.25/10} = 0.494$. The rank histogram is U-shaped, because the outcome often falls outside the ensemble.
6. Your generated weather maps have the correct distribution at every grid point, but the points are independent, while real maps are spatially correlated. Which scores notice?
   Answer

   The per-site CRPS does not: it depends only on marginals. The energy score notices only weakly. The variogram score does, because it compares $\mathbb E|X_i - X_j|^p$ with $|y_i - y_j|^p$ for pairs of sites. A plot of correlation against distance shows it too.

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

1. Give the bits-per-dimension formula for 8-bit data rescaled to $[0, 1]^D$.

   Show answer

2. Why does the dequantized continuous NLL bound the discrete NLL?

   Show answer

3. Write the Fréchet distance used by FID.

   Show answer

4. Why is the FID estimate biased upward at small $N$?

   Show answer

5. Why can KID be negative, and why is that fine?

   Show answer

6. What do precision and recall measure for a generative model?

   Show answer

7. What is a proper scoring rule?

   Show answer

8. State the kernel form of the CRPS.

   Show answer

9. What is the bias of the standard ensemble CRPS with $M$ members?

   Show answer

10. Why does the spread–skill ratio need the factor $\sqrt{(M+1)/M}$?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why can't you compare FID numbers from two papers?" — Different $N$ (FID is biased like $1/N$), different feature extractors and preprocessing. No unbiased estimator of FID exists.
• "Why is KID sometimes negative?" — It is the unbiased U-statistic for MMD²; for equal distributions its mean is 0.
• "Model A has a better likelihood but a worse FID than model B. Is that possible?" — Yes. The mixture argument (the likelihood moves by at most $\log 100$ while 99% of samples are bad), and the likelihood is dominated by low-level detail.
• "Derive the kernel form of the CRPS." — Theorem 4: $\int(\mathbb 1\{a \le x\} - \mathbb 1\{b \le x\})^2 dx = |a - b|$, then $F^2 = \mathbb E[\mathbb 1\{X \le x\}\mathbb 1\{X' \le x\}]$.
• "How do you check that a weather ensemble is calibrated?" — Rank histograms, the corrected spread–skill ratio, the fair CRPS, and for fields the variogram score.

10. Practice #

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

1. Bits per dimension of a Gaussian density on 8-bit data with uniform dequantization; check the bound of Theorem 1 against the exact discrete probabilities.
2. The Fréchet distance with the symmetric square root; measure its bias against the sample size.
3. Unbiased KID with the cubic kernel; compare it with FID on the same features at small $N$.
4. k-NN precision and recall on mode-collapsed and over-dispersed samples of the 8-Gaussians toy.
5. CRPS: the Gaussian closed form against the standard and fair ensemble estimators; rank histograms and the corrected spread–skill ratio.

11. Watch and read #

Cornell CS 6785: Deep Generative Models. Lecture 14: Evaluating Generative ModelsVolodymyr Kuleshov · 70 min

A full lecture on likelihood-based and sample-based evaluation, including FID and the problems of each metric.

FID Explained: The Fréchet Inception Distance, Derived and Computed by Hand | Generative AI 17MyCampus · 7 min

A short derivation and hand computation of FID; compare it with section 3.2.

Module 8: Verification of probabilistic forecastsRenewables in electricity markets · 19 min

Forecast verification as used in practice: reliability, sharpness, the CRPS and rank-based checks.

• Murphy, Probabilistic Machine Learning: Advanced Topics (2023), section 20.4 "Evaluating generative models" (likelihood, distances in feature space, precision and recall).
• Theis, van den Oord and Bethge, A note on the evaluation of generative models (2016): section 3.1 (log-likelihood and dequantization), the mixture argument.
• Heusel et al., GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium (2017): the definition of FID.
• Bińkowski, Sutherland, Arbel and Gretton, Demystifying MMD GANs (2018): section 4 (the bias of FID, KID) and appendix D (why FID has no unbiased estimator).
• Kynkäänniemi et al., Improved Precision and Recall Metric for Assessing Generative Models (2019); Naeem et al., Reliable Fidelity and Diversity Metrics for Generative Models (2020).
• Stein et al., Exposing flaws of generative model evaluation metrics and their unfair treatment of diffusion models (2023).
• Gneiting and Raftery, Strictly Proper Scoring Rules, Prediction, and Estimation (JASA 2007, author's PDF): section 2 (propriety), section 4.2 (CRPS), section 4.3 (energy score).
• Rasp et al., WeatherBench 2 (2023): the fair CRPS, the spread–skill ratio and rank histograms for data-driven weather ensembles.
• The scoringrules Python library: CRPS, energy and variogram scores with their ensemble estimators.