spb/qc-election

1# QC Élection Forecast — Plateforme de prévision électorale du Québec 20262# Auteur : Simon-Pierre Boucher3# Contact : contact@spboucher.ai4# https://www.qc-election.com5"""Modèle dynamique d'erreur des sondeurs (§6-§10) — v1 multi-élections.67Sépare explicitement, à partir des 6 élections rejouables (2007-2022) :89  * **erreur d'industrie** `industry_error[élection, parti]` : biais PARTAGÉ10    des sondages finaux (fenêtre 21 j) vs le résultat — 2022 n'était pas un11    accident : chaque élection a son erreur commune. σ_industrie devient une12    quantité ESTIMÉE (espace alr, LOEO : jamais sur l'élection testée) au lieu13    du 0,12 calibré à la main sur 2018/2022 (STATISTICAL_DEBT §industry).14  * **house effects hiérarchiques** : h[maison] global (rétréci par le nombre15    d'élections observées n/(n+k)) + déviation par élection (τ) — un sondeur16    vu 1 fois garde une GRANDE incertitude au lieu d'un simple malus.1718Un plancher conservateur (`floor`) protège contre la sous-estimation quand19l'historique est court (§60 : priors conservateurs si données faibles)."""20from __future__ import annotations2122from datetime import timedelta2324import numpy as np2526from ...config import settings27from ..compositions import alr, close2829FINAL_WINDOW_DAYS = 2130SIGMA_FLOOR = 0.08          # plancher conservateur (alr) — jamais en dessous31SIGMA_CAP = 0.25          # plafond de sécurité32SHRINK_K = 2.0              # h_global rétréci par n/(n+k)333435def _final_avg_shares(polls: list[dict], eday, parties: list[str]) -> np.ndarray | None:36    """Moyenne des sondages des 21 derniers jours (fractions, ordre `parties`)."""37    window = [p for p in polls38              if (eday - timedelta(days=FINAL_WINDOW_DAYS)) <= p["field_end"] <= eday]39    if len(window) < 2:40        return None41    rows = []42    for p in window:43        v = np.array([p["shares"].get(pt, 0.0) for pt in parties])44        if v[:len(parties) - 1].sum() > 60:45            rows.append(close(v / max(v.sum(), 1e-9)))46    return np.mean(rows, axis=0) if len(rows) >= 2 else None474849_CACHE: dict = {}505152def industry_errors() -> dict[int, dict]:53    """Erreur d'industrie par élection : alr(moyenne finale) − alr(résultat).54    (Mise en cache processus — relit les CSV une seule fois.)"""55    if "errs" in _CACHE:56        return _CACHE["errs"]57    from ...ingest.wikipedia_historical import ELECTIONS, load_polls, to_model_polls58    out: dict[int, dict] = {}59    for yr, meta in ELECTIONS.items():60        parties = meta["parties"]61        polls = to_model_polls(load_polls(yr))62        avg = _final_avg_shares(polls, meta["date"], parties)63        if avg is None:64            continue65        actual = close(np.array([meta["result"][p] for p in parties]) / 100.0)66        delta_alr = alr(avg) - alr(actual)67        delta_pp = (avg - actual) * 100.068        # composantes « principales » : le log-ratio d'un parti à 3 % explose69        # pour 1 pp d'erreur — l'estimation de σ ne retient que les composantes70        # dont le parti pèse ≥ 5 % (la référence alr est toujours principale).71        main = [float(delta_alr[i - 1]) for i in range(1, len(parties))72                if actual[i] >= 0.05]73        out[yr] = {"parties": parties,74                   "delta_alr": [round(float(v), 4) for v in delta_alr],75                   "main_delta_alr": [round(v, 4) for v in main],76                   "delta_pp": {p: round(float(delta_pp[i]), 2)77                                for i, p in enumerate(parties)},78                   "rms_alr_main": round(float(np.sqrt(np.mean(79                       np.square(main)))), 4) if main else None}80    _CACHE["errs"] = out81    return out828384def industry_sigma_loeo(test_year: int | None = None) -> float:85    """σ_industrie (alr) estimé sur les élections ≠ test_year : RMS des erreurs86    d'industrie composante par composante, borné [plancher, plafond]."""87    errs = industry_errors()88    vals = [v for yr, e in errs.items() if yr != test_year89            for v in e["main_delta_alr"]]90    if len(vals) < 6:91        return settings.industry_error_sd92    sigma = float(np.sqrt(np.mean(np.square(vals))))93    return float(np.clip(sigma, SIGMA_FLOOR, SIGMA_CAP))949596def hierarchical_house_effects() -> dict:97    """h[maison, parti] : moyenne inter-élections rétrécie n/(n+k), avec la98    dispersion inter-élections τ (les maisons vues 1 fois restent incertaines)."""99    from ...ingest.wikipedia_historical import ELECTIONS, load_polls, to_model_polls100    per: dict[str, dict[str, list[float]]] = {}101    for yr, meta in ELECTIONS.items():102        parties = [p for p in meta["parties"] if p != "AUT"]103        eday = meta["date"]104        for p in to_model_polls(load_polls(yr)):105            if not ((eday - timedelta(days=FINAL_WINDOW_DAYS))106                    <= p["field_end"] <= eday):107                continue108            for pt in parties:109                if pt in p["shares"] and pt in meta["result"]:110                    per.setdefault(p["pollster"], {}).setdefault(111                        f"{yr}:{pt}", []).append(112                        p["shares"][pt] - meta["result"][pt])113    out = {}114    for house, cells in per.items():115        by_party: dict[str, list[float]] = {}116        elections = set()117        for key, devs in cells.items():118            yr, pt = key.split(":")119            elections.add(yr)120            by_party.setdefault(pt, []).append(float(np.mean(devs)))121        n_e = len(elections)122        shrink = n_e / (n_e + SHRINK_K)123        h = {pt: round(float(np.mean(v)) * shrink, 2) for pt, v in by_party.items()}124        tau = round(float(np.mean([np.std(v) for v in by_party.values()125                                   if len(v) >= 2]) if any(126            len(v) >= 2 for v in by_party.values()) else 3.0), 2)127        out[house] = {"n_elections": n_e, "shrink": round(shrink, 2),128                      "house_effect_pp": h, "tau_pp": tau}129    return out130131132def report() -> dict:133    errs = industry_errors()134    return {135        "industry_errors": {str(k): v for k, v in errs.items()},136        "industry_sigma_alr": {137            "all": round(industry_sigma_loeo(None), 4),138            "loeo": {str(yr): round(industry_sigma_loeo(yr), 4) for yr in errs},139            "config_manual": settings.industry_error_sd,140            "floor": SIGMA_FLOOR, "cap": SIGMA_CAP},141        "house_effects_hierarchical": hierarchical_house_effects(),142        "note": ("Erreur d'industrie = biais PARTAGÉ des sondages finaux; "143                 "σ estimé LOEO remplace le 0,12 manuel dans le replay; la "144                 "production l'adopte via use_empirical_industry_error."),145    }146