Gaussian processes, scoring rules and calibration

Before this lesson: linear algebra lesson 07 (positive definite matrices and the Cholesky factor) and lesson 06 (Theorem 1 there: Gaussian conditioning; the Kalman update is the same formula). From lesson 01: the differential $d\log\det K = \operatorname{tr}(K^{-1}\,dK)$. From the probability track: the multivariate Gaussian. Helpful: kernel methods (for example kernel ridge regression), and time series lesson 10 for a first look at probabilistic forecasts. Time: about 100 minutes, plus the notebook.

---

Why this matters #

• Emulators of expensive simulators. A climate or ocean model run can cost thousands of node-hours. A GP trained on a few hundred runs is a cheap emulator: it predicts the model output for new parameter values, with an error bar. Emulators are used to calibrate model parameters (for example in the Calibrate–Emulate–Sample method, Cleary et al. 2021).
• Geostatistics and seismology. Kriging, the classic method to interpolate a field from scattered stations, is GP regression. Ground-motion models treat the residuals of nearby sites as a spatially correlated Gaussian field (Jayaram and Baker 2009); shake maps condition this field on the recorded stations.
• Bayesian optimisation. The next experiment or hyperparameter setting is chosen from the GP mean and variance.
• Honest uncertainty everywhere. Every probabilistic forecast must be checked. The CRPS, the PIT histogram and the spread–skill ratio score ensemble weather forecasts and diffusion models such as GenCast. Lesson 16 uses them to evaluate ML weather models.

1. Predict first #

2. The picture #

One data point #

Start with the smallest case. We observe one value $y_1$ at the input $x_1$, and we want to predict the function at a test input $x_*$. The model is

where $f$ is the unknown function and $\varepsilon$ is measurement noise. A kernel $k(x, x')$ gives the prior covariance of two function values: $\operatorname{Cov}(f(x), f(x')) = k(x, x')$. So the pair $(f(x_*), y_1)$ is a 2D Gaussian vector with mean zero and covariance

The noise adds $\sigma_n^2$ only to the variance of $y_1$: it does not change $f(x_*)$. Conditioning a 2D Gaussian (Theorem 1 of lesson 06) gives

Read it in three steps.

1. The mean copies a fraction of $y_1$. The fraction is large when $x_*$ is close to $x_1$ (large covariance $k(x_*, x_1)$) and when the noise is small.
2. The variance drops from the prior value $k(x_*, x_*)$ for the same two reasons.
3. Far from $x_1$, $k(x_*, x_1) \approx 0$: the mean goes back to 0 and the variance goes back to the prior variance. Far from the data, the GP knows only what the prior says.

This is the same formula as the scalar Kalman update of lesson 06, with the forecast replaced by the prior and the gain $k(x_*, x_1)/(k(x_1, x_1) + \sigma_n^2)$.

A distribution over functions #

A GP does this for all inputs at once. Pick any finite set of inputs $x_1, \dots, x_m$. Under a GP prior, the values $f(x_1), \dots, f(x_m)$ are jointly Gaussian, with mean 0 and covariance matrix $[k(x_i, x_j)]_{ij}$. Draw one sample of this Gaussian vector on a fine grid and connect the points: you get one random function. The most common kernel is the RBF (squared exponential) kernel

• The length scale $\ell$ is the distance over which function values stay strongly correlated. A small $\ell$ gives functions that wiggle fast; a large $\ell$ gives slow, smooth functions.
• The signal variance $\sigma_f^2 = k(x, x)$ is the prior variance of $f(x)$. It sets the typical size of the function: about 95 % of prior function values lie in $\pm 1.96\,\sigma_f$.
• The noise variance $\sigma_n^2$ is the variance of the measurement error.

Interactive figure: loading…

Scoring a forecast with one number #

Now turn the question around: the forecast is given, and the outcome arrives. How good was the forecast? A probabilistic forecast of a number is a cumulative distribution function (CDF) $F(x) = P(Y \le x)$. The observation $y$ is a single number. As a CDF, it is the step function $\mathbb 1\{x \ge y\}$: 0 to the left of $y$ and 1 from $y$ on. The continuous ranked probability score (CRPS) is the area of the squared gap between the two curves:

Lower is better. The gap is small when $F$ rises steeply (a sharp forecast) and when it rises at $y$ (a forecast centred on the outcome).

Interactive figure: loading…

3. Definitions and main results #

Gaussian processes and kernels #

A Gaussian process $f \sim \mathcal{GP}(m, k)$ is a collection of random variables $f(x)$, one for each input $x$, such that every finite subset is jointly Gaussian, with $\mathbb E f(x) = m(x)$ and $\operatorname{Cov}(f(x), f(x')) = k(x, x')$. In this lesson $m = 0$: centre the outputs first (subtract their mean), or the prior pulls predictions to 0.

A function $k$ is a valid kernel when every matrix $[k(x_i, x_j)]_{ij}$ (the Gram matrix) is symmetric positive semidefinite. Only then is it a possible covariance. With $r = |x - x'|$ (or $\|x - x'\|$ for vector inputs), two standard kernels are

The kernel sets the smoothness of the functions. RBF samples are infinitely differentiable. The second kernel belongs to the Matérn family with smoothness parameter $\nu$: a Matérn-$\nu$ process is $k$ times mean-square differentiable exactly when $\nu > k$. So Matérn-5/2 samples are twice differentiable. Matérn-1/2 is the exponential kernel $\sigma_f^2 e^{-r/\ell}$ (the Ornstein–Uhlenbeck process), whose samples are continuous but nowhere differentiable. Physical fields are often rougher than RBF assumes, so Matérn-5/2 is a common default (Rasmussen and Williams, §4.2.1).

Notation. Training inputs $X = (x_1, \dots, x_n)$, outputs $y \in \R^n$, test inputs $X_*$.

Symbol	Meaning	Size
$K_f = k(X, X)$	prior covariance of $f(X)$	$n \times n$
$K = K_f + \sigma_n^2 I$	covariance of the data $y$	$n \times n$
$K_* = k(X, X_*)$	covariance between $f(X)$ and $f_* = f(X_*)$	$n \times n_*$
$K_{**} = k(X_*, X_*)$	prior covariance of $f_*$	$n_* \times n_*$

Theorem 1 (the GP posterior) #

Write $m_*$ for the posterior mean and $\Sigma_*$ for the posterior covariance.

Proof.

1. Under the prior, $f(X) \sim \mathcal N(0, K_f)$, and the noise $\varepsilon \sim \mathcal N(0, \sigma_n^2 I)$ is independent of $f$.
2. So $(y, f_*)$ is jointly Gaussian with mean zero, $\operatorname{Cov}(y) = K_f + \sigma_n^2 I = K$, $\operatorname{Cov}(y, f_*) = K_*$ (the noise is independent of $f_*$, so it adds nothing here), and $\operatorname{Cov}(f_*) = K_{**}$.
3. Gaussian conditioning: for a zero-mean Gaussian vector $(a, b)$, $b \mid a \sim \mathcal N\big(\Sigma_{ba}\Sigma_{aa}^{-1} a,\; \Sigma_{bb} - \Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}\big)$ (the covariance is a Schur complement; linear algebra lesson 07).
4. Put $a = y$ and $b = f_*$. $\blacksquare$

