Proper scoring rules, calibration and the Murphy decomposition

"""M13 · Proper scoring rules, calibration and the Murphy decomposition. Read notes §13.7.

Binary events: forecasts p_i in [0, 1], outcomes y_i in {0, 1}, i = 1..N.

Run:  cd exercises && .venv/bin/python -m pytest modules/M13-forecasting/tests/test_m13_scoring.py
"""

from __future__ import annotations

import numpy as np
import pandas as pd


def brier(p, y) -> float:
    """Brier score BS = (1/N) sum_i (p_i - y_i)^2 (eq. 13.32). Lower is better; 0.25 for always saying 0.5."""
    raise NotImplementedError


def log_loss(p, y, eps: float = 1e-15) -> float:
    """Log score as a loss: -(1/N) sum_i [y_i ln p_i + (1 - y_i) ln(1 - p_i)].

    Natural log. Clip p to [eps, 1 - eps] first so a confident miss is large but finite.
    """
    raise NotImplementedError


def brier_skill_score(p, y, p_ref) -> float:
    """BSS = 1 - BS(p) / BS(p_ref). p_ref may be a scalar (0.5, or the base rate) or an array."""
    raise NotImplementedError


def expected_score(q: float, p, rule: str = "brier"):
    """Expected penalty, under true probability q, of reporting p (vectorised in p).

    rule="brier":  q (1 - p)^2 + (1 - q) p^2
    rule="log":    -q ln p - (1 - q) ln(1 - p)
    rule="linear": q (1 - p) + (1 - q) p        (the improper absolute-error score)
    Raise ValueError for any other rule.
    """
    raise NotImplementedError


def calibration_table(p, y, bins=10) -> pd.DataFrame:
    """Reliability table with one row per NON-EMPTY bin and columns n, mean_p, freq.

    bins=None: one bin per distinct forecast value.
    bins=K (int): equal-width bins on [0, 1]; bin k holds p in [k/K, (k+1)/K), with p = 1
    put in the last bin (np.digitize on the interior edges does exactly this).
    Index = bin id (0..K-1, or 0..#distinct-1 for bins=None).
    """
    raise NotImplementedError


def murphy_decomposition(p, y, bins=None) -> dict:
    """Murphy (1973) decomposition, eq. 13.36, using calibration_table(p, y, bins):

        REL = (1/N) sum_k n_k (pbar_k - ybar_k)^2
        RES = (1/N) sum_k n_k (ybar_k - ybar)^2
        UNC = ybar (1 - ybar)
    Return dict(reliability, resolution, uncertainty, brier). BS = REL - RES + UNC holds
    exactly when every forecast in a bin is the same number (always true for bins=None).
    """
    raise NotImplementedError