spb/traite-evaluation-qc
Public
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_reference(colonnes) -> dict:39 """A priori de référence (impropre, Jeffreys) : p(β, σ²) ∝ 1/σ². Cas limite de la NIG avec V0⁻¹ = 0, a0 = −k/2, b0 = 0 : la40 postérieure est NIG(β̂, (X'X)⁻¹, (n−k)/2, SCR/2) et la prédictive coïncide avec l'intervalle de prédiction fréquentiste classique."""41 k = len(colonnes)42 return {"m0": np.zeros(k), "V0inv": np.zeros((k, k)), "a0": -k / 2, "b0": 0.0, "nom": "référence (impropre)"}4344def prior_faible(colonnes, v0: float = 1e4, a0: float = 1e-3, b0: float = 1e-3) -> dict:45 """A priori faiblement informatif PROPRE : β | σ² ~ N(m0, σ² V0), σ² ~ IG(a0, b0). m0 = ln(médiane APCIQ de Hull) pour la constante,46 0 ailleurs ; V0 = v0·I (écart-type a priori ≈ 100 σ ≈ 8 en logarithme : pratiquement plat). Prix d'un a priori propre et vague :47 a_n = a0 + n/2 au lieu de (n−k)/2, donc E[σ²] ≈ SCR/(n−2) légèrement inférieur à s² = SCR/(n−k)."""48 k = len(colonnes); m0 = np.zeros(k)49 if "const" in colonnes: m0[list(colonnes).index("const")] = LN_MEDIANE_HULL50 return {"m0": m0, "V0inv": np.eye(k) / v0, "a0": a0, "b0": b0, "nom": "faiblement informatif (propre)"}5152def nig_posterieure(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:53 """Postérieure NIG : V_n = (V0⁻¹ + X'X)⁻¹ ; m_n = V_n (V0⁻¹ m0 + X'y) ; a_n = a0 + n/2 ; (prior : m0, V0inv, a0, b0)54 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)."""55 X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape56 V0inv = prior["V0inv"]; m0 = prior["m0"]57 Vn_inv = V0inv + X.T @ X; Vn = np.linalg.inv(Vn_inv)58 mn = Vn @ (V0inv @ m0 + X.T @ y)59 an = prior["a0"] + n / 260 bn = prior["b0"] + 0.5 * (y @ y + m0 @ V0inv @ m0 - mn @ Vn_inv @ mn)61 nu = 2 * an62 sd = np.sqrt(bn / an * np.diag(Vn)) * (math.sqrt(nu / (nu - 2)) if nu > 2 else np.nan) # écart-type de la marginale t63 echelle = np.sqrt(bn / an * np.diag(Vn))64 tq = stats.t.ppf(0.975, nu)65 return {"m_n": mn, "V_n": Vn, "a_n": an, "b_n": bn, "nu": nu, "sd": sd, "echelle": echelle,66 "ic95_bas": mn - tq * echelle, "ic95_haut": mn + tq * echelle,67 "sigma2_moy": bn / (an - 1), "sigma_moy": math.sqrt(bn / (an - 1)), "sigma2_mode": bn / (an + 1),68 "n": n, "k": k, "X": X, "y": y, "prior": prior}6970def nig_tirages(post: dict, S: int, rng: np.random.Generator) -> dict:71 """Tirages exacts de la postérieure NIG : σ² ~ IG(a_n, b_n) puis β | σ² ~ N(m_n, σ² V_n)."""72 sigma2 = 1.0 / rng.gamma(post["a_n"], 1.0 / post["b_n"], S)73 L = np.linalg.cholesky(post["V_n"])74 z = rng.standard_normal((S, post["k"]))75 beta = post["m_n"] + (z @ L.T) * np.sqrt(sigma2)[:, None]76 return {"beta": beta, "sigma2": sigma2}7778def predictive_t(post: dict, x: np.ndarray, niveaux=(0.90, 0.95)) -> dict:79 """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 $."""80 x = np.asarray(x, float); loc = float(x @ post["m_n"]); s2 = post["b_n"] / post["a_n"] * (1 + x @ post["V_n"] @ x)81 out = {"loc": loc, "echelle": math.sqrt(s2), "nu": post["nu"], "mediane": math.exp(loc),82 "part_coefficients": float((x @ post["V_n"] @ x) / (1 + x @ post["V_n"] @ x))}83 for p in niveaux:84 q = stats.t.ppf(0.5 + p / 2, post["nu"]) * math.sqrt(s2)85 out[f"bas{int(p*100)}"] = math.exp(loc - q); out[f"haut{int(p*100)}"] = math.exp(loc + q)86 return out8788def prior_marche_voisin(post_voisin: dict, colonnes_voisin, colonnes_cible, kappa: float, v0_const: float = 1e4) -> dict:89 """A priori informatif pour un petit marché : les prix implicites (pentes) sont centrés sur la postérieure du marché voisin,90 avec une matrice d'échelle V0 = κ·V_n(voisin) (κ = 1 : confiance entière ; κ = 4 : écart-type doublé) ; la constante (niveau de91 prix du marché cible) reste pratiquement plate. a0, b0 : on transmet aussi l'information sur σ² (a_n, b_n du voisin)."""92 cv = list(colonnes_voisin); cc = list(colonnes_cible); k = len(cc)93 m0 = np.zeros(k); V0 = np.zeros((k, k))94 idx = [cv.index(c) for c in cc if c != "const"]; pos = [i for i, c in enumerate(cc) if c != "const"]95 m0[pos] = post_voisin["m_n"][idx]96 V0[np.ix_(pos, pos)] = kappa * post_voisin["V_n"][np.ix_(idx, idx)]97 ic = cc.index("const"); m0[ic] = LN_MEDIANE_HULL; V0[ic, ic] = v0_const98 return {"m0": m0, "V0inv": np.linalg.inv(V0), "a0": post_voisin["a_n"], "b0": post_voisin["b_n"], "nom": f"marché voisin, κ = {kappa:g}"}99100def loo_exact_nig(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:101 """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)."""102 n = len(y); lp = np.empty(n)103 for i in range(n):104 m = np.ones(n, bool); m[i] = False105 p = nig_posterieure(X[m], y[m], prior)106 s2 = p["b_n"] / p["a_n"] * (1 + X[i] @ p["V_n"] @ X[i])107 lp[i] = stats.t.logpdf(y[i], p["nu"], loc=X[i] @ p["m_n"], scale=math.sqrt(s2))108 return {"elpd_loo": float(lp.sum()), "looic": float(-2 * lp.sum()), "pointwise": lp}109110# =============================================================================================================================111# 2. Modèle hiérarchique à effets de secteur aléatoires : échantillonneur de Gibbs112# =============================================================================================================================113114def matrice_hier(df: pd.DataFrame, moyennes: pd.Series | None = None):115 """Régresseurs du modèle hiérarchique : ceux de traite.cas1.matrice sans la constante ni les indicatrices de secteur, centrés116 (paramétrisation centrée : les niveaux α_s portent l'ordonnée à l'origine). Retourne (X centré, moyennes, indices de secteur)."""117 X = cas1.matrice(df).drop(columns=["const"] + [c for c in cas1.matrice(df).columns if c.startswith("sect_")])118 if moyennes is None: moyennes = X.mean()119 Xc = X - moyennes120 z = np.array([SECTEURS.index(s) for s in df["secteur"]])121 return Xc, moyennes, z122123def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5, n_iter: int = 3000, brulage: int = 1000,124 chaines: int = 4, graine: int = GRAINE, prior_tau: str = "uniforme", tau_fixe: float | None = None,125 ig_a: float = 1e-3, ig_b: float = 1e-3, bloc: bool = True) -> dict:126 """ln P_i = x_i'β + α_{s(i)} + ε_i, ε ~ N(0, σ²) ; α_s | α0, τ² ~ N(α0, τ²) ; p(β) ∝ 1, p(α0) ∝ 1, p(σ²) ∝ 1/σ² ;127 τ : p(τ) ∝ 1 (uniforme sur τ, BDA3 §5.7 ; « uniforme ») ou τ² ~ IG(ig_a, ig_b) (« ig ») ou τ fixé (bayésien empirique).128 Conditionnelles complètes (toutes conjuguées) :129 β | · ~ N((X'X)⁻¹X'(y − Zα), σ²(X'X)⁻¹)130 α_s | · ~ N([n_s r̄_s/σ² + α0/τ²] / [n_s/σ² + 1/τ²], 1/[n_s/σ² + 1/τ²]), r̄_s = moyenne de (y − Xβ) dans le secteur s131 α0 | · ~ N(moyenne des α_s, τ²/J)132 σ² | · ~ IG(n/2, ½ Σ (y − Xβ − Zα)²)133 τ² | · ~ IG((J−1)/2, ½ Σ (α_s − α0)²) [uniforme sur τ] ou IG(a + J/2, b + ½ Σ (α_s − α0)²) [IG]134 bloc=True : (β, α) sont tirés ENSEMBLE de leur conditionnelle jointe normale (Gibbs par blocs), ce qui supprime la lenteur de135 mélange due à la colinéarité entre la distance au centre et les niveaux des secteurs périphériques ; bloc=False : un à la fois.136 Retourne les tirages après brûlage, de forme (chaînes, n_iter − brulage, ·)."""137 y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape138 Z = np.zeros((n, J)); Z[np.arange(n), z] = 1.0139 W = np.hstack([X, Z]); WtW = W.T @ W; Wty = W.T @ y140 XtX_inv = np.linalg.inv(X.T @ X); L_beta = np.linalg.cholesky(XtX_inv)141 ns = np.bincount(z, minlength=J).astype(float)142 garder = n_iter - brulage143 out = {"beta": np.empty((chaines, garder, k)), "alpha": np.empty((chaines, garder, J)), "alpha0": np.empty((chaines, garder)),144 "sigma2": np.empty((chaines, garder)), "tau2": np.empty((chaines, garder))}145 for c in range(chaines):146 rng = np.random.default_rng([graine, c, 39])147 # valeurs initiales surdispersées (BDA3 : points de départ dispersés pour que R̂ ait un sens)148 beta = XtX_inv @ X.T @ y + rng.normal(0, 0.05, k)149 alpha = np.full(J, float(np.mean(y))) + rng.normal(0, 0.3, J); alpha0 = float(np.mean(alpha))150 sigma2 = float(np.exp(rng.normal(math.log(0.01), 0.5))); tau2 = float(np.exp(rng.normal(math.log(0.02), 1.0)))151 if tau_fixe is not None: tau2 = tau_fixe ** 2152 for it in range(n_iter):153 if bloc:154 # (β, α) | α0, σ², τ² : précision = W'W/σ² + diag(0_k, 1/τ² 1_J) ; moyenne = Préc⁻¹ (W'y/σ² + (0, α0/τ² 1_J))155 prec = WtW / sigma2; prec[k:, k:] += np.eye(J) / tau2156 b = Wty / sigma2; b[k:] += alpha0 / tau2157 Lp = np.linalg.cholesky(prec)158 moy = np.linalg.solve(Lp.T, np.linalg.solve(Lp, b))159 theta = moy + np.linalg.solve(Lp.T, rng.standard_normal(k + J))160 beta, alpha = theta[:k], theta[k:]161 e = y - X @ beta162 else:163 # β | α, σ²164 r = y - alpha[z]165 beta = XtX_inv @ X.T @ r + math.sqrt(sigma2) * (L_beta @ rng.standard_normal(k))166 # α_s | β, α0, σ², τ²167 e = y - X @ beta168 somme = np.bincount(z, weights=e, minlength=J)169 prec_a = ns / sigma2 + 1.0 / tau2170 moy_a = (somme / sigma2 + alpha0 / tau2) / prec_a171 alpha = moy_a + rng.standard_normal(J) / np.sqrt(prec_a)172 # α0 | α, τ²173 alpha0 = float(alpha.mean() + rng.standard_normal() * math.sqrt(tau2 / J))174 # σ² | β, α175 res = e - alpha[z]176 sigma2 = 1.0 / rng.gamma(n / 2, 2.0 / (res @ res))177 # τ² | α, α0178 ss = float(((alpha - alpha0) ** 2).sum())179 if tau_fixe is None:180 if prior_tau == "uniforme": tau2 = 1.0 / rng.gamma((J - 1) / 2, 2.0 / ss)181 elif prior_tau == "ig": tau2 = 1.0 / rng.gamma(ig_a + J / 2, 1.0 / (ig_b + ss / 2))182 else: raise ValueError(prior_tau)183 if it >= brulage:184 j = it - brulage185 out["beta"][c, j] = beta; out["alpha"][c, j] = alpha; out["alpha0"][c, j] = alpha0186 out["sigma2"][c, j] = sigma2; out["tau2"][c, j] = tau2187 out["meta"] = {"n": n, "k": k, "J": J, "ns": ns, "n_iter": n_iter, "brulage": brulage, "chaines": chaines, "prior_tau": prior_tau,188 "tau_fixe": tau_fixe, "bloc": bloc}189 return out190191def facteur_retrecissement(n_s: float, sigma2: float, tau2: float) -> float:192 """B_s = n_s τ² / (n_s τ² + σ²) : poids de la moyenne du secteur dans la moyenne conditionnelle de α_s (1 − B_s vers α0)."""193 return n_s * tau2 / (n_s * tau2 + sigma2)194195# ----------------------------------------------------------------------------------------------------------------- diagnostics196def rhat_scinde(ch: np.ndarray) -> float:197 """R̂ scindé (Gelman et Rubin 1992 ; BDA3 §11.4) : chaque chaîne est coupée en deux ; m demi-chaînes de longueur n ;198 B = n·Var(moyennes), W = moyenne des variances ; var⁺ = (n−1)/n W + B/n ; R̂ = sqrt(var⁺/W)."""199 ch = np.asarray(ch, float); c, T = ch.shape; T2 = T // 2200 halves = np.concatenate([ch[:, :T2], ch[:, T2:2 * T2]], axis=0) # (2c, T2)201 m, n = halves.shape202 moy = halves.mean(axis=1); B = n * moy.var(ddof=1); W = halves.var(axis=1, ddof=1).mean()203 var_plus = (n - 1) / n * W + B / n204 return float(math.sqrt(var_plus / W))205206def ess(ch: np.ndarray) -> float:207 """Taille d'échantillon effective (BDA3 §11.5) : ρ̂_t = 1 − V_t/(2 var⁺) avec V_t le variogramme empirique moyen des chaînes ;208 somme des autocorrélations arrêtée au premier couple (ρ̂_{2t'} + ρ̂_{2t'+1}) négatif (Geyer) ; ESS = m n / (1 + 2 Σ ρ̂_t)."""209 ch = np.asarray(ch, float); m, n = ch.shape210 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()211 var_plus = (n - 1) / n * W + B / n212 rho = np.empty(n)213 for t in range(n):214 if t == 0: Vt = 0.0215 else: Vt = float(((ch[:, t:] - ch[:, :-t]) ** 2).sum() / (m * (n - t)))216 rho[t] = 1.0 - Vt / (2 * var_plus)217 s = 0.0; t = 1218 while t + 1 < n:219 pair = rho[t] + rho[t + 1]220 if pair < 0: break221 s += pair; t += 2222 return float(m * n / (1 + 2 * s))223224def resume_tirages(x: np.ndarray, nom: str = "") -> dict:225 """Résumé d'un paramètre scalaire de forme (chaînes, tirages) : moyenne, écart-type, quantiles, R̂, ESS."""226 flat = x.reshape(-1)227 return {"nom": nom, "moyenne": float(flat.mean()), "sd": float(flat.std(ddof=1)), "mediane": float(np.median(flat)),228 "q025": float(np.quantile(flat, 0.025)), "q975": float(np.quantile(flat, 0.975)),229 "q05": float(np.quantile(flat, 0.05)), "q95": float(np.quantile(flat, 0.95)), "rhat": rhat_scinde(x), "ess": ess(x)}230231def resume_gibbs(g: dict, colonnes) -> pd.DataFrame:232 lignes = []233 for j, c in enumerate(colonnes): lignes.append(resume_tirages(g["beta"][:, :, j], c))234 for j, s in enumerate(SECTEURS): lignes.append(resume_tirages(g["alpha"][:, :, j], f"alpha_{s}"))235 lignes.append(resume_tirages(g["alpha0"], "alpha0"))236 lignes.append(resume_tirages(np.sqrt(g["sigma2"]), "sigma"))237 lignes.append(resume_tirages(np.sqrt(g["tau2"]), "tau"))238 return pd.DataFrame(lignes).set_index("nom")239240def effets_relatifs(g: dict, ref: str = "Hull") -> pd.DataFrame:241 """Effets de secteur relatifs à Hull (α_s − α_Hull), comparables aux coefficients des indicatrices du ch. 18."""242 i0 = SECTEURS.index(ref); lignes = []243 for j, s in enumerate(SECTEURS):244 if j == i0: continue245 lignes.append(resume_tirages(g["alpha"][:, :, j] - g["alpha"][:, :, i0], s))246 return pd.DataFrame(lignes).set_index("nom")247248# ----------------------------------------------------------------------------------------------------------------- critères249def log_vrais_pointwise_hier(g: dict, y: np.ndarray, X: np.ndarray, z: np.ndarray, S: int = 2000, graine: int = GRAINE) -> np.ndarray:250 """Matrice S × n des log-densités log p(y_i | θ_s) pour S tirages (sous-échantillonnés) du modèle hiérarchique."""251 rng = np.random.default_rng(graine)252 beta = g["beta"].reshape(-1, g["beta"].shape[-1]); alpha = g["alpha"].reshape(-1, g["alpha"].shape[-1]); s2 = g["sigma2"].reshape(-1)253 idx = rng.choice(len(s2), size=min(S, len(s2)), replace=False)254 mu = beta[idx] @ X.T + alpha[idx][:, z]255 return stats.norm.logpdf(y[None, :], loc=mu, scale=np.sqrt(s2[idx])[:, None])256257def log_vrais_pointwise_nig(post: dict, S: int = 2000, graine: int = GRAINE) -> np.ndarray:258 rng = np.random.default_rng(graine); t = nig_tirages(post, S, rng)259 mu = t["beta"] @ post["X"].T260 return stats.norm.logpdf(post["y"][None, :], loc=mu, scale=np.sqrt(t["sigma2"])[:, None])261262def waic(lp: np.ndarray) -> dict:263 """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) ;264 elpd_waic = lppd − p_waic ; WAIC = −2 elpd_waic."""265 S = lp.shape[0]266 lppd_i = np.logaddexp.reduce(lp, axis=0) - math.log(S)267 p_i = lp.var(axis=0, ddof=1)268 return {"lppd": float(lppd_i.sum()), "p_waic": float(p_i.sum()), "elpd_waic": float((lppd_i - p_i).sum()),269 "waic": float(-2 * (lppd_i - p_i).sum()), "se_elpd": float(math.sqrt(len(lppd_i)) * (lppd_i - p_i).std(ddof=1))}270271def loo_tronque(lp: np.ndarray) -> dict:272 """LOO approximatif par échantillonnage préférentiel tronqué (Ionides 2008 ; forme simplifiée de PSIS, Vehtari et al. 2017) :273 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."""274 S = lp.shape[0]275 lw = -lp; lw -= lw.max(axis=0, keepdims=True)276 w = np.exp(lw); w /= w.sum(axis=0, keepdims=True)277 w = np.minimum(w, math.sqrt(S) / S); w /= w.sum(axis=0, keepdims=True)278 lp_loo = np.log((w * np.exp(lp)).sum(axis=0))279 return {"elpd_loo": float(lp_loo.sum()), "looic": float(-2 * lp_loo.sum()), "pointwise": lp_loo}280281# =============================================================================================================================282# 3. Simulation à partir des lois du générateur (récupération des paramètres)283# =============================================================================================================================284285def simuler_depuis_generateur(v: pd.DataFrame, n: int, graine: int, sigma: float = VRAIS["sigma"], const: float = 13.2) -> pd.DataFrame:286 """Rééchantillonne n lignes d'attributs du jeu (bootstrap des covariables) et engendre ln P selon les lois déclarées287 (README synthétique) : ln P = const + Σ β_vrai x + prime de secteur + ε, ε ~ N(0, σ²). Les attributs sont ceux du jeu ; le prix288 est régénéré. Sert aux tests de récupération des paramètres."""289 rng = np.random.default_rng(graine)290 d = v.sample(n=n, replace=True, random_state=int(rng.integers(1 << 31))).reset_index(drop=True)291 X = cas1.matrice(d)292 beta = pd.Series({c: VRAIS.get(c, 0.0) for c in X.columns}); beta["const"] = const293 lnp = X.values @ beta.values + rng.normal(0, sigma, n)294 d["prix_vente"] = np.round(np.exp(lnp) / 500) * 500295 return d296297# =============================================================================================================================298# 4. Exemple intégral299# =============================================================================================================================300301def sous_echantillons(v: pd.DataFrame, s: dict) -> dict:302 import ch16_selection as c16303 vp = c16.preparer(v.copy(), s)304 micro = vp[(vp["mois"].between(0, 18)) & (vp["dist_km"] <= 3.0)].drop(columns=["mois", "dist_km", "prix_m2", "distance"])305 return {"tout": v, "hull": v[v["secteur"] == "Hull"], "voisin": v[v["secteur"] != "Hull"], "micro": micro}306307def colonnes_actives(X: pd.DataFrame) -> list:308 """Colonnes non constantes ; si le secteur de référence (Hull) est absent du sous-échantillon, les indicatrices restantes somment309 à 1 et l'une d'elles (la première) est retirée pour éviter la colinéarité avec la constante."""310 cols = [c for c in X.columns if c == "const" or X[c].std() > 0]311 sect = [c for c in cols if c.startswith("sect_")]312 if sect and np.allclose(X[sect].sum(axis=1), 1.0): cols.remove(sect[0])313 return cols314315def ajuster_nig(df: pd.DataFrame, prior_fn=None, colonnes=None) -> dict:316 X = cas1.matrice(df); cols = colonnes or colonnes_actives(X); X = X[cols]317 y = np.log(df["prix_vente"].astype(float)).values318 prior = prior_fn(cols) if prior_fn else prior_reference(cols)319 p = nig_posterieure(X.values, y, prior); p["colonnes"] = cols320 return p321322def x_sujet(s: dict, cols, secteur: str | None = None, dist_centre: float | None = None) -> np.ndarray:323 d = dict(s)324 if secteur: d["secteur"] = secteur325 if dist_centre is not None: d["dist_centre_km"] = dist_centre326 Xs = cas1.matrice(cas1.sujet_df(d))327 return Xs[cols].values[0]328329def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:330 """Étapes 1-2 : postérieure NIG sur 600 ventes (a priori faible) ; petit échantillon (40 ventes à 3 km / 18 mois) avec a priori331 faible vs informatif emprunté aux autres secteurs (κ = 1, 4, 16) ; prédictive du sujet."""332 se = sous_echantillons(v, s)333 p600 = ajuster_nig(se["tout"]); cols = p600["colonnes"] # a priori de référence334 p600_propre = ajuster_nig(se["tout"], prior_fn=prior_faible, colonnes=cols) # a priori propre vague (comparaison)335 xs = x_sujet(s, cols)336 tab = pd.DataFrame({"mco": ref["coef"]["coef"].reindex(cols), "se_hc1": ref["coef"]["se_hc1"].reindex(cols),337 "post_moy": p600["m_n"], "post_sd": p600["sd"], "ic_bas": p600["ic95_bas"], "ic_haut": p600["ic95_haut"],338 "vrai": [VRAIS.get(c, np.nan) for c in cols]}, index=cols)339 tab["couvre"] = (tab["vrai"] >= tab["ic_bas"]) & (tab["vrai"] <= tab["ic_haut"])340 pred600 = predictive_t(p600, xs)341 # cohérence séquentielle : 2023-2024 puis 2025-2026342 anc = v[v["date_vente"] < "2025-01-01"]; rec = v[v["date_vente"] >= "2025-01-01"]343 p_anc = ajuster_nig(anc, colonnes=cols)344 prior_seq = {"m0": p_anc["m_n"], "V0inv": np.linalg.inv(p_anc["V_n"]), "a0": p_anc["a_n"], "b0": p_anc["b_n"], "nom": "séquentiel"}345 p_seq = ajuster_nig(rec, prior_fn=lambda c: prior_seq, colonnes=cols)346 # petit échantillon : 40 ventes autour du sujet ; a priori emprunté aux autres secteurs347 p_vois = ajuster_nig(se["voisin"])348 micro = se["micro"]; Xm = cas1.matrice(micro); cols_m = colonnes_actives(Xm)349 petits = {}350 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)]:351 p = ajuster_nig(micro, prior_fn=pf, colonnes=cols_m)352 pr = predictive_t(p, x_sujet(s, cols_m))353 petits[nom] = {"post": p, "pred": pr, "largeur95_pct": (pr["haut95"] - pr["bas95"]) / pr["mediane"] * 100}354 # Hull seul (186), a priori faible355 p_hull = ajuster_nig(se["hull"]); pred_hull = predictive_t(p_hull, x_sujet(s, p_hull["colonnes"]))356 ih = p_hull["colonnes"].index("ln_S"); pred_hull["ln_S"] = float(p_hull["m_n"][ih]); pred_hull["ln_S_sd"] = float(p_hull["sd"][ih])357 pred_hull["largeur95_pct"] = (pred_hull["haut95"] - pred_hull["bas95"]) / pred_hull["mediane"] * 100358 return {"p600": p600, "p600_propre": p600_propre, "pred600_propre": predictive_t(p600_propre, xs), "cols": cols, "xs": xs, "tab600": tab, "pred600": pred600, "p_anc": p_anc, "p_rec_seq": p_seq, "n_anc": len(anc),359 "n_rec": len(rec), "p_vois": p_vois, "micro": micro, "cols_micro": cols_m, "petits": petits, "p_hull": p_hull,360 "pred_hull": pred_hull, "n_micro": len(micro)}361362def exemple_hierarchique(v: pd.DataFrame, s: dict, ref: dict, n_iter: int = 3000, brulage: int = 1000) -> dict:363 """Étapes 3-4 : Gibbs sur les 600 ventes ; effets relatifs à Hull vs effets fixes (ch. 18) ; prédictive du sujet ;364 expérience Masson-Angers réduit à 5, 10, 20 ventes ; sensibilité à l'a priori sur τ."""365 Xc, moy, z = matrice_hier(v); y = np.log(v["prix_vente"].astype(float)).values; cols = list(Xc.columns)366 g = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage) # Gibbs par blocs (résultats du chapitre)367 g_naif = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, bloc=False) # un paramètre à la fois (diagnostics)368 res = resume_gibbs(g, cols); rel = effets_relatifs(g); res_naif = resume_gibbs(g_naif, cols)369 diag = pd.DataFrame({"rhat_bloc": res["rhat"], "ess_bloc": res["ess"], "rhat_naif": res_naif["rhat"], "ess_naif": res_naif["ess"]})370 fixes = ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], ["coef", "se_hc1"]]; fixes.index = SECTEURS[1:]371 rel["fixe"] = fixes["coef"]; rel["fixe_se"] = fixes["se_hc1"]; rel["vrai"] = [VRAIS[f"sect_{x}"] for x in rel.index]372 # prédictive du sujet (Hull) : ln P* = (x_s − x̄)'β + α_Hull + ε373 xs = (cas1.matrice(cas1.sujet_df(s))[cols].iloc[0] - moy).values374 beta = g["beta"].reshape(-1, len(cols)); alpha = g["alpha"].reshape(-1, 5); s2 = g["sigma2"].reshape(-1)375 rng = np.random.default_rng(GRAINE)376 lnp = beta @ xs + alpha[:, 0] + rng.standard_normal(len(s2)) * np.sqrt(s2)377 P = np.exp(lnp)378 pred = {"tirages": P, "mediane": float(np.median(P)), "moyenne": float(P.mean()), "moyenne_log": float(lnp.mean()),379 "bas90": float(np.quantile(P, 0.05)), "haut90": float(np.quantile(P, 0.95)),380 "bas95": float(np.quantile(P, 0.025)), "haut95": float(np.quantile(P, 0.975)),381 "sd_log": float(lnp.std()), "sd_moyenne_log": float((beta @ xs + alpha[:, 0]).std())}382 # rétrécissement a posteriori : facteur B_s à la médiane de σ², τ²383 s2m = float(np.median(g["sigma2"])); t2m = float(np.median(g["tau2"]))384 B = {sct: facteur_retrecissement(g["meta"]["ns"][j], s2m, t2m) for j, sct in enumerate(SECTEURS)}385 # expérience Masson-Angers : 5, 10, 20, 47 ventes (sous-ensembles emboîtés)386 ma = v[v["secteur"] == "Masson-Angers"]; autres = v[v["secteur"] != "Masson-Angers"]387 perm = np.random.default_rng(GRAINE).permutation(len(ma))388 dist_ma = float(ma["dist_centre_km"].median())389 exp_rows = []390 for n_ma in (5, 10, 20, len(ma)):391 vv = pd.concat([autres, ma.iloc[perm[:n_ma]]]).reset_index(drop=True)392 # effets fixes (MCO, HC1) et prédictive conjuguée (plat) pour un sujet hypothétique à Masson-Angers393 pf = ajuster_nig(vv); i = pf["colonnes"].index("sect_Masson-Angers")394 xma = x_sujet(s, pf["colonnes"], secteur="Masson-Angers", dist_centre=dist_ma)395 pr_f = predictive_t(pf, xma)396 # hiérarchique397 Xv, moy_v, zv = matrice_hier(vv); yv = np.log(vv["prix_vente"].astype(float)).values398 gv = gibbs_hierarchique(yv, Xv.values, zv, n_iter=n_iter, brulage=brulage)399 rr = effets_relatifs(gv).loc["Masson-Angers"]400 xma_c = (cas1.matrice(cas1.sujet_df({**s, "secteur": "Masson-Angers", "dist_centre_km": dist_ma}))[cols].iloc[0] - moy_v).values401 bv = gv["beta"].reshape(-1, len(cols)); av = gv["alpha"].reshape(-1, 5); sv = gv["sigma2"].reshape(-1)402 lnv = bv @ xma_c + av[:, 4] + np.random.default_rng(GRAINE + n_ma).standard_normal(len(sv)) * np.sqrt(sv)403 Pv = np.exp(lnv)404 s2v = float(np.median(gv["sigma2"])); t2v = float(np.median(gv["tau2"]))405 # moyenne brute des résidus du secteur (sans rétrécissement) à la postérieure : r̄ − α0 relatif406 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],407 "hier": rr["moyenne"], "hier_sd": rr["sd"], "hier_bas": rr["q025"], "hier_haut": rr["q975"], "hier_rhat": rr["rhat"],408 "tau_med": math.sqrt(t2v), "B": facteur_retrecissement(n_ma, s2v, t2v),409 "pred_f_med": pr_f["mediane"], "pred_f_bas": pr_f["bas95"], "pred_f_haut": pr_f["haut95"],410 "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)),411 "vrai": VRAIS["sect_Masson-Angers"]})412 experience = pd.DataFrame(exp_rows).set_index("n_ma")413 experience["largeur_f_pct"] = (experience["pred_f_haut"] - experience["pred_f_bas"]) / experience["pred_f_med"] * 100414 experience["largeur_h_pct"] = (experience["pred_h_haut"] - experience["pred_h_bas"]) / experience["pred_h_med"] * 100415 # sensibilité à l'a priori sur τ (600 ventes et n_MA = 5)416 v5 = pd.concat([autres, ma.iloc[perm[:5]]]).reset_index(drop=True)417 X5, _, z5 = matrice_hier(v5); y5 = np.log(v5["prix_vente"].astype(float)).values418 tau_eb = float(np.std(ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], "coef"].tolist() + [0.0], ddof=1))419 sens_rows = []420 for nom, kw in [("uniforme sur τ", {"prior_tau": "uniforme"}), ("IG(0,001 ; 0,001) sur τ²", {"prior_tau": "ig"}),421 ("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})]:422 g600 = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, **kw)423 g5 = gibbs_hierarchique(y5, X5.values, z5, n_iter=n_iter, brulage=brulage, **kw)424 r600 = effets_relatifs(g600).loc["Masson-Angers"]; r5 = effets_relatifs(g5).loc["Masson-Angers"]425 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)),426 "tau_q95_600": float(np.quantile(np.sqrt(g600["tau2"]), 0.95)),427 "ma_600": r600["moyenne"], "ma_600_sd": r600["sd"], "tau_med_5": float(np.median(np.sqrt(g5["tau2"]))),428 "ma_5": r5["moyenne"], "ma_5_sd": r5["sd"], "ma_5_bas": r5["q025"], "ma_5_haut": r5["q975"]})429 sensibilite = pd.DataFrame(sens_rows).set_index("prior")430 return {"g": g, "g_naif": g_naif, "diag": diag, "cols": cols, "moyennes": moy, "z": z, "y": y, "Xc": Xc, "resume": res, "relatifs": rel, "pred": pred, "B": B,431 "sigma2_med": s2m, "tau2_med": t2m, "experience": experience, "dist_ma": dist_ma, "sensibilite": sensibilite, "tau_eb": tau_eb,432 "perm_ma": perm}433434def criteres(v: pd.DataFrame, hier: dict, ex: dict) -> pd.DataFrame:435 """Comparaison de trois modèles : pooling complet (sans secteur), effets fixes (indicatrices), hiérarchique."""436 y = hier["y"]437 # pooling complet : matrice de référence sans indicatrices de secteur438 X = cas1.matrice(v); cols_pool = [c for c in X.columns if not c.startswith("sect_")]439 p_pool = nig_posterieure(X[cols_pool].values, y, prior_reference(cols_pool)); p_pool["colonnes"] = cols_pool440 p_fixe = ex["p600"]441 lignes = []442 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"])),443 ("effets fixes de secteur (ch. 18)", log_vrais_pointwise_nig(p_fixe), loo_exact_nig(p_fixe["X"], y, p_fixe["prior"])),444 ("effets aléatoires de secteur (hiérarchique)", log_vrais_pointwise_hier(hier["g"], y, hier["Xc"].values, hier["z"]), None)]:445 w = waic(lp); lt = loo_tronque(lp)446 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"],447 "elpd_loo_tronque": lt["elpd_loo"], "elpd_loo_exact": loo["elpd_loo"] if loo else np.nan})448 t = pd.DataFrame(lignes).set_index("modele")449 t["d_elpd_waic"] = t["elpd_waic"] - t["elpd_waic"].max()450 return t451452# =============================================================================================================================453# 5. Exercices (chiffres des solutions)454# =============================================================================================================================455456def exercices(ex: dict, hier: dict) -> dict:457 r = {}458 # ex 1 : mise à jour normale-normale, variance connue (prix au m² d'un marché mince)459 m0, s0, sig, n, xbar = 3300.0, 200.0, 300.0, 4, 3520.0460 prec = 1 / s0 ** 2 + n / sig ** 2; mn = (m0 / s0 ** 2 + n * xbar / sig ** 2) / prec; sn = math.sqrt(1 / prec)461 w = (n / sig ** 2) / prec462 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),463 "n_equivalent_prior": sig ** 2 / s0 ** 2, "valeur_165": mn * 165, "sd_165": sn * 165}464 prec4 = 1 / 400.0 ** 2 + n / sig ** 2465 r["ex1"]["s0_400"] = {"m_n": (m0 / 400.0 ** 2 + n * xbar / sig ** 2) / prec4, "s_n": math.sqrt(1 / prec4), "poids_donnees": (n / sig ** 2) / prec4,466 "n_equivalent_prior": sig ** 2 / 400.0 ** 2}467 # ex 2 : lecture d'intervalles — probabilité a posteriori que la valeur dépasse 650 000 $ sous une postérieure log-normale468 mu, sd = ex["pred600"]["loc"], ex["pred600"]["echelle"]469 r["ex2"] = {"p_sup_650": float(1 - stats.t.cdf((math.log(650_000) - mu) / sd, ex["pred600"]["nu"])),470 "p_inf_550": float(stats.t.cdf((math.log(550_000) - mu) / sd, ex["pred600"]["nu"])),471 "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),472 "z_650": (math.log(650_000) - mu) / sd, "z_550": (math.log(550_000) - mu) / sd}473 # ex 3 : facteur de rétrécissement474 s2, t2 = 0.083 ** 2, 0.15 ** 2475 r["ex3"] = {n_: facteur_retrecissement(n_, s2, t2) for n_ in (5, 10, 20, 47)}476 r["ex3"]["alpha_5"] = facteur_retrecissement(5, s2, t2) * (-0.30) + (1 - facteur_retrecissement(5, s2, t2)) * (-0.08)477 r["ex3"]["tau_005"] = {n_: facteur_retrecissement(n_, s2, 0.05 ** 2) for n_ in (5, 10, 20, 47)}478 r["ex3"]["sd_cond_5"] = math.sqrt(1 / (5 / s2 + 1 / t2)); r["ex3"]["sd_fixe_5"] = math.sqrt(s2 / 5)479 r["ex3"]["n_demi"] = {"tau_015": s2 / t2, "tau_005": s2 / 0.05 ** 2}480 # ex 4 : prédictive Student-t481 p = ex["p600"]; xs = ex["xs"]; q = float(xs @ p["V_n"] @ xs)482 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"],483 "t975": float(stats.t.ppf(0.975, p["nu"])), "z975": 1.959964, "bas95": ex["pred600"]["bas95"], "haut95": ex["pred600"]["haut95"],484 "nu_petit": 20, "t975_20": float(stats.t.ppf(0.975, 20)), "mediane": ex["pred600"]["mediane"],485 "bas95_normal": math.exp(ex["pred600"]["loc"] - 1.959964 * ex["pred600"]["echelle"]), "haut95_normal": math.exp(ex["pred600"]["loc"] + 1.959964 * ex["pred600"]["echelle"]),486 "bas95_t20": math.exp(ex["pred600"]["loc"] - float(stats.t.ppf(0.975, 20)) * ex["pred600"]["echelle"]), "haut95_t20": math.exp(ex["pred600"]["loc"] + float(stats.t.ppf(0.975, 20)) * ex["pred600"]["echelle"]),487 "elargissement_pct": (float(stats.t.ppf(0.975, 20)) / float(stats.t.ppf(0.975, p["nu"])) - 1) * 100}488 # ex 5 : R̂ sur des moyennes/variances de chaînes données489 moy = np.array([-0.262, -0.258, -0.271, -0.249]); var = np.array([0.0049, 0.0052, 0.0047, 0.0055]); n_ = 1000490 B = n_ * moy.var(ddof=1); W = var.mean(); vp = (n_ - 1) / n_ * W + B / n_491 r["ex5"] = {"B": B, "W": W, "var_plus": vp, "rhat": math.sqrt(vp / W)}492 moy2 = np.array([-0.262, -0.258, -0.371, -0.249]); B2 = n_ * moy2.var(ddof=1); vp2 = (n_ - 1) / n_ * W + B2 / n_493 r["ex5"]["rhat_divergent"] = math.sqrt(vp2 / W); r["ex5"]["B2"] = B2494 r["ex5"]["ess_rho08"] = 4000 / (1 + 2 * 0.8 / (1 - 0.8))495 # ex 6 : script (κ) — valeurs du tableau des petits échantillons496 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")],497 "mediane": v_["pred"]["mediane"], "largeur95_pct": v_["largeur95_pct"]} for k, v_ in ex["petits"].items()}498 # ex 7 : conditionnelle de α0 avec a priori N(m, s²) — valeurs numériques à τ et α médians499 g = hier["g"]; alpha_bar = float(g["alpha"].reshape(-1, 5).mean()); tau2 = hier["tau2_med"]500 for (m_, s_) in [(13.0, 1.0), (13.0, 0.05)]:501 prec = 5 / tau2 + 1 / s_ ** 2502 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)}503 return r504505# =============================================================================================================================506# 6. Main507# =============================================================================================================================508509def main(verbose: bool = True, n_iter: int = 3000, brulage: int = 1000) -> dict:510 d = cas1.charger(); v = d["ventes"]; s = d["sujet"]511 ref = cas1.hedonique(v)512 ex = exemple_conjugue(v, s, ref)513 hier = exemple_hierarchique(v, s, ref, n_iter=n_iter, brulage=brulage)514 crit = criteres(v, hier, ex)515 exo = exercices(ex, hier)516 out = {"ventes": v, "sujet": s, "modele": ref, "conjugue": ex, "hier": hier, "criteres": crit, "exercices": exo}517 if verbose: imprimer(out)518 return out519520def _f(x, d=4): return f"{x:.{d}f}"521def _d(x): return f"{x:,.0f} $".replace(",", " ")522523def imprimer(r: dict) -> None:524 ex = r["conjugue"]; hi = r["hier"]525 print("=== CH39 — modèle conjugué, 600 ventes, a priori faiblement informatif ===")526 t = ex["tab600"]527 print(t.round(4).to_string())528 print(f"couverture des vrais paramètres (IC 95 %) : {int(t['couvre'].sum())}/{int(t['vrai'].notna().sum())}")529 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}")530 pp = ex["p600_propre"]; print(f"a priori propre vague : a_n = {pp['a_n']:.3f} ; ν = {pp['nu']:.0f} ; E[σ] ≈ {pp['sigma_moy']:.4f} ; max|m_n − β̂| = {np.abs(pp['m_n'] - r['modele']['beta']).max():.2e} ; sujet 95 % [{_d(ex['pred600_propre']['bas95'])} ; {_d(ex['pred600_propre']['haut95'])}]")531 pr = ex["pred600"]532 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} %")533 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'])}])")534 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}")535 print(f"\n=== petit échantillon : {ex['n_micro']} ventes à 3 km / 18 mois ; colonnes {ex['cols_micro']}")536 for nom, pe in ex["petits"].items():537 po = pe["post"]; i = ex["cols_micro"].index("ln_S"); iq = ex["cols_micro"].index("qualite"); ig = ex["cols_micro"].index("garage")538 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}) ; "539 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} %")540 ph = ex["pred_hull"]; print(f" Hull seul ({ex['p_hull']['n']}) : ln_S {ph['ln_S']:.3f} ({ph['ln_S_sd']:.3f}) ; médiane {_d(ph['mediane'])} [{_d(ph['bas95'])} ; {_d(ph['haut95'])}] largeur {ph['largeur95_pct']:.1f} %")541 print("\n=== modèle hiérarchique (Gibbs par blocs) : résumé ===")542 print(hi["resume"].round(4).to_string())543 print("\ndiagnostics : Gibbs par blocs vs un paramètre à la fois"); print(hi["diag"].round(3).to_string())544 print("\neffets relatifs à Hull (hiérarchique vs effets fixes ch. 18 vs vrai) :"); print(hi["relatifs"].round(4).to_string())545 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})")546 print("facteurs de rétrécissement B_s :", {k: round(v_, 4) for k, v_ in hi["B"].items()})547 pd_ = hi["pred"]548 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})")549 print(f"\n=== expérience Masson-Angers (dist. médiane {hi['dist_ma']:.1f} km) ===")550 print(hi["experience"].round(4).to_string())551 print("\n=== sensibilité à l'a priori sur τ ==="); print(hi["sensibilite"].round(4).to_string())552 print("\n=== critères ==="); print(r["criteres"].round(2).to_string())553 print("\n=== exercices ===")554 for k, v_ in r["exercices"].items(): print(k, v_)555556if __name__ == "__main__":557 main()558