Three remarks.

• The mean is a weighted sum of kernels. $m_* = \sum_i \alpha_i\, k(x_i, x_*)$ with $\alpha = K^{-1} y$. This is kernel ridge regression with ridge parameter $\lambda = \sigma_n^2$. It is also the Tikhonov / MAP solution of lesson 04, written with kernels.
• The covariance does not depend on $y$, only on where the data are.
• $f_*$ versus $y_*$. Theorem 1 is about the function value $f_*$. A new measurement $y_* = f_* + \varepsilon_*$ has one more independent noise term, so its variance is $\operatorname{diag}(\Sigma_*) + \sigma_n^2$.

Computing it: one Cholesky factor does everything #

$K$ is symmetric positive definite: its smallest eigenvalue is at least $\sigma_n^2 > 0$. So it has a Cholesky factor: $K = LL^\top$ with $L$ lower triangular and $L_{ii} > 0$. Rasmussen and Williams (Algorithm 2.1) compute everything from $L$:

1. $L = \operatorname{chol}(K)$. About $n^3/3$ flops: the only cubic step.
2. $\alpha = L^\top \backslash (L \backslash y)$. Two triangular solves, $O(n^2)$ each. ($A \backslash b$ means "solve $Ax = b$".)
3. Mean: $m_* = K_*^\top \alpha$.
4. $V = L \backslash K_*$, then $\Sigma_* = K_{**} - V^\top V$. This works because $K_*^\top K^{-1} K_* = K_*^\top L^{-\top} L^{-1} K_* = (L^{-1}K_*)^\top (L^{-1}K_*) = V^\top V$.
5. $\log\det K = 2\sum_i \log L_{ii}$, because $\det K = (\det L)^2 = \prod_i L_{ii}^2$.

Why not $K^{-1}$ and $\det K$?

• Forming $K^{-1}$ costs about three times more than the factor, and its rounding errors are larger. You never need the inverse itself, only products with it.
• $\det K$ is the product of $n$ eigenvalues. With many points, most of them are close to $\sigma_n^2 < 1$, so the product underflows to $0$ in float64. Section 5 shows this at $n = 400$, where $\log\det K \approx -860$ but det returns 0. The sum of logs has no problem.
• When $\sigma_n^2$ is tiny and $\ell$ is long, $K$ is nearly singular and the factorisation can fail. Then add a small jitter (for example $10^{-8}\sigma_f^2$) to the diagonal, and work in float64.

The cost is $O(n^3)$ time and $O(n^2)$ memory. At $n = 10^4$ the factor needs about $3 \times 10^{11}$ flops and $K$ needs 0.8 GB in float64. Larger problems need approximations (section 6).

Theorem 2 (the log marginal likelihood) #

The kernel has hyperparameters $\theta = (\sigma_f^2, \ell, \sigma_n^2)$. To choose them, look at the probability of the data under the model, with $f$ integrated out: the marginal likelihood (or evidence) $p(y \mid X, \theta)$.

