Source-linked AI summary
Evaluating probabilistic forecasts with scoringRules
Alexander Jordan, Fabian Krüger, Sebastian Lerch
TL;DR
The paper addresses how to evaluate and select among competing probabilistic forecasting models. It presents scoringRules, an R package that computes proper scoring rules across parametric and simulated forecast distributions, and illustrates its use in applications. The package supports comparative evaluation while leaving absolute forecast diagnostics outside its scope.
Problem
The proliferation of probabilistic models creates a need for principled comparative evaluation and selection of competing forecasts.
Method
The paper presents the scoringRules R package, which computes proper scoring rules for parametric and simulated forecast distributions and demonstrates their use in applications.
Results
scoringRules provides generally applicable, numerically efficient implementations covering a wide range of probabilistic forecasting situations.
Takeaways & Limitations
The package supports comparative evaluation of probabilistic models across classical and Bayesian forecasting methods.
Takeaways & Limitations
The package focuses on comparative forecast evaluation and does not cover absolute forecast evaluation diagnostics.
Abstract
from arXiv · showhide
Probabilistic forecasts in the form of probability distributions over future events have become popular in several fields including meteorology, hydrology, economics, and demography. In typical applications, many alternative statistical models and data sources can be used to produce probabilistic forecasts. Hence, evaluating and selecting among competing methods is an important task. The scoringRules package for R provides functionality for comparative evaluation of probabilistic models based on proper scoring rules, covering a wide range of situations in applied work. This paper discusses implementation and usage details, presents case studies from meteorology and economics, and points to the relevant background literature.
1. Introduction: Forecast evaluation
The paper addresses the need to evaluate competing probabilistic forecasts and introduces scoringRules as a general R framework for comparative evaluation with proper scoring rules.
- Probabilistic forecasting creates a need for principled tools to evaluate and compare models using diverse theories, data sources, and statistical specifications.
- scoringRules computes scoring rules for parametric distributions and simulated forecast samples, covering both frequentist and Bayesian forecasting methods.
- The package implements proper CRPS, logarithmic, multivariate energy, and variogram scores, with generally applicable formulas and statistically principled defaults.
- scoringRules is designed for comparative forecast evaluation, whereas absolute diagnostics such as checking predictive-distribution calibration are outside its scope.
- The package is positioned as a broadly applicable library, contrasting with existing tools that target narrower distributions, data structures, outcome types, or programming languages.
2. Theoretical background
The paper develops the theoretical basis for proper scoring rules and their computation for analytic and simulated forecast distributions. It emphasizes truthful comparative evaluation while addressing computational and density-estimation choices.
- Probabilistic forecasts are represented as parametric distributions or simulated samples, motivating separate computational approaches for analytic and non-analytic predictive distributions.
- Proper scoring rules rank competing forecasts through average scores, with lower scores preferred and truthful reporting optimized in expectation.
- CRPS and LogS are central proper scoring rules; CRPS uses the predictive CDF, whereas LogS requires a predictive density.
- Closed-form CRPS expressions support exact and efficient computation for many parametric families implemented in scoringRules.
- For simulated forecasts, the empirical CDF approximates the predictive CDF, and CRPS computation can improve from O(m^2) to average O(m log m) using sorted samples.
- Kernel density estimation for LogS can be fragile in forecast tails because the estimated score may depend strongly on bandwidth choices.
3. Package design and functionality
The package organizes scoring-rule computations across parametric and simulated predictive distributions, with parallel generic and family-specific interfaces. It also provides approximation methods for samples and illustrates their behavior through parametric and simulation examples.
- Core interfaces: Scoring functions follow the [score]_[suffix]() convention, including crps_norm() and generic crps() and logs() interfaces.Numeric-class wrappers add stricter checks and meaningful error messages for interactive use.
- Core interfaces: The package supports parametric families through computation functions and numeric-class methods, while sample-based functions are reserved for empirical predictive distributions.Sample functions such as [score]_sample() cannot be called through the numeric class method.
- Parametric predictive distributions: Fixed gamma-distribution examples show that LogS rises rapidly near zero and is minimized at the predictive mode, whereas CRPS is more symmetric and minimized at the predictive median.The example uses shape = 2 and scale = 1.5 through crps_gamma() and logs_gamma().
- Simulated predictive distributions: Sample-based functions accept observation vectors and simulated-sample matrices, with rows matched to observations; a vector is also allowed when y has length one.The package provides both univariate and multivariate sample functionality.
- Simulated predictive distributions: In the Gaussian toy experiment, score approximations become more accurate as the sample size increases, using samples of up to 5 000 and 500 independent replications.Figure 2 compares sampling approximations with analytically calculated CRPS and LogS values.
- Simulated predictive distributions: The package offers empirical-distribution and kernel-density methods for CRPS, while LogS uses kernel density estimation because it requires an estimated density.The default CRPS method is "edf"; the default LogS method is kernel density estimation with a Gaussian kernel.
4. Usage examples
The usage examples show how scoringRules supports probabilistic forecast evaluation and score-based parameter estimation in weather and economics. Post-processing improves precipitation forecasts over the raw ensemble, while alternative estimators produce broadly similar GDP forecast scores and parameter estimates.
- 4.1. Probabilistic weather forecasting via ensemble post-processing: Censored non-homogeneous regression models represent precipitation forecasts with a mixed discrete-continuous distribution, including a positive probability of zero precipitation.The distribution parameters depend on ensemble features, with location modeled from the ensemble mean and scale from the logarithm of ensemble spread.
- 4.1. Probabilistic weather forecasting via ensemble post-processing: All post-processing models substantially improve mean CRPS over the raw ensemble, with logistic and Student’s t models slightly outperforming the Gaussian model.CRPS is computed for the model outputs and for the raw ensemble using crps() and crps_sample().
- 4.2. Bayesian forecasts of US GDP growth rate: Bayesian Markov-switching autoregressive models generate US GDP growth forecast samples whose certainty decreases as the horizon extends from one to four quarters.Histograms show simulated forecast draws, while a black line gives the mixture-of-normals approximation.
- 4.2. Bayesian forecasts of US GDP growth rate: The two GDP forecast estimators have quite similar scores, while the mixture-of-parameters estimator is theoretically preferred for the logarithmic score over the KDE estimator.The similarity is attributed to the large number of 20,000 MCMC draws.
- 4.3. Parameter estimation with scoring rules: Proper scoring rules also support parameter estimation: minimum CRPS and maximum likelihood yield similar deviation distributions, although maximum likelihood has slightly lower variability for σ.The comparison uses 1,000 independent samples of size 500 from a normal distribution with μ = −1 and σ = 2.
5. Multivariate scoring rules
The paper extends proper scoring rules to multivariate forecast distributions, focusing on sample-based forecasts and implementations of the energy and variogram scores. It also illustrates joint evaluation across multiple forecast horizons.
- Multivariate proper scoring rules extend univariate LogS and CRPS to forecast distributions on R^d.The paper discusses multivariate generalizations including the energy score and variogram score.
- The scoringRules package implements the energy score and variogram score for multivariate forecasts represented by simulated samples.These implementations target forecasts given through discrete samples rather than requiring closed-form parametric expressions.
- Variogram score weights w_i,j can emphasize or downweight component pairs, while p controls the score order; typical values are 0.5 and 1.
- Multivariate scoring functions evaluate one forecast case at a time and return a single number, with vector observations and d × m sample matrices as inputs.Multiple cases require sequential application through apply() functions or loops.
- The economic case study jointly evaluates forecast performance over four horizons using a four-variate simulated sample.This contrasts with evaluating each horizon separately.
- Empirical analyses typically average scores across several forecast cases and compare two or more models.
6. Summary and discussion
The paper presents scoringRules as a broad R library for proper scoring-rule computation across parametric and simulated forecasts. It emphasizes extensible implementations while noting that score selection remains non-trivial.
- The scoringRules package computes proper scoring rules for both parametric and simulated forecast distributions across many applied forecasting settings.The implementations are designed to be generally applicable and numerically efficient.
- The package’s main crps() and logs() functions are extensible S3 generics that can support additional forecast-object classes.Users can define S3 methods for classes beyond numeric, including fitted model objects.
- The package includes CRPS, logarithmic score, multivariate energy score, and variogram score, with statistically principled defaults where alternatives exist.
- Choosing an appropriate proper scoring rule for model evaluation or parameter estimation remains a non-trivial task.The paper identifies future extensions as including additional scoring rules and parametric families when suitable expressions are available.
A. Formulas for the CRPS
Appendix A.1 defines notation and special functions used in the paper’s CRPS formulas.
- A.1. Notation: γ denotes the Euler–Mascheroni constant.
- A.1. Notation: ⌊x⌋ denotes the floor function.
- A.1. Notation: sgn(x) denotes the sign function.
- A.1. Notation: Ei(x) denotes the exponential integral, while ϕ(x) and Φ(x) denote the standard Gaussian density and distribution functions.
- A.1. Notation: Γ(a) denotes the gamma function, with Γ_l(a, x) and Γ_u(a, x) denoting lower and upper incomplete gamma functions.
- A.1. Notation: B(a, b) denotes the beta function, and I(a, b, x) denotes the regularized incomplete beta function.
- A.1. Notation: I_m(x) denotes the modified Bessel function of the first kind.
- A.1. Notation: 2F1(a, b; c; x) denotes the hypergeometric function.
B. Computation of multivariate scores for multiple forecast cases
The appendix demonstrates how to evaluate multivariate scores across many forecast cases by generating repeated observations and simulated forecasts, organizing them in alternative data structures, and applying scoring functions sequentially.
- The example generates observations from a ten-dimensional multivariate normal distribution with unit variances and pairwise covariance 0.2.
- The simulation is independently repeated 1,000 times to create multiple forecast cases.
- Observations and forecasts are stored both as list elements and as arrays whose final dimension indexes the forecast case.
- All 1,000 forecast cases can be evaluated sequentially with sapply() over list elements or the last array dimension.
- The resulting ES vectors are computed equivalently from list-based and array-based inputs using es_sample().