Scoring a Forecast Distribution
A distributional regression model predicts a conditional law, not just a single response value. For a covariate value \(x\), let \(F_x\) be the forecast distribution of the response and let \(y\) be the outcome that is later observed. A scoring rule assigns a numerical loss to that forecast and outcome; lower scores are better. A useful score should reward a forecast distribution that matches the conditional law of the response.
Two such scores are the continuous ranked probability score (CRPS), designed for real-valued outcomes, and the energy score, which extends distance-based scoring to vector outcomes. Their distance representations make both their calculations and their propriety properties especially clear.
The integral definition compares the forecast CDF with the step function associated with the observed value. The pairwise form has a different interpretation: it rewards forecast values that are close to the observation, while subtracting half the forecast’s expected internal pairwise distance. The subtraction matters. Using only \(E[|X-y|]\) would generally favor an overly concentrated forecast.
Why the CRPS Has a Pairwise-Distance Form
For any real numbers \(a\) and \(b\), the distance \(|a-b|\) is the length of the set of thresholds \(z\) at which the indicators \(\mathbf{1}\{a\le z\}\) and \(\mathbf{1}\{b\le z\}\) differ. Thus $$ |a-b|=\int_{-\infty}^{\infty}\left|\mathbf{1}\{a\le z\}-\mathbf{1}\{b\le z\}\right|\,dz. $$ Applying this identity to \(X,y\), and to two independent draws \(X,X'\), gives an integral expression for each expectation in the pairwise formula.
At a fixed threshold \(z\), write \(A=\mathbf{1}\{X\le z\}\), \(A'=\mathbf{1}\{X'\le z\}\), and \(c=\mathbf{1}\{y\le z\}\). Since \(A,A'\), and \(c\) are binary, $$ E[|A-c|]=|F(z)-c| \quad\text{and}\quad \frac{1}{2}E[|A-A'|]=F(z)(1-F(z)). $$ If \(c=0\), their difference is \(F(z)-F(z)(1-F(z))=F(z)^2\). If \(c=1\), it is \(1-F(z)-F(z)(1-F(z))=(1-F(z))^2\). In either case the difference is \((F(z)-c)^2\). Integrating over \(z\) proves the pairwise representation. The finite first moment assumption ensures the distance expectations are finite, so the representation is well-defined.
Worked Example: A Bernoulli Forecast
Suppose \(F\) assigns probability \(p\) to 1 and \(1-p\) to 0. Calculate the CRPS for each possible observation, then compare the forecast with a true Bernoulli distribution having success probability \(q\).
If \(y=1\), then \(E[|X-y|]=1-p\). Two independent draws differ with probability \(2p(1-p)\), and their distance when they differ is 1. Therefore \(E[|X-X'|]=2p(1-p)\), and $$ \operatorname{CRPS}(F,1)=(1-p)-\frac{1}{2}\,2p(1-p)=(1-p)^2. $$
If \(y=0\), then \(E[|X-y|]=p\), while the pairwise term is unchanged. Hence $$ \operatorname{CRPS}(F,0)=p-p(1-p)=p^2. $$
Under the true Bernoulli law, the expected score of the forecast is \(q(1-p)^2+(1-q)p^2\). The true forecast, with \(p=q\), has expected score \(q(1-q)\). Subtracting gives $$ q(1-p)^2+(1-q)p^2-q(1-q)=(p-q)^2. $$ The expected-score penalty is nonnegative and vanishes exactly when \(p=q\).
Propriety of the CRPS
A scoring rule is proper if the expected score under the true distribution is minimized by forecasting that distribution. It is strictly proper if no different forecast achieves the same minimum. For the CRPS, the precise regret has a useful CDF expression.
To prove the result, fix \(z\). The indicator \(\mathbf{1}\{Y\le z\}\) is Bernoulli with mean \(G(z)\). Expanding its squared error around that mean gives $$ E\left[(F(z)-\mathbf{1}\{Y\le z\})^2\right] =(F(z)-G(z))^2+G(z)(1-G(z)). $$ For forecast \(G\), the corresponding expected squared error at \(z\) is \(G(z)(1-G(z))\). Subtracting the two expressions gives \((F(z)-G(z))^2\).
Integrating this identity proves the regret formula. The interchange of expectation and integration is valid for the nonnegative squared-error integrands by Tonelli’s theorem. The expected CRPS values are finite under the stated first-moment assumptions: the pairwise representation bounds them by sums of expected absolute distances. Thus the difference of the two expected scores is well-defined. The right side is nonnegative, proving propriety.
If the integral is zero, then \(F(z)=G(z)\) for almost every \(z\) with respect to Lebesgue measure. CDFs are right-continuous, so this equality almost everywhere implies equality everywhere: for any \(z\), choose a sequence \(z_k>z\) approaching \(z\) from the full-measure set where the CDFs agree, and take right-hand limits. Therefore \(F=G\), proving strict propriety. \(\square\)
Worked Example: CRPS for a Standard Normal Forecast at Zero
Let \(F=N(0,1)\) and observe \(y=0\). For a standard normal random variable, direct integration of its density gives \(E[|X|]=\sqrt{2/\pi}\). Since \(X-X'\sim N(0,2)\), the same absolute-normal calculation gives \(E[|X-X'|]=\sqrt{2}\sqrt{2/\pi}=2/\sqrt{\pi}\). Therefore $$ \operatorname{CRPS}(F,0)=\sqrt{\frac{2}{\pi}}-\frac{1}{2}\frac{2}{\sqrt{\pi}} =\frac{\sqrt{2}-1}{\sqrt{\pi}}\approx 0.234. $$
The score is positive even though the observation is exactly at the forecast mean: it assesses the whole forecast distribution, not just how far the observation is from a point prediction.
The Multivariate Energy Score
For a response in \(\mathbb{R}^d\), the energy score uses Euclidean distance. Let \(X,X'\) be independent draws from the forecast distribution \(F\), and let \(y\in\mathbb{R}^d\). The finite first moment condition is \(E[\|X\|]<\infty\), where \(\|\cdot\|\) is the Euclidean norm.
The energy score is not obtained by putting a multivariate CDF into the one-dimensional CRPS formula. It is a distance-based multivariate score. Its propriety follows from the fact that Euclidean distance distinguishes probability distributions through a nonnegative energy distance.
The equality follows directly from the score definition: the forecast’s expected distance to \(Y\) is \(E[\|X-Y\|]\), and the truth’s expected score is \(\tfrac{1}{2}E[\|Y-Y'\|]\). For nonnegativity and strictness, define the expression in parentheses as \(D_E(F,G)\). A standard integral representation of the Euclidean norm gives, for some constant \(c_d>0\), $$ \|u\|=c_d\int_{\mathbb{R}^d}\frac{1-\cos(t^\top u)}{\|t\|^{d+1}}\,dt. $$ The integral is finite: near zero, \(1-\cos(t^\top u)\le \|t\|^2\|u\|^2/2\), and at infinity the numerator is bounded while the radial tail is proportional to \(\int_1^\infty r^{-2}\,dr\). Rotation and scaling show that the integral is a positive constant times \(\|u\|\), which establishes the representation.
Substitute this representation into the three terms defining \(D_E\). At each \(t\), the constant terms cancel, and the cosine terms combine to give $$ D_E(F,G)=c_d\int_{\mathbb{R}^d} \frac{|\varphi_F(t)-\varphi_G(t)|^2}{\|t\|^{d+1}}\,dt, $$ where \(\varphi_F(t)=E[e^{it^\top X}]\) and \(\varphi_G(t)=E[e^{it^\top Y}]\) are the characteristic functions. The finite first moments ensure the original distance expectations are finite; the norm representation can be integrated against each distance expectation by Tonelli’s theorem. The displayed integrand is nonnegative, so \(D_E(F,G)\ge 0\). If \(D_E(F,G)=0\), the characteristic functions agree almost everywhere away from zero. They are continuous, so they agree everywhere. Uniqueness of characteristic functions then gives \(F=G\). Conversely, if \(F=G\), the expression is zero. This proves the theorem. \(\square
Worked Example: A Multivariate Forecast Sensitive to Dependence
Consider two distributions on \(\mathbb{R}^2\). Let \(G\) assign probability \(1/2\) to each of \((1,1)\) and \((-1,-1)\), and let \(F\) assign probability \(1/2\) to each of \((1,-1)\) and \((-1,1)\). Each coordinate has the same marginal distribution under \(F\) and \(G\), but the dependence differs.
Within either distribution, two independent draws differ with probability \(1/2\); when they differ, their distance is \(2\sqrt{2}\). Thus \(E[\|X-X'\|]=E[\|Y-Y'\|]=\sqrt{2}\). Every point in the support of \(F\) is distance \(2\) from each point in the support of \(G\), so \(E[\|X-Y\|]=2\). The energy-distance regret of forecasting \(F\) when \(G\) is true is therefore $$ \frac{1}{2}\left(2(2)-\sqrt{2}-\sqrt{2}\right)=2-\sqrt{2}\approx 0.586. $$
Thus the multivariate energy score can distinguish these joint distributions even though their coordinatewise marginal distributions are identical. It assesses the vector forecast as a whole.
Using the Scores in Distributional Regression
At a fixed covariate value \(x\), apply the definitions to the forecast conditional distribution \(F_x\) and the observed response \(y\). If the true conditional law is \(G_x\) and both laws have finite first moments, the propriety results apply conditionally at that \(x\): the expected score is minimized by forecasting \(G_x\). For scalar responses this is the CRPS; for vector responses the energy score uses Euclidean distance.
This distinction is useful when choosing a score for a distributional forecast. CRPS directly compares univariate CDFs and has an exact integrated squared-CDF regret. The energy score provides a natural multivariate distance-based rule and, under its stated moment condition, strictly rewards the true joint distribution.
Common Pitfall
Wrong reasoning: “The energy score uses Euclidean distance, so if two forecasts have the same coordinate marginals, they must receive the same score and cannot differ in dependence.” The worked example shows why this conclusion fails: cross-distances between vector draws depend on how the coordinates occur together, not only on their separate marginal laws. A particular observation can happen to give equal scores for two forecasts, but propriety concerns expected scores under the true joint distribution. Compare the expected scores, not just one realized score or a collection of marginal scores.
Check Your Understanding
Use the definitions and regret results to answer the following questions.
- In the CRPS pairwise representation, what role does the subtracted half-distance term play?
- If the CRPS expected regret is zero, why does equality of the CDFs almost everywhere imply equality of the distributions?
- What finite-moment condition is used for the energy score on \(\mathbb{R}^d\)?
- Why is the energy score for a vector outcome not simply a one-dimensional CDF formula applied coordinate by coordinate?
- For a conditional forecast at a fixed covariate value, which distribution minimizes expected CRPS or energy score under the stated moment conditions?