Dawid-Sebastiani Score

scoringrules/__init__.py

@@ -43,6 +43,7 @@
     twenergy_score,
     vrenergy_score,
 )
+from scoringrules._dawid_sebastiani import dawid_sebastiani_score
 from scoringrules._error_spread import error_spread_score
 from scoringrules._interval import interval_score, weighted_interval_score
 from scoringrules._logs import (
@@ -113,6 +114,7 @@
     "crps_quantile",
     "crps_t",
     "crps_uniform",
+    "dawid_sebastiani_score",
     "owcrps_ensemble",
     "twcrps_ensemble",
     "vrcrps_ensemble",

scoringrules/_dawid_sebastiani.py

@@ -0,0 +1,49 @@
+import typing as tp
+
+from scoringrules.core import dss
+from scoringrules.core.utils import multivariate_array_check
+
+if tp.TYPE_CHECKING:
+    from scoringrules.core.typing import Array, Backend
+
+
+def dawid_sebastiani_score(
+    observations: "Array",
+    forecasts: "Array",
+    /,
+    m_axis: int = -2,
+    v_axis: int = -1,
+    *,
+    backend: "Backend" = None,
+) -> "Array":
+    r"""Compute the Dawid-Sebastiani-Score for a finite multivariate ensemble.
+
+    ## TO-DO ##
+
+    Parameters
+    ----------
+    forecasts: Array
+        The predicted forecast ensemble, where the ensemble dimension is by default
+        represented by the second last axis and the variables dimension by the last axis.
+    observations: Array
+        The observed values, where the variables dimension is by default the last axis.
+    m_axis: int
+        The axis corresponding to the ensemble dimension. Defaults to -2.
+    v_axis: int
+        The axis corresponding to the variables dimension. Defaults to -1.
+    backend: str
+        The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'.
+
+    Returns
+    -------
+    ds_score: Array
+        The computed Dawid-Sebastiani-Score.
+    """
+    observations, forecasts = multivariate_array_check(
+        observations, forecasts, m_axis, v_axis, backend=backend
+    )
+
+    if backend == "numba":
+        return dss._dss_gufunc(observations, forecasts)
+
+    return dss._dss(observations, forecasts, backend=backend)

scoringrules/backend/numpy.py

@@ -224,6 +224,20 @@ def where(self, condition: "NDArray", x1: "NDArray", x2: "NDArray") -> "NDArray"
     def size(self, x: "NDArray") -> int:
         return x.size
 
+    def inv(self, x: "NDArray") -> "NDArray":
+        return np.linalg.inv(x)
+
+    def cov(self, x: "NDArray") -> "NDArray":
+        # TODO: This deserves some axis parameter probably?
+        if len(x.shape) == 2:
+            return np.cov(x)
+        else:
+            centered = x - x.mean(axis=-2, keepdims=True)
+            return (centered.transpose(0, 2, 1) @ centered) / (x.shape[-2] - 1)
+
+    def det(self, x: "NDArray") -> "NDArray":
+        return np.linalg.det(x)
+
 
 class NumbaBackend(NumpyBackend):
     """Numba backend."""

scoringrules/core/dss/__init__.py

@@ -0,0 +1,10 @@
+try:
+    from ._gufuncs import (
+        _dss_gufunc,
+    )
+except ImportError:
+    _dss_gufunc = None
+
+from ._score import _ds_score as _dss
+
+__all__ = ["_dss_gufunc", "_dss"]

scoringrules/core/dss/_gufuncs.py

@@ -0,0 +1,19 @@
+import numpy as np
+from numba import guvectorize
+
+
+@guvectorize(
+    [
+        "void(float32[:], float32[:,:], float32[:])",
+        "void(float64[:], float64[:,:], float64[:])",
+    ],
+    "(d),(m,d)->()",
+)
+def _dss_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray):
+    M = fct.shape[0]
+    bias = obs - (np.sum(fct, axis=0) / M)
+    cov = np.cov(fct, rowvar=False).astype(bias.dtype)
+    prec = np.linalg.inv(cov).astype(bias.dtype)
+    log_det = np.log(np.linalg.det(cov))
+    bias_precision = bias.T @ prec @ bias
+    out[0] = log_det + bias_precision

scoringrules/core/dss/_score.py

@@ -0,0 +1,21 @@
+import typing as tp
+
+from scoringrules.backend import backends
+
+if tp.TYPE_CHECKING:
+    from scoringrules.core.typing import Array, Backend
+
+
+def _ds_score(
+    obs: "Array",  # (... D)
+    fct: "Array",  # (... M D)
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the Dawid Sebastiani score Score for a multivariate finite ensemble."""
+    B = backends.active if backend is None else backends[backend]
+    cov = B.cov(fct)
+    precision = B.inv(cov)
+    log_det = B.log(B.det(cov))
+    bias = obs - B.mean(fct, axis=-2)
+    bias_precision = B.squeeze(bias[:, None, :] @ precision @ bias[:, :, None])
+    return log_det + bias_precision