Proof.

1. $f(X) \sim \mathcal N(0, K_f)$ and $\varepsilon \sim \mathcal N(0, \sigma_n^2 I)$ are independent.
2. A sum of independent Gaussian vectors is Gaussian, and the covariances add: $y = f(X) + \varepsilon \sim \mathcal N(0, K_f + \sigma_n^2 I) = \mathcal N(0, K)$. (This does the integral $\int p(y \mid f)\,p(f)\,df$ in one line.)
3. Take the log of the Gaussian density $(2\pi)^{-n/2} (\det K)^{-1/2} \exp(-\tfrac12 y^\top K^{-1} y)$. $\blacksquare$

With the factor: $\log p = -\tfrac12\, y^\top\alpha - \sum_i \log L_{ii} - \tfrac n2\log 2\pi$.

The two terms pull in opposite directions.

• A short length scale lets the GP follow every data point. The data-fit term is then close to 0 (good). But such a model could produce many different data sets: $K$ is close to $(\sigma_f^2 + \sigma_n^2) I$, its determinant is large, and the complexity penalty is large.
• A long length scale gives a simple model with a small $\det K$ (the values are strongly correlated, so the model "expects" only smooth data sets). But it fits wiggly data badly.

The maximum lies in between. The marginal likelihood prefers the simplest model that explains the data: an automatic Occam's razor. Choosing $\theta$ by maximising it is called type-II maximum likelihood (or empirical Bayes).

Interactive figure: loading…

Two local maxima are common. Rasmussen and Williams (§5.4.1, Fig. 5.5) show 7 data points with a "wiggly, low noise" optimum and a "smooth, high noise" optimum. The objective is not concave, so restart the optimiser from several initial values.

Theorem 3 (the gradient of the evidence) #

To maximise $\log p$ with a gradient method, we need its derivatives. For any hyperparameter $\theta_j$,

Proof.

1. Differentiate $K K^{-1} = I$: $dK\, K^{-1} + K\, d(K^{-1}) = 0$, so $d(K^{-1}) = -K^{-1}\, dK\, K^{-1}$.
2. Jacobi's formula (lesson 01): $d \log\det K = \operatorname{tr}(K^{-1} dK)$.
3. Data-fit term: by step 1, $d(-\tfrac12 y^\top K^{-1} y) = \tfrac12 y^\top K^{-1} dK\, K^{-1} y = \tfrac12 \alpha^\top dK\, \alpha = \tfrac12 \operatorname{tr}(\alpha\alpha^\top dK)$.
4. Complexity term: by step 2, $d(-\tfrac12 \log\det K) = -\tfrac12 \operatorname{tr}(K^{-1} dK)$.
5. Add steps 3 and 4, and divide by $d\theta_j$. $\blacksquare$

Cost. Get $K^{-1}$ from the factor once ($O(n^3)$). Then each hyperparameter costs only $O(n^2)$, because $\operatorname{tr}(AB) = \sum_{ij} A_{ij} B_{ji}$ needs no matrix product.

Use log hyperparameters. Optimise $\theta = (\log\sigma_f^2, \log\ell, \log\sigma_n^2)$. Every value is then positive without constraints, and a step changes a value by a factor, not by an amount. The chain rule gives $\partial K / \partial \log\eta = \eta\, \partial K / \partial\eta$. For the RBF kernel, with $R^2$ the matrix of squared distances $(x_i - x_j)^2$ and $\circ$ the elementwise product:

The Matérn-5/2 derivative

Write $a = \sqrt5\, r/\ell$, so $k_{5/2} = \sigma_f^2 (1 + a + a^2/3)\,e^{-a}$. Then $dk/da = \sigma_f^2\big[(1 + 2a/3) - (1 + a + a^2/3)\big]e^{-a} = -\sigma_f^2\,\tfrac{a}{3}(1 + a)\,e^{-a}$, and $da/d\log\ell = -a$. So

$$\frac{\partial k_{5/2}}{\partial \log\ell} = \sigma_f^2\, \frac{a^2}{3}\,(1 + a)\, e^{-a} = \sigma_f^2\, \frac{5r^2}{3\ell^2}\Big(1 + \frac{\sqrt5\, r}{\ell}\Big) e^{-\sqrt5 r/\ell}.$$

The derivative in $\log\sigma_f^2$ is $k_{5/2}$ itself, as for every kernel of the form $\sigma_f^2 \times (\text{correlation})$. Check such formulas against finite differences (the notebook does this for the RBF kernel).

Proper scoring rules #

Now the second half: how to judge a probabilistic forecast. A scoring rule $S(F, y)$ gives a number to the forecast $F$ when the outcome $y$ happens. Here it is negatively oriented: lower is better. If the outcome is really drawn from a distribution $G$, the expected score is $S(F, G) = \mathbb E_{Y \sim G}\, S(F, Y)$.

Definition. $S$ is proper if $S(G, G) \le S(F, G)$ for all $F$ and $G$. It is strictly proper if, in addition, equality holds only for $F = G$.

