1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Aggregation of the committee's final estimates. Everything is computed in log space and converted back.34Point rules (configs/treatments.yaml ``aggregation_rules``):5 mean, median, trimmed_mean_20 (20 % trimmed from each tail, scipy convention), inverse_variance (weights6 = 1 / var of each role's validation log errors), stacking (non-negative least squares on log points,7 weights normalised to sum 1), robust_bayes (Student-t, 4 df, location by IRLS; scale from validation).8Interval rules: median of each quantile across roles (default) and the equal-weight mixture quantile.910Weights are estimated ONCE on validation transcripts (``fit_weights``), saved to11``results/aggregation_weights.json`` and applied unchanged to test/holdout. When no weights are supplied,12inverse_variance and stacking fall back to equal weights and robust_bayes uses a MAD scale from the13committee itself; the runner records ``weights_source`` in the manifest so that this is never silent.14"""15from __future__ import annotations1617import json18import math19from pathlib import Path20from typing import Any2122import numpy as np23from scipy.optimize import nnls24from scipy.stats import trim_mean2526from ..config import RESULTS27from ..schemas import DebateTranscript2829RULES: tuple[str, ...] = ("mean", "median", "trimmed_mean_20", "inverse_variance", "stacking", "robust_bayes")30QUANTILE_LEVELS: tuple[float, ...] = (0.05, 0.25, 0.50, 0.75, 0.95)31T_DF = 432WEIGHTS_PATH = RESULTS / "aggregation_weights.json"333435# ---------------------------------------------------------------------------------------------------36# point rules37# ---------------------------------------------------------------------------------------------------38def _logs(points: dict[str, float]) -> tuple[list[str], np.ndarray]:39 roles = [r for r, v in points.items() if v is not None and v > 0]40 return roles, np.array([math.log(points[r]) for r in roles], dtype=float)414243def _normalise(w: dict[str, float], roles: list[str]) -> np.ndarray:44 arr = np.array([max(float(w.get(r, 0.0)), 0.0) for r in roles], dtype=float)45 if arr.sum() <= 0:46 arr = np.ones(len(roles))47 return arr / arr.sum()484950def student_t_location(x: np.ndarray, scale: float, df: int = T_DF, iters: int = 100, tol: float = 1e-10) -> float:51 """IRLS location estimate under a Student-t likelihood with known scale (robust Bayes posterior mode under a flat prior)."""52 if len(x) == 0:53 return float("nan")54 mu = float(np.median(x))55 scale = max(float(scale), 1e-6)56 for _ in range(iters):57 z2 = ((x - mu) / scale) ** 258 w = (df + 1.0) / (df + z2)59 new = float(np.sum(w * x) / np.sum(w))60 if abs(new - mu) < tol:61 mu = new62 break63 mu = new64 return mu656667def mad_scale(x: np.ndarray, floor: float = 0.02) -> float:68 if len(x) < 2:69 return floor70 return max(1.4826 * float(np.median(np.abs(x - np.median(x)))), floor)717273def aggregate_points(points: dict[str, float], weights: dict[str, Any] | None = None) -> dict[str, float]:74 """rule -> aggregated value (CAD). Empty dict when no valid point."""75 roles, lx = _logs(points)76 if len(roles) == 0:77 return {}78 out: dict[str, float] = {}79 out["mean"] = math.exp(float(lx.mean()))80 out["median"] = math.exp(float(np.median(lx)))81 out["trimmed_mean_20"] = math.exp(float(trim_mean(lx, 0.2))) if len(lx) >= 3 else out["mean"]82 w_iv = _normalise((weights or {}).get("inverse_variance", {}), roles)83 out["inverse_variance"] = math.exp(float(np.dot(w_iv, lx)))84 w_st = _normalise((weights or {}).get("stacking", {}), roles)85 out["stacking"] = math.exp(float(np.dot(w_st, lx)))86 rb = (weights or {}).get("robust_bayes", {}) or {}87 scale = float(rb.get("scale")) if rb.get("scale") else mad_scale(lx)88 out["robust_bayes"] = math.exp(student_t_location(lx, scale=scale, df=int(rb.get("df", T_DF))))89 return out909192def aggregate(transcript: DebateTranscript, weights: dict[str, Any] | None = None) -> dict[str, float]:93 return aggregate_points(transcript.final_estimates, weights)949596# ---------------------------------------------------------------------------------------------------97# interval rules98# ---------------------------------------------------------------------------------------------------99def median_of_quantiles(quantiles: dict[str, list[float]]) -> list[float]:100 """Element-wise median across roles of [q05, q25, q50, q75, q95] (in log space)."""101 rows = [q for q in quantiles.values() if q and len(q) == 5 and all(v > 0 for v in q)]102 if not rows:103 return []104 arr = np.log(np.array(rows, dtype=float))105 return [float(v) for v in np.exp(np.median(arr, axis=0))]106107108def mixture_quantiles(quantiles: dict[str, list[float]], levels: tuple[float, ...] = QUANTILE_LEVELS, grid: int = 2001) -> list[float]:109 """Quantiles of the equal-weight mixture of the roles' distributions. Each role's distribution is110 represented by linear interpolation of its quantile function in log space (tails extended linearly111 to 0.5 % / 99.5 %)."""112 rows = [np.log(np.array(q, dtype=float)) for q in quantiles.values() if q and len(q) == 5 and all(v > 0 for v in q)]113 if not rows:114 return []115 p = np.array(QUANTILE_LEVELS)116 ext_p = np.concatenate([[0.005], p, [0.995]])117 cdfs = []118 for r in rows:119 lo = r[0] - (r[1] - r[0]) * (0.05 - 0.005) / 0.20 if r[1] > r[0] else r[0] - 0.05120 hi = r[4] + (r[4] - r[3]) * (0.995 - 0.95) / 0.20 if r[4] > r[3] else r[4] + 0.05121 cdfs.append((np.concatenate([[lo], r, [hi]]), ext_p))122 xmin = min(c[0][0] for c in cdfs)123 xmax = max(c[0][-1] for c in cdfs)124 xs = np.linspace(xmin, xmax, grid)125 mix = np.zeros_like(xs)126 for xq, pq in cdfs:127 xq_m = np.maximum.accumulate(xq) # guard against flat/non-monotone quantiles128 mix += np.interp(xs, xq_m, pq, left=0.0, right=1.0)129 mix /= len(cdfs)130 out = []131 for lv in levels:132 idx = int(np.searchsorted(mix, lv, side="left"))133 idx = min(max(idx, 0), grid - 1)134 out.append(float(math.exp(xs[idx])))135 return out136137138def aggregate_intervals(transcript_or_quantiles: DebateTranscript | dict[str, list[float]]) -> dict[str, list[float]]:139 """rule -> [q05, q25, q50, q75, q95]."""140 q = transcript_or_quantiles.final_quantiles if isinstance(transcript_or_quantiles, DebateTranscript) else transcript_or_quantiles141 out: dict[str, list[float]] = {}142 m = median_of_quantiles(q)143 if m:144 out["median_of_quantiles"] = m145 mx = mixture_quantiles(q)146 if mx:147 out["mixture"] = mx148 return out149150151def interval_90(iv: dict[str, list[float]]) -> dict[str, list[float]]:152 """DebateTranscript.aggregated_interval contract: rule -> [lo90, hi90]."""153 return {k: [v[0], v[4]] for k, v in iv.items() if len(v) == 5}154155156# ---------------------------------------------------------------------------------------------------157# weights estimated on validation158# ---------------------------------------------------------------------------------------------------159def fit_weights(valid_transcripts: list[DebateTranscript], truth: dict[str, float], *, min_obs: int = 3, df: int = T_DF) -> dict[str, Any]:160 """Estimate aggregation weights on validation transcripts (truth = realised sale price by tx_id).161162 inverse_variance: w_r ∝ 1 / Var(log point_r − log P) over the properties where role r answered.163 stacking: NNLS of log P on the roles' log points (properties where every role answered; missing164 roles imputed by the row median), weights normalised to sum 1.165 robust_bayes: scale = robust pooled SD (1.4826·MAD) of the role log errors; df fixed at 4.166 """167 roles: list[str] = []168 for t in valid_transcripts:169 for r in t.final_estimates:170 if r not in roles:171 roles.append(r)172 errs: dict[str, list[float]] = {r: [] for r in roles}173 rows: list[list[float]] = []174 ys: list[float] = []175 for t in valid_transcripts:176 p = truth.get(t.tx_id)177 if p is None or p <= 0 or not t.final_estimates:178 continue179 lp = math.log(p)180 row = []181 for r in roles:182 v = t.final_estimates.get(r)183 if v is not None and v > 0:184 errs[r].append(math.log(v) - lp)185 row.append(math.log(v))186 else:187 row.append(float("nan"))188 if any(not math.isnan(x) for x in row):189 med = float(np.nanmedian(row))190 rows.append([med if math.isnan(x) else x for x in row])191 ys.append(lp)192 pooled = [e for r in roles for e in errs[r]]193 pooled_var = float(np.var(pooled, ddof=1)) if len(pooled) > 1 else 0.05194 inv_var: dict[str, float] = {}195 for r in roles:196 e = errs[r]197 var = float(np.var(e, ddof=1)) if len(e) >= min_obs else pooled_var198 inv_var[r] = 1.0 / max(var, 1e-6)199 s = sum(inv_var.values())200 inv_var = {r: v / s for r, v in inv_var.items()} if s > 0 else {r: 1.0 / len(roles) for r in roles}201 stacking: dict[str, float]202 if len(rows) >= max(len(roles), 2):203 X = np.array(rows, dtype=float)204 y = np.array(ys, dtype=float)205 w, _ = nnls(X, y)206 if w.sum() <= 0:207 w = np.ones(len(roles))208 w = w / w.sum()209 stacking = {r: float(v) for r, v in zip(roles, w)}210 else:211 stacking = {r: 1.0 / len(roles) for r in roles} if roles else {}212 scale = mad_scale(np.array(pooled)) if pooled else 0.15213 return {214 "roles": roles,215 "n_transcripts": len(rows),216 "inverse_variance": inv_var,217 "stacking": stacking,218 "robust_bayes": {"df": df, "scale": float(scale)},219 "role_error_sd": {r: (float(np.std(errs[r], ddof=1)) if len(errs[r]) > 1 else None) for r in roles},220 "role_n": {r: len(errs[r]) for r in roles},221 }222223224def save_weights(weights: dict[str, Any], path: Path | None = None) -> Path:225 path = path or WEIGHTS_PATH226 path.parent.mkdir(parents=True, exist_ok=True)227 path.write_text(json.dumps(weights, indent=1, ensure_ascii=False))228 return path229230231def load_weights(path: Path | None = None) -> dict[str, Any] | None:232 path = path or WEIGHTS_PATH233 if not path.exists():234 return None235 return json.loads(path.read_text())236