1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Per-record errors and summary metrics (definitions in docs/evaluation_protocol.md §1–§3).34All errors are computed on the sale price P_i (a noisy proxy of value, never "truth" in the economic sense):5 e_i = log(point_i) − log(P_i).6Interval metrics use the five quantiles q05..q95 (CAD). The 50 % and 90 % intervals are read directly from the7quantiles ([q25, q75] and [q05, q95]); the 80 % and 95 % intervals are *interpolated* from a normal distribution8fitted in log space to the five quantiles (least squares on the probit grid). Both raw and fitted versions are reported.9"""10from __future__ import annotations1112from collections.abc import Sequence13from typing import Any1415import numpy as np16import pandas as pd17from scipy import stats1819from ._common import QUANTILE_COLUMNS, TAUS, join_truth2021LOG_125 = float(np.log(1.25))22Z_TAUS = stats.norm.ppf(TAUS) # probit grid of the five quantiles23Z_SUMSQ = float(np.sum(Z_TAUS**2))24NOMINAL_LEVELS = (0.50, 0.80, 0.90, 0.95)2526PRIMARY_METRIC = "mean_abs_log_error"27METRIC_LABELS: dict[str, str] = {28 "n": "N",29 "mean_abs_log_error": r"Mean $|e|$",30 "mae_cad": "MAE (CAD)",31 "mdape": "MdAPE",32 "rmse_log": "RMSE (log)",33 "p90_abs_log_error": r"P90 $|e|$",34 "p95_abs_log_error": r"P95 $|e|$",35 "bias_log": r"Bias $\bar e$",36 "severe_rate": r"$|e|>\log 1.25$",37 "over_rate": "Over +25\\%",38 "under_rate": "Under $-$25\\%",39 "coverage_50": "Cov. 50\\%",40 "coverage_80": "Cov. 80\\%",41 "coverage_90": "Cov. 90\\%",42 "coverage_95": "Cov. 95\\%",43 "coverage_90_normal_fit": "Cov. 90\\% (fit)",44 "mean_log_width_90": "Width 90\\% (log)",45 "mean_log_width_50": "Width 50\\% (log)",46 "pinball_log": "Pinball (log)",47 "crps_log": "CRPS (log)",48 "crps_log_normal_fit": "CRPS (log, fit)",49 "cost_usd_mean": "Cost (USD)",50 "latency_s_mean": "Latency (s)",51 "tokens_mean": "Tokens",52 "cost_usd_total": "Total cost (USD)",53}545556# ----------------------------------------------------------------------------------------------- per-record5758def fit_lognormal_from_quantiles(logq: np.ndarray) -> tuple[np.ndarray, np.ndarray]:59 """Least-squares fit of log q_τ = μ + σ z_τ over the five quantiles (rows = records). Returns (μ, σ)."""60 logq = np.asarray(logq, dtype=float)61 mu = np.full(logq.shape[0], np.nan)62 sigma = np.full(logq.shape[0], np.nan)63 ok = np.all(np.isfinite(logq), axis=1)64 if ok.any():65 mu[ok] = logq[ok].mean(axis=1)66 sigma[ok] = np.maximum((logq[ok] - mu[ok, None]) @ Z_TAUS / Z_SUMSQ, 0.0)67 return mu, sigma686970def _pinball(y: np.ndarray, q: np.ndarray, tau: float) -> np.ndarray:71 d = y - q72 return np.where(d >= 0, tau * d, (tau - 1) * d)737475def per_record(preds: Any, truth: Any) -> pd.DataFrame:76 """Join predictions to truth and add per-record error / interval columns.7778 Added columns: log_price, log_point, e, abs_e, abs_err_cad, ape, severe, over, under,79 cov_50, cov_90 (raw), cov_80, cov_95 (normal fit), cov_90_fit, width_90, width_50 (log), mu_fit, sigma_fit,80 pinball_log (mean over the 5 quantiles), crps_log (quantile approximation), crps_log_fit (closed form of fitted normal).81 Interval columns are NaN when the quantiles are missing.82 """83 d = join_truth(preds, truth)84 y = np.log(d["price"].to_numpy(dtype=float))85 yhat = np.log(d["point"].to_numpy(dtype=float))86 e = yhat - y87 d["log_price"], d["log_point"], d["e"], d["abs_e"] = y, yhat, e, np.abs(e)88 d["abs_err_cad"] = np.abs(d["point"] - d["price"])89 d["ape"] = np.abs(d["point"] / d["price"] - 1.0)90 d["severe"] = (d["abs_e"] > LOG_125).astype(float)91 d["over"] = (d["point"] / d["price"] - 1.0 > 0.25).astype(float)92 d["under"] = (d["point"] / d["price"] - 1.0 < -0.25).astype(float)9394 q = d[QUANTILE_COLUMNS].to_numpy(dtype=float)95 has_q = np.all(np.isfinite(q), axis=1) & np.all(q > 0, axis=1)96 logq = np.where(has_q[:, None], np.log(np.where(q > 0, q, np.nan)), np.nan)97 nan = np.full(len(d), np.nan)9899 def masked(values: np.ndarray) -> np.ndarray:100 return np.where(has_q, values, nan)101102 price = d["price"].to_numpy(dtype=float)103 d["cov_50"] = masked(((q[:, 1] <= price) & (price <= q[:, 3])).astype(float))104 d["cov_90"] = masked(((q[:, 0] <= price) & (price <= q[:, 4])).astype(float))105 mu, sigma = fit_lognormal_from_quantiles(logq)106 d["mu_fit"], d["sigma_fit"] = masked(mu), masked(sigma)107 for level, name in ((0.80, "cov_80"), (0.95, "cov_95"), (0.90, "cov_90_fit")):108 z = stats.norm.ppf(0.5 + level / 2)109 d[name] = masked(((mu - z * sigma <= y) & (y <= mu + z * sigma)).astype(float))110 d["width_90"] = masked(logq[:, 4] - logq[:, 0])111 d["width_50"] = masked(logq[:, 3] - logq[:, 1])112 pb = np.column_stack([_pinball(y, logq[:, k], float(TAUS[k])) for k in range(5)])113 d["pinball_log"] = masked(pb.mean(axis=1))114 # CRPS ≈ 2 ∫ ρ_τ dτ, trapezoidal rule on the τ-grid 0.05..0.95 (the two outer 5 % bands are not observed and are115 # dropped, so this is a lower-biased approximation that is consistent across arms; see protocol §3.4)116 d["crps_log"] = masked(2.0 * np.trapezoid(pb, TAUS, axis=1))117 # closed-form CRPS of the fitted normal (Gneiting & Raftery 2007, eq. 21)118 zz = np.where(sigma > 0, (y - mu) / np.where(sigma > 0, sigma, 1.0), 0.0)119 crps_n = sigma * (zz * (2 * stats.norm.cdf(zz) - 1) + 2 * stats.norm.pdf(zz) - 1 / np.sqrt(np.pi))120 d["crps_log_fit"] = masked(np.where(sigma > 0, crps_n, np.abs(y - mu)))121 d["has_interval"] = has_q122 return d123124125# ----------------------------------------------------------------------------------------------- summaries126127def _nanmean(x: pd.Series) -> float:128 v = x.to_numpy(dtype=float)129 return float(np.nanmean(v)) if np.isfinite(v).any() else float("nan")130131132def summarise_records(d: pd.DataFrame) -> dict[str, float]:133 """Summary metrics from a per-record frame (output of `per_record`)."""134 abs_e = d["abs_e"].to_numpy(dtype=float)135 out: dict[str, float] = {136 "n": float(len(d)),137 "mean_abs_log_error": float(np.mean(abs_e)) if len(d) else float("nan"),138 "mae_cad": float(d["abs_err_cad"].mean()),139 "mdape": float(d["ape"].median()),140 "rmse_log": float(np.sqrt(np.mean(d["e"].to_numpy(dtype=float) ** 2))) if len(d) else float("nan"),141 "p90_abs_log_error": float(np.percentile(abs_e, 90)) if len(d) else float("nan"),142 "p95_abs_log_error": float(np.percentile(abs_e, 95)) if len(d) else float("nan"),143 "bias_log": float(d["e"].mean()),144 "severe_rate": float(d["severe"].mean()),145 "over_rate": float(d["over"].mean()),146 "under_rate": float(d["under"].mean()),147 "coverage_50": _nanmean(d["cov_50"]),148 "coverage_80": _nanmean(d["cov_80"]),149 "coverage_90": _nanmean(d["cov_90"]),150 "coverage_95": _nanmean(d["cov_95"]),151 "coverage_90_normal_fit": _nanmean(d["cov_90_fit"]),152 "mean_log_width_90": _nanmean(d["width_90"]),153 "mean_log_width_50": _nanmean(d["width_50"]),154 "pinball_log": _nanmean(d["pinball_log"]),155 "crps_log": _nanmean(d["crps_log"]),156 "crps_log_normal_fit": _nanmean(d["crps_log_fit"]),157 "n_with_interval": float(d["has_interval"].sum()),158 "cost_usd_mean": _nanmean(d["cost_usd"]),159 "cost_usd_total": float(np.nansum(d["cost_usd"].to_numpy(dtype=float))) if d["cost_usd"].notna().any() else float("nan"),160 "latency_s_mean": _nanmean(d["latency_s"]),161 "tokens_mean": _nanmean(d["tokens"]),162 }163 return out164165166def summarise(preds: Any, truth: Any) -> dict[str, float]:167 """All summary metrics for one prediction table (one arm) joined to truth."""168 return summarise_records(per_record(preds, truth))169170171def summarise_by(preds: Any, truth: Any, by: Sequence[str] = ("arm", "committee", "aggregation")) -> pd.DataFrame:172 """One row of summary metrics per group (default: per arm × committee × aggregation)."""173 d = per_record(preds, truth)174 by = list(by)175 rows = []176 for key, g in d.groupby(by, sort=True, dropna=False):177 key = key if isinstance(key, tuple) else (key,)178 row = dict(zip(by, key))179 row.update(summarise_records(g))180 rows.append(row)181 return pd.DataFrame(rows)182183184# per-record loss columns usable as bootstrap statistics (metric name -> (column, reducer))185RECORD_METRICS: dict[str, tuple[str, str]] = {186 "mean_abs_log_error": ("abs_e", "mean"),187 "mae_cad": ("abs_err_cad", "mean"),188 "mdape": ("ape", "median"),189 "rmse_log": ("e", "rms"),190 "p90_abs_log_error": ("abs_e", "p90"),191 "p95_abs_log_error": ("abs_e", "p95"),192 "bias_log": ("e", "mean"),193 "severe_rate": ("severe", "mean"),194 "over_rate": ("over", "mean"),195 "under_rate": ("under", "mean"),196 "coverage_50": ("cov_50", "nanmean"),197 "coverage_80": ("cov_80", "nanmean"),198 "coverage_90": ("cov_90", "nanmean"),199 "coverage_95": ("cov_95", "nanmean"),200 "mean_log_width_90": ("width_90", "nanmean"),201 "pinball_log": ("pinball_log", "nanmean"),202 "crps_log": ("crps_log", "nanmean"),203}204205206def metric_statistic(name: str):207 """Return a callable DataFrame -> float computing the named summary metric from a per-record frame."""208 col, red = RECORD_METRICS[name]209 reducers = {210 "mean": lambda v: float(np.mean(v)),211 "nanmean": lambda v: float(np.nanmean(v)) if np.isfinite(v).any() else float("nan"),212 "median": lambda v: float(np.median(v)),213 "rms": lambda v: float(np.sqrt(np.mean(v**2))),214 "p90": lambda v: float(np.percentile(v, 90)),215 "p95": lambda v: float(np.percentile(v, 95)),216 }217 f = reducers[red]218 return lambda df: f(df[col].to_numpy(dtype=float))219