Details on the metrics implemented in scoringutils

The crps is a proper scoring rule that generalises the absolute error to probabilistic forecasts. It measures the ‘distance’ of the predictive distribution to the observed data-generating distribution. The CRPS is given as \[\text{CRPS}(F, y) = \int_{-\infty}^\infty \left( F(x) - 1(x \geq y) \right)^2 dx,\] where y is the true observed value and F the CDF of predictive distribution. Often An alternative representation is used: \[ \text{CRPS}(F, y) = \frac{1}{2} \mathbb{E}_{F} |X - X'| - \mathbb{E}_P |X - y|,\] where \(X\) and \(X'\) are independent realisations from the predictive distributions \(F\) with finite first moment and \(y\) is the true value. In this represenation we can simply replace \(X\) and \(X'\) by samples sum over all possible combinations to obtain the CRPS. For integer-valued forecasts, the RPS is given as \[ \text{RPS}(F, y) = \sum_{x = 0}^\infty (F(x) - 1(x \geq y))^2. \]

The Log score is a proper scoring rule that is simply compuated as the log of the predictive density evaluated at the true observed value. It is given as \[ \text{log score} = \log f(y), \] where \(f\) is the predictive density function and y is the true value. For integer-valued forecasts, the log score can be computed as \[ \text{log score} = \log p_y, \] where \(p_y\) is the probability assigned to outcome p by the forecast F.

The (weighted) interval score is a proper scoring rule for quantile forecasts that converges to the crps for an increasing number of intervals. The score can be decomposed into a sharpness (uncertainty) component and penalties for over- and underprediction. For a single interval, the score is computed as \[IS_\alpha(F,y) = (u-l) + \frac{2}{\alpha} \cdot (l-y) \cdot 1(y \leq l) + \frac{2}{\alpha} \cdot (y-u) \cdot 1(y \geq u), \] where \(1()\) is the indicator function, \(y\) is the true value, and \(l\) and \(u\) are the \(\frac{\alpha}{2}\) and \(1 - \frac{\alpha}{2}\) quantiles of the predictive distribution \(F\), i.e. the lower and upper bound of a single prediction interval. For a set of \(K\) prediction intervals and the median \(m\), the score is computed as a weighted sum, \[WIS = \frac{1}{K + 0.5} \cdot (w_0 \cdot |y - m| + \sum_{k = 1}^{K} w_k \cdot IS_{\alpha}(F, y)),\] where \(w_k\) is a weight for every interval. Usually, \(w_k = \frac{\alpha_k}{2}\) and \(w_0 = 0.5\).

The Dawid-Sebastiani-Score is a proper scoring rule proposed that only relies on the first moments of the predictive distribution and is therefore easy to compute. It is given as

\[\text{dss}(F, y) = \left( \frac{y - \mu}{\sigma} \right)^2 + 2 \cdot \log \sigma,\] where \(F\) is the predictive distribution with mean \(\mu\) and standard deviation \(\sigma\) and \(y\) is the true observed value.

Proper scoring rule for binary forecasts. The Brier score is computed as \[\text{Brier Score} = \frac{1}{N} \sum_{n = 1}^{N} (f_n - y_n),\] where \(f_n\), with \(n = 1, \dots, N\) are the predicted probablities that the corresponding events, \(y_n \in (0, 1)\) will be equal to one.)

Interval coverage measures the proportion of observed values that fall in a given prediction interval range. Interval coverage for a single prediction interval range can be calculated as \[IC_{\alpha} = \text{nominal coverage} - \text{empirical coverage},\] where nominal coverage is \(1 - \alpha\) and empirical coverage is the proportion of true values actually covered by all \(1 - \alpha\) prediction intervals.

To summarise interval coverage over different over multiple interval ranges, we can compute coverage deviation defined as the mean interval coverage over all \(K\) interval ranges \(\alpha_k\) with \(k = 1, \dots, K\): \[\text{Coverage deviation} = \frac{1}{K} \sum_{k = 1}^{K} \text{IC}_{\alpha_k}\]

