1 Introduction
In many real-world decision problems — from reinsurance loss models to weather predictions — what matters is not a single point forecast but a full probabilistic forecast that quantifies the uncertainty around future outcomes. A probabilistic forecast issues a predictive distribution \(F\) for a future quantity \(Y\), and the natural question becomes: how do we evaluate such forecasts?
Point forecasts are commonly assessed by the mean absolute error (MAE) or root mean squared error (RMSE). But these metrics are blind to the shape, spread, and calibration of a predictive distribution. A climate model that consistently underestimates extreme rainfall events may have a low RMSE but is clearly unfit for flood risk management. Conversely, a model that issues overly wide prediction intervals may be well-calibrated but uninformative.
Scoring rules provide a framework to address these problems. A scoring rule \(S(F, y)\) assigns a numerical penalty to a forecast distribution \(F\) when the outcome \(y\) materializes. The key requirement for a proper scoring rule is that it incentivizes honesty: the forecaster’s expected score is optimized when the issued forecast \(F\) coincides with the true data-generating distribution \(G\). This property is fundamental for principled forecast evaluation and comparison (Gneiting and Raftery 2007; Winkler et al. 1996).
In this post we develop the theory of proper scoring rules from the ground up, survey the most important examples — Brier Score, Logarithmic Score, CRPS, Energy Score, and Variogram Score — and discuss essential diagnostic tools: PIT histograms for calibration assessment and the Diebold-Mariano test for formal model comparison.
2 Proper Scoring Rules
2.1 Setup and Core Definitions
Let \((\Omega, \mathcal{F})\) be a measurable space, where \(\Omega\) is the sample space and \(\mathcal{F}\) is a \(\sigma\)-algebra of events. Let \(\mathcal{P}\) be a class of probability distributions on this space. We write \(F,G \in \mathcal{P}\) for distributions; when \(\Omega \subseteq \mathbb{R}\)^d, d, \(F(z)\) also denotes the CDF of \(F\), and \(f\) denotes its density when one exists.
A scoring rule is a function \[ S: \mathcal{P} \times \Omega \to \bar{\mathbb{R}}, \quad \bar{\mathbb{R}} := \mathbb{R} \cup \{-\infty, +\infty\}, \] where \(S(F, y)\) quantifies the quality of the predictive distribution \(F\) when the outcome \(y\) is observed.
By convention we adopt the negatively oriented framework: lower scores indicate better forecasts. The expected score under the true distribution \(G\) is \[ S(F, G) := \int_{\Omega} S(F,\omega) \text{d}G(\omega) = \mathbb{E}_{Y \sim G}[S(F, Y)]. \]
A scoring rule \(S\) is proper relative to \(\mathcal{P}\) if \[ S(G, G) \leq S(F, G), \qquad \forall\, F,G \in \mathcal{P}, \] i.e. \(G \in \arg\min_{F \in \mathcal{P}} S(F,G)\) when \(G\) is the data-generating distribution.
It is strictly proper if equality holds only when \(F = G\), i.e. when the minimizer is unique.
The divergence \(d(F, G)\) measures the penalty for misspecification: the excess expected loss incurred by reporting \(F\) when the truth is \(G\). For the Logarithmic Score, this divergence is the Kullback–Leibler divergence; for the CRPS, it is related to the integrated squared difference between CDFs.
2.2 Univariate and Categorical Scoring Rules
The following scores cover binary, categorical, and continuous univariate forecasts.
For a binary event with forecast probability \(p \in [0,1]\) and outcome \(y \in \{0, 1\}\), the Brier Score is \[ \operatorname{BS}(p, y) = (p - y)^2. \] More generally, for \(K\) mutually exclusive categories with forecast probabilities \(\mathbf{p} = (p_1, \ldots, p_K)\) and one-hot outcome \(\mathbf{y}\): \[ \operatorname{BS}(\mathbf{p}, \mathbf{y}) = \sum_{k=1}^{K} (p_k - y_k)^2. \]
The Brier Score is strictly proper and can be decomposed into reliability, resolution, and uncertainty components.
For a distribution \(F\) with predictive density \(f\) and outcome \(y\), the Logarithmic Score (negative log-likelihood) is \[ \operatorname{LogS}(F, y) = -\log f(y). \]
The Logarithmic Score is strictly proper and, up to equivalence, the only proper score that is local: it depends on the forecast only through \(f(y)\) (Gneiting and Raftery 2007). It strongly penalizes forecasts that assign very low density to the realized outcome, making it sensitive to both tail behavior and numerical density estimates.
The Continuous Ranked Probability Score (CRPS) generalizes absolute error to continuous distributional forecasts.
For a predictive CDF \(F\) and observation \(y \in \mathbb{R}\), the Continuous Ranked Probability Score is \[ \operatorname{CRPS}(F, y) = \int_{-\infty}^{\infty} \bigl(F(z) - \mathbf{1}\{y \leq z\}\bigr)^2 \, dz. \tag{1}\] Equivalently, using the representation in terms of expectations: \[ \operatorname{CRPS}(F, y) = \mathbb{E}_F|X - y| - \tfrac{1}{2}\mathbb{E}_F|X - X'|, \tag{2}\] where \(X, X' \stackrel{\text{iid}}{\sim} F\).
The CRPS is strictly proper for distributions with a finite first moment. Its divergence is \[ d(F, G) = \int_{-\infty}^{\infty} \bigl(F(z) - G(z)\bigr)^2 \, dz, \] the integrated squared CDF distance (Gneiting and Raftery 2007). It has the same unit as the observation and reduces to absolute error when \(F\) is a point mass: \(\operatorname{CRPS}(\delta_{\hat y},y)=|\hat y-y|\).
For an ensemble \(\{x^{(1)},\ldots,x^{(m)}\}\), let \(F_m\) denote its empirical distribution. Then \[ \widehat{\operatorname{CRPS}}(F_m,y) = \frac{1}{m}\sum_{k=1}^{m}|x^{(k)}-y| - \frac{1}{2m^2}\sum_{k=1}^{m}\sum_{\ell=1}^{m}|x^{(k)}-x^{(\ell)}|. \]
The fair version replaces the second term by \(\{2m(m-1)\}^{-1}\sum_{k\ne\ell}|x^{(k)}-x^{(\ell)}|\). Closed forms are also available for many parametric families (Jordan, Krüger, and Lerch 2019).
The threshold-weighted CRPS emphasizes selected outcome regions, such as unusually high electricity prices.
For a non-negative weight function \(w: \mathbb{R} \to [0, \infty)\), the threshold-weighted CRPS is \[ \operatorname{twCRPS}(F, y) = \int_{-\infty}^{\infty} \bigl(F(z) - \mathbf{1}\{y \leq z\}\bigr)^2 \, w(z) \, dz. \] This is a proper scoring rule for any non-negative weight function \(w\).
A convenient sample representation uses a non-decreasing chaining function \(v\) satisfying \[ v(z)-v(z')=\int_{z'}^z w(t)\,dt \] (equivalently, \(v'(z)=w(z)\) wherever the derivative exists). For the empirical distribution \(F_m\), \[ \widehat{\operatorname{twCRPS}}_w(F_m,y) =\frac{1}{m}\sum_{k=1}^{m}\left|v\!\left(x^{(k)}\right)-v(y)\right| -\frac{1}{2m^2}\sum_{k=1}^{m}\sum_{\ell=1}^{m} \left|v\!\left(x^{(k)}\right)-v\!\left(x^{(\ell)}\right)\right|. \] Thus, the usual ensemble CRPS formula applies after transforming both forecasts and the observation by \(v\). The fair version again replaces \(2m^2\) by \(2m(m-1)\) and sums over \(k\ne\ell\).
Common choices are:
| Focus | Weight \(w(z)\) | Chaining function \(v(z)\) |
|---|---|---|
| No weighting | \(1\) | \(z\) |
| Upper tail above \(c\) | \(\mathbf{1}\{z\geq c\}\) | \((z-c)_+\) |
| Lower tail below \(c\) | \(\mathbf{1}\{z\leq c\}\) | \(\min(z,c)\) |
| Interval \([a,b]\) | \(\mathbf{1}\{a\leq z\leq b\}\) | \(\min\{\max(z,a),b\}\) |
| Smooth upper tail | \(\{1+\exp[-(z-c)/s]\}^{-1}\) | \(s\log\{1+\exp[(z-c)/s]\}\) |
Here \((u)_+=\max(u,0)\) and \(s>0\) controls the smoothness of the transition. Adding a constant to \(v\) does not change the score.
2.3 Multivariate Scoring Rules
For a joint forecast \(F\) on \(\mathbb{R}^d\) and observation \(\mathbf{y}=(y_1,\ldots,y_d)\), the Logarithmic, Energy, and Variogram Scores capture different aspects of forecast quality.
If \(F\) has joint density \(f\), its multivariate Logarithmic Score is \[ \operatorname{LogS}(F,\mathbf{y})=-\log f(\mathbf{y}). \]
It is strictly proper and directly assesses the full joint density, including dependence, but requires that density to be available and can be dominated by a single near-zero value at the observation.
The Energy Score is a multivariate generalization of the CRPS. For a predictive distribution \(F\) on \(\mathbb{R}^d\) and observation \(\mathbf{y} \in \mathbb{R}^d\): \[ \operatorname{ES}(F, \mathbf{y}) = \mathbb{E}_F \|\mathbf{X} - \mathbf{y}\|^\beta - \frac{1}{2} \mathbb{E}_F \|\mathbf{X} - \mathbf{X}'\|^\beta, \] with \(0 < \beta < 2\) (typically \(\beta = 1\)), where \(\mathbf{X}, \mathbf{X}' \stackrel{\text{iid}}{\sim} F\).
The Energy Score is strictly proper for \(0<\beta<2\) (assuming the required moment exists); for \(d=1\) and \(\beta=1\), it equals the CRPS. Given an ensemble \(\{\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(m)}\}\), let \(F_m\) denote its empirical distribution. Then \[ \widehat{\operatorname{ES}}_\beta(F_m,\mathbf{y}) =\frac{1}{m}\sum_{k=1}^{m}\|\mathbf{x}^{(k)}-\mathbf{y}\|^\beta -\frac{1}{2m^2}\sum_{k=1}^{m}\sum_{\ell=1}^{m}\|\mathbf{x}^{(k)}-\mathbf{x}^{(\ell)}\|^\beta. \]
Its fair version replaces \(2m^2\) by \(2m(m-1)\) and restricts the double sum to \(k\ne\ell\). Computing all pairwise distances costs \(\mathcal{O}(m^2d)\). The score can also be relatively insensitive to misspecified dependence, which motivates the Variogram Score.
The Variogram Score of order \(p > 0\) for a predictive distribution \(F\) on \(\mathbb{R}^d\) and observation \(\mathbf{y} = (y_1, \ldots, y_d)\) is \[ \operatorname{VS}_p(F, \mathbf{y}) = \sum_{i=1}^{d}\sum_{j=1}^{d} w_{ij}\Bigl(|y_i - y_j|^p - \mathbb{E}_F|X_i - X_j|^p\Bigr)^2, \] where \(w_{ij} \geq 0\) are non-negative weights (typically \(w_{ij} = 1\)) and \(p = 0.5\) or \(p = 1\) are common choices.
For \(p\in(0,2)\) and non-negative weights, the Variogram Score is proper, though generally not strictly proper. It is sensitive to dependence through component-wise differences. Its ensemble estimate is \[ \widehat{\operatorname{VS}}_p(F_m,\mathbf{y}) =\sum_{i=1}^{d}\sum_{j=1}^{d}w_{ij}\left(|y_i-y_j|^p -\frac{1}{m}\sum_{k=1}^{m}|x_i^{(k)}-x_j^{(k)}|^p\right)^2. \]
In practice, the Energy and Variogram Scores are complementary: the former captures overall distributional accuracy, while the latter more directly probes dependence.