spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 39 — Modèles bayésiens et intervalles : exemple intégral sur les 600 ventes du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES,3graine 20261031). Numpy/scipy uniquement (pymc et arviz absents) : tout est implémenté ici.45Contenu :6  * régression linéaire bayésienne conjuguée (a priori normale–inverse-gamma, NIG) : postérieure, prédictive Student-t, a priori7    faiblement informatif vs a priori informatif « emprunté » au marché voisin (autres secteurs) pour un petit échantillon ;8  * modèle hiérarchique à effets de secteur aléatoires, échantillonneur de Gibbs (paramétrisation centrée de Gelfand, Sahu et Carlin9    1995), rétrécissement, expérience « Masson-Angers réduit à 5, 10, 20 ventes » ;10  * diagnostics MCMC : R̂ de Gelman–Rubin (version scindée, BDA3 §11.4), taille d'échantillon effective (autocorrélations, Geyer) ;11  * distribution prédictive du sujet (intervalles de crédibilité 90/95 %) vs intervalle de prédiction fréquentiste du ch. 18 ;12  * critères : WAIC (Watanabe 2010 ; Vehtari, Gelman et Gabry 2017), LOO exact pour le modèle conjugué, LOO par échantillonnage13    préférentiel tronqué pour le modèle hiérarchique ;14  * simulation à partir des lois du générateur (data/synthetic/README.md) pour la récupération des paramètres (tests).15Le modèle de référence (régresseurs) est `traite.cas1.matrice` ; il n'est pas redéfini.16Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch39_bayes.py"""17from __future__ import annotations18import sys, math, pathlib19import numpy as np, pandas as pd20from scipy import stats21sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))22sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))23from traite import cas124from traite import hedonique as h25import ch18_hedonique as c182627ROOT = pathlib.Path(__file__).resolve().parents[3]28GRAINE = 2026103129VRAIS = c18.VRAIS30SECTEURS = cas1.SECTEURS                      # Hull, Gatineau, Aylmer, Buckingham, Masson-Angers31IC18 = {"point": 637_953, "bas95": 539_682, "haut95": 749_087}   # chapitre 18 (traite.cas1.intervalle_prediction)32LN_MEDIANE_HULL = math.log(475_000)           # médiane APCIQ T3 2026 secteur Hull : centre de l'a priori sur la constante3334# =============================================================================================================================35# 1. Régression linéaire bayésienne conjuguée (normale–inverse-gamma)36# =============================================================================================================================3738def prior_faible(colonnes, v0: float = 1e4, a0: float = 1e-3, b0: float = 1e-3) -> dict:39    """A priori faiblement informatif : β | σ² ~ N(m0, σ² V0), σ² ~ IG(a0, b0). m0 = ln(médiane APCIQ de Hull) pour la constante,40    0 ailleurs ; V0 = v0·I (écart-type a priori ≈ 100 σ ≈ 8 en logarithme : pratiquement plat)."""41    k = len(colonnes); m0 = np.zeros(k)42    if "const" in colonnes: m0[list(colonnes).index("const")] = LN_MEDIANE_HULL43    return {"m0": m0, "V0": np.eye(k) * v0, "a0": a0, "b0": b0, "nom": "faiblement informatif"}4445def nig_posterieure(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:46    """Postérieure NIG : V_n = (V0⁻¹ + X'X)⁻¹ ; m_n = V_n (V0⁻¹ m0 + X'y) ; a_n = a0 + n/2 ;47    b_n = b0 + ½ (y'y + m0' V0⁻¹ m0 − m_n' V_n⁻¹ m_n). Marginale de β : Student-t à 2a_n ddl, échelle² = (b_n/a_n) diag(V_n)."""48    X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape49    V0inv = np.linalg.inv(prior["V0"]); m0 = prior["m0"]50    Vn_inv = V0inv + X.T @ X; Vn = np.linalg.inv(Vn_inv)51    mn = Vn @ (V0inv @ m0 + X.T @ y)52    an = prior["a0"] + n / 253    bn = prior["b0"] + 0.5 * (y @ y + m0 @ V0inv @ m0 - mn @ Vn_inv @ mn)54    nu = 2 * an55    sd = np.sqrt(bn / an * np.diag(Vn)) * (math.sqrt(nu / (nu - 2)) if nu > 2 else np.nan)   # écart-type de la marginale t56    echelle = np.sqrt(bn / an * np.diag(Vn))57    tq = stats.t.ppf(0.975, nu)58    return {"m_n": mn, "V_n": Vn, "a_n": an, "b_n": bn, "nu": nu, "sd": sd, "echelle": echelle,59            "ic95_bas": mn - tq * echelle, "ic95_haut": mn + tq * echelle,60            "sigma2_moy": bn / (an - 1), "sigma_moy": math.sqrt(bn / (an - 1)), "sigma2_mode": bn / (an + 1),61            "n": n, "k": k, "X": X, "y": y, "prior": prior}6263def nig_tirages(post: dict, S: int, rng: np.random.Generator) -> dict:64    """Tirages exacts de la postérieure NIG : σ² ~ IG(a_n, b_n) puis β | σ² ~ N(m_n, σ² V_n)."""65    sigma2 = 1.0 / rng.gamma(post["a_n"], 1.0 / post["b_n"], S)66    L = np.linalg.cholesky(post["V_n"])67    z = rng.standard_normal((S, post["k"]))68    beta = post["m_n"] + (z @ L.T) * np.sqrt(sigma2)[:, None]69    return {"beta": beta, "sigma2": sigma2}7071def predictive_t(post: dict, x: np.ndarray, niveaux=(0.90, 0.95)) -> dict:72    """Prédictive a posteriori d'une nouvelle observation : y* ~ t_{2a_n}(x'm_n, (b_n/a_n)(1 + x'V_n x)) ; quantiles retransformés en $."""73    x = np.asarray(x, float); loc = float(x @ post["m_n"]); s2 = post["b_n"] / post["a_n"] * (1 + x @ post["V_n"] @ x)74    out = {"loc": loc, "echelle": math.sqrt(s2), "nu": post["nu"], "mediane": math.exp(loc),75           "part_coefficients": float((x @ post["V_n"] @ x) / (1 + x @ post["V_n"] @ x))}76    for p in niveaux:77        q = stats.t.ppf(0.5 + p / 2, post["nu"]) * math.sqrt(s2)78        out[f"bas{int(p*100)}"] = math.exp(loc - q); out[f"haut{int(p*100)}"] = math.exp(loc + q)79    return out8081def prior_marche_voisin(post_voisin: dict, colonnes_voisin, colonnes_cible, kappa: float, v0_const: float = 1e4) -> dict:82    """A priori informatif pour un petit marché : les prix implicites (pentes) sont centrés sur la postérieure du marché voisin,83    avec une matrice d'échelle V0 = κ·V_n(voisin) (κ = 1 : confiance entière ; κ = 4 : écart-type doublé) ; la constante (niveau de84    prix du marché cible) reste pratiquement plate. a0, b0 : on transmet aussi l'information sur σ² (a_n, b_n du voisin)."""85    cv = list(colonnes_voisin); cc = list(colonnes_cible); k = len(cc)86    m0 = np.zeros(k); V0 = np.zeros((k, k))87    idx = [cv.index(c) for c in cc if c != "const"]; pos = [i for i, c in enumerate(cc) if c != "const"]88    m0[pos] = post_voisin["m_n"][idx]89    V0[np.ix_(pos, pos)] = kappa * post_voisin["V_n"][np.ix_(idx, idx)]90    ic = cc.index("const"); m0[ic] = LN_MEDIANE_HULL; V0[ic, ic] = v0_const91    return {"m0": m0, "V0": V0, "a0": post_voisin["a_n"], "b0": post_voisin["b_n"], "nom": f"marché voisin, κ = {kappa:g}"}9293def loo_exact_nig(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:94    """LOO exact du modèle conjugué : Σ_i log p(y_i | y_{−i}) par n réestimations en forme fermée (prédictive Student-t)."""95    n = len(y); lp = np.empty(n)96    for i in range(n):97        m = np.ones(n, bool); m[i] = False98        p = nig_posterieure(X[m], y[m], prior)99        s2 = p["b_n"] / p["a_n"] * (1 + X[i] @ p["V_n"] @ X[i])100        lp[i] = stats.t.logpdf(y[i], p["nu"], loc=X[i] @ p["m_n"], scale=math.sqrt(s2))101    return {"elpd_loo": float(lp.sum()), "looic": float(-2 * lp.sum()), "pointwise": lp}102103# =============================================================================================================================104# 2. Modèle hiérarchique à effets de secteur aléatoires : échantillonneur de Gibbs105# =============================================================================================================================106107def matrice_hier(df: pd.DataFrame, moyennes: pd.Series | None = None):108    """Régresseurs du modèle hiérarchique : ceux de traite.cas1.matrice sans la constante ni les indicatrices de secteur, centrés109    (paramétrisation centrée : les niveaux α_s portent l'ordonnée à l'origine). Retourne (X centré, moyennes, indices de secteur)."""110    X = cas1.matrice(df).drop(columns=["const"] + [c for c in cas1.matrice(df).columns if c.startswith("sect_")])111    if moyennes is None: moyennes = X.mean()112    Xc = X - moyennes113    z = np.array([SECTEURS.index(s) for s in df["secteur"]])114    return Xc, moyennes, z115116def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5, n_iter: int = 3000, brulage: int = 1000,117                       chaines: int = 4, graine: int = GRAINE, prior_tau: str = "uniforme", tau_fixe: float | None = None,118                       ig_a: float = 1e-3, ig_b: float = 1e-3) -> dict:119    """ln P_i = x_i'β + α_{s(i)} + ε_i, ε ~ N(0, σ²) ; α_s | α0, τ² ~ N(α0, τ²) ; p(β) ∝ 1, p(α0) ∝ 1, p(σ²) ∝ 1/σ² ;120    τ : p(τ) ∝ 1 (uniforme sur τ, BDA3 §5.7 ; « uniforme ») ou τ² ~ IG(ig_a, ig_b) (« ig ») ou τ fixé (« fixe », bayésien empirique).121    Conditionnelles complètes (toutes conjuguées) :122      β | ·  ~ N((X'X)⁻¹X'(y − Zα), σ²(X'X)⁻¹)123      α_s | · ~ N([n_s r̄_s/σ² + α0/τ²] / [n_s/σ² + 1/τ²], 1/[n_s/σ² + 1/τ²]),  r̄_s = moyenne de (y − Xβ) dans le secteur s124      α0 | · ~ N(moyenne des α_s, τ²/J)125      σ² | · ~ IG(n/2, ½ Σ (y − Xβ − Zα)²)126      τ² | · ~ IG((J−1)/2, ½ Σ (α_s − α0)²)   [uniforme sur τ]   ou   IG(a + J/2, b + ½ Σ (α_s − α0)²)   [IG]127    Retourne les tirages après brûlage, de forme (chaînes, n_iter − brulage, ·)."""128    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape129    XtX_inv = np.linalg.inv(X.T @ X); L_beta = np.linalg.cholesky(XtX_inv)130    ns = np.bincount(z, minlength=J).astype(float)131    garder = n_iter - brulage132    out = {"beta": np.empty((chaines, garder, k)), "alpha": np.empty((chaines, garder, J)), "alpha0": np.empty((chaines, garder)),133           "sigma2": np.empty((chaines, garder)), "tau2": np.empty((chaines, garder))}134    for c in range(chaines):135        rng = np.random.default_rng([graine, c, 39])136        # valeurs initiales dispersées (BDA3 : points de départ surdispersés pour que R̂ ait un sens)137        beta = XtX_inv @ X.T @ y + rng.normal(0, 0.05, k)138        alpha = np.full(J, float(np.mean(y))) + rng.normal(0, 0.3, J); alpha0 = float(np.mean(alpha))139        sigma2 = float(np.exp(rng.normal(math.log(0.01), 0.5))); tau2 = float(np.exp(rng.normal(math.log(0.02), 1.0)))140        if tau_fixe is not None: tau2 = tau_fixe ** 2141        for it in range(n_iter):142            # β | α, σ²143            r = y - alpha[z]144            beta = XtX_inv @ X.T @ r + math.sqrt(sigma2) * (L_beta @ rng.standard_normal(k))145            # α_s | β, α0, σ², τ²146            e = y - X @ beta147            somme = np.bincount(z, weights=e, minlength=J)148            prec = ns / sigma2 + 1.0 / tau2149            moy = (somme / sigma2 + alpha0 / tau2) / prec150            alpha = moy + rng.standard_normal(J) / np.sqrt(prec)151            # α0 | α, τ²152            alpha0 = float(alpha.mean() + rng.standard_normal() * math.sqrt(tau2 / J))153            # σ² | β, α154            res = e - alpha[z]155            sigma2 = 1.0 / rng.gamma(n / 2, 2.0 / (res @ res))156            # τ² | α, α0157            ss = float(((alpha - alpha0) ** 2).sum())158            if tau_fixe is None:159                if prior_tau == "uniforme": tau2 = 1.0 / rng.gamma((J - 1) / 2, 2.0 / ss)160                elif prior_tau == "ig": tau2 = 1.0 / rng.gamma(ig_a + J / 2, 1.0 / (ig_b + ss / 2))161                else: raise ValueError(prior_tau)162            if it >= brulage:163                j = it - brulage164                out["beta"][c, j] = beta; out["alpha"][c, j] = alpha; out["alpha0"][c, j] = alpha0165                out["sigma2"][c, j] = sigma2; out["tau2"][c, j] = tau2166    out["meta"] = {"n": n, "k": k, "J": J, "ns": ns, "n_iter": n_iter, "brulage": brulage, "chaines": chaines, "prior_tau": prior_tau,167                   "tau_fixe": tau_fixe}168    return out169170def facteur_retrecissement(n_s: float, sigma2: float, tau2: float) -> float:171    """B_s = n_s τ² / (n_s τ² + σ²) : poids de la moyenne du secteur dans la moyenne conditionnelle de α_s (1 − B_s vers α0)."""172    return n_s * tau2 / (n_s * tau2 + sigma2)173174# ----------------------------------------------------------------------------------------------------------------- diagnostics175def rhat_scinde(ch: np.ndarray) -> float:176    """R̂ scindé (Gelman et Rubin 1992 ; BDA3 §11.4) : chaque chaîne est coupée en deux ; m demi-chaînes de longueur n ;177    B = n·Var(moyennes), W = moyenne des variances ; var⁺ = (n−1)/n W + B/n ; R̂ = sqrt(var⁺/W)."""178    ch = np.asarray(ch, float); c, T = ch.shape; T2 = T // 2179    halves = np.concatenate([ch[:, :T2], ch[:, T2:2 * T2]], axis=0)   # (2c, T2)180    m, n = halves.shape181    moy = halves.mean(axis=1); B = n * moy.var(ddof=1); W = halves.var(axis=1, ddof=1).mean()182    var_plus = (n - 1) / n * W + B / n183    return float(math.sqrt(var_plus / W))184185def ess(ch: np.ndarray) -> float:186    """Taille d'échantillon effective (BDA3 §11.5) : ρ̂_t = 1 − V_t/(2 var⁺) avec V_t le variogramme empirique moyen des chaînes ;187    somme des autocorrélations arrêtée au premier couple (ρ̂_{2t'} + ρ̂_{2t'+1}) négatif (Geyer) ; ESS = m n / (1 + 2 Σ ρ̂_t)."""188    ch = np.asarray(ch, float); m, n = ch.shape189    moy = ch.mean(axis=1); B = n * moy.var(ddof=1) if m > 1 else 0.0; W = ch.var(axis=1, ddof=1).mean()190    var_plus = (n - 1) / n * W + B / n191    rho = np.empty(n)192    for t in range(n):193        if t == 0: Vt = 0.0194        else: Vt = float(((ch[:, t:] - ch[:, :-t]) ** 2).sum() / (m * (n - t)))195        rho[t] = 1.0 - Vt / (2 * var_plus)196    s = 0.0; t = 1197    while t + 1 < n:198        pair = rho[t] + rho[t + 1]199        if pair < 0: break200        s += pair; t += 2201    return float(m * n / (1 + 2 * s))202203def resume_tirages(x: np.ndarray, nom: str = "") -> dict:204    """Résumé d'un paramètre scalaire de forme (chaînes, tirages) : moyenne, écart-type, quantiles, R̂, ESS."""205    flat = x.reshape(-1)206    return {"nom": nom, "moyenne": float(flat.mean()), "sd": float(flat.std(ddof=1)), "mediane": float(np.median(flat)),207            "q025": float(np.quantile(flat, 0.025)), "q975": float(np.quantile(flat, 0.975)),208            "q05": float(np.quantile(flat, 0.05)), "q95": float(np.quantile(flat, 0.95)), "rhat": rhat_scinde(x), "ess": ess(x)}209210def resume_gibbs(g: dict, colonnes) -> pd.DataFrame:211    lignes = []212    for j, c in enumerate(colonnes): lignes.append(resume_tirages(g["beta"][:, :, j], c))213    for j, s in enumerate(SECTEURS): lignes.append(resume_tirages(g["alpha"][:, :, j], f"alpha_{s}"))214    lignes.append(resume_tirages(g["alpha0"], "alpha0"))215    lignes.append(resume_tirages(np.sqrt(g["sigma2"]), "sigma"))216    lignes.append(resume_tirages(np.sqrt(g["tau2"]), "tau"))217    return pd.DataFrame(lignes).set_index("nom")218219def effets_relatifs(g: dict, ref: str = "Hull") -> pd.DataFrame:220    """Effets de secteur relatifs à Hull (α_s − α_Hull), comparables aux coefficients des indicatrices du ch. 18."""221    i0 = SECTEURS.index(ref); lignes = []222    for j, s in enumerate(SECTEURS):223        if j == i0: continue224        lignes.append(resume_tirages(g["alpha"][:, :, j] - g["alpha"][:, :, i0], s))225    return pd.DataFrame(lignes).set_index("nom")226227# ----------------------------------------------------------------------------------------------------------------- critères228def log_vrais_pointwise_hier(g: dict, y: np.ndarray, X: np.ndarray, z: np.ndarray, S: int = 2000, graine: int = GRAINE) -> np.ndarray:229    """Matrice S × n des log-densités log p(y_i | θ_s) pour S tirages (sous-échantillonnés) du modèle hiérarchique."""230    rng = np.random.default_rng(graine)231    beta = g["beta"].reshape(-1, g["beta"].shape[-1]); alpha = g["alpha"].reshape(-1, g["alpha"].shape[-1]); s2 = g["sigma2"].reshape(-1)232    idx = rng.choice(len(s2), size=min(S, len(s2)), replace=False)233    mu = beta[idx] @ X.T + alpha[idx][:, z]234    return stats.norm.logpdf(y[None, :], loc=mu, scale=np.sqrt(s2[idx])[:, None])235236def log_vrais_pointwise_nig(post: dict, S: int = 2000, graine: int = GRAINE) -> np.ndarray:237    rng = np.random.default_rng(graine); t = nig_tirages(post, S, rng)238    mu = t["beta"] @ post["X"].T239    return stats.norm.logpdf(post["y"][None, :], loc=mu, scale=np.sqrt(t["sigma2"])[:, None])240241def waic(lp: np.ndarray) -> dict:242    """WAIC (Watanabe 2010 ; Vehtari, Gelman et Gabry 2017) : lppd = Σ_i log moyenne_s p(y_i|θ_s) ; p_waic = Σ_i Var_s log p(y_i|θ_s) ;243    elpd_waic = lppd − p_waic ; WAIC = −2 elpd_waic."""244    S = lp.shape[0]245    lppd_i = np.logaddexp.reduce(lp, axis=0) - math.log(S)246    p_i = lp.var(axis=0, ddof=1)247    return {"lppd": float(lppd_i.sum()), "p_waic": float(p_i.sum()), "elpd_waic": float((lppd_i - p_i).sum()),248            "waic": float(-2 * (lppd_i - p_i).sum()), "se_elpd": float(math.sqrt(len(lppd_i)) * (lppd_i - p_i).std(ddof=1))}249250def loo_tronque(lp: np.ndarray) -> dict:251    """LOO approximatif par échantillonnage préférentiel tronqué (Ionides 2008 ; forme simplifiée de PSIS, Vehtari et al. 2017) :252    poids w_is ∝ 1/p(y_i|θ_s), tronqués à sqrt(S)·w̄_i ; p(y_i|y_{−i}) ≈ Σ_s w_is p(y_i|θ_s) / Σ_s w_is."""253    S = lp.shape[0]254    lw = -lp; lw -= lw.max(axis=0, keepdims=True)255    w = np.exp(lw); w /= w.sum(axis=0, keepdims=True)256    w = np.minimum(w, math.sqrt(S) / S); w /= w.sum(axis=0, keepdims=True)257    lp_loo = np.log((w * np.exp(lp)).sum(axis=0))258    return {"elpd_loo": float(lp_loo.sum()), "looic": float(-2 * lp_loo.sum()), "pointwise": lp_loo}259260# =============================================================================================================================261# 3. Simulation à partir des lois du générateur (récupération des paramètres)262# =============================================================================================================================263264def simuler_depuis_generateur(v: pd.DataFrame, n: int, graine: int, sigma: float = VRAIS["sigma"], const: float = 13.2) -> pd.DataFrame:265    """Rééchantillonne n lignes d'attributs du jeu (bootstrap des covariables) et engendre ln P selon les lois déclarées266    (README synthétique) : ln P = const + Σ β_vrai x + prime de secteur + ε, ε ~ N(0, σ²). Les attributs sont ceux du jeu ; le prix267    est régénéré. Sert aux tests de récupération des paramètres."""268    rng = np.random.default_rng(graine)269    d = v.sample(n=n, replace=True, random_state=int(rng.integers(1 << 31))).reset_index(drop=True)270    X = cas1.matrice(d)271    beta = pd.Series({c: VRAIS.get(c, 0.0) for c in X.columns}); beta["const"] = const272    lnp = X.values @ beta.values + rng.normal(0, sigma, n)273    d["prix_vente"] = np.round(np.exp(lnp) / 500) * 500274    return d275276# =============================================================================================================================277# 4. Exemple intégral278# =============================================================================================================================279280def sous_echantillons(v: pd.DataFrame, s: dict) -> dict:281    import ch16_selection as c16282    vp = c16.preparer(v.copy(), s)283    micro = vp[(vp["mois"].between(0, 18)) & (vp["dist_km"] <= 3.0)].drop(columns=["mois", "dist_km", "prix_m2", "distance"])284    return {"tout": v, "hull": v[v["secteur"] == "Hull"], "voisin": v[v["secteur"] != "Hull"], "micro": micro}285286def colonnes_actives(X: pd.DataFrame) -> list:287    return [c for c in X.columns if c == "const" or X[c].std() > 0]288289def ajuster_nig(df: pd.DataFrame, prior_fn=None, colonnes=None) -> dict:290    X = cas1.matrice(df); cols = colonnes or colonnes_actives(X); X = X[cols]291    y = np.log(df["prix_vente"].astype(float)).values292    prior = prior_fn(cols) if prior_fn else prior_faible(cols)293    p = nig_posterieure(X.values, y, prior); p["colonnes"] = cols294    return p295296def x_sujet(s: dict, cols, secteur: str | None = None, dist_centre: float | None = None) -> np.ndarray:297    d = dict(s)298    if secteur: d["secteur"] = secteur299    if dist_centre is not None: d["dist_centre_km"] = dist_centre300    Xs = cas1.matrice(cas1.sujet_df(d))301    return Xs[cols].values[0]302303def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:304    """Étapes 1-2 : postérieure NIG sur 600 ventes (a priori faible) ; petit échantillon (40 ventes à 3 km / 18 mois) avec a priori305    faible vs informatif emprunté aux autres secteurs (κ = 1, 4, 16) ; prédictive du sujet."""306    se = sous_echantillons(v, s)307    p600 = ajuster_nig(se["tout"]); cols = p600["colonnes"]308    xs = x_sujet(s, cols)309    tab = pd.DataFrame({"mco": ref["coef"]["coef"].reindex(cols), "se_hc1": ref["coef"]["se_hc1"].reindex(cols),310                        "post_moy": p600["m_n"], "post_sd": p600["sd"], "ic_bas": p600["ic95_bas"], "ic_haut": p600["ic95_haut"],311                        "vrai": [VRAIS.get(c, np.nan) for c in cols]}, index=cols)312    tab["couvre"] = (tab["vrai"] >= tab["ic_bas"]) & (tab["vrai"] <= tab["ic_haut"])313    pred600 = predictive_t(p600, xs)314    # cohérence séquentielle : 2023-2024 puis 2025-2026315    anc = v[v["date_vente"] < "2025-01-01"]; rec = v[v["date_vente"] >= "2025-01-01"]316    p_anc = ajuster_nig(anc, colonnes=cols)317    prior_seq = {"m0": p_anc["m_n"], "V0": p_anc["V_n"], "a0": p_anc["a_n"], "b0": p_anc["b_n"], "nom": "séquentiel"}318    p_seq = ajuster_nig(rec, prior_fn=lambda c: prior_seq, colonnes=cols)319    # petit échantillon : 40 ventes autour du sujet ; a priori emprunté aux autres secteurs320    p_vois = ajuster_nig(se["voisin"])321    micro = se["micro"]; Xm = cas1.matrice(micro); cols_m = colonnes_actives(Xm)322    petits = {}323    for nom, pf in [("faible", None)] + [(f"kappa_{k}", (lambda k: (lambda c: prior_marche_voisin(p_vois, p_vois["colonnes"], c, k)))(k)) for k in (1, 4, 16)]:324        p = ajuster_nig(micro, prior_fn=pf, colonnes=cols_m)325        pr = predictive_t(p, x_sujet(s, cols_m))326        petits[nom] = {"post": p, "pred": pr, "largeur95_pct": (pr["haut95"] - pr["bas95"]) / pr["mediane"] * 100}327    # Hull seul (186), a priori faible328    p_hull = ajuster_nig(se["hull"]); pred_hull = predictive_t(p_hull, x_sujet(s, p_hull["colonnes"]))329    return {"p600": p600, "cols": cols, "xs": xs, "tab600": tab, "pred600": pred600, "p_anc": p_anc, "p_rec_seq": p_seq, "n_anc": len(anc),330            "n_rec": len(rec), "p_vois": p_vois, "micro": micro, "cols_micro": cols_m, "petits": petits, "p_hull": p_hull,331            "pred_hull": pred_hull, "n_micro": len(micro)}332333def exemple_hierarchique(v: pd.DataFrame, s: dict, ref: dict, n_iter: int = 3000, brulage: int = 1000) -> dict:334    """Étapes 3-4 : Gibbs sur les 600 ventes ; effets relatifs à Hull vs effets fixes (ch. 18) ; prédictive du sujet ;335    expérience Masson-Angers réduit à 5, 10, 20 ventes ; sensibilité à l'a priori sur τ."""336    Xc, moy, z = matrice_hier(v); y = np.log(v["prix_vente"].astype(float)).values; cols = list(Xc.columns)337    g = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage)338    res = resume_gibbs(g, cols); rel = effets_relatifs(g)339    fixes = ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], ["coef", "se_hc1"]]; fixes.index = SECTEURS[1:]340    rel["fixe"] = fixes["coef"]; rel["fixe_se"] = fixes["se_hc1"]; rel["vrai"] = [VRAIS[f"sect_{x}"] for x in rel.index]341    # prédictive du sujet (Hull) : ln P* = (x_s − x̄)'β + α_Hull + ε342    xs = (cas1.matrice(cas1.sujet_df(s))[cols].iloc[0] - moy).values343    beta = g["beta"].reshape(-1, len(cols)); alpha = g["alpha"].reshape(-1, 5); s2 = g["sigma2"].reshape(-1)344    rng = np.random.default_rng(GRAINE)345    lnp = beta @ xs + alpha[:, 0] + rng.standard_normal(len(s2)) * np.sqrt(s2)346    P = np.exp(lnp)347    pred = {"tirages": P, "mediane": float(np.median(P)), "moyenne": float(P.mean()), "moyenne_log": float(lnp.mean()),348            "bas90": float(np.quantile(P, 0.05)), "haut90": float(np.quantile(P, 0.95)),349            "bas95": float(np.quantile(P, 0.025)), "haut95": float(np.quantile(P, 0.975)),350            "sd_log": float(lnp.std()), "sd_moyenne_log": float((beta @ xs + alpha[:, 0]).std())}351    # rétrécissement a posteriori : facteur B_s à la médiane de σ², τ²352    s2m = float(np.median(g["sigma2"])); t2m = float(np.median(g["tau2"]))353    B = {sct: facteur_retrecissement(g["meta"]["ns"][j], s2m, t2m) for j, sct in enumerate(SECTEURS)}354    # expérience Masson-Angers : 5, 10, 20, 47 ventes (sous-ensembles emboîtés)355    ma = v[v["secteur"] == "Masson-Angers"]; autres = v[v["secteur"] != "Masson-Angers"]356    perm = np.random.default_rng(GRAINE).permutation(len(ma))357    dist_ma = float(ma["dist_centre_km"].median())358    exp_rows = []359    for n_ma in (5, 10, 20, len(ma)):360        vv = pd.concat([autres, ma.iloc[perm[:n_ma]]]).reset_index(drop=True)361        # effets fixes (MCO, HC1) et prédictive conjuguée (plat) pour un sujet hypothétique à Masson-Angers362        pf = ajuster_nig(vv); i = pf["colonnes"].index("sect_Masson-Angers")363        xma = x_sujet(s, pf["colonnes"], secteur="Masson-Angers", dist_centre=dist_ma)364        pr_f = predictive_t(pf, xma)365        # hiérarchique366        Xv, moy_v, zv = matrice_hier(vv); yv = np.log(vv["prix_vente"].astype(float)).values367        gv = gibbs_hierarchique(yv, Xv.values, zv, n_iter=n_iter, brulage=brulage)368        rr = effets_relatifs(gv).loc["Masson-Angers"]369        xma_c = (cas1.matrice(cas1.sujet_df({**s, "secteur": "Masson-Angers", "dist_centre_km": dist_ma}))[cols].iloc[0] - moy_v).values370        bv = gv["beta"].reshape(-1, len(cols)); av = gv["alpha"].reshape(-1, 5); sv = gv["sigma2"].reshape(-1)371        lnv = bv @ xma_c + av[:, 4] + np.random.default_rng(GRAINE + n_ma).standard_normal(len(sv)) * np.sqrt(sv)372        Pv = np.exp(lnv)373        s2v = float(np.median(gv["sigma2"])); t2v = float(np.median(gv["tau2"]))374        # moyenne brute des résidus du secteur (sans rétrécissement) à la postérieure : r̄ − α0 relatif375        exp_rows.append({"n_ma": n_ma, "fixe": pf["m_n"][i], "fixe_sd": pf["sd"][i], "fixe_bas": pf["ic95_bas"][i], "fixe_haut": pf["ic95_haut"][i],376                         "hier": rr["moyenne"], "hier_sd": rr["sd"], "hier_bas": rr["q025"], "hier_haut": rr["q975"], "hier_rhat": rr["rhat"],377                         "tau_med": math.sqrt(t2v), "B": facteur_retrecissement(n_ma, s2v, t2v),378                         "pred_f_med": pr_f["mediane"], "pred_f_bas": pr_f["bas95"], "pred_f_haut": pr_f["haut95"],379                         "pred_h_med": float(np.median(Pv)), "pred_h_bas": float(np.quantile(Pv, 0.025)), "pred_h_haut": float(np.quantile(Pv, 0.975)),380                         "vrai": VRAIS["sect_Masson-Angers"]})381    experience = pd.DataFrame(exp_rows).set_index("n_ma")382    experience["largeur_f_pct"] = (experience["pred_f_haut"] - experience["pred_f_bas"]) / experience["pred_f_med"] * 100383    experience["largeur_h_pct"] = (experience["pred_h_haut"] - experience["pred_h_bas"]) / experience["pred_h_med"] * 100384    # sensibilité à l'a priori sur τ (600 ventes et n_MA = 5)385    v5 = pd.concat([autres, ma.iloc[perm[:5]]]).reset_index(drop=True)386    X5, _, z5 = matrice_hier(v5); y5 = np.log(v5["prix_vente"].astype(float)).values387    tau_eb = float(np.std(ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], "coef"].tolist() + [0.0], ddof=1))388    sens_rows = []389    for nom, kw in [("uniforme sur τ", {"prior_tau": "uniforme"}), ("IG(0,001 ; 0,001) sur τ²", {"prior_tau": "ig"}),390                    ("IG(1 ; 0,01) sur τ²", {"prior_tau": "ig", "ig_a": 1.0, "ig_b": 0.01}), (f"τ fixé à {tau_eb:.3f} (bayésien empirique)", {"tau_fixe": tau_eb})]:391        g600 = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, **kw)392        g5 = gibbs_hierarchique(y5, X5.values, z5, n_iter=n_iter, brulage=brulage, **kw)393        r600 = effets_relatifs(g600).loc["Masson-Angers"]; r5 = effets_relatifs(g5).loc["Masson-Angers"]394        sens_rows.append({"prior": nom, "tau_med_600": float(np.median(np.sqrt(g600["tau2"]))), "tau_q05_600": float(np.quantile(np.sqrt(g600["tau2"]), 0.05)),395                          "tau_q95_600": float(np.quantile(np.sqrt(g600["tau2"]), 0.95)),396                          "ma_600": r600["moyenne"], "ma_600_sd": r600["sd"], "tau_med_5": float(np.median(np.sqrt(g5["tau2"]))),397                          "ma_5": r5["moyenne"], "ma_5_sd": r5["sd"], "ma_5_bas": r5["q025"], "ma_5_haut": r5["q975"]})398    sensibilite = pd.DataFrame(sens_rows).set_index("prior")399    return {"g": g, "cols": cols, "moyennes": moy, "z": z, "y": y, "Xc": Xc, "resume": res, "relatifs": rel, "pred": pred, "B": B,400            "sigma2_med": s2m, "tau2_med": t2m, "experience": experience, "dist_ma": dist_ma, "sensibilite": sensibilite, "tau_eb": tau_eb,401            "perm_ma": perm}402403def criteres(v: pd.DataFrame, hier: dict, ex: dict) -> pd.DataFrame:404    """Comparaison de trois modèles : pooling complet (sans secteur), effets fixes (indicatrices), hiérarchique."""405    y = hier["y"]406    # pooling complet : matrice de référence sans indicatrices de secteur407    X = cas1.matrice(v); cols_pool = [c for c in X.columns if not c.startswith("sect_")]408    p_pool = nig_posterieure(X[cols_pool].values, y, prior_faible(cols_pool)); p_pool["colonnes"] = cols_pool409    p_fixe = ex["p600"]410    lignes = []411    for nom, lp, loo in [("pooling complet (sans secteur)", log_vrais_pointwise_nig(p_pool), loo_exact_nig(p_pool["X"], y, p_pool["prior"])),412                         ("effets fixes de secteur (ch. 18)", log_vrais_pointwise_nig(p_fixe), loo_exact_nig(p_fixe["X"], y, p_fixe["prior"])),413                         ("effets aléatoires de secteur (hiérarchique)", log_vrais_pointwise_hier(hier["g"], y, hier["Xc"].values, hier["z"]), None)]:414        w = waic(lp); lt = loo_tronque(lp)415        lignes.append({"modele": nom, "lppd": w["lppd"], "p_waic": w["p_waic"], "elpd_waic": w["elpd_waic"], "waic": w["waic"], "se_elpd": w["se_elpd"],416                       "elpd_loo_tronque": lt["elpd_loo"], "elpd_loo_exact": loo["elpd_loo"] if loo else np.nan})417    t = pd.DataFrame(lignes).set_index("modele")418    t["d_elpd_waic"] = t["elpd_waic"] - t["elpd_waic"].max()419    return t420421# =============================================================================================================================422# 5. Exercices (chiffres des solutions)423# =============================================================================================================================424425def exercices(ex: dict, hier: dict) -> dict:426    r = {}427    # ex 1 : mise à jour normale-normale, variance connue (prix au m² d'un marché mince)428    m0, s0, sig, n, xbar = 3300.0, 200.0, 300.0, 4, 3520.0429    prec = 1 / s0 ** 2 + n / sig ** 2; mn = (m0 / s0 ** 2 + n * xbar / sig ** 2) / prec; sn = math.sqrt(1 / prec)430    w = (n / sig ** 2) / prec431    r["ex1"] = {"m_n": mn, "s_n": sn, "poids_donnees": w, "ic95": (mn - 1.96 * sn, mn + 1.96 * sn), "ic90": (mn - 1.645 * sn, mn + 1.645 * sn),432                "n_equivalent_prior": sig ** 2 / s0 ** 2, "valeur_165": mn * 165, "sd_165": sn * 165}433    # ex 2 : lecture d'intervalles — probabilité a posteriori que la valeur dépasse 650 000 $ sous une postérieure log-normale434    mu, sd = ex["pred600"]["loc"], ex["pred600"]["echelle"]435    r["ex2"] = {"p_sup_650": float(1 - stats.t.cdf((math.log(650_000) - mu) / sd, ex["pred600"]["nu"])),436                "p_inf_550": float(stats.t.cdf((math.log(550_000) - mu) / sd, ex["pred600"]["nu"])),437                "q10": math.exp(mu + stats.t.ppf(0.10, ex["pred600"]["nu"]) * sd), "q90": math.exp(mu + stats.t.ppf(0.90, ex["pred600"]["nu"]) * sd)}438    # ex 3 : facteur de rétrécissement439    s2, t2 = 0.083 ** 2, 0.15 ** 2440    r["ex3"] = {n_: facteur_retrecissement(n_, s2, t2) for n_ in (5, 10, 20, 47)}441    r["ex3"]["alpha_5"] = facteur_retrecissement(5, s2, t2) * (-0.30) + (1 - facteur_retrecissement(5, s2, t2)) * (-0.08)442    r["ex3"]["tau_005"] = {n_: facteur_retrecissement(n_, s2, 0.05 ** 2) for n_ in (5, 10, 20, 47)}443    r["ex3"]["sd_cond_5"] = math.sqrt(1 / (5 / s2 + 1 / t2)); r["ex3"]["sd_fixe_5"] = math.sqrt(s2 / 5)444    # ex 4 : prédictive Student-t445    p = ex["p600"]; xs = ex["xs"]; q = float(xs @ p["V_n"] @ xs)446    r["ex4"] = {"a_n": p["a_n"], "b_n": p["b_n"], "nu": p["nu"], "xVx": q, "b_sur_a": p["b_n"] / p["a_n"], "echelle": ex["pred600"]["echelle"],447                "t975": float(stats.t.ppf(0.975, p["nu"])), "z975": 1.959964, "bas95": ex["pred600"]["bas95"], "haut95": ex["pred600"]["haut95"],448                "nu_petit": 20, "t975_20": float(stats.t.ppf(0.975, 20)), "mediane": ex["pred600"]["mediane"]}449    # ex 5 : R̂ sur des moyennes/variances de chaînes données450    moy = np.array([-0.262, -0.258, -0.271, -0.249]); var = np.array([0.0049, 0.0052, 0.0047, 0.0055]); n_ = 1000451    B = n_ * moy.var(ddof=1); W = var.mean(); vp = (n_ - 1) / n_ * W + B / n_452    r["ex5"] = {"B": B, "W": W, "var_plus": vp, "rhat": math.sqrt(vp / W)}453    moy2 = np.array([-0.262, -0.258, -0.371, -0.249]); B2 = n_ * moy2.var(ddof=1); vp2 = (n_ - 1) / n_ * W + B2 / n_454    r["ex5"]["rhat_divergent"] = math.sqrt(vp2 / W); r["ex5"]["B2"] = B2455    # ex 6 : script (κ) — valeurs du tableau des petits échantillons456    r["ex6"] = {k: {"ln_S": v_["post"]["m_n"][ex["cols_micro"].index("ln_S")], "ln_S_sd": v_["post"]["sd"][ex["cols_micro"].index("ln_S")],457                    "mediane": v_["pred"]["mediane"], "largeur95_pct": v_["largeur95_pct"]} for k, v_ in ex["petits"].items()}458    # ex 7 : conditionnelle de α0 avec a priori N(m, s²) — valeurs numériques à τ et α médians459    g = hier["g"]; alpha_bar = float(g["alpha"].reshape(-1, 5).mean()); tau2 = hier["tau2_med"]460    for (m_, s_) in [(13.0, 1.0), (13.0, 0.05)]:461        prec = 5 / tau2 + 1 / s_ ** 2462        r.setdefault("ex7", {})[f"m{m_}_s{s_}"] = {"moy": (5 * alpha_bar / tau2 + m_ / s_ ** 2) / prec, "sd": math.sqrt(1 / prec), "alpha_bar": alpha_bar, "tau": math.sqrt(tau2)}463    return r464465# =============================================================================================================================466# 6. Main467# =============================================================================================================================468469def main(verbose: bool = True, n_iter: int = 3000, brulage: int = 1000) -> dict:470    d = cas1.charger(); v = d["ventes"]; s = d["sujet"]471    ref = cas1.hedonique(v)472    ex = exemple_conjugue(v, s, ref)473    hier = exemple_hierarchique(v, s, ref, n_iter=n_iter, brulage=brulage)474    crit = criteres(v, hier, ex)475    exo = exercices(ex, hier)476    out = {"ventes": v, "sujet": s, "modele": ref, "conjugue": ex, "hier": hier, "criteres": crit, "exercices": exo}477    if verbose: imprimer(out)478    return out479480def _f(x, d=4): return f"{x:.{d}f}"481def _d(x): return f"{x:,.0f} $".replace(",", " ")482483def imprimer(r: dict) -> None:484    ex = r["conjugue"]; hi = r["hier"]485    print("=== CH39 — modèle conjugué, 600 ventes, a priori faiblement informatif ===")486    t = ex["tab600"]487    print(t.round(4).to_string())488    print(f"couverture des vrais paramètres (IC 95 %) : {int(t['couvre'].sum())}/{int(t['vrai'].notna().sum())}")489    p = ex["p600"]; print(f"a_n = {p['a_n']:.1f} ; b_n = {p['b_n']:.4f} ; ν = {p['nu']:.0f} ; E[σ] ≈ {p['sigma_moy']:.4f} ; MCO σ = {r['modele']['sigma']:.4f}")490    pr = ex["pred600"]491    print(f"prédictive du sujet : loc {pr['loc']:.4f}, échelle {pr['echelle']:.5f}, médiane {_d(pr['mediane'])}, part coefficients {pr['part_coefficients']*100:.2f} %")492    print(f"  crédibilité 90 % : [{_d(pr['bas90'])} ; {_d(pr['haut90'])}]   95 % : [{_d(pr['bas95'])} ; {_d(pr['haut95'])}]   (ch. 18 : [{_d(IC18['bas95'])} ; {_d(IC18['haut95'])}])")493    print(f"cohérence séquentielle : n_anc={ex['n_anc']}, n_rec={ex['n_rec']} ; max|m_seq − m_600| = {np.abs(ex['p_rec_seq']['m_n'] - p['m_n']).max():.2e} ; a_n seq {ex['p_rec_seq']['a_n']:.1f} vs {p['a_n']:.1f} ; b_n {ex['p_rec_seq']['b_n']:.6f} vs {p['b_n']:.6f}")494    print(f"\n=== petit échantillon : {ex['n_micro']} ventes à 3 km / 18 mois ; colonnes {ex['cols_micro']}")495    for nom, pe in ex["petits"].items():496        po = pe["post"]; i = ex["cols_micro"].index("ln_S"); iq = ex["cols_micro"].index("qualite"); ig = ex["cols_micro"].index("garage")497        print(f"  {nom:9s}: ln_S {po['m_n'][i]:.3f} ({po['sd'][i]:.3f}) ; qualité {po['m_n'][iq]:.3f} ({po['sd'][iq]:.3f}) ; garage {po['m_n'][ig]:.3f} ({po['sd'][ig]:.3f}) ; "498              f"E[σ] {po['sigma_moy']:.4f} ; sujet médiane {_d(pe['pred']['mediane'])} [{_d(pe['pred']['bas95'])} ; {_d(pe['pred']['haut95'])}] largeur {pe['largeur95_pct']:.1f} %")499    ph = ex["pred_hull"]; print(f"  Hull seul ({ex['p_hull']['n']}) : médiane {_d(ph['mediane'])} [{_d(ph['bas95'])} ; {_d(ph['haut95'])}]")500    print("\n=== modèle hiérarchique (Gibbs) : résumé ===")501    print(hi["resume"].round(4).to_string())502    print("\neffets relatifs à Hull (hiérarchique vs effets fixes ch. 18 vs vrai) :"); print(hi["relatifs"].round(4).to_string())503    print(f"σ² médian {hi['sigma2_med']:.5f} (σ {math.sqrt(hi['sigma2_med']):.4f}) ; τ² médian {hi['tau2_med']:.5f} (τ {math.sqrt(hi['tau2_med']):.4f})")504    print("facteurs de rétrécissement B_s :", {k: round(v_, 4) for k, v_ in hi["B"].items()})505    pd_ = hi["pred"]506    print(f"prédictive hiérarchique du sujet : médiane {_d(pd_['mediane'])}, moyenne {_d(pd_['moyenne'])}, 90 % [{_d(pd_['bas90'])} ; {_d(pd_['haut90'])}], 95 % [{_d(pd_['bas95'])} ; {_d(pd_['haut95'])}] ; sd log {pd_['sd_log']:.4f} (dont moyenne {pd_['sd_moyenne_log']:.4f})")507    print(f"\n=== expérience Masson-Angers (dist. médiane {hi['dist_ma']:.1f} km) ===")508    print(hi["experience"].round(4).to_string())509    print("\n=== sensibilité à l'a priori sur τ ==="); print(hi["sensibilite"].round(4).to_string())510    print("\n=== critères ==="); print(r["criteres"].round(2).to_string())511    print("\n=== exercices ===")512    for k, v_ in r["exercices"].items(): print(k, v_)513514if __name__ == "__main__":515    main()516