Add CRPS for Student's t distribution

docs/api/crps.md

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

scoringrules/__init__.py

@@ -76,6 +76,7 @@
     "crps_normal",
     "crps_poisson",
     "crps_quantile",
+    "crps_t",
     "crps_uniform",
     "owcrps_ensemble",
     "twcrps_ensemble",

scoringrules/_crps.py

@@ -1384,6 +1384,51 @@ def crps_poisson(
     return crps.poisson(observation, mean, backend=backend)
 
 
+def crps_t(
+    observation: "ArrayLike",
+    df: "ArrayLike",
+    /,
+    location: "ArrayLike" = 0.0,
+    scale: "ArrayLike" = 1.0,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the student's t distribution.
+
+    It is based on the following formulation from
+    [Jordan et al. (2019)](https://www.jstatsoft.org/article/view/v090i12):
+
+    $$ \mathrm{CRPS}(F, y) = \sigma \left\{ \omega (2 F_{\nu} (\omega) - 1) + 2 f_{\nu} \left( \frac{\nu + \omega^{2}}{\nu - 1} \right) - \frac{2 \sqrt{\nu}}{\nu - 1} \frac{B(\frac{1}{2}, \nu - \frac{1}{2})}{B(\frac{1}{2}, \frac{\nu}{2}^{2})}  \right},$$
+
+    where $\omega = (y - \mu)/\sigma$, where $\nu > 1, \mu$, and $\sigma > 0$ are the 
+    degrees of freedom, location, and scale parameters respectively of the Student's t 
+    distribution, and $f_{\nu}$ and $F_{\nu}$ are the PDF and CDF of the standard Student's
+    t distribution with $\nu$ degrees of freedom.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    df: ArrayLike
+        Degrees of freedom parameter of the forecast t distribution.
+    location: ArrayLike
+        Location parameter of the forecast t distribution.
+    sigma: ArrayLike
+        Scale parameter of the forecast t distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between t(df, location, scale) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.crps_t(0.0, 0.1, 0.4, 0.1)
+    """
+    return crps.t(observation, df, location, scale, backend=backend)
+
+
 def crps_uniform(
     observation: "ArrayLike",
     min: "ArrayLike",

scoringrules/backend/jax.py

@@ -208,7 +208,8 @@ def factorial(self, n: "ArrayLike") -> "ArrayLike":
         return jsp.special.factorial(n)
 
     def hypergeometric(self, a: "Array", b: "Array", c: "Array", z: "Array"):
-        return jsp.special.hyp2f1(a, b, c, z)
+        #return jsp.special.hyp2f1(a, b, c, z)
+        raise NotImplementedError 
 
     def comb(self, n: "ArrayLike", k: "ArrayLike") -> "ArrayLike":
         return jsp.special.factorial(n) // (

scoringrules/backend/torch.py

@@ -225,7 +225,7 @@ def factorial(self, n: "TensorLike") -> "TensorLike":
     def hypergeometric(
         self, a: "Tensor", b: "Tensor", c: "Tensor", z: "Tensor"
     ) -> "Tensor":
-        return None
+        return NotImplementedError
 
     def comb(self, n: "Tensor", k: "Tensor") -> "Tensor":
         return self.factorial(n) // (self.factorial(k) * self.factorial(n - k))

scoringrules/core/crps/__init__.py

@@ -17,6 +17,7 @@
     lognormal,
     normal,
     poisson,
+    t,
     uniform,
 )
 from ._gufuncs import estimator_gufuncs, quantile_pinball_gufunc
@@ -42,6 +43,7 @@
     "lognormal",
     "normal",
     "poisson",
+    "t",
     "uniform",
     "estimator_gufuncs",
     "quantile_pinball",

scoringrules/core/crps/_closed.py

@@ -15,6 +15,8 @@
     _norm_pdf,
     _pois_cdf,
     _pois_pdf,
+    _t_cdf,
+    _t_pdf,
 )
 
 if tp.TYPE_CHECKING:
@@ -517,6 +519,22 @@ def poisson(
     return s
 
 
+def t(
+    obs: "ArrayLike", df: "ArrayLike", location: "ArrayLike", scale: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the CRPS for the t distribution."""
+    B = backends.active if backend is None else backends[backend]
+    df, mu, sigma, obs = map(B.asarray, (df, location, scale, obs))
+    z = (obs - mu) / sigma
+    F_z = _t_cdf(z, df, backend=backend)
+    f_z = _t_pdf(z, df, backend=backend)
+    G_z = (df + z**2) / (df - 1)    
+    s1 = z * (2 * F_z - 1)
+    s2 = 2 * f_z * G_z
+    s3 = (2 * B.sqrt(df) / (df - 1)) * B.beta(1 / 2, df - 1 / 2) / (B.beta(1 / 2, df / 2)**2)
+    return sigma * (s1 + s2 - s3)
+
+
 def uniform(
     obs: "ArrayLike", min: "ArrayLike", max: "ArrayLike", lmass: "ArrayLike", umass: "ArrayLike", backend: "Backend" = None
 ) -> "Array":

scoringrules/core/stats.py

@@ -64,7 +64,7 @@ def _t_pdf(x: "ArrayLike", df: "ArrayLike", backend: "Backend" = None) -> "Array
 def _t_cdf(x: "ArrayLike", df: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Cumulative distribution function for the standard Student's t distribution."""
     B = backends.active if backend is None else backends[backend]
-    return 1 / 2 + x * B.hypergeo(1 / 2, (df + 1) / 2, 3 / 2, -(x**2) / df) / (
+    return 1 / 2 + x * B.hypergeometric(1 / 2, (df + 1) / 2, 3 / 2, -(x**2) / df) / (
         B.sqrt(df) * B.beta(1 / 2, df / 2)
     )
 

tests/test_crps.py

@@ -433,6 +433,22 @@ def test_poisson(backend):
     assert np.isclose(res, expected)
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_t(backend):
+    if backend in ["jax", "torch", "tensorflow"]:
+        pytest.skip("Not implemented in torch or tensorflow backends")
+
+    obs, df, mu, sigma = 11.1, 5.2, 13.8, 2.3
+    expected = 1.658226
+    res = _crps.crps_t(obs, df, mu, sigma, backend=backend)
+    assert np.isclose(res, expected)
+
+    obs, df = 0.7, 4.0
+    expected = 0.4387929
+    res = _crps.crps_t(obs, df, backend=backend)
+    assert np.isclose(res, expected)
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_uniform(backend):
     obs, min, max, lmass, umass = 0.3, -1.0, 2.1, 0.3, 0.1