Implement SCRPS

examples/time_series_learning.py

@@ -15,7 +15,7 @@
 import torch
 from matplotlib.figure import Figure
 
-from torch_crps import crps_analytical
+from torch_crps.analytical import crps_analytical
 
 EXAMPLES_DIR = pathlib.Path(pathlib.Path(__file__).parent)
 

pyproject.toml

@@ -40,7 +40,7 @@ classifiers = [
 dependencies = [
   "torch>=2.7",
 ]
-description = "Implementations of the CRPS using PyTorch"
+description = "PyTorch-based implementations of the Continuously-Ranked Probability Score (CRPS) as well as its locally scale-invariant version (SCRPS)"
 dynamic = ["version"]
 license = "CC-BY-4.0"
 maintainers = [{name = "Fabio Muratore", email = "accounts@famura.net"}]
@@ -140,6 +140,7 @@ ignore = [
   "PLC0415", # import should be at the top-level of a file
   "PLR", # pylint refactor
   "RUF001", # allow for greek letters
+  "RUF003", # allow for greek letters
   "UP035", # allow importing Callable from typing
 ]
 preview = true

readme.md

@@ -13,18 +13,58 @@
 [![Ruff][ruff-badge]][ruff]
 [![uv][uv-badge]][uv]
 
-Implementations of the Continuously-Ranked Probability Score (CRPS) using PyTorch
+PyTorch-based implementations of the Continuously-Ranked Probability Score (CRPS) as well as its locally scale-invariant
+version (SCRPS)
 
 ## Background
 
-The Continuously-Ranked Probability Score (CRPS) is a strictly proper scoring rule.
-It assesses how well a distribution with the cumulative distribution function $F$ is explaining an observation $y$
+### Continuously-Ranked Probability Score (CRPS)
 
-$$ \text{CRPS}(F,y) = \int _{\mathbb {R} }(F(x)-\mathbb {1} (x\geq y))^{2}dx \qquad (\text{integral formulation}) $$
+The CRPS is a strictly proper scoring rule.
+It assesses how well a distribution with the cumulative distribution function $F(X)$ of the estimate $X$ (a random
+variable) is explaining an observation $y$
+
+$$
+\text{CRPS}(F,y) = \int _{\mathbb {R}} \left( F(x)-\mathbb {1} (x\geq y) \right)^{2} dx
+$$
 
 where $1$ denoted the indicator function.
 
-In Section 2 of this [paper][crps-folumations] Zamo & Naveau list 3 different formulations of the CRPS.
+In Section 2 of this [paper][crps-folumations] Zamo & Naveau list 3 different formulations of the CRPS. One of them is
+
+$$
+\text{CRPS}(F, y) = E[|X - y|] - 0.5 E[|X - X'|] = E[|X - y|] + E[X] - 2 E[X F(X)]
+$$
+
+which can be shortened to
+
+$$
+\text{CRPS}(F, y) = A - 0.5 D
+$$
+
+where $A$ is called the accuracy term and $D$ is called the disperion term (at least I do it in this repo).
+
+### Scaled Continuously-Ranked Probability Score (SCRPS)
+
+The SCRPS is a locally scale-invariant version of the CRPS.
+In their [paper][scrps-paper], Bolling & Wallin define it in a positively-oriented, i.e., higher is better.
+In contrast, I implement the SCRPS in this repo negatively-oriented, just like a loss function.
+
+Oversimplifying the notation, the (negatively-oriented) SCRPS can be written as
+
+$$
+\text{SCRPS}(F, y) = -\frac{E[|X - y|]}{E[|X - X'|]} - 0.5 \log \left( E[|X - X'|] \right)
+$$
+
+which can be shortened to
+
+$$
+\text{SCRPS}(F, y) = \frac{A}{D} + 0.5 \log(D)
+$$
+
+The scale-invariance, i.e., the SCRPS value does not depend on the magnitude of $D$, comes from the division by $D$.
+
+Note that the SCRPS can, in contrast to the CRPS, yield negative values.
 
 ### Incomplete list of sources that I came across while researching about the CRPS
 
@@ -33,15 +73,16 @@ In Section 2 of this [paper][crps-folumations] Zamo & Naveau list 3 different fo
 - Gneiting & Raftery; "Strictly Proper Scoring Rules, Prediction, and Estimation"; 2007
 - Zamo & Naveau; "Estimation of the Continuous Ranked Probability Score with Limited Information and Applications to Ensemble Weather Forecasts"; 2018
 - Jordan et al.; "Evaluating Probabilistic Forecasts with scoringRules"; 2019
+- Bollin & Wallin; "Local scale invariance and robustness of proper scoring rules"; 2029
 - Olivares & Négiar & Ma et al; "CLOVER: Probabilistic Forecasting with Coherent Learning Objective Reparameterization"; 2023
 - Vermorel & Tikhonov; "Continuously-Ranked Probability Score (CRPS)" [blog post][Lokad-post]; 2024
 - Nvidia; "PhysicsNeMo Framework" [source code][nvidia-crps-implementation]; 2025
 - Zheng & Sun; "MVG-CRPS: A Robust Loss Function for Multivariate Probabilistic Forecasting"; 2025
 
 ## Application to Machine Learning
 
-The CRPS can be used as a loss function in machine learning, just like the well-known negative log-likelihood loss which
-is the log scoring rule.
+The CRPS, as well as the SCRPS, can be used as a loss function in machine learning, just like the well-known negative
+log-likelihood loss which is the log scoring rule.
 
 The parametrized model outputs a distribution $q(x)$. The CRPS loss evaluates how good $q(x)$ is explaining the
 observation $y$.
@@ -54,7 +95,15 @@ There is [work on multi-variate CRPS estimation][multivariate-crps], but it is n
 
 ## Implementation
 
-The integral formulation is infeasible to naively evaluate on a computer due to the infinite integration over $x$.
+The direct implementation of the integral formulation is not suited to evaluate on a computer due to the infinite
+integration over the domain of the random variable $X$.
+Nevertheless, this repository includes such an implementation to verify the others.
+
+The normalization-by-observation variants are improper solutions to normalize the CPRS values. The goal is to use the
+CPRS as a loss function in machine learning tasks. For that, it is highly beneficial if the loss does not depend on
+the scale of the problem.
+However, deviding by the absolute maximum of the observations is a bad proxy for doing this.
+I plan on removing these methods once I gained trust in my SCRPS implementation.
 
 I found [Nvidia's implementation][nvidia-crps-implementation] of the CRPS for ensemble preductions in $M log(M)$ time
 inspiring to read.
@@ -89,6 +138,7 @@ inspiring to read.
 [uv]: https://docs.astral.sh/uv
 <!-- Paper URLS-->
 [crps-folumations]: https://link.springer.com/article/10.1007/s11004-017-9709-7
+[scrps-paper]: https://arxiv.org/abs/1912.05642
 [Lokad-post]: https://www.lokad.com/continuous-ranked-probability-score/
 [multivariate-crps]: https://arxiv.org/pdf/2410.09133
 [nvidia-crps-implementation]: https://docs.nvidia.com/physicsnemo/25.11/_modules/physicsnemo/metrics/general/crps.html

tests/conftest.py

@@ -2,6 +2,9 @@
 
 import pytest
 import torch
+from torch.distributions import Normal, StudentT
+
+from torch_crps.analytical.studentt import standardized_studentt_cdf_via_scipy
 
 results_dir = Path(__file__).parent / "results"
 results_dir.mkdir(parents=True, exist_ok=True)
@@ -49,3 +52,97 @@ def case_batched_3d():
         "y": torch.randn(2, 3) * 10 + 50,
         "expected_shape": torch.Size([2, 3]),
     }
+
+
+def crps_analytical_normal_gneiting(
+    q: Normal,
+    y: torch.Tensor,
+) -> torch.Tensor:
+    """Compute the analytical CRPS assuming a normal distribution.
+
+    See Also:
+        Gneiting & Raftery; "Strictly Proper Scoring Rules, Prediction, and Estimation"; 2007
+        Equation (5) for the analytical formula for CRPS of Normal distribution.
+
+    Args:
+        q: A PyTorch Normal distribution object, typically a model's output distribution.
+        y: Observed values, of shape (num_samples,).
+
+    Returns:
+        CRPS values for each observation, of shape (num_samples,).
+    """
+    # Compute standard normal CDF and PDF.
+    z = (y - q.loc) / q.scale
+    standard_normal = torch.distributions.Normal(0, 1)
+    phi_z = standard_normal.cdf(z)  # Φ(z)
+    pdf_z = torch.exp(standard_normal.log_prob(z))  # φ(z)
+
+    # Analytical CRPS formula.
+    crps = q.scale * (z * (2 * phi_z - 1) + 2 * pdf_z - 1 / torch.sqrt(torch.tensor(torch.pi)))
+
+    return crps
+
+
+def crps_analytical_studentt_jordan(
+    q: StudentT,
+    y: torch.Tensor,
+) -> torch.Tensor:
+    r"""Compute the (negatively-oriented) CRPS in closed-form assuming a StudentT distribution.
+
+    This is the previous implementation of the analytical CRPS for StudentT distributions.. It is provided here for
+    testing and comparison purposes.
+
+    See Also:
+        Jordan et al.; "Evaluating Probabilistic Forecasts with scoringRules"; 2019.
+
+    Args:
+        q: A PyTorch StudentT distribution object, typically a model's output distribution.
+        y: Observed values, of shape (num_samples,).
+
+    Returns:
+        CRPS values for each observation, of shape (num_samples,).
+    """
+    # Extract degrees of freedom ν, location μ, and scale σ.
+    nu, mu, sigma = q.df, q.loc, q.scale
+    if torch.any(nu <= 1):
+        raise ValueError("StudentT CRPS requires degrees of freedom > 1")
+
+    # Standardize, and create standard StudentT distribution for CDF and PDF.
+    z = (y - mu) / sigma
+    standard_t = torch.distributions.StudentT(nu, loc=0, scale=1)
+
+    # Compute standardized CDF F_ν(z) and PDF f_ν(z).
+    cdf_z = standardized_studentt_cdf_via_scipy(z, nu)
+    pdf_z = torch.exp(standard_t.log_prob(z))
+
+    # Compute the beta function ratio: B(1/2, ν - 1/2) / B(1/2, ν/2)^2
+    # Using the relationship: B(a,b) = Gamma(a) * Gamma(b) / Gamma(a+b)
+    # B(1/2, ν - 1/2) / B(1/2, ν/2)^2 = ( Gamma(1/2) * Gamma(ν-1/2) / Gamma(ν) ) /
+    #                                   ( Gamma(1/2) * Gamma(ν/2) / Gamma(ν/2 + 1/2) )^2
+    # Simplifying to Gamma(ν - 1/2) Gamma(ν/2 + 1/2)^2 / ( Gamma(ν)Gamma(ν/2)^2 )
+    # For numerical stability, we compute in log space.
+    log_gamma_half = torch.lgamma(torch.tensor(0.5, dtype=nu.dtype, device=nu.device))
+    log_gamma_df_minus_half = torch.lgamma(nu - 0.5)
+    log_gamma_df_half = torch.lgamma(nu / 2)
+    log_gamma_df_half_plus_half = torch.lgamma(nu / 2 + 0.5)
+
+    # log[B(1/2, ν-1/2)] = log Gamma(1/2) + log Gamma(ν-1/2) - log Gamma(ν)
+    # log[B(1/2, ν/2)] = log Gamma(1/2) + log Gamma(ν/2) - log Gamma(ν/2 + 1/2)
+    # log[B(1/2, ν-1/2) / B(1/2, ν/2)^2] = log B(1/2, ν-1/2) - 2*log B(1/2, ν/2)
+    log_beta_ratio = (
+        log_gamma_half
+        + log_gamma_df_minus_half
+        - torch.lgamma(nu)
+        - 2 * (log_gamma_half + log_gamma_df_half - log_gamma_df_half_plus_half)
+    )
+    beta_frac = torch.exp(log_beta_ratio)
+
+    # Compute the CRPS for standardized values.
+    crps_standard = (
+        z * (2 * cdf_z - 1) + 2 * pdf_z * (nu + z**2) / (nu - 1) - (2 * torch.sqrt(nu) / (nu - 1)) * beta_frac
+    )
+
+    # Apply location-scale transformation CRPS(F_{ν,μ,σ}, y) = σ * CRPS(F_{ν}, z) with z = (y - μ) / σ.
+    crps = sigma * crps_standard
+
+    return crps