From e081b4d95277699b89a4ec885a356ebb4212edb8 Mon Sep 17 00:00:00 2001 From: tobiasbiegert Date: Fri, 25 Apr 2025 17:06:08 +0200 Subject: [PATCH 1/2] Add crps_csg0 for 0-censored shifted Gamma; update API, docs, and tests --- docs/api/crps.md | 2 + scoringrules/__init__.py | 2 + scoringrules/_crps.py | 70 ++++++++++++++++++++++++++++++ scoringrules/core/crps/__init__.py | 2 + scoringrules/core/crps/_closed.py | 25 +++++++++++ tests/test_crps.py | 26 +++++++++++ 6 files changed, 127 insertions(+) diff --git a/docs/api/crps.md b/docs/api/crps.md index 9e8895c..5bda4d2 100644 --- a/docs/api/crps.md +++ b/docs/api/crps.md @@ -25,6 +25,8 @@ also be viewed as the Brier score integrated over all real-valued thresholds. ::: scoringrules.crps_gamma +::: scoringrules.crps_csg0 + ::: scoringrules.crps_gev ::: scoringrules.crps_gpd diff --git a/scoringrules/__init__.py b/scoringrules/__init__.py index 4e2725b..47cf351 100644 --- a/scoringrules/__init__.py +++ b/scoringrules/__init__.py @@ -14,6 +14,7 @@ crps_exponentialM, crps_2pexponential, crps_gamma, + crps_csg0, crps_gev, crps_gpd, crps_gtclogistic, @@ -113,6 +114,7 @@ "crps_exponentialM", "crps_2pexponential", "crps_gamma", + "crps_csg0", "crps_gev", "crps_gpd", "crps_gtclogistic", diff --git a/scoringrules/_crps.py b/scoringrules/_crps.py index c02ed61..70f1d39 100644 --- a/scoringrules/_crps.py +++ b/scoringrules/_crps.py @@ -664,6 +664,75 @@ def crps_gamma( return crps.gamma(observation, shape, rate, backend=backend) +def crps_csg0( + observation: "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 Appendix B of + [Scheuerer and Hamill (2015)](https://doi.org/10.1175/MWR-D-15-0061.1): + + $$ + \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 $F^0_{\alpha,\beta,\delta}$ is the censored, shifted gamma distribution function with + shape parameter $\alpha > 0$, rate parameter $\beta > 0$ (equivalently, with scale parameter $1/\beta$) + and shift parameter $\delta > 0$. + + Parameters + ---------- + observation: + The observed values. + shape: + Shape parameter of the forecast CSG distribution. + rate: + Rate parameter of the forecast CSG distribution. + scale: + Scale parameter of the forecast CSG distribution, where `scale = 1 / rate`. + shift: + Shift parameter of the forecast CSG distribution. + backend: + The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + + Returns + ------- + score: + The CRPS between obs and CSG(shape, rate, shift). + + 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(observation, shape, rate, shift, backend=backend) + + def crps_gev( observation: "ArrayLike", shape: "ArrayLike", @@ -1957,6 +2026,7 @@ def crps_uniform( "crps_exponentialM", "crps_2pexponential", "crps_gamma", + "crps_csg0", "crps_gev", "crps_gpd", "crps_gtclogistic", diff --git a/scoringrules/core/crps/__init__.py b/scoringrules/core/crps/__init__.py index 74f563a..c54b179 100644 --- a/scoringrules/core/crps/__init__.py +++ b/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", diff --git a/scoringrules/core/crps/_closed.py b/scoringrules/core/crps/_closed.py index 7edc13e..024f740 100644 --- a/scoringrules/core/crps/_closed.py +++ b/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 diff --git a/tests/test_crps.py b/tests/test_crps.py index 5fb7447..ff49400 100644 --- a/tests/test_crps.py +++ b/tests/test_crps.py @@ -229,6 +229,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": From d7ebeea58c99aad8f9a56047625f0760c1fc84d3 Mon Sep 17 00:00:00 2001 From: tobiasbiegert Date: Thu, 2 Oct 2025 16:42:21 +0200 Subject: [PATCH 2/2] Update docstring using RST syntax --- scoringrules/_crps.py | 58 ++++++++++++++++++++++++------------------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/scoringrules/_crps.py b/scoringrules/_crps.py index 35e2e91..62e806e 100644 --- a/scoringrules/_crps.py +++ b/scoringrules/_crps.py @@ -847,7 +847,7 @@ def crps_gamma( def crps_csg0( - observation: "ArrayLike", + obs: "ArrayLike", /, shape: "ArrayLike", rate: "ArrayLike | None" = None, @@ -858,41 +858,49 @@ def crps_csg0( ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the censored, shifted gamma distribution. - It is based on the following formulation from Appendix B of - [Scheuerer and Hamill (2015)](https://doi.org/10.1175/MWR-D-15-0061.1): + It is based on the following formulation from [1]_: - $$ - \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} - $$ + .. 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 $F^0_{\alpha,\beta,\delta}$ is the censored, shifted gamma distribution function with - shape parameter $\alpha > 0$, rate parameter $\beta > 0$ (equivalently, with scale parameter $1/\beta$) - and shift parameter $\delta > 0$. + 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 ---------- - observation: + obs : array_like The observed values. - shape: + shape : array_like Shape parameter of the forecast CSG distribution. - rate: + rate : array_like, optional Rate parameter of the forecast CSG distribution. - scale: - Scale parameter of the forecast CSG distribution, where `scale = 1 / rate`. - shift: + 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: - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + backend : str, optional + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. Returns ------- - score: + 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 @@ -902,17 +910,17 @@ def crps_csg0( Raises ------ ValueError - If both `rate` and `scale` are provided, or if neither is provided. + 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." + "Either ``rate`` or ``scale`` must be provided, but not both or neither." ) if rate is None: rate = 1.0 / scale - return crps.csg0(observation, shape, rate, shift, backend=backend) + return crps.csg0(obs, shape, rate, shift, backend=backend) def crps_gev(