Add utility functions and refactor score function interface

scoringrules/_brier.py

@@ -1,7 +1,7 @@
 import typing as tp
 
-from scoringrules.backend import backends
 from scoringrules.core import brier
+from scoringrules.core.utils import binary_array_check, categorical_array_check
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, ArrayLike, Backend
@@ -14,30 +14,43 @@ def brier_score(
     backend: "Backend" = None,
 ) -> "Array":
     r"""
-    Compute the Brier Score (BS).
+    Brier Score
 
-    The BS is formulated as
+    The Brier Score is defined as
 
     .. math::
-       BS(f, y) = (f - y)^2,
+       \text{BS}(F, y) = (F - y)^2,
 
-    where :math:`f \in [0, 1]` is the predicted probability of an event and :math:`y \in \{0, 1\}` the actual outcome.
+    where :math:`F \in [0, 1]` is the predicted probability of an event and :math:`y \in \{0, 1\}`
+    is the outcome [1]_.
 
     Parameters
     ----------
     obs : array_like
-        Observed outcome, either 0 or 1.
+        Observed outcomes, either 0 or 1.
     fct : array_like
-        Forecasted probabilities between 0 and 1.
+        Forecast probabilities, between 0 and 1.
     backend : str
-        The name of the backend used for computations. Defaults to 'numpy'.
+        The name of the backend used for computations. Default is 'numpy'.
 
     Returns
     -------
-    brier_score : array_like
+    score : array_like
         The computed Brier Score.
 
+    References
+    ----------
+    .. [1] Brier, G. W. (1950).
+        Verification of forecasts expressed in terms of probability.
+        Monthly Weather Review, 78, 1-3.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.brier_score(1, 0.2)
+    0.64000
     """
+    obs, fct = binary_array_check(obs, fct, backend=backend)
     return brier.brier_score(obs=obs, fct=fct, backend=backend)
 
 
@@ -46,45 +59,68 @@ def rps_score(
     fct: "ArrayLike",
     k_axis: int = -1,
     *,
+    onehot: bool = False,
     backend: "Backend" = None,
 ) -> "Array":
     r"""
-    Compute the (Discrete) Ranked Probability Score (RPS).
+    (Discrete) Ranked Probability Score (RPS)
 
     Suppose the outcome corresponds to one of :math:`K` ordered categories. The RPS is defined as
 
     .. math::
-       RPS(f, y) = \sum_{k=1}^{K}(\tilde{f}_{k} - \tilde{y}_{k})^2,
+       \text{RPS}(F, y) = \sum_{k=1}^{K}(\tilde{F}_{k} - \tilde{y}_{k})^2,
 
-    where :math:`f \in [0, 1]^{K}` is a vector of length :math:`K` containing forecast probabilities
-    that each of the :math:`K` categories will occur, and :math:`y \in \{0, 1\}^{K}` is a vector of
-    length :math:`K`, with the :math:`k`-th element equal to one if the :math:`k`-th category occurs. We
-    have :math:`\sum_{k=1}^{K} y_{k} = \sum_{k=1}^{K} f_{k} = 1`, and, for :math:`k = 1, \dots, K`,
-    :math:`\tilde{y}_{k} = \sum_{i=1}^{k} y_{i}` and :math:`\tilde{f}_{k} = \sum_{i=1}^{k} f_{i}`.
+    where :math:`F \in [0, 1]^{K}` is a vector of length :math:`K`, containing forecast probabilities
+    that each of the :math:`K` categories will occur, with :math:`\sum_{k=1}^{K} F_{k} = 1` and
+    :math:`\tilde{F}_{k} = \sum_{i=1}^{k} F_{i}` for all :math:`k = 1, \dots, K`, and where
+    :math:`y \in \{1, \dots, K\}` is the category that occurs, with :math:`\tilde{y}_{k} = 1\{y \le i\}`
+    for all :math:`k = 1, \dots, K` [1]_.
+
+    The outcome can alternatively be interpreted as a vector :math:`y \in \{0, 1\}^K` of length :math:`K`, with the
+    :math:`k`-th element equal to one if the :math:`k`-th category occurs, and zero otherwise.
+    Using this one-hot encoding, the RPS is defined analogously to as above, but with
+    :math:`\tilde{y}_{k} = \sum_{i=1}^{k} y_{i}`.
 
     Parameters
     ----------
     obs : array_like
-        Array of 0's and 1's corresponding to unobserved and observed categories
-    forecasts :
+        Category that occurs. Or array of 0's and 1's corresponding to unobserved and
+        observed categories if `onehot=True`.
+    fct : array
         Array of forecast probabilities for each category.
     k_axis: int
-        The axis corresponding to the categories. Default is the last axis.
+        The axis of `obs` and `fct` corresponding to the categories. Default is the last axis.
+    onehot: bool
+        Boolean indicating whether the observation is the category that occurs or a onehot
+        encoded vector of 0's and 1's. Default is False.
     backend : str
-        The name of the backend used for computations. Defaults to 'numpy'.
+        The name of the backend used for computations. Default is 'numpy'.
 
     Returns
     -------
-    score:
+    score: array_like
         The computed Ranked Probability Score.
 
+    References
+    ----------
+    .. [1] Epstein, E. S. (1969).
+        A scoring system for probability forecasts of ranked categories.
+        Journal of Applied Meteorology, 8, 985-987.
+        https://www.jstor.org/stable/26174707.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> import numpy as np
+    >>> fct = np.array([0.1, 0.2, 0.3, 0.4])
+    >>> obs = 3
+    >>> sr.rps_score(obs, fct)
+    0.25999999999999995
+    >>> obs = np.array([0, 0, 1, 0])
+    >>> sr.rps_score(obs, fct, onehot=True)
+    0.25999999999999995
     """
