Add CRPS for the generalised truncated and censored normal distribution

docs/api/crps.md

@@ -27,6 +27,12 @@ When the true forecast CDF is not fully known, but represented by a finite ensem
 
 ::: scoringrules.crps_gpd
 
+::: scoringrules.crps_gtcnormal
+
+::: scoringrules.crps_tnormal
+
+::: scoringrules.crps_cnormal
+
 ::: scoringrules.crps_hypergeometric
 
 ::: scoringrules.crps_laplace

scoringrules/__init__.py

@@ -10,6 +10,9 @@
     crps_gamma,
     crps_gev,
     crps_gpd,
+    crps_gtcnormal,
+    crps_tnormal,
+    crps_cnormal,
     crps_hypergeometric,
     crps_laplace,
     crps_logistic,
@@ -48,18 +51,21 @@
     "crps_ensemble",
     "crps_beta",
     "crps_binomial",
-    "crps_normal",
     "crps_exponential",
     "crps_exponentialM",
     "crps_gamma",
     "crps_gev",
     "crps_gpd",
+    "crps_gtcnormal",
+    "crps_tnormal",
+    "crps_cnormal",
     "crps_hypergeometric",
     "crps_laplace",
     "crps_loglaplace",
     "crps_loglogistic",
     "crps_logistic",
     "crps_lognormal",
+    "crps_normal",
     "crps_quantile",
     "owcrps_ensemble",
     "twcrps_ensemble",

scoringrules/_crps.py

@@ -1,7 +1,7 @@
 import typing as tp
 
 from scoringrules.backend import backends