In words: if you believe $G$, the best thing to report is $G$ itself. No hedging (making the forecast wider or narrower than you believe) can improve your expected score.

• The log score $-\log p_F(y)$ (the negative log-likelihood) is strictly proper. It needs a density, and it punishes a too-narrow forecast very hard when the outcome falls in its tail.
• The squared error of the forecast mean, $(\mu_F - y)^2$, is proper but not strictly proper. Its expected value $(\mu_F - \mu_G)^2 + \operatorname{Var}_G(Y)$ is the same for every $F$ with the right mean, whatever its spread. So it cannot reward a correct spread.
• The CRPS has two equal forms (Gneiting and Raftery 2007, eqs. 20–21):

The second form (the kernel form) needs only samples of $F$. So it scores an ensemble directly (lesson 16). Its first term rewards being close to $y$; its second term rewards spread, so a forecast cannot win by being narrow alone.

The Gaussian CRPS. For $F = \mathcal N(\mu, \sigma^2)$, with $z = (y - \mu)/\sigma$ and $\varphi, \Phi$ the standard normal density and CDF,

The CRPS has the units of $y$. For a point forecast ($\sigma \to 0$) it becomes $|y - \mu|$, the absolute error: the CRPS generalises the absolute error to probabilistic forecasts.

Where the Gaussian formula comes from

Use the kernel form. For $X \sim \mathcal N(\mu, \sigma^2)$, $X - y \sim \mathcal N(\mu - y, \sigma^2)$, and the mean of $|W|$ for $W \sim \mathcal N(m, \sigma^2)$ is $\sigma\,[2\varphi(m/\sigma) + (m/\sigma)(2\Phi(m/\sigma) - 1)]$. With $m/\sigma = -z$ and the symmetry $\varphi(-z) = \varphi(z)$, $2\Phi(-z) - 1 = -(2\Phi(z) - 1)$, this gives $\mathbb E|X - y| = \sigma\,[2\varphi(z) + z(2\Phi(z) - 1)]$. For the second term, $X - X' \sim \mathcal N(0, 2\sigma^2)$, so $\mathbb E|X - X'| = \sqrt2\,\sigma \cdot \sqrt{2/\pi} = 2\sigma/\sqrt\pi$, and half of it is $\sigma/\sqrt\pi$. Subtract. $\blacksquare$

Theorem 4 (the CRPS is strictly proper) #

If $Y \sim G$, and $F$ and $G$ have finite first moments, then

The second term does not depend on $F$. The first term is never negative and is zero only for $F = G$. So the honest forecast $F = G$ is the unique best one.

Why it is true, in one picture. Fix one point $x$. Whether the outcome falls left or right of $x$ is a coin toss: heads ($Y \le x$) with probability $G(x)$. The forecast says "heads with probability $F(x)$". The squared error of a probability for a coin toss (the Brier score) is smallest on average when the stated probability is the true one. The CRPS adds up these Brier scores over all thresholds $x$.

Proof of Theorem 4

1. Fix $x$. Then $B = \mathbb 1\{x \ge Y\}$ is a Bernoulli variable with $P(B = 1) = P(Y \le x) = G(x)$.
2. Since $B^2 = B$: $\mathbb E(F(x) - B)^2 = F(x)^2 - 2F(x)G(x) + G(x) = (F(x) - G(x))^2 + G(x)(1 - G(x))$.
3. Integrate over $x$ and swap the integral and the expectation. Both integrands are nonnegative, so this is allowed (Tonelli's theorem).
4. The first term is $\ge 0$. It is $0$ only when $F(x) = G(x)$ for almost every $x$, and two right-continuous CDFs that agree almost everywhere agree everywhere. So $F = G$ is the unique minimiser: the CRPS is strictly proper. $\blacksquare$

Propriety is a statement about the expected score. On a single observation, a wrong forecast can score better than the right one; section 4 shows such a case. Always average a score over many cases before you compare forecasts.

Calibration and sharpness #

A proper score gives one number. To see what is wrong with a forecast, look at calibration.

• PIT. The probability integral transform of an observation is $u = F(y)$: the forecast probability that the outcome is below the value that happened. If $y$ really comes from the continuous distribution $F$, then $u$ is uniform on $(0, 1)$. A forecasting system is probabilistically calibrated when its PIT values, over many forecasts, are uniform. Draw them as a histogram:
  • flat: calibrated;
  • U shape: too many outcomes in the tails, so the forecasts are too narrow (overconfident, under-dispersed);
  • hump in the middle: the forecasts are too wide (underconfident, over-dispersed);
  • slope: a bias; the outcomes fall more often on one side.
• Coverage. The central interval at level $p$ contains $y$ exactly when $|u - \tfrac12| \le p/2$. So coverage at 50 %, 80 % and 95 % checks the PIT distribution at three places. One level is not enough: a Gaussian forecast for heavy-tailed data can cover 95 % correctly and still be wrong at 50 %.
• Sampling error. Observed coverage on $n$ independent test cases has standard error $\sqrt{p(1 - p)/n}$: $0.0069$ for $p = 0.95$ and $n = 1000$. Test cases that share one weather situation or one sample path are correlated, so the effective $n$ is smaller.
• Sharpness is how concentrated the forecasts are, for example the mean predictive standard deviation. It depends on the forecasts alone, not on the outcomes.

The rule of Gneiting, Balabdaoui and Raftery (2007) is: maximise sharpness subject to calibration. A proper score such as the CRPS rewards both at once.

Interactive figure: loading…

Answers to section 1.

1. The prior: mean $0$, and variance $\sigma_f^2$ for $f$ (plus $\sigma_n^2$ for a new measurement). With $\ell = 1$, $k(100, x_i) \le \sigma_f^2 e^{-90^2/2} \approx 0$ for every data point, so the data have no effect. This is why you centre $y$ (or add a mean function).
2. No. Both are calibrated, so coverage cannot separate them, but B is sharper, and a proper score can see this. In the example of Gneiting, Balabdaoui and Raftery, $y \sim \mathcal N(\mu, 1)$ with $\mu \sim \mathcal N(0, 1)$. B knows $\mu$ and issues $\mathcal N(\mu, 1)$. A issues the climatology $\mathcal N(0, 2)$. For a calibrated Gaussian forecast the expected CRPS is $\tfrac12\,\mathbb E|X - X'| = \sigma/\sqrt\pi$ (Theorem 4 with $F = G$). So B scores $1/\sqrt\pi = 0.564$ and A scores $\sqrt2/\sqrt\pi = 0.798$ on average. (The paper's simulation, Table 5, reports 0.56 and 0.78.)
3. Fewer. The band of $f$ leaves out the noise variance $\sigma_n^2$, so it is too narrow for new measurements. Section 4 computes 82 % instead of 95 % at one test point; section 5 finds 40 % on held-out data with many training points.

