Add closed-form CRPS for exponential distribution

docs/api/crps.md

@@ -13,10 +13,11 @@ When the true forecast CDF is not fully known, but represented by a finite ensem
 
 <h2>Analytical formulations</h2>
 
-::: scoringrules.crps_normal
+::: scoringrules.crps_exponential
 
 ::: scoringrules.crps_lognormal
 
+::: scoringrules.crps_normal
 
 <h2>Ensemble-based estimators</h2>
 

scoringrules/__init__.py

@@ -3,6 +3,7 @@
 from scoringrules._brier import brier_score
 from scoringrules._crps import (
     crps_ensemble,
+    crps_exponential,
     crps_logistic,
     crps_lognormal,
     crps_normal,
@@ -35,6 +36,7 @@
     "backends",
     "crps_ensemble",
     "crps_normal",
+    "crps_exponential",
     "crps_lognormal",
     "crps_logistic",
     "crps_quantile",

scoringrules/_crps.py

@@ -335,6 +335,48 @@ def vrcrps_ensemble(
     )
 
 
+def crps_exponential(
+    observation: "ArrayLike",
+    rate: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the exponential distribution.
+
+    It is based on the following formulation from
+    [Jordan et al. (2019)](https://www.jstatsoft.org/article/view/v090i12):
+
+    $$\mathrm{CRPS}(F_{\lambda}, y) = |y| - \frac{2F_{\lambda}(y)}{\lambda} + \frac{1}{2 \lambda},$$
+
+    where $F_{\lambda}$ is exponential distribution function with rate parameter $\lambda > 0$.
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    rate:
+        Rate parameter of the forecast exponential distribution.
+
+    Returns
+    -------
+    score:
+        The CRPS between Exp(rate) and obs.
+
+    Examples
+    --------
+    ```pycon
+    >>> import scoringrules as sr
+    >>> import numpy as np
+    >>> sr.crps_exponential(0.8, 3.0)
+    0.360478635526275
+    >>> sr.crps_exponential(np.array([0.8, 0.9]), np.array([3.0, 2.0]))
+    array([0.36047864, 0.24071795])
+    ```
+    """
+    return crps.exponential(observation, rate, backend=backend)
+
+
 def crps_normal(
     observation: "ArrayLike",
     mu: "ArrayLike",
@@ -398,11 +440,11 @@ def crps_lognormal(
 
     Parameters
     ----------
-    observations: ArrayLike
+    observation:
         The observed values.
-    mulog: ArrayLike
+    mulog:
         Mean of the normal underlying distribution.
-    sigmalog: ArrayLike
+    sigmalog:
         Standard deviation of the underlying normal distribution.
 
     Returns

scoringrules/backend/base.py

@@ -39,6 +39,12 @@ def mean(
     ) -> "Array":
         """Calculate the arithmetic mean of the input array ``x``."""
 
+    @abc.abstractmethod
+    def max(
+        self, x: "Array", axis: int | tuple[int, ...] | None, keepdims: bool = False
+    ) -> "Array":
+        """Return the maximum value of an input array ``x``."""
+
     @abc.abstractmethod
     def moveaxis(
         self,

scoringrules/backend/jax.py

@@ -44,6 +44,11 @@ def mean(
     ) -> "Array":
         return jnp.mean(x, axis=axis, keepdims=keepdims)
 
+    def max(
+        self, x: "Array", axis: int | tuple[int, ...] | None, keepdims: bool = False
+    ) -> "Array":
+        return jnp.max(x, axis=axis, keepdims=keepdims)
+
     def moveaxis(
         self,
         x: "Array",

scoringrules/backend/numpy.py

@@ -43,6 +43,11 @@ def mean(
     ) -> "NDArray":
         return np.mean(x, axis=axis, keepdims=keepdims)
 
+    def max(
+        self, x: "NDArray", axis: int | tuple[int, ...] | None, keepdims: bool = False
+    ) -> "NDArray":
+        return np.max(x, axis=axis, keepdims=keepdims)
+
     def moveaxis(
         self,
         x: "NDArray",

scoringrules/backend/tensorflow.py

@@ -49,6 +49,11 @@ def mean(
     ) -> "Tensor":
         return tf.math.reduce_mean(x, axis=axis, keepdims=keepdims)
 
+    def max(
+        self, x: "Tensor", axis: int | tuple[int, ...] | None, keepdims: bool = False
+    ) -> "Tensor":
+        return tf.math.reduce_max(x, axis=axis, keepdims=keepdims)
+
     def moveaxis(
         self,
         x: "Tensor",

scoringrules/backend/torch.py

@@ -44,6 +44,11 @@ def mean(
     ) -> "Tensor":
         return torch.mean(x, axis=axis, keepdim=keepdims)
 
+    def max(
+        self, x: "Tensor", axis: int | tuple[int, ...] | None, keepdims: bool = False
+    ) -> "Tensor":
+        return torch.max(x, axis=axis, keepdim=keepdims)
+
     def moveaxis(
         self,
         x: "Tensor",

scoringrules/core/crps/__init__.py

@@ -1,14 +1,15 @@
 from ._approx import ensemble, ow_ensemble, quantile_pinball, vr_ensemble
-from ._closed import logistic, lognormal, normal
+from ._closed import exponential, logistic, lognormal, normal
 from ._gufuncs import estimator_gufuncs, quantile_pinball_gufunc
 
 __all__ = [
     "ensemble",
     "ow_ensemble",
     "vr_ensemble",
-    "normal",
-    "lognormal",
+    "exponential",
     "logistic",
+    "lognormal",
+    "normal",
     "estimator_gufuncs",
     "quantile_pinball",
     "quantile_pinball_gufunc",

scoringrules/core/crps/_closed.py

@@ -1,12 +1,22 @@
 import typing as tp
 
 from scoringrules.backend import backends
-from scoringrules.core.stats import _logis_cdf, _norm_cdf, _norm_pdf
+from scoringrules.core.stats import _exp_cdf, _logis_cdf, _norm_cdf, _norm_pdf
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, ArrayLike, Backend
 
 
+def exponential(
+    obs: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the CRPS for the exponential distribution."""
+    B = backends.active if backend is None else backends[backend]
+    rate, obs = map(B.asarray, (rate, obs))
+    s = B.abs(obs) - (2 * _exp_cdf(obs, rate, backend=backend) / rate) + 1 / (2 * rate)
+    return s
+
+
 def normal(
     obs: "ArrayLike", mu: "ArrayLike", sigma: "ArrayLike", backend: "Backend" = None
 ) -> "Array":