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 bayésien des intentions de vote (state space / Kalman).67État latent x_t ∈ R^5 : log-ratios des parts de vote (référence = CAQ).8Dynamique : marche aléatoire gaussienne x_t = x_{t-1} + w_t, w_t ~ N(0, q·Δt·I).9Observations : chaque sondage est une observation bruitée y = alr(parts ajustées)10avec covariance issue de la multinomiale (méthode delta) + erreur excédentaire.1112La variance quotidienne q est estimée par maximum de vraisemblance (recherche13sur grille sur la log-vraisemblance des innovations), bornée par configuration.14"""15from __future__ import annotations1617from dataclasses import dataclass18from datetime import date, timedelta1920import numpy as np2122from ..config import settings23from .compositions import alr, alr_obs_cov, close, inv_alr24from .house_effects import PollsterProfile, adjust_shares, profile_for25from .weights import effective_n2627M = len(settings.parties) - 1  # dimension de l'état (K-1) — élection cible2829# Le cœur est paramétré par la liste de partis : le replay historique rejoue30# des ères différentes (ADQ 2007-2008, pas de PCQ avant 2022…) avec le même code.313233def shares_vector(shares: dict[str, float],34                  parties: list[str] | None = None) -> np.ndarray:35    """dict {parti: %} → vecteur ordonné (fractions, somme 1)."""36    ps = parties or settings.parties37    v = np.array([shares.get(p, 0.0) for p in ps], dtype=float)38    if "AUT" in ps and v[ps.index("AUT")] == 0.0:39        v[ps.index("AUT")] = max(0.5, 100.0 - v.sum())40    return close(v / 100.0)414243@dataclass44class PreparedPoll:45    field_end: date46    pollster: str47    y: np.ndarray          # observation alr (5,)48    R: np.ndarray          # covariance d'observation (5,5)49    p_adj: np.ndarray      # parts ajustées (6,)50    n_eff: float51    multiplier: float525354@dataclass55class TrendResult:56    dates: list[date]57    share_mean: np.ndarray   # (T, 6) parts lissées58    share_lo: np.ndarray     # (T, 6) borne 2.5 %59    share_hi: np.ndarray     # (T, 6) borne 97.5 %60    x: np.ndarray            # état filtré à as_of (5,)61    P: np.ndarray            # covariance filtrée (5,5)62    q: float                 # variance quotidienne retenue63    loglik: float64    n_polls: int656667def prepare_polls(polls: list[dict], profiles: dict[str, PollsterProfile],68                  parties: list[str] | None = None) -> list[PreparedPoll]:69    ps = parties or settings.parties70    min_main = max(3, len(ps) - 2)71    out = []72    for poll in polls:73        prof = profile_for(profiles, poll["pollster"])74        shares = {p: v for p, v in poll["shares"].items() if p in ps}75        main = {p: v for p, v in shares.items() if p != "AUT"}76        if len(main) < min_main:77            continue78        adj = adjust_shares(shares, prof)79        p_adj = shares_vector(adj, ps)80        n_eff = effective_n(poll.get("sample_size"), poll.get("mode", "unknown"))81        # la cote du sondeur module la précision effective (jamais > ±40 %)82        n_scaled = n_eff * prof.weight_multiplier83        # maison sans historique : erreur excédentaire accrue; les observations84        # spéciales (ex. partielles) fixent leur propre erreur excédentaire85        excess = poll.get("excess_pp") or (86            settings.excess_poll_sd_pp * (1.0 if prof.mae_pp is not None else 1.25))87        R = alr_obs_cov(p_adj, n_scaled, excess_pp=excess)88        out.append(PreparedPoll(89            field_end=poll["field_end"], pollster=poll["pollster"], y=alr(p_adj),90            R=R, p_adj=p_adj, n_eff=n_eff, multiplier=prof.weight_multiplier))91    out.sort(key=lambda p: p.field_end)92    return out939495def _filter(prepared: list[PreparedPoll], grid: list[date], q: float,96            x0: np.ndarray, P0: np.ndarray, smooth: bool = False):97    """Filtre de Kalman sur grille quotidienne; lisseur RTS optionnel."""98    T = len(grid)99    Mx = len(x0)100    by_day: dict[date, list[PreparedPoll]] = {}101    for p in prepared:102        by_day.setdefault(p.field_end, []).append(p)103104    x, P = x0.copy(), P0.copy()105    loglik = 0.0106    xs_pred = np.zeros((T, Mx)); Ps_pred = np.zeros((T, Mx, Mx))107    xs_filt = np.zeros((T, Mx)); Ps_filt = np.zeros((T, Mx, Mx))108    I = np.eye(Mx)109    for t, d in enumerate(grid):110        if t > 0:111            P = P + q * I112        xs_pred[t], Ps_pred[t] = x, P113        for obs in by_day.get(d, []):114            S = P + obs.R115            Sinv = np.linalg.inv(S)116            innov = obs.y - x117            K = P @ Sinv118            x = x + K @ innov119            P = (I - K) @ P120            P = 0.5 * (P + P.T)121            sign, logdet = np.linalg.slogdet(S)122            loglik += -0.5 * (Mx * np.log(2 * np.pi) + logdet + innov @ Sinv @ innov)123        xs_filt[t], Ps_filt[t] = x, P124125    if not smooth:126        return x, P, loglik, None, None127    # Lisseur de Rauch–Tung–Striebel128    xs_s = xs_filt.copy(); Ps_s = Ps_filt.copy()129    for t in range(T - 2, -1, -1):130        G = Ps_filt[t] @ np.linalg.inv(Ps_pred[t + 1])131        xs_s[t] = xs_filt[t] + G @ (xs_s[t + 1] - xs_pred[t + 1])132        Ps_s[t] = Ps_filt[t] + G @ (Ps_s[t + 1] - Ps_pred[t + 1]) @ G.T133    return x, P, loglik, xs_s, Ps_s134135136def fit_trend(polls: list[dict], profiles: dict[str, PollsterProfile],137              as_of: date, rng: np.random.Generator | None = None,138              parties: list[str] | None = None) -> TrendResult | None:139    ps = parties or settings.parties140    Mx = len(ps) - 1141    prepared = [p for p in prepare_polls(polls, profiles, ps) if p.field_end <= as_of]142    if len(prepared) < 3:143        return None144    rng = rng or np.random.default_rng(20261005)145    start = prepared[0].field_end146    grid = [start + timedelta(days=i) for i in range((as_of - start).days + 1)]147148    x0 = prepared[0].y.copy()149    P0 = np.eye(Mx) * 0.09  # a priori large (~±6 pp)150151    # Estimation de q par vraisemblance profilée (grille log-uniforme)152    qs = np.geomspace(settings.rw_daily_sd_min ** 2, settings.rw_daily_sd_max ** 2, 12)153    lls = []154    for q in qs:155        _, _, ll, _, _ = _filter(prepared, grid, float(q), x0, P0)156        lls.append(ll)157    q_best = float(qs[int(np.argmax(lls))])158159    x, P, ll, xs_s, Ps_s = _filter(prepared, grid, q_best, x0, P0, smooth=True)160161    # Parts lissées + intervalle crédible 95 % par échantillonnage162    T = len(grid)163    share_mean = np.zeros((T, Mx + 1)); share_lo = np.zeros((T, Mx + 1)); share_hi = np.zeros((T, Mx + 1))164    n_draw = 400165    for t in range(T):166        L = np.linalg.cholesky(Ps_s[t] + 1e-10 * np.eye(Mx))167        draws = xs_s[t] + rng.standard_normal((n_draw, Mx)) @ L.T168        s = inv_alr(draws)169        share_mean[t] = inv_alr(xs_s[t])170        share_lo[t] = np.percentile(s, 2.5, axis=0)171        share_hi[t] = np.percentile(s, 97.5, axis=0)172173    return TrendResult(dates=grid, share_mean=share_mean, share_lo=share_lo,174                       share_hi=share_hi, x=x, P=P, q=q_best, loglik=float(ll),175                       n_polls=len(prepared))176