4. Worked example #

Two training points $x = (0, 1)$ with outputs $y = (1, 0)$. RBF kernel with $\sigma_f^2 = 1$ and $\ell = 1$; noise $\sigma_n^2 = 0.1$. Test input $x_* = 0.5$.

Step 1: the matrix $K$. $k(0, 0) = k(1, 1) = 1$ and $k(0, 1) = e^{-1/2} = 0.6065$. Add the noise to the diagonal:

Step 2: the Cholesky factor. $L_{11} = \sqrt{1.1} = 1.0488$, $L_{21} = 0.6065/1.0488 = 0.5783$, $L_{22} = \sqrt{1.1 - 0.5783^2} = \sqrt{0.7656} = 0.8750$:

Step 3: $\alpha = K^{-1}y$ by two triangular solves.

• Forward, $Lw = y$: $w_1 = 1/1.0488 = 0.9535$; $w_2 = (0 - 0.5783 \times 0.9535)/0.8750 = -0.6302$.
• Backward, $L^\top \alpha = w$: $\alpha_2 = -0.6302/0.8750 = -0.7202$; $\alpha_1 = (0.9535 - 0.5783 \times (-0.7202))/1.0488 = 1.3062$.

Step 4: the posterior at $x_* = 0.5$. With one test point, $K_*$ is a vector: $k_* = (k(0, 0.5), k(1, 0.5)) = (e^{-1/8}, e^{-1/8}) = (0.8825, 0.8825)$.

• Mean: $m_* = k_*^\top \alpha = 0.8825 \times (1.3062 - 0.7202) = 0.8825 \times 0.5860 = 0.5171$.
• $v = L^{-1} k_*$: $v_1 = 0.8825/1.0488 = 0.8414$; $v_2 = (0.8825 - 0.5783 \times 0.8414)/0.8750 = 0.4525$.
• $\operatorname{Var}[f_*] = k(x_*, x_*) - \|v\|^2 = 1 - (0.7080 + 0.2047) = 0.0873$.
• For a new measurement: $\operatorname{Var}[y_*] = 0.0873 + 0.1 = 0.1873$, standard deviation $0.4327$.

Step 5: the evidence.

• Data fit: $-\tfrac12 y^\top \alpha = -\tfrac12 (1 \times 1.3062 + 0 \times (-0.7202)) = -0.6531$.
• Complexity: $-\sum_i \log L_{ii} = -(\log 1.0488 + \log 0.8750) = -(0.0477 - 0.1336) = +0.0859$. It is positive here because $\det K = 0.8421 < 1$.
• Constant: $-\tfrac22 \log 2\pi = -1.8379$.
• Total: $\log p(y) = -0.6531 + 0.0859 - 1.8379 = -2.4051$.

Step 6: one gradient. $\partial K/\partial\log\sigma_n^2 = \sigma_n^2 I$, so Theorem 3 gives $\tfrac12 \sigma_n^2\,(\|\alpha\|^2 - \operatorname{tr} K^{-1}) = 0.05 \times (2.2250 - 2.6125) = -0.0194$. The same formula gives $-0.3275$ for $\log\sigma_f^2$ and $-0.1338$ for $\log\ell$; section 5 and the notebook check all three against finite differences. All three are negative: at this $\theta$, the evidence would grow with a smaller signal variance, a shorter length scale and less noise.

Step 7: score one observation. Suppose $y_* = 0.8$ is measured at $x_* = 0.5$.

