Add crps_csg0 for 0-censored shifted gamma distribution

docs/reference.md

@@ -58,6 +58,7 @@ Parametric distributions forecasts
     crps_exponentialM
     crps_2pexponential
     crps_gamma
+    crps_csg0
     crps_gev
     crps_gpd
     crps_gtclogistic

scoringrules/__init__.py

@@ -16,6 +16,7 @@
     crps_exponentialM,
     crps_2pexponential,
     crps_gamma,
+    crps_csg0,
     crps_gev,
     crps_gpd,
     crps_gtclogistic,
@@ -139,6 +140,7 @@ def wrapper(*args, **kwargs):
     "crps_exponentialM",
     "crps_2pexponential",
     "crps_gamma",
+    "crps_csg0",
     "crps_gev",
     "crps_gpd",
     "crps_gtclogistic",

scoringrules/_crps.py

@@ -846,6 +846,83 @@ def crps_gamma(
     return crps.gamma(obs, shape, rate, backend=backend)
 
 
+def crps_csg0(
+    obs: "ArrayLike",
+    /,
+    shape: "ArrayLike",
+    rate: "ArrayLike | None" = None,
+    *,
+    scale: "ArrayLike | None" = None,
+    shift: "ArrayLike" = 0.0,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the censored, shifted gamma distribution.
+
+    It is based on the following formulation from [1]_:
+
+    .. math::
+        \mathrm{CRPS}\bigl(F^0_{\alpha,\beta,\delta},\,y\bigr) =
+        (y+\delta)\,\bigl(2F_{\alpha,\beta}(y+\delta)-1\bigr)
+        - \frac{\alpha}{\beta\pi}\,B\!\Bigl(\tfrac12,\alpha+\tfrac12\Bigr)\,\bigl(1 - F_{2\alpha,\beta}(2\delta)\bigr)
+        + \frac{\alpha}{\beta}\,\Bigl(1 + 2\,F_{\alpha,\beta}(\delta)\,F_{\alpha+1,\beta}(\delta) - F_{\alpha,\beta}(\delta)^2 - 2\,F_{\alpha+1,\beta}(y+\delta)\Bigr)
+        - \delta\,\bigl(F_{\alpha,\beta}(\delta)\bigr)^{2},
+
+    where :math:`F^0_{\alpha,\beta,\delta}` is the censored, shifted gamma distribution function with
+    shape parameter :math:`\alpha > 0`, rate parameter :math:`\beta > 0` (equivalently, with scale parameter :math:`1/\beta`)
+    and shift parameter :math:`\delta > 0`.
+
+    Parameters
+    ----------
+    obs : array_like
+        The observed values.
+    shape : array_like
+        Shape parameter of the forecast CSG distribution.
+    rate : array_like, optional
+        Rate parameter of the forecast CSG distribution.
+        Either ``rate`` or ``scale`` must be provided.
+    scale : array_like, optional
+        Scale parameter of the forecast CSG distribution, where ``scale = 1 / rate``.
+        Either ``rate`` or ``scale`` must be provided.
+    shift : array_like
+        Shift parameter of the forecast CSG distribution.
+    backend : str, optional
+        The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``.
+
+    Returns
+    -------
+    crps : array_like
+        The CRPS between obs and CSG(shape, rate, shift).
+
+    References
+    ----------
+    .. [1] Scheuerer, M., & Hamill, T. M. (2015).
+        Statistical postprocessing of ensemble precipitation forecasts
+        by fitting censored, shifted gamma distributions.
+        Monthly Weather Review, 143(11), 4578-4596.
+        https://doi.org/10.1175/MWR-D-15-0061.1
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.crps_csg0(0.7, shape=0.5, rate=2.0, shift=0.3)
+    0.5411044348806484
+
+    Raises
+    ------
+    ValueError
+        If both ``rate`` and ``scale`` are provided, or if neither is provided.
+    """
+    if (scale is None and rate is None) or (scale is not None and rate is not None):
+        raise ValueError(
+            "Either ``rate`` or ``scale`` must be provided, but not both or neither."
+        )
+
+    if rate is None:
+        rate = 1.0 / scale
+
+    return crps.csg0(obs, shape, rate, shift, backend=backend)
+
+
 def crps_gev(
     obs: "ArrayLike",
     shape: "ArrayLike",
@@ -2282,6 +2359,7 @@ def crps_uniform(
     "crps_exponentialM",
     "crps_2pexponential",
     "crps_gamma",
+    "crps_csg0",
     "crps_gev",
     "crps_gpd",
     "crps_gtclogistic",

scoringrules/core/crps/__init__.py

@@ -6,6 +6,7 @@
     exponentialM,
     twopexponential,
     gamma,
+    csg0,
     gev,
     gpd,
     gtclogistic,
@@ -41,6 +42,7 @@
     "exponentialM",
     "twopexponential",
     "gamma",
+    "csg0",
     "gev",
     "gpd",
     "gtclogistic",

scoringrules/core/crps/_closed.py

@@ -196,6 +196,31 @@ def gamma(
     return s
 
 
+def csg0(
+    obs: "ArrayLike",
+    shape: "ArrayLike",
+    rate: "ArrayLike",
+    shift: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the CRPS for the censored, shifted gamma distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, shape, rate, shift = map(B.asarray, (obs, shape, rate, shift))
+    obs_shifted = obs + shift
+    F_ab_shifted = _gamma_cdf(obs_shifted, shape, rate, backend=backend)
+    F_2ab_2d = _gamma_cdf(2 * shift, 2 * shape, rate, backend=backend)
+    F_ab_d = _gamma_cdf(shift, shape, rate, backend=backend)
+    F_ab1_d = _gamma_cdf(shift, shape + 1, rate, backend=backend)
+    F_ab1_shifted = _gamma_cdf(obs_shifted, shape + 1, rate, backend=backend)
+    s = (
+        obs_shifted * (2 * F_ab_shifted - 1)
+        - (shape / (rate * B.pi)) * B.beta(B.asarray(0.5), shape + 0.5) * (1 - F_2ab_2d)
+        + shape / rate * (1 + 2 * F_ab_d * F_ab1_d - F_ab_d**2 - 2 * F_ab1_shifted)
+        - shift * F_ab_d**2
+    )
+    return s
+
+
 EULERMASCHERONI = 0.57721566490153286060651209008240243
 
 

tests/test_crps.py

@@ -243,6 +243,32 @@ def test_crps_gamma(backend):
         return
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_crps_csg0(backend):
+    obs, shape, rate, shift = 0.7, 0.5, 2.0, 0.3
+    expected = 0.5411044
+
+    expected_gamma = sr.crps_gamma(obs, shape, rate, backend=backend)
+    res_gamma = sr.crps_csg0(obs, shape=shape, rate=rate, shift=0.0, backend=backend)
+    assert np.isclose(res_gamma, expected_gamma)
+
+    res = sr.crps_csg0(obs, shape=shape, rate=rate, shift=shift, backend=backend)
+    assert np.isclose(res, expected)
+
+    res = sr.crps_csg0(obs, shape=shape, scale=1.0 / rate, shift=shift, backend=backend)
+    assert np.isclose(res, expected)
+
+    with pytest.raises(ValueError):
+        sr.crps_csg0(
+            obs, shape=shape, rate=rate, scale=1.0 / rate, shift=shift, backend=backend
+        )
+        return
+
+    with pytest.raises(ValueError):
+        sr.crps_csg0(obs, shape=shape, shift=shift, backend=backend)
+        return
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_crps_gev(backend):
     if backend == "torch":