-from scoringrules.core import crps
+from scoringrules.core import crps, stats
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, ArrayLike, Backend
@@ -759,6 +759,147 @@ def crps_gpd(
     return crps.gpd(observation, shape, location, scale, mass, backend=backend)
 
 
+def crps_gtcnormal(
+    observation: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    /,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    lmass: "ArrayLike" = 0.0,
+    umass: "ArrayLike" = 0.0,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the generalised truncated and censored normal distribution.
+
+    It is based on the following formulation from
+    [Jordan et al. (2019)](https://www.jstatsoft.org/article/view/v090i12):
+
+    $$ \mathrm{CRPS}(F_{l, L}^{u, U}, y) = |y - z| + uU^{2} - lL^{2} + \left( \frac{1 - L - U}{\Phi(u) - \Phi(l)} \right) z \left( 2 \Phi(z) - \frac{(1 - 2L) \Phi(u) + (1 - 2U) \Phi(l)}{1 - L - U} \right) + \left( \frac{1 - L - U}{\Phi(u) - \Phi(l)} \right) \left( 2 \phi(z) - 2 \phi(u)U - 2 \phi(l)L \right) - \left( \frac{1 - L - U}{\Phi(u) - \Phi(l)} \right)^{2} \left( \frac{1}{\sqrt{\pi}} \right) \left( \Phi(u \sqrt{2}) - \Phi(l \sqrt{2}) \right), $$
+
+    $$ \mathrm{CRPS}(F_{l, L, \mu, \sigma}^{u, U}, y) = \sigma \mathrm{CRPS}(F_{(l - \mu)/\sigma, L}^{(u - \mu)/\sigma, U}, \frac{y - \mu}{\sigma}), $$
+
+    where $\Phi$ and $\phi$ are respectively the CDF and PDF of the standard normal
+    distribution, $F_{l, L, \mu, \sigma}^{u, U}$ is the CDF of the normal distribution 
+    truncated below at $l$ and above at $u$, with point masses $L, U > 0$ at the lower and upper
+    boundaries, respectively, and location and scale parameters $\mu$ and $\sigma > 0$.
+    $F_{l, L}^{u, U} = F_{l, L, 0, 1}^{u, U}$.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+    lmass: ArrayLike
+        Point mass assigned to the lower boundary of the forecast distribution.
+    umass: ArrayLike
+        Point mass assigned to the upper boundary of the forecast distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between gtcNormal(location, scale, lower, upper, lmass, umass) and obs.
+
+    Examples
+    --------
+    >>> from scoringrules import crps
+    >>> crps.gtcnormal(0.0, 0.1, 0.4, -1.0, 1.0, 0.1, 0.1)
+    """
+    return crps.gtcnormal(observation, location, scale, lower, upper, lmass, umass, backend=backend)
+
+
+def crps_tnormal(
+    observation: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    /,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the truncated normal distribution.
+
+    It is based on the formulation for the generalised truncated and censored normal distribution with 
+    lmass and umass set to zero.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between tNormal(location, scale, lower, upper) and obs.
+
+    Examples
+    --------
+    >>> from scoringrules import crps
+    >>> crps.tnormal(0.0, 0.1, 0.4, -1.0, 1.0)
+    """
+    return crps.gtcnormal(observation, location, scale, lower, upper, 0.0, 0.0, backend=backend)
+
+
+def crps_cnormal(
+    observation: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    /,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the censored normal distribution.
+
+    It is based on the formulation for the generalised truncated and censored normal distribution with 
+    lmass and umass set to the tail probabilities of the predictive distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between cNormal(location, scale, lower, upper) and obs.
+
+    Examples
+    --------
+    >>> from scoringrules import crps
+    >>> crps.cnormal(0.0, 0.1, 0.4, -1.0, 1.0)
+    """
+    lmass = stats._norm_cdf((lower - location) / scale)
+    umass = 1 - stats._norm_cdf((upper - location) / scale)
+    return crps.gtcnormal(observation, location, scale, lower, upper, lmass, umass, backend=backend)
+
+
 def crps_hypergeometric(
     observation: "ArrayLike",
     m: "ArrayLike",
@@ -1100,6 +1241,9 @@ def crps_loglogistic(
     "crps_gamma",
     "crps_gev",
     "crps_gpd",
+    "crps_gtcnormal",
+    "crps_tnormal",
+    "crps_cnormal",
     "crps_hypergeometric",
     "crps_logistic",
     "crps_lognormal",

scoringrules/core/crps/__init__.py

@@ -7,6 +7,7 @@
     gamma,
     gev,
     gpd,
+    gtcnormal,
     hypergeometric,
     laplace,
     logistic,
@@ -28,6 +29,7 @@
     "gamma",
     "gev",
     "gpd",
+    "gtcnormal",
     "hypergeometric",
     "laplace",
     "loglaplace",

scoringrules/core/crps/_closed.py

@@ -249,6 +249,42 @@ def gpd(
     return scale * s
 
 
+def gtcnormal(
+    obs: "ArrayLike", location: "ArrayLike", scale: "ArrayLike", lower: "ArrayLike", upper: "ArrayLike", lmass: "ArrayLike", umass: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the CRPS for the generalised truncated and censored normal distribution."""
+    B = backends.active if backend is None else backends[backend]
+    mu, sigma, lower, upper, lmass, umass, obs = map(B.asarray, (location, scale, lower, upper, lmass, umass, obs))
+    ω = (obs - mu) / sigma
+    u = (upper - mu) / sigma
+    l = (lower - mu) / sigma
+    z = B.minimum(B.maximum(ω, l), u)
+    F_u = _norm_cdf(u, backend=backend)
+    F_l = _norm_cdf(l, backend=backend)
+    F_z = _norm_cdf(z, backend=backend)
+    F_u2 = _norm_cdf(u * B.sqrt(2), backend=backend)
+    F_l2 = _norm_cdf(l * B.sqrt(2), backend=backend)
+    f_u = _norm_pdf(u, backend=backend)
+    f_l = _norm_pdf(l, backend=backend)
+    f_z = _norm_pdf(z, backend=backend)
+
+    u_inf = u == float("inf")
+    l_inf = l == float("-inf")
+
+    u = B.where(u_inf, B.nan, u)
+    l = B.where(l_inf, B.nan, l)
+    s1_u = B.where(u_inf and umass == 0.0, 0.0, u * umass**2)
+    s1_l = B.where(l_inf and lmass == 0.0, 0.0, l * lmass**2)
+
+    c = (1 - lmass - umass) / (F_u - F_l)
+
+    s1 = B.abs(ω - z) + s1_u - s1_l
+    s2 = c * z * (2 * F_z - ((1 - 2 * lmass) * F_u + (1 - 2 * umass) * F_l) / (1 - lmass - umass))
+    s3 = c * (2 * f_z - 2 * f_u * umass - 2 * f_l * lmass)
+    s4 = c**2 * (F_u2 - F_l2) / B.sqrt(B.pi)
+    return sigma * (s1 + s2 + s3 - s4)
+
+
 def hypergeometric(
     obs: "ArrayLike",
     m: "ArrayLike",

tests/test_crps.py

@@ -276,6 +276,50 @@ def test_gpd(backend):
     )
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_gtcnormal(backend):
+    obs, location, scale, lower, upper, lmass, umass = 0.9, -2.3, 4.1, -7.3, 1.7, 0.0, 0.21
+    expected = 1.422805
+    res = _crps.crps_gtcnormal(obs, location, scale, lower, upper, lmass, umass, backend=backend)
+    assert np.isclose(res, expected)
+
+    # aligns with crps_normal
+    res0 = _crps.crps_normal(obs, location, scale, backend=backend)
+    res = _crps.crps_gtcnormal(obs, location, scale, backend=backend)
+    assert np.isclose(res, res0)
+
+    # aligns with crps_tnormal
+    res0 = _crps.crps_tnormal(obs, location, scale, lower, upper, backend=backend)
+    res = _crps.crps_gtcnormal(obs, location, scale, lower, upper, backend=backend)
+    assert np.isclose(res, res0)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_tnormal(backend):
+    obs, location, scale, lower, upper = -1.0, 2.9, 2.2, 1.5, 17.3
+    expected = 3.982434
+    res = _crps.crps_tnormal(obs, location, scale, lower, upper, backend=backend)
+    assert np.isclose(res, expected)
+
+    # aligns with crps_normal
+    res0 = _crps.crps_normal(obs, location, scale, backend=backend)
+    res = _crps.crps_tnormal(obs, location, scale, backend=backend)
+    assert np.isclose(res, res0)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_cnormal(backend):
+    obs, location, scale, lower, upper = 1.8, 0.4, 1.1, 0.0, 2.0
+    expected = 0.8296078
+    res = _crps.crps_cnormal(obs, location, scale, lower, upper, backend=backend)
+    assert np.isclose(res, expected)
+
+    # aligns with crps_normal
+    res0 = _crps.crps_normal(obs, location, scale, backend=backend)
+    res = _crps.crps_cnormal(obs, location, scale, backend=backend)
+    assert np.isclose(res, res0)
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_hypergeometric(backend):
     res = _crps.crps_hypergeometric(5 * np.ones((2, 2)), 7, 13, 12, backend=backend)