• Correct forecast $\mathcal N(0.5171, 0.4327^2)$: $z = (0.8 - 0.5171)/0.4327 = 0.6537$, CRPS $= 0.1724$.
• Latent-only forecast $\mathcal N(0.5171, 0.2954^2)$ (it forgets $\sigma_n^2$): $z = 0.9575$, CRPS $= 0.1695$.

The wrong forecast wins on this one observation. The log score agrees: $0.158$ for the wrong one against $0.295$ for the correct one. This is not a contradiction: propriety is about the expected score. If the correct forecast is the truth, the latent-only 95 % interval $0.5171 \pm 1.96 \times 0.2954$ contains $y_*$ with probability $2\Phi(1.96 \times 0.2954/0.4327) - 1 = 0.82$, not $0.95$. Over many test points this shows up, as the next section shows.

5. In code #

NumPy and SciPy check the worked example. The notebook for this lesson has these cells, the scikit-learn cross-check and the exercises.

import numpy as np
from scipy.linalg import cho_solve, solve_triangular
from scipy.stats import norm

def rbf(a, b, sf2=1.0, ell=1.0):
    return sf2 * np.exp(-(a[:, None] - b[None, :]) ** 2 / (2 * ell**2))

# (a) The worked example with one Cholesky factor
X, y, xs, sn2 = np.array([0.0, 1.0]), np.array([1.0, 0.0]), np.array([0.5]), 0.1
K = rbf(X, X) + sn2 * np.eye(2)
L = np.linalg.cholesky(K)                       # K = L L^T
alpha = cho_solve((L, True), y)                 # K^{-1} y by two triangular solves
ks = rbf(X, xs)
v = solve_triangular(L, ks, lower=True)
mean = ks.T @ alpha
var_f = rbf(xs, xs) - v.T @ v
print(alpha.round(4), mean.round(4), var_f.round(4), (var_f + sn2).round(4))
# [ 1.3062 -0.7202] [0.5171] [[0.0873]] [[0.1873]]

# (b) The log marginal likelihood: data fit, complexity, constant
fit = -0.5 * y @ alpha
comp = -np.sum(np.log(np.diag(L)))              # -1/2 log det K
const = -len(y) / 2 * np.log(2 * np.pi)
print(round(fit, 4), round(comp, 4), round(fit + comp + const, 4))
# -0.6531 0.0859 -2.4051

# (c) The CRPS of the two forecasts for y* = 0.8
def crps_gauss(mu, s, y):
    z = (y - mu) / s
    return s * (z * (2 * norm.cdf(z) - 1) + 2 * norm.pdf(z) - 1 / np.sqrt(np.pi))
sd_y, sd_f = np.sqrt(var_f + sn2).item(), np.sqrt(var_f).item()
print(round(crps_gauss(mean.item(), sd_y, 0.8), 4), round(crps_gauss(mean.item(), sd_f, 0.8), 4))
# 0.1724 0.1695

scikit-learn gives the same numbers and the gradient. Its kernel.theta holds $(\log\sigma_f^2, \log\ell, \log\sigma_n^2)$, the log hyperparameters of Theorem 3. The last lines show why $\log\det K$ must come from a factor:

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel as C, RBF, WhiteKernel

kernel = C(1.0) * RBF(1.0) + WhiteKernel(0.1)
gp = GaussianProcessRegressor(kernel, optimizer=None).fit(X[:, None], y)
print(gp.log_marginal_likelihood_value_)               # -2.4051
th, h = kernel.theta, 1e-6
_, grad = gp.log_marginal_likelihood(th, eval_gradient=True)
fd = [(gp.log_marginal_likelihood(th + h * e) - gp.log_marginal_likelihood(th - h * e)) / (2 * h)
      for e in np.eye(3)]
print(grad, np.array(fd))                              # [-0.3275 -0.1338 -0.0194] twice

rng = np.random.default_rng(0)
Kb = kernel(np.sort(rng.uniform(0, 10, 400))[:, None])
print(np.linalg.det(Kb), np.linalg.slogdet(Kb)[1])     # 0.0 and about -860

Calibration on held-out data. Data come from the model itself ($\sigma_f^2 = 1$, $\ell = 1$, $\sigma_n^2 = 0.09$): four independent functions, each with 150 training and 1000 test points. The GP uses the true hyperparameters. PIT values come from norm.cdf; coverage at level $p$ is the fraction with $|u - \tfrac12| \le p/2$:

rng = np.random.default_rng(1)
z_y, z_f = [], []
for _ in range(4):
    Xa = rng.uniform(0, 10, 1150)[:, None]
    Kf = (C(1.0) * RBF(1.0))(Xa)
    ya = np.linalg.cholesky(Kf + 1e-8 * np.eye(1150)) @ rng.standard_normal(1150)
    ya += 0.3 * rng.standard_normal(1150)
    gp = GaussianProcessRegressor(C(1.0) * RBF(1.0) + WhiteKernel(0.09), optimizer=None)
    gp.fit(Xa[:150], ya[:150])
    mu, sd_y = gp.predict(Xa[150:], return_std=True)    # sd for new measurements y_*
    sd_f = np.sqrt(sd_y**2 - 0.09)                     # sd of f_* only: the wrong one for y_*
    z_y.append((ya[150:] - mu) / sd_y)
    z_f.append((ya[150:] - mu) / sd_f)

