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

Aperta
  • Qualcosa non va?
#140 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

@lucaluo925 ci sta già lavorando.

Dal 29/9/2026.

  • #141 di @lucaluo925 — aperta

Valutazione

Difficoltà
2/5
Tempo stimato
1-3 ore
Idoneità per principianti
35/100
Tipo di issue
Bug
Chiarezza
Specificata chiaramente
Stato di attività
Ferma
Stack tecnologico
numpy, python
Ambito
data

Direzione di ricerca

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.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

bug
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

Lingua principale
Python
Stelle
108
Fork
15
Merge medio
1h 42m
PR unite (30g)
9

Preparare l'ambiente

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di frazane/scoringrules

Tutte le issue di frazane/scoringrules

Issue simili

Altre issue su Python

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.