1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Interval calibration: PIT, reliability of the five quantiles, conditional coverage, sharpness and a3split-conformal rescaling of the log-interval (docs/evaluation_protocol.md §5)."""4from __future__ import annotations56from collections.abc import Sequence7from dataclasses import dataclass, field8from typing import Any910import numpy as np11import pandas as pd12from scipy import stats1314from ._common import QUANTILE_COLUMNS, TAUS, to_pandas, validate_predictions15from .metrics import NOMINAL_LEVELS, fit_lognormal_from_quantiles, per_record161718# ----------------------------------------------------------------------------------------------- PIT1920def pit_from_quantiles(records: pd.DataFrame) -> np.ndarray:21 """Approximate probability integral transform F_i(log P_i) from the five quantiles.2223 The predictive CDF is interpolated in *probit* space: z = Φ⁻¹(F) is piecewise linear in log price between the24 knots (log q_τ, Φ⁻¹(τ)), τ ∈ {.05,.25,.5,.75,.95}. This is exact when the forecast is log-normal and reduces to25 the usual piecewise-linear CDF otherwise up to the probit map. Below q05 / above q95 the tails of the normal fitted to26 the five quantiles are used, so PIT ∈ (0, 1). NaN when the record has no interval.27 """28 q = records[QUANTILE_COLUMNS].to_numpy(dtype=float)29 y = np.log(records["price"].to_numpy(dtype=float))30 out = np.full(len(records), np.nan)31 ok = np.all(np.isfinite(q), axis=1) & np.all(q > 0, axis=1)32 if not ok.any():33 return out34 lq = np.log(q[ok])35 yy = y[ok]36 mu, sigma = fit_lognormal_from_quantiles(lq)37 z_knots = stats.norm.ppf(TAUS)38 z = np.empty(len(yy))39 for i in range(len(yy)):40 xs = lq[i]41 if yy[i] < xs[0] or yy[i] > xs[-1] or not np.all(np.diff(xs) > 0):42 s = sigma[i] if sigma[i] > 0 else 1e-643 z[i] = (yy[i] - mu[i]) / s44 else:45 z[i] = np.interp(yy[i], xs, z_knots)46 out[ok] = np.clip(stats.norm.cdf(z), 1e-9, 1 - 1e-9)47 return out484950def pit_histogram(pit: np.ndarray, bins: int = 10) -> pd.DataFrame:51 """Histogram of PIT values with the uniform reference (`expected` = 1/bins)."""52 p = np.asarray(pit, dtype=float)53 p = p[np.isfinite(p)]54 counts, edges = np.histogram(p, bins=bins, range=(0.0, 1.0))55 share = counts / counts.sum() if counts.sum() else np.full(bins, np.nan)56 return pd.DataFrame({"bin_lo": edges[:-1], "bin_hi": edges[1:], "count": counts, "share": share, "expected": 1.0 / bins})575859# ----------------------------------------------------------------------------------------------- reliability6061def reliability(records: pd.DataFrame) -> pd.DataFrame:62 """Empirical frequency P(price ≤ q_τ) for each nominal τ, and interval coverage at each nominal level."""63 q = records[QUANTILE_COLUMNS].to_numpy(dtype=float)64 y = records["price"].to_numpy(dtype=float)65 ok = np.all(np.isfinite(q), axis=1)66 rows = []67 for k, tau in enumerate(TAUS):68 rows.append({"kind": "quantile", "nominal": float(tau), "empirical": float(np.mean(y[ok] <= q[ok, k])) if ok.any() else np.nan, "n": int(ok.sum())})69 for level in NOMINAL_LEVELS:70 col = {0.5: "cov_50", 0.8: "cov_80", 0.9: "cov_90", 0.95: "cov_95"}[level]71 v = records[col].to_numpy(dtype=float) if col in records.columns else np.full(len(records), np.nan)72 rows.append({"kind": "interval", "nominal": level, "empirical": float(np.nanmean(v)) if np.isfinite(v).any() else np.nan, "n": int(np.isfinite(v).sum())})73 df = pd.DataFrame(rows)74 df["gap"] = df["empirical"] - df["nominal"]75 return df767778def coverage_by(records: pd.DataFrame, by: Sequence[str], levels: Sequence[float] = NOMINAL_LEVELS, min_n: int = 30) -> pd.DataFrame:79 """Conditional coverage per subgroup (long table). Cells with n < min_n are suppressed (NaN, suppressed=True)."""80 by = list(by)81 cols = {0.5: "cov_50", 0.8: "cov_80", 0.9: "cov_90", 0.95: "cov_95"}82 rows = []83 for key, g in records.groupby(by, sort=True, dropna=False):84 key = key if isinstance(key, tuple) else (key,)85 base = dict(zip(by, key))86 for level in levels:87 v = g[cols[level]].to_numpy(dtype=float)88 n = int(np.isfinite(v).sum())89 row = dict(base)90 row.update({"nominal": level, "n": n, "coverage": float(np.nanmean(v)) if n >= min_n else np.nan,91 "mean_log_width": float(np.nanmean(g["width_90" if level == 0.9 else "width_50"])) if (n >= min_n and level in (0.5, 0.9)) else np.nan,92 "suppressed": n < min_n})93 rows.append(row)94 return pd.DataFrame(rows)959697def sharpness(records: pd.DataFrame) -> dict[str, float]:98 """Width statistics of the intervals (log scale and relative to the median)."""99 w90, w50 = records["width_90"].to_numpy(dtype=float), records["width_50"].to_numpy(dtype=float)100 q05, q95, q50 = (records[c].to_numpy(dtype=float) for c in ("q05", "q95", "q50"))101 rel = (q95 - q05) / q50102 f = lambda v: float(np.nanmean(v)) if np.isfinite(v).any() else float("nan")103 return {"mean_log_width_90": f(w90), "median_log_width_90": float(np.nanmedian(w90)) if np.isfinite(w90).any() else float("nan"),104 "mean_log_width_50": f(w50), "mean_relative_width_90": f(rel), "n": int(np.isfinite(w90).sum())}105106107# ----------------------------------------------------------------------------------------------- conformal108109@dataclass(frozen=True)110class ConformalScales:111 """Multiplicative factors applied to the log half-widths around log q50 (one per interval)."""112113 scale_90: float114 scale_50: float115 n_valid: int116 level_90: float = 0.90117 level_50: float = 0.50118 achieved_valid_90: float = float("nan")119 achieved_valid_50: float = float("nan")120 note: str = field(default="split conformal, asymmetric scaled residual score, centre = log q50")121122 def as_dict(self) -> dict[str, float]:123 return {"scale_90": self.scale_90, "scale_50": self.scale_50, "n_valid": self.n_valid,124 "achieved_valid_90": self.achieved_valid_90, "achieved_valid_50": self.achieved_valid_50}125126127def _scores(logy: np.ndarray, c: np.ndarray, lo: np.ndarray, hi: np.ndarray) -> np.ndarray:128 eps = 1e-6129 lo, hi = np.maximum(lo, eps), np.maximum(hi, eps)130 return np.maximum((c - logy) / lo, (logy - c) / hi)131132133def _conformal_quantile(s: np.ndarray, level: float) -> float:134 n = len(s)135 k = int(np.ceil((n + 1) * level))136 if k > n:137 return float(np.max(s))138 return float(np.sort(s)[k - 1])139140141def conformalise(valid_preds: Any, valid_truth: Any, level_90: float = 0.90, level_50: float = 0.50) -> ConformalScales:142 """Split-conformal scale factors estimated on the validation table (pure function).143144 Score s_i = max((c_i − y_i)/ℓ_i, (y_i − c_i)/u_i) with c = log q50, ℓ = log q50 − log q05, u = log q95 − log q50145 (and the analogous q25/q75 pair for the 50 % interval). The factor is the ⌈(n+1)(1−α)⌉-th smallest score,146 which guarantees marginal coverage ≥ 1−α on exchangeable test data (Vovk et al. 2005; Romano et al. 2019).147 Factor > 1 widens, < 1 narrows.148 """149 d = per_record(valid_preds, valid_truth)150 d = d[d["has_interval"]]151 if d.empty:152 raise ValueError("no validation records with intervals")153 y = np.log(d["price"].to_numpy(dtype=float))154 lq = np.log(d[QUANTILE_COLUMNS].to_numpy(dtype=float))155 c = lq[:, 2]156 s90 = _scores(y, c, c - lq[:, 0], lq[:, 4] - c)157 s50 = _scores(y, c, c - lq[:, 1], lq[:, 3] - c)158 k90, k50 = _conformal_quantile(s90, level_90), _conformal_quantile(s50, level_50)159 return ConformalScales(scale_90=k90, scale_50=k50, n_valid=len(d),160 achieved_valid_90=float(np.mean(s90 <= k90)), achieved_valid_50=float(np.mean(s50 <= k50)))161162163def apply_conformal(preds: Any, scales: ConformalScales) -> pd.DataFrame:164 """Return a copy of a prediction table with q05/q95 (and q25/q75) rescaled around log q50; q50 and point unchanged."""165 p = validate_predictions(to_pandas(preds))166 lq = np.log(p[QUANTILE_COLUMNS].to_numpy(dtype=float))167 c = lq[:, 2]168 new = lq.copy()169 new[:, 0] = c - scales.scale_90 * np.maximum(c - lq[:, 0], 1e-6)170 new[:, 4] = c + scales.scale_90 * np.maximum(lq[:, 4] - c, 1e-6)171 new[:, 1] = c - scales.scale_50 * np.maximum(c - lq[:, 1], 1e-6)172 new[:, 3] = c + scales.scale_50 * np.maximum(lq[:, 3] - c, 1e-6)173 # keep the point inside [q05, q95] and the quantiles ordered174 new = np.sort(new, axis=1)175 out = p.copy()176 for k, col in enumerate(QUANTILE_COLUMNS):177 out[col] = np.exp(new[:, k])178 lp = np.log(out["point"].to_numpy(dtype=float))179 out["q05"] = np.minimum(out["q05"], np.exp(lp))180 out["q95"] = np.maximum(out["q95"], np.exp(lp))181 return out182