for name, z in [("y_* (with sn2)", z_y), ("f_* only      ", z_f)]:
    pit = norm.cdf(np.concatenate(z))                 # uniform on (0, 1) if calibrated
    cover = [np.mean(np.abs(pit - 0.5) <= p / 2) for p in (0.5, 0.8, 0.95)]
    print(name, np.round(cover, 3), np.histogram(pit, bins=5, range=(0, 1))[0])
# y_* (with sn2) [0.504 0.812 0.953] [799 831 783 809 778]   flat: calibrated
# f_* only       [0.147 0.28  0.404] [1659 256 213 250 1622]  U-shaped: far too narrow

With 150 points on $[0, 10]$ the variance of $f_*$ is much smaller than $\sigma_n^2 = 0.09$, so the latent-only "95 % interval" covers only 40 % of new measurements.

6. Where it shows up #

• Machine learning.
  • Bayesian optimisation picks the next hyperparameter setting or experiment from the GP mean and variance (for example with the expected-improvement rule).
  • Kernel ridge regression is the GP mean.
  • Infinitely wide neural networks are GPs (Neal 1996; the NNGP and NTK limits).
  • Deep ensembles and MC dropout are judged with the same PIT, coverage and proper-score tools.
• ML for science.
  • Emulators of expensive simulators: GP surrogates for the calibration of climate-model parameters and for history matching; the Calibrate–Emulate–Sample method (Cleary et al. 2021).
  • Geostatistics: kriging is GP regression, used for mineral deposits, soil properties and rainfall maps.
  • Seismology: ground-motion models treat the within-event residuals at nearby sites as a spatially correlated Gaussian field, often with an exponential (Matérn-1/2) correlation (Jayaram and Baker 2009). Conditioning this field on the recorded stations gives shake maps and regional risk estimates.
  • Weather: ensemble and diffusion forecasts (GenCast) are compared by the CRPS and the spread–skill ratio; WeatherBench 2 reports both (lesson 16).
  • Data assimilation: the Kalman update is the same Gaussian conditioning as Theorem 1, and background-error covariances are kernels (lesson 06).
• Large data. Exact GPs cost $O(n^3)$ time and $O(n^2)$ memory. Beyond about $10^4$ points, use inducing points (sparse variational GPs), structured kernels, or conjugate gradients on GPUs (GPyTorch). Standardise inputs and outputs, and keep the factorisation in float64.

7. Common mistakes #

• Intervals for new measurements from the variance of $f_*$. Add $\sigma_n^2$ for $y_*$. Without it the intervals are far too narrow (40 % coverage in section 5).
• Forming $K^{-1}$ or $\det K$. Use the Cholesky factor: cho_solve, solve_triangular and $2\sum_i \log L_{ii}$.
• Not centring $y$. The zero-mean prior pulls predictions to 0 away from the data.
• Raw hyperparameters and one start. Optimise log hyperparameters, and restart from several points: the evidence can have several local maxima.
• Calibration at one level, on training data, or without an error bar. Check several levels (or the PIT histogram) on held-out data, with $\sqrt{p(1 - p)/n}$ in mind.
• Coverage as the only metric. Climatology has perfect coverage and no skill. Report a proper score as well.
• Judging from one observation. Propriety holds in expectation (section 4).
• Sign conventions. Gneiting and Raftery write scores as rewards (higher is better); their CRPS is the negative of the one here.

8. Check yourself #

1. One training point $x = 0$, $y = 1$; RBF kernel with $\sigma_f^2 = 1$, $\ell = 1$; $\sigma_n^2 = 0.25$. Give the posterior mean and variance of $f$ at $x_* = 0$ and at $x_* = 1$.
   Answer

   At $x_* = 0$: mean $1/1.25 = 0.8$, variance $1 - 1/1.25 = 0.2$. At $x_* = 1$: $k = e^{-1/2} = 0.6065$, mean $0.6065/1.25 = 0.4852$, variance $1 - e^{-1}/1.25 = 0.7057$.
2. Why compute $\Sigma_* = K_{**} - V^\top V$ with $V = L^{-1} K_*$, and not $K_{**} - K_*^\top (K^{-1} K_*)$?
   Answer

   It reuses the one factor $L$ (one triangular solve per column, no inverse). Also, $V^\top V$ is symmetric positive semidefinite by construction, so $\Sigma_*$ stays symmetric and its diagonal is at most the prior variance. Forming $K^{-1}$ costs more and is less accurate.
3. Derive $\partial K/\partial \log\ell$ for the RBF kernel. Why is its diagonal zero?
   Answer

   $\partial k/\partial\ell = k\, r^2/\ell^3$, so $\partial k/\partial\log\ell = \ell\,\partial k/\partial\ell = k\, r^2/\ell^2$: the matrix $K_f \circ R^2/\ell^2$. On the diagonal $r = 0$: the prior variance $k(x, x) = \sigma_f^2$ does not depend on $\ell$, and neither does the noise term.
