Add censored and conditional likelihood scores

docs/api/logarithmic.md

@@ -49,3 +49,7 @@
 ::: scoringrules.logs_tt
 
 ::: scoringrules.logs_uniform
+
+## Conditional and Censored Likelihood Score
+
+::: scoringrules.clogs_ensemble

scoringrules/__init__.py

@@ -70,6 +70,7 @@
     logs_tnormal,
     logs_tt,
     logs_uniform,
+    clogs_ensemble,
 )
 from scoringrules._variogram import (
     owvariogram_score,
@@ -148,6 +149,7 @@
     "logs_tnormal",
     "logs_tt",
     "logs_uniform",
+    "clogs_ensemble",
     "brier_score",
     "gksuv_ensemble",
     "gksmv_ensemble",

scoringrules/_logs.py

@@ -72,6 +72,88 @@ def logs_ensemble(
     return logarithmic.mixnorm(observations, forecasts, bw, w, backend=backend)
 
 
+def clogs_ensemble(
+    observations: "ArrayLike",
+    forecasts: "Array",
+    /,
+    a: "ArrayLike" = float("-inf"),
+    b: "ArrayLike" = float("inf"),
+    axis: int = -1,
+    *,
+    bw: "ArrayLike" = None,
+    cens: bool = True,
+    backend: "Backend" = None,
+) -> "Array":
+    r"""Estimate the conditional and censored likelihood score for an ensemble forecast.
+
+    The conditional and censored likelihood scores are introduced by
+    [Diks et al. (2011)](https://doi.org/10.1016/j.jeconom.2011.04.001):
+
+    The weight function is an indicator function of the form $w(z) = 1\{a < z < b\}$.
+
+    The ensemble forecast is converted to a mixture of normal distributions using Gaussian
+    kernel density estimation. The score is then calculated for this smoothed distribution.
+
+    Parameters
+    ----------
+    observations: ArrayLike
+        The observed values.
+    forecasts: ArrayLike
+        The predicted forecast ensemble, where the ensemble dimension is by default
+        represented by the last axis.
+    a: ArrayLike
+        The lower bound in the weight function.
+    b: ArrayLike
+        The upper bound in the weight function.
+    axis: int
+        The axis corresponding to the ensemble. Default is the last axis.
+    bw : ArrayLike
+        The bandwidth parameter for each forecast ensemble. If not given, estimated using
+        Silverman's rule of thumb.
+    cens : Boolean
+        Boolean specifying whether to return the conditional ('cens = False') or the censored
+        likelihood score ('cens = True').
+    backend: str
+        The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'.
+
+    Returns
+    -------
+    score:
+        The CoLS or CeLS between the forecast ensemble and obs for the chosen weight parameters.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.clogs_ensemble(obs, pred, -1.0, 1.0)
+    """
+    B = backends.active if backend is None else backends[backend]
+    forecasts = B.asarray(forecasts)
+
+    if axis != -1:
+        forecasts = B.moveaxis(forecasts, axis, -1)
+
+    M = forecasts.shape[-1]
+
+    # Silverman's rule of thumb for estimating the bandwidth parameter
+    if bw is None:
+        sigmahat = B.std(forecasts, axis=-1)
+        q75 = B.quantile(forecasts, 0.75, axis=-1)
+        q25 = B.quantile(forecasts, 0.25, axis=-1)
+        iqr = q75 - q25
+        bw = 1.06 * B.minimum(sigmahat, iqr / 1.34) * (M ** (-1 / 5))
+    bw = B.stack([bw] * M, axis=-1)
+
+    return logarithmic.clogs_ensemble(
+        observations,
+        forecasts,
+        a=a,
+        b=b,
+        bw=bw,
+        cens=cens,
+        backend=backend,
+    )
+
+
 def logs_beta(
     observation: "ArrayLike",
     a: "ArrayLike",
@@ -1010,6 +1092,7 @@ def logs_uniform(
 __all__ = [
     "logs_beta",
     "logs_binomial",
+    "logs_ensemble",
     "logs_exponential",
     "logs_exponential2",
     "logs_2pexponential",
@@ -1028,4 +1111,5 @@ def logs_uniform(
     "logs_poisson",
     "logs_t",
     "logs_uniform",
+    "clogs_ensemble",
 ]

scoringrules/core/logarithmic.py

@@ -23,6 +23,44 @@
     from scoringrules.core.typing import Array, ArrayLike, Backend
 
 
+def clogs_ensemble(
+    obs: "ArrayLike",
+    fct: "Array",
+    a: "ArrayLike",
+    b: "ArrayLike",
+    bw: "ArrayLike",
+    cens: bool = True,
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the conditional or censored likelihood score for an ensemble forecast."""
+    B = backends.active if backend is None else backends[backend]
+    a, b, bw, fct, obs = map(B.asarray, (a, b, bw, fct, obs))
+
+    z = (obs[..., None] - fct) / bw
+    a_z = (a - fct) / bw
+    b_z = (b - fct) / bw
+
+    dens = _norm_pdf(z, backend=backend) / bw
+    dens = B.mean(dens, axis=-1)
+    s1 = -B.log(dens)
+
+    F_a = _norm_cdf(a_z, backend=backend)
+    F_a = B.mean(F_a, axis=-1)
+    F_b = _norm_cdf(b_z, backend=backend)
+    F_b = B.mean(F_b, axis=-1)
+    cons = F_b - F_a
+
+    w_ind = (obs > a) & (obs < b)
+    if cens:
+        s2 = -B.log(1 - cons)
+        s = B.where(w_ind, s1, s2)
+    else:
+        s2 = B.log(cons)
+        s = B.where(w_ind, s1 + s2, 0.0)
+
+    return s
+
+
 def beta(
     obs: "ArrayLike",
     a: "ArrayLike",

tests/test_logs.py

@@ -33,6 +33,43 @@ def test_ensemble(backend):
     assert np.isclose(res, expected)
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_clogs(backend):
+    obs = np.random.randn(N)
+    mu = obs + np.random.randn(N) * 0.1
+    sigma = abs(np.random.randn(N))
+    fct = np.random.randn(N, ENSEMBLE_SIZE) * sigma[..., None] + mu[..., None]
+
+    res0 = _logs.logs_ensemble(obs, fct, axis=-1, backend=backend)
+    res = _logs.clogs_ensemble(obs, fct, axis=-1, backend=backend)
+    res_co = _logs.clogs_ensemble(obs, fct, axis=-1, cens=False, backend=backend)
+    assert res.shape == (N,)
+    assert res_co.shape == (N,)
+    assert np.allclose(res, res0, atol=1e-7)
+    assert np.allclose(res_co, res0, atol=1e-7)
+
+    fct = fct.T
+    res0 = _logs.clogs_ensemble(obs, fct, axis=0, backend=backend)
+    assert np.allclose(res, res0, atol=1e-7)
+
+    obs, fct = 6.2, [4.2, 5.1, 6.1, 7.6, 8.3, 9.5]
+    a = 8.0
+    res = _logs.clogs_ensemble(obs, fct, a=a, backend=backend)
+    expected = 0.3929448
+    assert np.isclose(res, expected)
+
+    a = 5.0
+    res = _logs.clogs_ensemble(obs, fct, a=a, cens=False, backend=backend)
+    expected = 1.646852
+    assert np.isclose(res, expected)
+
+    fct = [-142.7, -160.3, -179.4, -184.5]
+    b = -150.0
+    res = _logs.clogs_ensemble(obs, fct, b=b, backend=backend)
+    expected = 1.420427
+    assert np.isclose(res, expected)
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_beta(backend):
     if backend == "torch":