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"""Simulateur Monte Carlo de l'élection — erreurs corrélées en couches.67Chaque simulation tire :8  1. le vote national (posterior du forecast) — erreur GLOBALE partagée;9  2. une erreur RÉGIONALE par région de modélisation (partagée par les10     circonscriptions de la région — partial pooling implicite);11  3. une erreur PROPRE à chaque circonscription.1213Le vote d'une circonscription = softmax( alr(référence 2022 locale)14        + swing national (alr) + ε_région + ε_circonscription − pénalités ).15Swing proportionnel en espace log-ratio (pas de swing uniforme naïf).16"""17from __future__ import annotations1819from dataclasses import dataclass2021import numpy as np2223from ..config import settings24from .compositions import alr, close, inv_alr25from .trend import M2627CATEGORY_BOUNDS = [(0.95, "Solide"), (0.80, "Probable"), (0.60, "Favorisé")]282930def category_for(p: float) -> str:31    for bound, label in CATEGORY_BOUNDS:32        if p >= bound:33            return label34    return "Chaudement disputé"353637@dataclass38class SimulationInput:39    x_mean: np.ndarray            # état national (5,)40    P: np.ndarray                 # covariance nationale (5,5)41    baseline_national: np.ndarray  # parts nationales élection de référence (6,)42    district_names: list[str]43    district_baselines: np.ndarray  # (D, 6) parts locales élection de référence44    district_regions: list[str]45    retirement_flags: np.ndarray    # (D,) 1 si le sortant ne se représente pas46    incumbent_party_idx: np.ndarray  # (D,) index parti sortant, -1 si aucun47    majority_seats: int48    n_sims: int = 0495051@dataclass52class SimulationResult:53    n_sims: int54    seats: dict          # {party: {mean, p05..p95, hist, prob_most, prob_majority, prob_minority_win}}55    national_vote: dict  # {party: mean} moyenne des tirages56    districts: list      # [{district, region, favorite, category, win_probs, expected}]57    seat_matrix_summary: dict58    raw_winners: np.ndarray | None = None   # (S,D) indices — pour analyses pivots59    raw_shares: np.ndarray | None = None    # (S,D,K)606162def tipping_points(winners: np.ndarray, shares: np.ndarray,63                   district_names: list[str], majority: int,64                   parties: list[str]) -> dict:65    """Circonscriptions pivots (v3 « 127 », §24) — méthode des simulations :66    dans chaque simulation où le parti de tête atteint la majorité, on classe67    SES circonscriptions gagnées par marge décroissante; la `majority`-ième est68    LA circonscription qui fait basculer la majorité. P(pivot) = fréquence."""69    S, D = winners.shape70    # marge du vainqueur dans chaque (sim, circ.)71    part = np.partition(shares, -2, axis=-1)72    margins = part[..., -1] - part[..., -2]                    # (S,D)73    seat_counts = np.zeros((S, len(parties)), dtype=np.int32)74    for k in range(len(parties)):75        seat_counts[:, k] = (winners == k).sum(axis=1)76    top = seat_counts.argmax(axis=1)77    maj_tips = np.zeros(D)78    n_maj = 079    for s in range(S):80        k = top[s]81        if seat_counts[s, k] < majority:82            continue83        n_maj += 184        won = np.flatnonzero(winners[s] == k)85        order = won[np.argsort(-margins[s, won])]86        maj_tips[order[majority - 1]] += 187    out = []88    if n_maj:89        for d in np.argsort(-maj_tips)[:20]:90            if maj_tips[d] == 0:91                break92            out.append({"district": district_names[d],93                        "prob_majority_tipping": round(float(maj_tips[d] / n_maj), 4)})94    return {"majority_tipping": out,95            "n_sims_with_majority": int(n_maj), "n_sims": int(S)}969798def _scenario_alr_delta(base_shares: np.ndarray, pp_delta: dict[str, float]) -> np.ndarray:99    """Delta en pp sur les parts → delta additif en espace alr."""100    p0 = close(base_shares)101    p1 = p0.copy()102    for party, dpp in pp_delta.items():103        if party in settings.parties:104            p1[settings.parties.index(party)] += dpp / 100.0105    p1 = close(np.clip(p1, 0.001, None))106    return alr(p1) - alr(p0)107108109def run_simulation(inp: SimulationInput, n_sims: int | None = None,110                   scenario: dict | None = None,111                   rng: np.random.Generator | None = None,112                   keep_raw: bool = False) -> SimulationResult:113    """scenario (optionnel) :114      national_pp   : {party: ±pp}                — choc national115      regional_pp   : {region: {party: ±pp}}      — choc régional116      turnout_mult  : {party: facteur}            — participation différentielle117      poll_error_pp : {party: ±pp}                — les sondages se trompent de X118    """119    rng = rng or np.random.default_rng(20261005)120    S = n_sims or settings.n_simulations121    D = len(inp.district_names)122    K = M + 1123    scenario = scenario or {}124125    x_mean = inp.x_mean.copy()126    for key in ("national_pp", "poll_error_pp"):127        if scenario.get(key):128            x_mean = x_mean + _scenario_alr_delta(inv_alr(inp.x_mean), scenario[key])129130    # 1. erreur globale : tirages du vote national131    L = np.linalg.cholesky(inp.P + 1e-10 * np.eye(M))132    x_nat = x_mean + rng.standard_normal((S, M)).astype(np.float64) @ L.T   # (S,5)133    nat_shares = inv_alr(x_nat)                                             # (S,6)134135    # swing national en alr par rapport à l'élection de référence136    delta_nat = x_nat - alr(inp.baseline_national)                          # (S,5)137138    # 2. erreurs régionales corrélées139    regions = sorted(set(inp.district_regions))140    reg_idx = np.array([regions.index(r) for r in inp.district_regions])141    eps_reg = rng.standard_normal((S, len(regions), M)) * settings.regional_error_sd142143    # chocs régionaux de scénario144    reg_scenario = np.zeros((len(regions), M))145    for region, pp in (scenario.get("regional_pp") or {}).items():146        if region in regions:147            mask = np.array([r == region for r in inp.district_regions])148            base_reg = close(inp.district_baselines[mask].mean(axis=0))149            reg_scenario[regions.index(region)] = _scenario_alr_delta(base_reg, pp)150151    # 3. erreurs propres aux circonscriptions152    eps_riding = (rng.standard_normal((S, D, M)) * settings.riding_error_sd).astype(np.float32)153154    base_alr = alr(inp.district_baselines)                                  # (D,5)155156    # pénalité de retraite du sortant (appliquée au parti sortant, espace alr)157    retire_pen = np.zeros((D, M))158    for d in range(D):159        pi = inp.incumbent_party_idx[d]160        if inp.retirement_flags[d] and 0 <= pi < K:161            if pi == 0:  # parti de référence alr : pénaliser = bonifier les autres162                retire_pen[d, :] += settings.incumbent_retirement_penalty163            else:164                retire_pen[d, pi - 1] -= settings.incumbent_retirement_penalty165166    y = (base_alr[None, :, :] + retire_pen[None, :, :]167         + delta_nat[:, None, :] + eps_reg[:, reg_idx, :]168         + reg_scenario[None, reg_idx, :]169         + eps_riding)                                                       # (S,D,5)170    shares = inv_alr(y)                                                      # (S,D,6)171172    if scenario.get("turnout_mult"):173        mult = np.array([scenario["turnout_mult"].get(p, 1.0) for p in settings.parties])174        shares = shares * mult[None, None, :]175        shares = shares / shares.sum(axis=-1, keepdims=True)176177    winners = shares.argmax(axis=-1)                                         # (S,D)178    seat_counts = np.zeros((S, K), dtype=np.int32)179    for k in range(K):180        seat_counts[:, k] = (winners == k).sum(axis=1)181182    most = seat_counts.argmax(axis=1)                                        # (S,)183    seats = {}184    for k, party in enumerate(settings.parties):185        sc = seat_counts[:, k]186        hist = np.bincount(sc, minlength=D + 1)187        seats[party] = {188            "mean": round(float(sc.mean()), 1),189            "p05": int(np.percentile(sc, 5)), "p25": int(np.percentile(sc, 25)),190            "p50": int(np.percentile(sc, 50)), "p75": int(np.percentile(sc, 75)),191            "p95": int(np.percentile(sc, 95)),192            "prob_most": round(float((most == k).mean()), 4),193            "prob_majority": round(float((sc >= inp.majority_seats).mean()), 4),194            "prob_minority_win": round(float(((most == k) & (sc < inp.majority_seats)).mean()), 4),195            "hist": {str(i): int(c) for i, c in enumerate(hist) if c > 0},196        }197198    districts = []199    for d in range(D):200        wp = {settings.parties[k]: round(float((winners[:, d] == k).mean()), 4)201              for k in range(K)}202        exp = {settings.parties[k]: round(float(shares[:, d, k].mean() * 100), 2)203               for k in range(K)}204        lo = {settings.parties[k]: round(float(np.percentile(shares[:, d, k], 2.5) * 100), 2)205              for k in range(K)}206        hi = {settings.parties[k]: round(float(np.percentile(shares[:, d, k], 97.5) * 100), 2)207              for k in range(K)}208        fav = max(wp, key=wp.get)209        districts.append({210            "district": inp.district_names[d], "region": inp.district_regions[d],211            "favorite": fav, "category": category_for(wp[fav]),212            "win_probs": wp, "expected": exp, "lo95": lo, "hi95": hi,213        })214215    nat = {settings.parties[k]: round(float(nat_shares[:, k].mean() * 100), 2)216           for k in range(K)}217    hung = float((seat_counts.max(axis=1) < inp.majority_seats).mean())218    return SimulationResult(219        n_sims=S, seats=seats, national_vote=nat, districts=districts,220        seat_matrix_summary={"prob_no_majority": round(hung, 4),221                             "majority_threshold": inp.majority_seats,222                             "total_seats": D},223        raw_winners=winners if keep_raw else None,224        raw_shares=shares if keep_raw else None)225