4. Show from the kernel form that the CRPS of a point forecast at $\mu$ is $|y - \mu|$.
   Answer

   For a point mass at $\mu$, $X = X' = \mu$, so $\mathbb E|X - X'| = 0$ and $\mathbb E|X - y| = |\mu - y|$.
5. On 400 independent held-out points, your 95 % intervals cover 0.91. Is this consistent with calibration?
   Answer

   The standard error is $\sqrt{0.95 \times 0.05/400} = 0.0109$, and $0.95 - 0.91 = 0.04$ is about 3.7 standard errors. Not consistent: the intervals are too narrow (overconfident). If the test points are correlated, the effective $n$ is smaller and the evidence is weaker.
6. Why is the squared error of the forecast mean not enough to judge a probabilistic forecast?
   Answer

   Its expected value $(\mu_F - \mu_G)^2 + \operatorname{Var}_G$ depends only on the mean of $F$. Two forecasts with the same mean and very different spreads get the same expected score, so it is proper but not strictly proper. The CRPS and the log score also reward the right spread.

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 GP posterior mean and covariance of $f_*$.

   Show answer

2. What is the predictive variance of a new measurement $y_*$, compared with that of $f_*$?

   Show answer

3. How do you get $\log\det K$ from the Cholesky factor $L$?

   Show answer

4. Write the GP log marginal likelihood and name its terms.

   Show answer

5. What is the gradient of the log marginal likelihood in a hyperparameter $\theta_j$?

   Show answer

6. What does the length scale $\ell$ of an RBF kernel control?

   Show answer

7. When is a scoring rule proper?

   Show answer

8. Give the CRPS of a Gaussian forecast $\mathcal N(\mu, \sigma^2)$.

   Show answer

9. What is the PIT value of an observation, and how is it distributed for a calibrated forecast?

   Show answer

10. When does the central interval at level $p$ contain the observation, in terms of the PIT value $u$?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What does an exact GP cost, and how do you scale it?" — $O(n^3)$ time and $O(n^2)$ memory for the Cholesky factor. For large $n$: inducing points (sparse variational GPs), structured kernels, or iterative solvers (conjugate gradients) on GPUs.
• "Why maximise the marginal likelihood and not the training fit?" — The complexity penalty $-\tfrac12\log\det K$ guards against overfitting (Occam's razor). It can still fail: local optima, very few points, a wrong kernel, or many hyperparameters with little data.
• "Is the negative log-likelihood a proper score? Why use the CRPS?" — Both are strictly proper. The CRPS needs no density (it works on ensembles), has the units of $y$, reduces to the absolute error for a point forecast, and is less sensitive to a single outlier.
• "Your 95 % intervals cover 95 %. Is your model calibrated?" — Not proven. Check other levels or the PIT histogram, check by subgroup (lead time, region, input range), and check sharpness with a proper score.
• "GP or neural surrogate?" — A GP is data-efficient and its variance grows away from the data. A neural surrogate scales to large $n$ and to field outputs, but its uncertainty must be checked and often recalibrated.

10. Practice #

The notebook for this lesson runs in Colab and in the browser. It has five exercises, each with a check cell and a folded solution:

1. Write the GP posterior with one Cholesky factor, and check it against scikit-learn.
2. Compute the log marginal likelihood stably for $n = 400$ points, where det underflows.
3. Write the gradient in the log hyperparameters, and check it against finite differences.
4. Write the Gaussian CRPS in closed form, and check it against the integral of the squared gap and against the kernel form with samples.
5. Measure coverage at 50, 80 and 95 % on held-out data, for intervals of $f_*$ and of $y_*$.

11. Watch and read #

Machine Learning Lecture 26 "Gaussian Processes" -Cornell CS4780 SP17Kilian Weinberger · Cornell CS4780 · 53 min

A clear lecture on GP regression as Bayesian linear regression with kernels, with the posterior formulas of Theorem 1.

Introduction to Gaussian processesGaussian Process Summer School · 100 min

A longer introduction to kernels, prior samples, the posterior and the marginal likelihood.

Forecast verificationICTP Earth System Physics · Andrew Robertson · 68 min

How probabilistic weather and climate forecasts are verified: reliability, sharpness and skill scores.

• Rasmussen and Williams, Gaussian Processes for Machine Learning (free PDF): §2.2 (the function-space view, eqs. 2.22–2.26; Algorithm 2.1), §4.2.1 (the Matérn class), §5.4.1 (the marginal likelihood, its three terms and two local optima), App. A.4 (Cholesky).
• Gneiting and Raftery, Strictly Proper Scoring Rules, Prediction, and Estimation, JASA 2007: §2.1 (proper and strictly proper), §4.2 (the CRPS and its Gaussian closed form).
• Gneiting, Balabdaoui and Raftery, Probabilistic Forecasts, Calibration and Sharpness, JRSS-B 2007: §2 (modes of calibration, the climatological forecaster, the sharpness principle), §3.1 (PIT histograms).
• The scikit-learn user guide on Gaussian processes and GPyTorch for large-scale GPs.