Add crps for the uniform distribution (#50)

* add crps for uniform distribution

* fix bug in crps uniform for torch and tensorflow

docs/api/crps.md

@@ -45,6 +45,8 @@ When the true forecast CDF is not fully known, but represented by a finite ensem
 
 ::: scoringrules.crps_normal
 
+::: scoringrules.crps_uniform
+
 ## Ensemble-based estimators
 
 ::: scoringrules.crps_ensemble

scoringrules/__init__.py

@@ -20,6 +20,7 @@
     crps_loglogistic,
     crps_lognormal,
     crps_normal,
+    crps_uniform,
     crps_quantile,
     owcrps_ensemble,
     twcrps_ensemble,
@@ -67,6 +68,7 @@
     "crps_lognormal",
     "crps_normal",
     "crps_quantile",
+    "crps_uniform",
     "owcrps_ensemble",
     "twcrps_ensemble",
     "vrcrps_ensemble",

scoringrules/_crps.py

@@ -1228,6 +1228,53 @@ def crps_loglogistic(
     return crps.loglogistic(observation, mulog, sigmalog, backend=backend)
 
 
+def crps_uniform(
+    observation: "ArrayLike",
+    min: "ArrayLike",
+    max: "ArrayLike",
+    /,
+    lmass: "ArrayLike" = 0.0,
+    umass: "ArrayLike" = 0.0,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the uniform distribution.
+    
+    It is based on the following formulation from
+    [Jordan et al. (2019)](https://www.jstatsoft.org/article/view/v090i12):
+
+    $$ \mathrm{CRPS}(\mathcal{U}_{L}^{U}(l, u), y) = (u - l) \left\{ | \frac{y - l}{u - l} - F \left( \frac{y - l}{u - l} \right) | + F \left( \frac{y - l}{u - l} \right)^{2} (1 - L - U) - F \left( \frac{y - l}{u - l} \right) (1 - 2L) + \frac{(1 - L - U)^{2}}{3} + (1 - L)U \right\},$$
+
+    where $\mathcal{U}_{L}^{U}(l, u)$ is the uniform distribution with lower bound $l$, 
+    upper bound $u > l$, point mass $L$ on the lower bound, and point mass $U$ on the upper bound.
+    We must have that $L, U \ge 0, L + U < 1$. 
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    min: ArrayLike
+        Lower bound of the forecast uniform distribution.
+    max: ArrayLike
+        Upper bound of the forecast uniform distribution.
+    lmass: ArrayLike
+        Point mass on the lower bound of the forecast uniform distribution.
+    umass: ArrayLike
+        Point mass on the upper bound of the forecast uniform distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between U(min, max, lmass, umass) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.crps_uniform(0.4, 0.0, 1.0, 0.0, 0.0)
+    """
+    return crps.uniform(observation, min, max, lmass, umass, backend=backend)
+
+
 __all__ = [
     "crps_ensemble",
     "twcrps_ensemble",
@@ -1248,4 +1295,5 @@ def crps_loglogistic(
     "crps_logistic",
     "crps_lognormal",
     "crps_normal",
+    "crps_uniform",
 ]

scoringrules/core/crps/__init__.py

@@ -15,6 +15,7 @@
     loglogistic,
     lognormal,
     normal,
+    uniform,
 )
 from ._gufuncs import estimator_gufuncs, quantile_pinball_gufunc
 
@@ -37,6 +38,7 @@
     "logistic",
     "lognormal",
     "normal",
+    "uniform",
     "estimator_gufuncs",
     "quantile_pinball",
     "quantile_pinball_gufunc",

scoringrules/core/crps/_closed.py

@@ -458,6 +458,18 @@ def loglogistic(
     return s
 
 
+def uniform(
+    obs: "ArrayLike", min: "ArrayLike", max: "ArrayLike", lmass: "ArrayLike", umass: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the CRPS for the uniform distribution."""
+    B = backends.active if backend is None else backends[backend]
+    min, max, lmass, umass, obs = map(B.asarray, (min, max, lmass, umass, obs))
+    ω = (obs - min) / (max - min)
+    F_ω = B.minimum(B.maximum(ω, B.asarray(0)), B.asarray(1))
+    s = B.abs(ω - F_ω) + (F_ω**2) * (1 - lmass - umass)  - F_ω * (1 - 2 * lmass) + ((1 - lmass - umass)**2) / 3 + (1 - lmass) * umass
+    return (max - min) * s
+
+
 def _is_scalar_value(x, value):
     if x.size != 1:
         return False

tests/test_crps.py

@@ -66,50 +66,6 @@ def test_quantile_pinball(backend):
         return
 
 
-@pytest.mark.parametrize("backend", BACKENDS)
-def test_normal(backend):
-    obs = np.random.randn(N)
-    mu = obs + np.random.randn(N) * 0.1
-    sigma = abs(np.random.randn(N)) * 0.3
-
-    # non-negative values
-    res = _crps.crps_normal(obs, mu, sigma, backend=backend)
-    res = np.asarray(res)
-    assert not np.any(np.isnan(res))
-    assert not np.any(res < 0.0)
-
-    # approx zero when perfect forecast
-    mu = obs + np.random.randn(N) * 1e-6
-    sigma = abs(np.random.randn(N)) * 1e-6
-    res = _crps.crps_normal(obs, mu, sigma, backend=backend)
-    res = np.asarray(res)
-
-    assert not np.any(np.isnan(res))
-    assert not np.any(res - 0.0 > 0.0001)
-
-
-@pytest.mark.parametrize("backend", BACKENDS)
-def test_lognormal(backend):
-    obs = np.exp(np.random.randn(N))
-    mulog = np.log(obs) + np.random.randn(N) * 0.1
-    sigmalog = abs(np.random.randn(N)) * 0.3
-
-    # non-negative values
-    res = _crps.crps_lognormal(obs, mulog, sigmalog, backend=backend)
-    res = np.asarray(res)
-    assert not np.any(np.isnan(res))
-    assert not np.any(res < 0.0)
-
-    # approx zero when perfect forecast
-    mulog = np.log(obs) + np.random.randn(N) * 1e-6
-    sigmalog = abs(np.random.randn(N)) * 1e-6
-    res = _crps.crps_lognormal(obs, mulog, sigmalog, backend=backend)
-    res = np.asarray(res)
-
-    assert not np.any(np.isnan(res))
-    assert not np.any(res - 0.0 > 0.0001)
-
-
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_beta(backend):
     if backend == "torch":
@@ -357,3 +313,65 @@ def test_loglogistic(backend):
     assert np.isclose(
         _crps.crps_loglogistic(3.0, 0.1, 0.9, backend=backend), 1.13295277, atol=1e-4
     )
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_lognormal(backend):
+    obs = np.exp(np.random.randn(N))
+    mulog = np.log(obs) + np.random.randn(N) * 0.1
+    sigmalog = abs(np.random.randn(N)) * 0.3
+
+    # non-negative values
+    res = _crps.crps_lognormal(obs, mulog, sigmalog, backend=backend)
+    res = np.asarray(res)
+    assert not np.any(np.isnan(res))
+    assert not np.any(res < 0.0)
+
+    # approx zero when perfect forecast
+    mulog = np.log(obs) + np.random.randn(N) * 1e-6
+    sigmalog = abs(np.random.randn(N)) * 1e-6
+    res = _crps.crps_lognormal(obs, mulog, sigmalog, backend=backend)
+    res = np.asarray(res)
+
+    assert not np.any(np.isnan(res))
+    assert not np.any(res - 0.0 > 0.0001)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_normal(backend):
+    obs = np.random.randn(N)
+    mu = obs + np.random.randn(N) * 0.1
+    sigma = abs(np.random.randn(N)) * 0.3
+
+    # non-negative values
+    res = _crps.crps_normal(obs, mu, sigma, backend=backend)
+    res = np.asarray(res)
+    assert not np.any(np.isnan(res))
+    assert not np.any(res < 0.0)
+
+    # approx zero when perfect forecast
+    mu = obs + np.random.randn(N) * 1e-6
+    sigma = abs(np.random.randn(N)) * 1e-6
+    res = _crps.crps_normal(obs, mu, sigma, backend=backend)
+    res = np.asarray(res)
+
+    assert not np.any(np.isnan(res))
+    assert not np.any(res - 0.0 > 0.0001)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_uniform(backend):
+    obs, min, max, lmass, umass = 0.3, -1.0, 2.1, 0.3, 0.1
+    res = _crps.crps_uniform(obs, min, max, lmass, umass, backend=backend)
+    expected = 0.3960968
+    assert np.isclose(res, expected)
+
+    obs, min, max, lmass = -17.9, -15.2, -8.7, 0.2
+    res = _crps.crps_uniform(obs, min, max, lmass, backend=backend)
+    expected = 4.086667
+    assert np.isclose(res, expected)
+
+    obs, min, max = 2.2, 0.1, 3.1
+    res = _crps.crps_uniform(obs, min, max, backend=backend)
+    expected = 0.37
+    assert np.isclose(res, expected)