-    B = backends.active if backend is None else backends[backend]
-    fct = B.asarray(fct)
-
-    if k_axis != -1:
-        fct = B.moveaxis(fct, k_axis, -1)
-
+    obs, fct = categorical_array_check(obs, fct, k_axis, onehot, backend=backend)
     return brier.rps_score(obs=obs, fct=fct, backend=backend)
 
 
@@ -119,6 +155,7 @@ def log_score(
         The computed Log Score.
 
     """
+    obs, fct = binary_array_check(obs, fct, backend=backend)
     return brier.log_score(obs=obs, fct=fct, backend=backend)
 
 
@@ -127,6 +164,7 @@ def rls_score(
     fct: "ArrayLike",
     k_axis: int = -1,
     *,
+    onehot: bool = False,
     backend: "Backend" = None,
 ) -> "Array":
     r"""
@@ -151,6 +189,9 @@ def rls_score(
         Forecasted probabilities between 0 and 1.
     k_axis: int
         The axis corresponding to the categories. Default is the last axis.
+    onehot: bool
+        Boolean indicating whether the observation is the category that occurs or a onehot
+        encoded vector of 0's and 1's. Default is False.
     backend : str
         The name of the backend used for computations. Defaults to 'numpy'.
 
@@ -160,12 +201,7 @@ def rls_score(
         The computed Ranked Logarithmic Score.
 
     """
-    B = backends.active if backend is None else backends[backend]
-    fct = B.asarray(fct)
-
-    if k_axis != -1:
-        fct = B.moveaxis(fct, k_axis, -1)
-
+    obs, fct = categorical_array_check(obs, fct, k_axis, onehot, backend=backend)
     return brier.rls_score(obs=obs, fct=fct, backend=backend)
 
 

scoringrules/_crps.py

@@ -2,6 +2,13 @@
 
 from scoringrules.backend import backends
 from scoringrules.core import crps, stats
+from scoringrules.core.utils import (
+    univariate_array_check,
+    estimator_check,
+    uv_weighted_score_weights,
+    uv_weighted_score_chain,
+    univariate_sort_ens,
+)
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, ArrayLike, Backend
@@ -95,29 +102,13 @@ def crps_ensemble(
     >>> sr.crps_ensemble(obs, pred)
     array([0.69605316, 0.32865417, 0.39048665])
     """
-    B = backends.active if backend is None else backends[backend]
-    obs, fct = map(B.asarray, (obs, fct))
-
-    if m_axis != -1:
-        fct = B.moveaxis(fct, m_axis, -1)
-
-    if not sorted_ensemble and estimator not in [
-        "nrg",
-        "akr",
-        "akr_circperm",
-        "fair",
-    ]:
-        fct = B.sort(fct, axis=-1)
-
+    obs, fct = univariate_array_check(obs, fct, m_axis, backend=backend)
+    fct = univariate_sort_ens(fct, estimator, sorted_ensemble, backend=backend)
     if backend == "numba":
-        if estimator not in crps.estimator_gufuncs:
-            raise ValueError(
-                f"{estimator} is not a valid estimator. "
-                f"Must be one of {crps.estimator_gufuncs.keys()}"
-            )
+        estimator_check(estimator, crps.estimator_gufuncs)
         return crps.estimator_gufuncs[estimator](obs, fct)
-
-    return crps.ensemble(obs, fct, estimator, backend=backend)
+    else:
+        return crps.ensemble(obs, fct, estimator, backend=backend)
 
 
 def twcrps_ensemble(
@@ -202,14 +193,9 @@ def twcrps_ensemble(
     >>> sr.twcrps_ensemble(obs, fct, v_func=v_func)
     array([0.69605316, 0.32865417, 0.39048665])
     """
-    if v_func is None:
-        B = backends.active if backend is None else backends[backend]
-        a, b, obs, fct = map(B.asarray, (a, b, obs, fct))
-
-        def v_func(x):
-            return B.minimum(B.maximum(x, a), b)
-
-    obs, fct = map(v_func, (obs, fct))
+    obs, fct = uv_weighted_score_chain(
+        obs, fct, a=a, b=b, v_func=v_func, backend=backend
+    )
     return crps_ensemble(
         obs,
         fct,
@@ -300,25 +286,12 @@ def owcrps_ensemble(
     >>> sr.owcrps_ensemble(obs, fct, w_func=w_func)
     array([0.91103733, 0.45212402, 0.35686667])
     """
-
-    B = backends.active if backend is None else backends[backend]
-    obs, fct = map(B.asarray, (obs, fct))
-
-    if m_axis != -1:
-        fct = B.moveaxis(fct, m_axis, -1)
-
-    if w_func is None:
-
-        def w_func(x):
-            return ((a <= x) & (x <= b)) * 1.0
-
-    obs_weights, fct_weights = map(w_func, (obs, fct))
-    obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights))
-
+    obs, fct = univariate_array_check(obs, fct, m_axis, backend=backend)
+    obs_w, fct_w = uv_weighted_score_weights(obs, fct, a, b, w_func, backend=backend)
     if backend == "numba":
-        return crps.estimator_gufuncs["ownrg"](obs, fct, obs_weights, fct_weights)
-
-    return crps.ow_ensemble(obs, fct, obs_weights, fct_weights, backend=backend)
+        return crps.estimator_gufuncs["ownrg"](obs, fct, obs_w, fct_w)
+    else:
+        return crps.ow_ensemble(obs, fct, obs_w, fct_w, backend=backend)
 
 
 def vrcrps_ensemble(
@@ -399,24 +372,12 @@ def vrcrps_ensemble(
     >>> sr.vrcrps_ensemble(obs, fct, w_func)
     array([0.90036433, 0.41515255, 0.41653833])
     """
-    B = backends.active if backend is None else backends[backend]
-    obs, fct = map(B.asarray, (obs, fct))
-
-    if m_axis != -1:
-        fct = B.moveaxis(fct, m_axis, -1)
-
-    if w_func is None:
-
-        def w_func(x):
-            return ((a <= x) & (x <= b)) * 1.0
-
-    obs_weights, fct_weights = map(w_func, (obs, fct))
-    obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights))
-
+    obs, fct = univariate_array_check(obs, fct, m_axis, backend=backend)
+    obs_w, fct_w = uv_weighted_score_weights(obs, fct, a, b, w_func, backend=backend)
     if backend == "numba":
-        return crps.estimator_gufuncs["vrnrg"](obs, fct, obs_weights, fct_weights)
-
-    return crps.vr_ensemble(obs, fct, obs_weights, fct_weights, backend=backend)
+        return crps.estimator_gufuncs["vrnrg"](obs, fct, obs_w, fct_w)
+    else:
+        return crps.vr_ensemble(obs, fct, obs_w, fct_w, backend=backend)
 
 
 def crps_quantile(
@@ -473,10 +434,11 @@ def crps_quantile(
     # TODO: add example
     """
     B = backends.active if backend is None else backends[backend]
-    obs, fct, alpha = map(B.asarray, (obs, fct, alpha))
+    obs, fct = univariate_array_check(obs, fct, m_axis, backend=backend)
 
-    if m_axis != -1:
-        fct = B.moveaxis(fct, m_axis, -1)
+    alpha = B.asarray(alpha)
+    if B.any(alpha <= 0) or B.any(alpha >= 1):
+        raise ValueError("`alpha` contains entries that are not between 0 and 1.")
 
     if not fct.shape[-1] == alpha.shape[-1]:
         raise ValueError("Expected matching length of `fct` and `alpha` values.")

scoringrules/_dss.py

@@ -1,8 +1,7 @@
 import typing as tp
 
-from scoringrules.backend import backends
 from scoringrules.core import dss
-from scoringrules.core.utils import multivariate_array_check
+from scoringrules.core.utils import univariate_array_check, multivariate_array_check
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, Backend
@@ -45,16 +44,11 @@ def dssuv_ensemble(
     score: Array
         The computed Dawid-Sebastiani Score.
     """
-    B = backends.active if backend is None else backends[backend]
-    obs, fct = map(B.asarray, (obs, fct))
-
-    if m_axis != -1:
-        fct = B.moveaxis(fct, m_axis, -1)
-
+    obs, fct = univariate_array_check(obs, fct, m_axis, backend=backend)
     if backend == "numba":
         return dss._dss_uv_gufunc(obs, fct, bias)
-
-    return dss.ds_score_uv(obs, fct, bias, backend=backend)
+    else:
+        return dss.ds_score_uv(obs, fct, bias, backend=backend)
 
 
 def dssmv_ensemble(
@@ -99,8 +93,7 @@ def dssmv_ensemble(
         The computed Dawid-Sebastiani Score.
     """
     obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend)
-
     if backend == "numba":
         return dss._dss_mv_gufunc(obs, fct, bias)
-
-    return dss.ds_score_mv(obs, fct, bias, backend=backend)
+    else:
+        return dss.ds_score_mv(obs, fct, bias, backend=backend)