Add quantile score (#78)

* implement quantile score

* [skip ci] fix docstring

Author: frazane

docs/api/interval.md

@@ -4,23 +4,29 @@
 
 For a prediction interval (PI), the interval or Winkler score is given by:
 
-$\text{IS} = \begin{cases}
+$$
+\text{IS} = \begin{cases}
     (u - l) + \frac{2}{\alpha}(l - y)  & \text{for } y < l \\
     (u - l)                             & \text{for } l \leq y \leq u \\
     (u - l) + \frac{2}{\alpha}(y - u)   & \text{for } y > u. \\
-\end{cases}$
+\end{cases}
+$$
 
 for an $(1 - \alpha)$PI of $[l, u]$ and the true value $y$ [@gneiting_strictly_2007, @bracher2021evaluating @winkler1972decision].
 
 ## Weighted Interval Score
 
 The weighted interval score (WIS) is defined as
 
-$\text{WIS}_{\alpha_{0:K}}(F, y) = \frac{1}{K+0.5}(w_0 \times |y - m| + \sum_{k=1}^K (w_k \times IS_{\alpha_k}(F, y)))$
+$$
+\text{WIS}_{\alpha_{0:K}}(F, y) = \frac{1}{K+0.5}(w_0 \times |y - m| + \sum_{k=1}^K (w_k \times IS_{\alpha_k}(F, y)))
+$$
 
 where $m$ denotes the median prediction, $w_0$ denotes the weight of the median prediction, $IS_{\alpha_k}(F, y)$ denotes the interval score for the $1 - \alpha$ prediction interval and $w_k$ is the according weight. The WIS is calculated for a set of (central) PIs and the predictive median [@bracher2021evaluating]. The weights are an optional parameter and default weight is the canonical weight $w_k = \frac{2}{\alpha_k}$ and $w_0 = 0.5$. For these weights, it holds that:
 
-$\text{WIS}_{\alpha_{0:K}}(F, y) \approx \text{CRPS}(F, y).$
+$$
+\text{WIS}_{\alpha_{0:K}}(F, y) \approx \text{CRPS}(F, y).
+$$
 
 
 ::: scoringrules.interval_score

docs/api/quantile.md

@@ -0,0 +1,3 @@
+# Quantile score
+
+::: scoringrules.quantile_score

mkdocs.yaml

@@ -15,6 +15,7 @@ nav:
     - Energy Score: api/energy.md
     - Variogram Score: api/variogram.md
     - Interval Score: api/interval.md
+    - Quantile Score: api/quantile.md
     - Visualization: api/visualization.md
 
 

scoringrules/__init__.py

@@ -96,6 +96,8 @@
     vrgksmv_ensemble,
 )
 
+from scoringrules._quantile import quantile_score
+
 from scoringrules.backend import backends, register_backend
 
 __version__ = version("scoringrules")
@@ -188,4 +190,5 @@
     "vrgksmv_ensemble",
     "interval_score",
     "weighted_interval_score",
+    "quantile_score",
 ]

scoringrules/_quantile.py

@@ -0,0 +1,58 @@
+import typing as tp
+
+from scoringrules.backend import backends
+
+if tp.TYPE_CHECKING:
+    from scoringrules.core.typing import Array, ArrayLike, Backend
+
+
+def quantile_score(
+    obs: "ArrayLike",
+    fct: "ArrayLike",
+    alpha: "ArrayLike",
+    backend: "Backend | None" = None,
+) -> "Array":
+    r"""Compute the quantile score for a given quantile level.
+
+    The quantile score (Koenker, R. and G. Bassett, 1978) is defined as
+
+    $$
+        S_{\alpha}(q_{\alpha}, y) = \begin{cases}
+        (1 - \alpha) (q_{\alpha} - y), & \text{if } y \leq q_{\alpha}, \\
+        \alpha (y - q_{\alpha}), & \text{if } y > q_{\alpha}.
+        \end{cases}
+    $$
+
+    where $y$ is the observed value and $q_{\alpha}$ is the predicted value at the
+    $\alpha$ quantile level.
+
+    Parameters
+    ----------
+    obs:
+        The observed values.
+    fct:
+        The forecast values.
+    alpha:
+        The quantile level.
+
+    Returns
+    -------
+    score:
+        The quantile score.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.quantile_score(0.3, 0.5, 0.2)
+    0.16
+
+    Raises
+    ------
+    ValueError
+        If the quantile level is not between 0 and 1.
+    """
+    B = backends.active if backend is None else backends[backend]
+    obs, fct, alpha = map(B.asarray, (obs, fct, alpha))
+    if B.any(alpha < 0) or B.any(alpha > 1):
+        raise ValueError("The quantile level must be between 0 and 1.")
+    return (B.asarray((obs < fct), dtype=alpha.dtype) - alpha) * (fct - obs)

tests/test_quantile.py

@@ -0,0 +1,23 @@
+import numpy as np
+import pytest
+
+import scoringrules as sr
+
+from .conftest import BACKENDS
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_quantile_score(backend):
+    obs = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
+    fct = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
+    alpha = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
+    res = sr.quantile_score(obs, fct, alpha, backend=backend)
+    res.__class__.__module__ == backend
+
+    # correctness
+    np.isclose(sr.quantile_score(0.3, 0.5, 0.2, backend=backend), 0.16)
+    np.isclose(sr.quantile_score(0.3, 0.1, 0.2, backend=backend), 0.04)
+
+    # value checks
+    with pytest.raises(ValueError):
+        sr.quantile_score(0.3, 0.5, -0.2, backend=backend)