Quantile coverage for a given quantile level is the proportion of true values smaller than the predictions corresponding to that quantile level.

The probability integral transform (PIT, Dawid 1984) represents a succinct way to visualise deviations between the predictive distribution \(F\) and the true data-generating distribution \(G\). The idea is to transform the observed values such that agreement between forecasts and data can then be examined by observing whether or not the transformed values follow a uniform distribution. The PIT is given by \[u = F (y),\] where \(u\) is the transformed variable and \(F(y)\) is the predictive distribution \(F\) evaluated at the true observed value \(y\). If \(F = G\), then \(u\) follows a uniform distribution.

For integer outcomes, the PIT is no longer uniform even when forecasts are ideal. Instead, a randomised PIT can be used: \[u = P(y) + v \cdot (P(y) - P(y - 1) ),\] where \(y\) is again the observed value \(P()\) is the cumulative probability assigned to all values smaller or equal to \(y\) (where \(P(-1) = 0\) by definition, and \(v\) is a standard uniform variable independent of \(y\). If \(P\) is equal to the true data-generating distribution function, then \(u\) is standard uniform. also propose a non-randomised version of the PIT for count data that could be used alternatively.

Sharpness is the ability to produce narrow forecasts and is a feature of the forecasts only and does not depend on the observations. Sharpness is therefore only of interest conditional on calibration: a very precise forecast is not useful if it is clearly wrong.

As suggested by Funk et al. (2019), we measure sharpness for continuous and integer forecasts represented by predictive samples as the normalised median absolute deviation about the median (MADN) ), i.e. \[ S(F) = \frac{1}{0.675} \cdot \text{median}(|x - \text{median(x)}|), \] where \(x\) is the vector of all predictive samples and \(\frac{1}{0.675}\) is a normalising constant. If the predictive distribution \(F\) is the CDF of a normal distribution, then sharpness will equal the standard deviation of \(F\).

Bias is a measure of the tendency of a forecaster to over- or underpredict. For continuous forecasts, bias is given as \[B(F, y) = 1 - 2 \cdot (F (y)), \] where \(F\) is the CDF of the predictive distribution and \(y\) is the observed value.

For integer-valued forecasts, bias can be calculated as \[B(P, y) = 1 - (P(y) + P(y + 1)), \] where \(P(y)\) is the cumulative probability assigned to all outcomes smaller or equal to \(y\).

For quantile forecasts, Bias can be calculated as the maximum percentile rank for which the prediction is smaller than \(y\), if the true value is smaller than the median of the predictive distribution. If the true value is above the median of the predictive distribution, then bias is the minimum percentile rank for which the corresponding quantile is still larger than the true value. If the true value is exactly the median, bias is zero. For a large enough number of quantiles, the percentile rank will equal the proportion of predictive samples below the observed true value, and this metric coincides with the one for continuous forecasts.

The mean score ratio is used to compare two models on the overlapping set of forecast targets for which both models have made a prediction. The mean score ratio is calculated as the mean score achieved by the first model over the mean score achieved by the second model. More precisely, for two models \(i, j\), we determine the set of overlapping forecasts, denoted by \(\mathcal{A}_{ij}\) and compute the mean score ratio \(\theta_{ij}\) as \[\theta_{ij} =\frac{\text{mean score model } i \text{ on } \mathcal{A}_{ij}}{\text{mean score model } j \text{ on } \mathcal{A}_{ij}}.\] The mean score ratio can in principle be computed for any arbitrary score.

Relative skill scores can be used to obtain a ranking of models based on pairwise comparisons between all models. To compute the relative skill \(\theta_i\) of model \(i\), we take the geometric mean of all mean score ratios that involve model \(i\), i.e. \[ \theta_{i} = \left(\prod_{m = 1}^M \theta_{im}\right)^{1/M}, \] where M is the number of models.