weighted_interval_score adds w_median * median instead of w_median * |obs - median| on all non-numba backends

The issue identifies scoringrules/core/interval/_score.py and tests/test_interval.py. Read the weighted interval score implementation and its existing test, then use a median different from the observation in a regression case. Done means the non-numba backends match the stated formula and the regression test catches the previous result; a linked open pull request is already addressing this issue.

What happened?

weighted_interval_score returns incorrect values under the numpy backend. The numba backend produces correct results, meaning the same function yields different outputs across backends.

Line 39 in scoringrules/core/interval/_score.py reads WIS += w_median * median. This adds the median itself instead of the absolute deviation between the observation and the median. Per equation (3) in Bracher et al. (2021), this term should be w_median * |obs - median|. The numba implementation in core/interval/_gufunc.py correctly implements weight_median * np.abs(obs - median).

Since _score.py is the shared implementation for all non-numba backends, numpy, jax and torch are all affected. I have verified numpy and jax: both return 0.458, while numba returns 0.338.

The issue is more than a numerical discrepancy: WIS is a weighted average of non-negative terms and mathematically cannot be negative. However, the numpy backend returns negative WIS values when the median is negative. The third example in my reproduction script returns -19.82.

Why existing tests did not catch this: test_weighted_interval_score in tests/test_interval.py calls sr.weighted_interval_score(obs, obs, lower, upper, alpha), passing obs as both the observation and the median. In this case median = obs and |obs - median| = 0, so both formulas evaluate to zero and the faulty term vanishes. As a result, CI passes for all four backends. This also explains why the issue persisted even after #81 and #82 improved backend test coverage.

Fix: change line 39 to WIS += w_median * B.abs(obs - median) (using the absolute value function of the active backend) to match the numba implementation. For tests, it is recommended to use a median value different from obs; otherwise the bug will remain untestable after the fix.

There is also a mismatch between documentation and code: the docstring states the default weights are w_k = 2/alpha_k, but the code uses alpha/2. The code is correct (consistent with the paper), so the docstring needs updating.

I also compared all 97 exported functions one-by-one between numpy and numba backends, and only this function shows inconsistency. I independently verified all closed-form CRPS variants, energy score, variogram score, interval score and pinball score via numerical integration and reference paper formulas, and all are correct. This is therefore an isolated bug.

If this assessment holds, I can open a PR to fix this single line and add a regression test covering cases where median != obs.

AI note: This bug was discovered with assistance from an AI assistant. The reproduction script and validations above were run by the assistant, and I have manually verified every line prior to submission.

Minimal reproducer

import numpy as np, scoringrules as sr

# Bracher et al. (2021), Eq. (3):
#   WIS = (1/(K+1/2)) * ( w0*|y - m| + sum_k w_k * IS_k ),  w0 = 1/2, w_k = alpha_k/2
def wis_reference(y, m, lower, upper, alpha):
    IS = (upper - lower) \
         + 2/alpha*(lower - y)*(y < lower) \
         + 2/alpha*(y - upper)*(y > upper)
    return (0.5*abs(y - m) + np.sum(alpha/2 * IS)) / (len(alpha) + 0.5)

cases = {
    "median = obs (the case the test suite uses)":
        (0.0, 0.0, np.array([-1.2, -0.7]), np.array([1.2, 0.7])),
    "median != obs":
        (1.4, 1.0, np.array([0.5, 0.8]), np.array([1.2, 1.1])),
    "negative location":
        (-100.0, -100.0, np.array([-101.0, -100.5]), np.array([-99.0, -99.5])),
}
alpha = np.array([0.2, 0.5])
for name, (y, m, lo, hi) in cases.items():
    ref = wis_reference(y, m, lo, hi, alpha)
    npv = float(sr.weighted_interval_score(y, m, lo, hi, alpha, backend="numpy"))
    nbv = float(sr.weighted_interval_score(y, m, lo, hi, alpha, backend="numba"))
    print(f"{name}\n  reference={ref:12.6f}  numpy={npv:12.6f}  numba={nbv:12.6f}")

# Output on scoringrules 0.10.0 / numpy 2.4.6 / numba 0.67.0 / Python 3.11:
#
# median = obs (the case the test suite uses)
#   reference=    0.236000  numpy=    0.236000  numba=    0.236000
# median != obs
#   reference=    0.338000  numpy=    0.458000  numba=    0.338000
# negative location
#   reference=    0.180000  numpy=  -19.820000  numba=    0.180000

Backend

numpy, jax, torch

Python version

3.11

scoringrules version

0.10.0