Fix WIS: use w_median * |obs - median| on the array backends

scoringrules/_interval.py

@@ -132,7 +132,7 @@ def weighted_interval_score(
     :math:`1 - \alpha` prediction interval and :math:`w_k` is the according weight.
     The WIS is calculated for a set of (central) PIs and the predictive median.
     The weights are an optional parameter and default weight is the canonical
-    weight :math:`w_k = \frac{2}{\alpha_k}` and :math:`w_0 = 0.5`.
+    weight :math:`w_k = \frac{\alpha_k}{2}` and :math:`w_0 = 0.5`.
     For these weights, it holds that:
 
     .. math::

scoringrules/core/interval/_score.py

@@ -36,6 +36,6 @@ def weighted_interval_score(
     K = w_alpha.shape[0]
     IS = interval_score(obs[..., None], lower, upper, alpha)
     WIS = B.sum(IS * w_alpha, axis=-1)
-    WIS += w_median * median
+    WIS += w_median * B.abs(obs - median)
     WIS /= K + 1 / 2
     return B.squeeze(WIS)

tests/test_interval.py

@@ -56,13 +56,36 @@ def test_interval_score(backend):
 
 ## We use Bracher et al (2021) Eq. (3) to test the WIS
 def test_weighted_interval_score(backend):
-    obs = np.zeros(N)
+    # `obs` differs from the predictive median, so the median term contributes
+    # `w_median * |obs - median|` and cannot be confused with `w_median * median`.
+    obs = np.full(N, 0.4)
+    median = np.zeros(N)
     alpha = np.linspace(0.01, 0.99, 99)
     upper = st.norm(0, 1).ppf(np.tile(1 - alpha / 2, (N, 1)))
     lower = st.norm(0, 1).ppf(np.tile(alpha / 2, (N, 1)))
 
-    WIS = sr.weighted_interval_score(obs, obs, lower, upper, alpha, backend=backend)
+    WIS = sr.weighted_interval_score(obs, median, lower, upper, alpha, backend=backend)
     CRPS = sr.crps_normal(obs, 0, 1, backend=backend)
     WIS, CRPS = map(np.asarray, (WIS, CRPS))
 
     assert np.all(1 - WIS / CRPS <= 0.001 * CRPS)
+
+
+def test_weighted_interval_score_median_term(backend):
+    # Closed form of Bracher et al (2021) Eq. (3) for a single set of intervals.
+    obs, median = 1.4, 1.0
+    alpha = np.array([0.2, 0.5])
+    lower = np.array([0.5, 0.8])
+    upper = np.array([1.2, 1.1])
+
+    interval_scores = (
+        (upper - lower)
+        + 2 / alpha * (lower - obs) * (obs < lower)
+        + 2 / alpha * (obs - upper) * (obs > upper)
+    )
+    expected = (0.5 * abs(obs - median) + np.sum(alpha / 2 * interval_scores)) / (
+        len(alpha) + 0.5
+    )
+
+    res = sr.weighted_interval_score(obs, median, lower, upper, alpha, backend=backend)
+    assert np.isclose(np.asarray(res), expected)