#!/usr/bin/env python3
"""Chapitre 39 — Modèles bayésiens et intervalles : exemple intégral sur les 600 ventes du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES,
graine 20261031). Numpy/scipy uniquement (pymc et arviz absents) : tout est implémenté ici.

Contenu :
  * régression linéaire bayésienne conjuguée (a priori normale–inverse-gamma, NIG) : postérieure, prédictive Student-t, a priori
    faiblement informatif vs a priori informatif « emprunté » au marché voisin (autres secteurs) pour un petit échantillon ;
  * modèle hiérarchique à effets de secteur aléatoires, échantillonneur de Gibbs (paramétrisation centrée de Gelfand, Sahu et Carlin
    1995), rétrécissement, expérience « Masson-Angers réduit à 5, 10, 20 ventes » ;
  * diagnostics MCMC : R̂ de Gelman–Rubin (version scindée, BDA3 §11.4), taille d'échantillon effective (autocorrélations, Geyer) ;
  * distribution prédictive du sujet (intervalles de crédibilité 90/95 %) vs intervalle de prédiction fréquentiste du ch. 18 ;
  * critères : WAIC (Watanabe 2010 ; Vehtari, Gelman et Gabry 2017), LOO exact pour le modèle conjugué, LOO par échantillonnage
    préférentiel tronqué pour le modèle hiérarchique ;
  * simulation à partir des lois du générateur (data/synthetic/README.md) pour la récupération des paramètres (tests).
Le modèle de référence (régresseurs) est `traite.cas1.matrice` ; il n'est pas redéfini.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch39_bayes.py"""
from __future__ import annotations
import sys, math, pathlib
import numpy as np, pandas as pd
from scipy import stats
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
from traite import cas1
from traite import hedonique as h
import ch18_hedonique as c18

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031
VRAIS = c18.VRAIS
SECTEURS = cas1.SECTEURS                      # Hull, Gatineau, Aylmer, Buckingham, Masson-Angers
IC18 = {"point": 637_953, "bas95": 539_682, "haut95": 749_087}   # chapitre 18 (traite.cas1.intervalle_prediction)
LN_MEDIANE_HULL = math.log(475_000)           # médiane APCIQ T3 2026 secteur Hull : centre de l'a priori sur la constante

# =============================================================================================================================
# 1. Régression linéaire bayésienne conjuguée (normale–inverse-gamma)
# =============================================================================================================================

def prior_reference(colonnes) -> dict:
    """A priori de référence (impropre, Jeffreys) : p(β, σ²) ∝ 1/σ². Cas limite de la NIG avec V0⁻¹ = 0, a0 = −k/2, b0 = 0 : la
    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."""
    k = len(colonnes)
    return {"m0": np.zeros(k), "V0inv": np.zeros((k, k)), "a0": -k / 2, "b0": 0.0, "nom": "référence (impropre)"}

def prior_faible(colonnes, v0: float = 1e4, a0: float = 1e-3, b0: float = 1e-3) -> dict:
    """A priori faiblement informatif PROPRE : β | σ² ~ N(m0, σ² V0), σ² ~ IG(a0, b0). m0 = ln(médiane APCIQ de Hull) pour la constante,
    0 ailleurs ; V0 = v0·I (écart-type a priori ≈ 100 σ ≈ 8 en logarithme : pratiquement plat). Prix d'un a priori propre et vague :
    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)."""
    k = len(colonnes); m0 = np.zeros(k)
    if "const" in colonnes: m0[list(colonnes).index("const")] = LN_MEDIANE_HULL
    return {"m0": m0, "V0inv": np.eye(k) / v0, "a0": a0, "b0": b0, "nom": "faiblement informatif (propre)"}

def nig_posterieure(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:
    """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)
    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)."""
    X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape
    V0inv = prior["V0inv"]; m0 = prior["m0"]
    Vn_inv = V0inv + X.T @ X; Vn = np.linalg.inv(Vn_inv)
    mn = Vn @ (V0inv @ m0 + X.T @ y)
    an = prior["a0"] + n / 2
    bn = prior["b0"] + 0.5 * (y @ y + m0 @ V0inv @ m0 - mn @ Vn_inv @ mn)
    nu = 2 * an
    sd = np.sqrt(bn / an * np.diag(Vn)) * (math.sqrt(nu / (nu - 2)) if nu > 2 else np.nan)   # écart-type de la marginale t
    echelle = np.sqrt(bn / an * np.diag(Vn))
    tq = stats.t.ppf(0.975, nu)
    return {"m_n": mn, "V_n": Vn, "a_n": an, "b_n": bn, "nu": nu, "sd": sd, "echelle": echelle,
            "ic95_bas": mn - tq * echelle, "ic95_haut": mn + tq * echelle,
            "sigma2_moy": bn / (an - 1), "sigma_moy": math.sqrt(bn / (an - 1)), "sigma2_mode": bn / (an + 1),
            "n": n, "k": k, "X": X, "y": y, "prior": prior}

def nig_tirages(post: dict, S: int, rng: np.random.Generator) -> dict:
    """Tirages exacts de la postérieure NIG : σ² ~ IG(a_n, b_n) puis β | σ² ~ N(m_n, σ² V_n)."""
    sigma2 = 1.0 / rng.gamma(post["a_n"], 1.0 / post["b_n"], S)
    L = np.linalg.cholesky(post["V_n"])
    z = rng.standard_normal((S, post["k"]))
    beta = post["m_n"] + (z @ L.T) * np.sqrt(sigma2)[:, None]
    return {"beta": beta, "sigma2": sigma2}

def predictive_t(post: dict, x: np.ndarray, niveaux=(0.90, 0.95)) -> dict:
    """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 $."""
    x = np.asarray(x, float); loc = float(x @ post["m_n"]); s2 = post["b_n"] / post["a_n"] * (1 + x @ post["V_n"] @ x)
    out = {"loc": loc, "echelle": math.sqrt(s2), "nu": post["nu"], "mediane": math.exp(loc),
           "part_coefficients": float((x @ post["V_n"] @ x) / (1 + x @ post["V_n"] @ x))}
    for p in niveaux:
        q = stats.t.ppf(0.5 + p / 2, post["nu"]) * math.sqrt(s2)
        out[f"bas{int(p*100)}"] = math.exp(loc - q); out[f"haut{int(p*100)}"] = math.exp(loc + q)
    return out

def prior_marche_voisin(post_voisin: dict, colonnes_voisin, colonnes_cible, kappa: float, v0_const: float = 1e4) -> dict:
    """A priori informatif pour un petit marché : les prix implicites (pentes) sont centrés sur la postérieure du marché voisin,
    avec une matrice d'échelle V0 = κ·V_n(voisin) (κ = 1 : confiance entière ; κ = 4 : écart-type doublé) ; la constante (niveau de
    prix du marché cible) reste pratiquement plate. a0, b0 : on transmet aussi l'information sur σ² (a_n, b_n du voisin)."""
    cv = list(colonnes_voisin); cc = list(colonnes_cible); k = len(cc)
    m0 = np.zeros(k); V0 = np.zeros((k, k))
    idx = [cv.index(c) for c in cc if c != "const"]; pos = [i for i, c in enumerate(cc) if c != "const"]
    m0[pos] = post_voisin["m_n"][idx]
    V0[np.ix_(pos, pos)] = kappa * post_voisin["V_n"][np.ix_(idx, idx)]
    ic = cc.index("const"); m0[ic] = LN_MEDIANE_HULL; V0[ic, ic] = v0_const
    return {"m0": m0, "V0inv": np.linalg.inv(V0), "a0": post_voisin["a_n"], "b0": post_voisin["b_n"], "nom": f"marché voisin, κ = {kappa:g}"}

def loo_exact_nig(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:
    """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)."""
    n = len(y); lp = np.empty(n)
    for i in range(n):
        m = np.ones(n, bool); m[i] = False
        p = nig_posterieure(X[m], y[m], prior)
        s2 = p["b_n"] / p["a_n"] * (1 + X[i] @ p["V_n"] @ X[i])
        lp[i] = stats.t.logpdf(y[i], p["nu"], loc=X[i] @ p["m_n"], scale=math.sqrt(s2))
    return {"elpd_loo": float(lp.sum()), "looic": float(-2 * lp.sum()), "pointwise": lp}

# =============================================================================================================================
# 2. Modèle hiérarchique à effets de secteur aléatoires : échantillonneur de Gibbs
# =============================================================================================================================

def matrice_hier(df: pd.DataFrame, moyennes: pd.Series | None = None):
    """Régresseurs du modèle hiérarchique : ceux de traite.cas1.matrice sans la constante ni les indicatrices de secteur, centrés
    (paramétrisation centrée : les niveaux α_s portent l'ordonnée à l'origine). Retourne (X centré, moyennes, indices de secteur)."""
    X = cas1.matrice(df).drop(columns=["const"] + [c for c in cas1.matrice(df).columns if c.startswith("sect_")])
    if moyennes is None: moyennes = X.mean()
    Xc = X - moyennes
    z = np.array([SECTEURS.index(s) for s in df["secteur"]])
    return Xc, moyennes, z

def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5, n_iter: int = 3000, brulage: int = 1000,
                       chaines: int = 4, graine: int = GRAINE, prior_tau: str = "uniforme", tau_fixe: float | None = None,
                       ig_a: float = 1e-3, ig_b: float = 1e-3, bloc: bool = True) -> dict:
    """ln P_i = x_i'β + α_{s(i)} + ε_i, ε ~ N(0, σ²) ; α_s | α0, τ² ~ N(α0, τ²) ; p(β) ∝ 1, p(α0) ∝ 1, p(σ²) ∝ 1/σ² ;
    τ : p(τ) ∝ 1 (uniforme sur τ, BDA3 §5.7 ; « uniforme ») ou τ² ~ IG(ig_a, ig_b) (« ig ») ou τ fixé (bayésien empirique).
    Conditionnelles complètes (toutes conjuguées) :
      β | ·  ~ N((X'X)⁻¹X'(y − Zα), σ²(X'X)⁻¹)
      α_s | · ~ N([n_s r̄_s/σ² + α0/τ²] / [n_s/σ² + 1/τ²], 1/[n_s/σ² + 1/τ²]),  r̄_s = moyenne de (y − Xβ) dans le secteur s
      α0 | · ~ N(moyenne des α_s, τ²/J)
      σ² | · ~ IG(n/2, ½ Σ (y − Xβ − Zα)²)
      τ² | · ~ IG((J−1)/2, ½ Σ (α_s − α0)²)   [uniforme sur τ]   ou   IG(a + J/2, b + ½ Σ (α_s − α0)²)   [IG]
    bloc=True : (β, α) sont tirés ENSEMBLE de leur conditionnelle jointe normale (Gibbs par blocs), ce qui supprime la lenteur de
    mélange due à la colinéarité entre la distance au centre et les niveaux des secteurs périphériques ; bloc=False : un à la fois.
    Retourne les tirages après brûlage, de forme (chaînes, n_iter − brulage, ·)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    Z = np.zeros((n, J)); Z[np.arange(n), z] = 1.0
    W = np.hstack([X, Z]); WtW = W.T @ W; Wty = W.T @ y
    XtX_inv = np.linalg.inv(X.T @ X); L_beta = np.linalg.cholesky(XtX_inv)
    ns = np.bincount(z, minlength=J).astype(float)
    garder = n_iter - brulage
    out = {"beta": np.empty((chaines, garder, k)), "alpha": np.empty((chaines, garder, J)), "alpha0": np.empty((chaines, garder)),
           "sigma2": np.empty((chaines, garder)), "tau2": np.empty((chaines, garder))}
    for c in range(chaines):
        rng = np.random.default_rng([graine, c, 39])
        # valeurs initiales surdispersées (BDA3 : points de départ dispersés pour que R̂ ait un sens)
        beta = XtX_inv @ X.T @ y + rng.normal(0, 0.05, k)
        alpha = np.full(J, float(np.mean(y))) + rng.normal(0, 0.3, J); alpha0 = float(np.mean(alpha))
        sigma2 = float(np.exp(rng.normal(math.log(0.01), 0.5))); tau2 = float(np.exp(rng.normal(math.log(0.02), 1.0)))
        if tau_fixe is not None: tau2 = tau_fixe ** 2
        for it in range(n_iter):
            if bloc:
                # (β, α) | α0, σ², τ² : précision = W'W/σ² + diag(0_k, 1/τ² 1_J) ; moyenne = Préc⁻¹ (W'y/σ² + (0, α0/τ² 1_J))
                prec = WtW / sigma2; prec[k:, k:] += np.eye(J) / tau2
                b = Wty / sigma2; b[k:] += alpha0 / tau2
                Lp = np.linalg.cholesky(prec)
                moy = np.linalg.solve(Lp.T, np.linalg.solve(Lp, b))
                theta = moy + np.linalg.solve(Lp.T, rng.standard_normal(k + J))
                beta, alpha = theta[:k], theta[k:]
                e = y - X @ beta
            else:
                # β | α, σ²
                r = y - alpha[z]
                beta = XtX_inv @ X.T @ r + math.sqrt(sigma2) * (L_beta @ rng.standard_normal(k))
                # α_s | β, α0, σ², τ²
                e = y - X @ beta
                somme = np.bincount(z, weights=e, minlength=J)
                prec_a = ns / sigma2 + 1.0 / tau2
                moy_a = (somme / sigma2 + alpha0 / tau2) / prec_a
                alpha = moy_a + rng.standard_normal(J) / np.sqrt(prec_a)
            # α0 | α, τ²
            alpha0 = float(alpha.mean() + rng.standard_normal() * math.sqrt(tau2 / J))
            # σ² | β, α
            res = e - alpha[z]
            sigma2 = 1.0 / rng.gamma(n / 2, 2.0 / (res @ res))
            # τ² | α, α0
            ss = float(((alpha - alpha0) ** 2).sum())
            if tau_fixe is None:
                if prior_tau == "uniforme": tau2 = 1.0 / rng.gamma((J - 1) / 2, 2.0 / ss)
                elif prior_tau == "ig": tau2 = 1.0 / rng.gamma(ig_a + J / 2, 1.0 / (ig_b + ss / 2))
                else: raise ValueError(prior_tau)
            if it >= brulage:
                j = it - brulage
                out["beta"][c, j] = beta; out["alpha"][c, j] = alpha; out["alpha0"][c, j] = alpha0
                out["sigma2"][c, j] = sigma2; out["tau2"][c, j] = tau2
    out["meta"] = {"n": n, "k": k, "J": J, "ns": ns, "n_iter": n_iter, "brulage": brulage, "chaines": chaines, "prior_tau": prior_tau,
                   "tau_fixe": tau_fixe, "bloc": bloc}
    return out

def facteur_retrecissement(n_s: float, sigma2: float, tau2: float) -> float:
    """B_s = n_s τ² / (n_s τ² + σ²) : poids de la moyenne du secteur dans la moyenne conditionnelle de α_s (1 − B_s vers α0)."""
    return n_s * tau2 / (n_s * tau2 + sigma2)

# ----------------------------------------------------------------------------------------------------------------- diagnostics
def rhat_scinde(ch: np.ndarray) -> float:
    """R̂ scindé (Gelman et Rubin 1992 ; BDA3 §11.4) : chaque chaîne est coupée en deux ; m demi-chaînes de longueur n ;
    B = n·Var(moyennes), W = moyenne des variances ; var⁺ = (n−1)/n W + B/n ; R̂ = sqrt(var⁺/W)."""
    ch = np.asarray(ch, float); c, T = ch.shape; T2 = T // 2
    halves = np.concatenate([ch[:, :T2], ch[:, T2:2 * T2]], axis=0)   # (2c, T2)
    m, n = halves.shape
    moy = halves.mean(axis=1); B = n * moy.var(ddof=1); W = halves.var(axis=1, ddof=1).mean()
    var_plus = (n - 1) / n * W + B / n
    return float(math.sqrt(var_plus / W))

def ess(ch: np.ndarray) -> float:
    """Taille d'échantillon effective (BDA3 §11.5) : ρ̂_t = 1 − V_t/(2 var⁺) avec V_t le variogramme empirique moyen des chaînes ;
    somme des autocorrélations arrêtée au premier couple (ρ̂_{2t'} + ρ̂_{2t'+1}) négatif (Geyer) ; ESS = m n / (1 + 2 Σ ρ̂_t)."""
    ch = np.asarray(ch, float); m, n = ch.shape
    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()
    var_plus = (n - 1) / n * W + B / n
    rho = np.empty(n)
    for t in range(n):
        if t == 0: Vt = 0.0
        else: Vt = float(((ch[:, t:] - ch[:, :-t]) ** 2).sum() / (m * (n - t)))
        rho[t] = 1.0 - Vt / (2 * var_plus)
    s = 0.0; t = 1
    while t + 1 < n:
        pair = rho[t] + rho[t + 1]
        if pair < 0: break
        s += pair; t += 2
    return float(m * n / (1 + 2 * s))

def resume_tirages(x: np.ndarray, nom: str = "") -> dict:
    """Résumé d'un paramètre scalaire de forme (chaînes, tirages) : moyenne, écart-type, quantiles, R̂, ESS."""
    flat = x.reshape(-1)
    return {"nom": nom, "moyenne": float(flat.mean()), "sd": float(flat.std(ddof=1)), "mediane": float(np.median(flat)),
            "q025": float(np.quantile(flat, 0.025)), "q975": float(np.quantile(flat, 0.975)),
            "q05": float(np.quantile(flat, 0.05)), "q95": float(np.quantile(flat, 0.95)), "rhat": rhat_scinde(x), "ess": ess(x)}

def resume_gibbs(g: dict, colonnes) -> pd.DataFrame:
    lignes = []
    for j, c in enumerate(colonnes): lignes.append(resume_tirages(g["beta"][:, :, j], c))
    for j, s in enumerate(SECTEURS): lignes.append(resume_tirages(g["alpha"][:, :, j], f"alpha_{s}"))
    lignes.append(resume_tirages(g["alpha0"], "alpha0"))
    lignes.append(resume_tirages(np.sqrt(g["sigma2"]), "sigma"))
    lignes.append(resume_tirages(np.sqrt(g["tau2"]), "tau"))
    return pd.DataFrame(lignes).set_index("nom")

def effets_relatifs(g: dict, ref: str = "Hull") -> pd.DataFrame:
    """Effets de secteur relatifs à Hull (α_s − α_Hull), comparables aux coefficients des indicatrices du ch. 18."""
    i0 = SECTEURS.index(ref); lignes = []
    for j, s in enumerate(SECTEURS):
        if j == i0: continue
        lignes.append(resume_tirages(g["alpha"][:, :, j] - g["alpha"][:, :, i0], s))
    return pd.DataFrame(lignes).set_index("nom")

# ----------------------------------------------------------------------------------------------------------------- critères
def log_vrais_pointwise_hier(g: dict, y: np.ndarray, X: np.ndarray, z: np.ndarray, S: int = 2000, graine: int = GRAINE) -> np.ndarray:
    """Matrice S × n des log-densités log p(y_i | θ_s) pour S tirages (sous-échantillonnés) du modèle hiérarchique."""
    rng = np.random.default_rng(graine)
    beta = g["beta"].reshape(-1, g["beta"].shape[-1]); alpha = g["alpha"].reshape(-1, g["alpha"].shape[-1]); s2 = g["sigma2"].reshape(-1)
    idx = rng.choice(len(s2), size=min(S, len(s2)), replace=False)
    mu = beta[idx] @ X.T + alpha[idx][:, z]
    return stats.norm.logpdf(y[None, :], loc=mu, scale=np.sqrt(s2[idx])[:, None])

def log_vrais_pointwise_nig(post: dict, S: int = 2000, graine: int = GRAINE) -> np.ndarray:
    rng = np.random.default_rng(graine); t = nig_tirages(post, S, rng)
    mu = t["beta"] @ post["X"].T
    return stats.norm.logpdf(post["y"][None, :], loc=mu, scale=np.sqrt(t["sigma2"])[:, None])

def waic(lp: np.ndarray) -> dict:
    """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) ;
    elpd_waic = lppd − p_waic ; WAIC = −2 elpd_waic."""
    S = lp.shape[0]
    lppd_i = np.logaddexp.reduce(lp, axis=0) - math.log(S)
    p_i = lp.var(axis=0, ddof=1)
    return {"lppd": float(lppd_i.sum()), "p_waic": float(p_i.sum()), "elpd_waic": float((lppd_i - p_i).sum()),
            "waic": float(-2 * (lppd_i - p_i).sum()), "se_elpd": float(math.sqrt(len(lppd_i)) * (lppd_i - p_i).std(ddof=1))}

def loo_tronque(lp: np.ndarray) -> dict:
    """LOO approximatif par échantillonnage préférentiel tronqué (Ionides 2008 ; forme simplifiée de PSIS, Vehtari et al. 2017) :
    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."""
    S = lp.shape[0]
    lw = -lp; lw -= lw.max(axis=0, keepdims=True)
    w = np.exp(lw); w /= w.sum(axis=0, keepdims=True)
    w = np.minimum(w, math.sqrt(S) / S); w /= w.sum(axis=0, keepdims=True)
    lp_loo = np.log((w * np.exp(lp)).sum(axis=0))
    return {"elpd_loo": float(lp_loo.sum()), "looic": float(-2 * lp_loo.sum()), "pointwise": lp_loo}

# =============================================================================================================================
# 3. Simulation à partir des lois du générateur (récupération des paramètres)
# =============================================================================================================================

def simuler_depuis_generateur(v: pd.DataFrame, n: int, graine: int, sigma: float = VRAIS["sigma"], const: float = 13.2) -> pd.DataFrame:
    """Rééchantillonne n lignes d'attributs du jeu (bootstrap des covariables) et engendre ln P selon les lois déclarées
    (README synthétique) : ln P = const + Σ β_vrai x + prime de secteur + ε, ε ~ N(0, σ²). Les attributs sont ceux du jeu ; le prix
    est régénéré. Sert aux tests de récupération des paramètres."""
    rng = np.random.default_rng(graine)
    d = v.sample(n=n, replace=True, random_state=int(rng.integers(1 << 31))).reset_index(drop=True)
    X = cas1.matrice(d)
    beta = pd.Series({c: VRAIS.get(c, 0.0) for c in X.columns}); beta["const"] = const
    lnp = X.values @ beta.values + rng.normal(0, sigma, n)
    d["prix_vente"] = np.round(np.exp(lnp) / 500) * 500
    return d

# =============================================================================================================================
# 4. Exemple intégral
# =============================================================================================================================

def sous_echantillons(v: pd.DataFrame, s: dict) -> dict:
    import ch16_selection as c16
    vp = c16.preparer(v.copy(), s)
    micro = vp[(vp["mois"].between(0, 18)) & (vp["dist_km"] <= 3.0)].drop(columns=["mois", "dist_km", "prix_m2", "distance"])
    return {"tout": v, "hull": v[v["secteur"] == "Hull"], "voisin": v[v["secteur"] != "Hull"], "micro": micro}

def colonnes_actives(X: pd.DataFrame) -> list:
    """Colonnes non constantes ; si le secteur de référence (Hull) est absent du sous-échantillon, les indicatrices restantes somment
    à 1 et l'une d'elles (la première) est retirée pour éviter la colinéarité avec la constante."""
    cols = [c for c in X.columns if c == "const" or X[c].std() > 0]
    sect = [c for c in cols if c.startswith("sect_")]
    if sect and np.allclose(X[sect].sum(axis=1), 1.0): cols.remove(sect[0])
    return cols

def ajuster_nig(df: pd.DataFrame, prior_fn=None, colonnes=None) -> dict:
    X = cas1.matrice(df); cols = colonnes or colonnes_actives(X); X = X[cols]
    y = np.log(df["prix_vente"].astype(float)).values
    prior = prior_fn(cols) if prior_fn else prior_reference(cols)
    p = nig_posterieure(X.values, y, prior); p["colonnes"] = cols
    return p

def x_sujet(s: dict, cols, secteur: str | None = None, dist_centre: float | None = None) -> np.ndarray:
    d = dict(s)
    if secteur: d["secteur"] = secteur
    if dist_centre is not None: d["dist_centre_km"] = dist_centre
    Xs = cas1.matrice(cas1.sujet_df(d))
    return Xs[cols].values[0]

def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:
    """Étapes 1-2 : postérieure NIG sur 600 ventes (a priori faible) ; petit échantillon (40 ventes à 3 km / 18 mois) avec a priori
    faible vs informatif emprunté aux autres secteurs (κ = 1, 4, 16) ; prédictive du sujet."""
    se = sous_echantillons(v, s)
    p600 = ajuster_nig(se["tout"]); cols = p600["colonnes"]          # a priori de référence
    p600_propre = ajuster_nig(se["tout"], prior_fn=prior_faible, colonnes=cols)   # a priori propre vague (comparaison)
    xs = x_sujet(s, cols)
    tab = pd.DataFrame({"mco": ref["coef"]["coef"].reindex(cols), "se_hc1": ref["coef"]["se_hc1"].reindex(cols),
                        "post_moy": p600["m_n"], "post_sd": p600["sd"], "ic_bas": p600["ic95_bas"], "ic_haut": p600["ic95_haut"],
                        "vrai": [VRAIS.get(c, np.nan) for c in cols]}, index=cols)
    tab["couvre"] = (tab["vrai"] >= tab["ic_bas"]) & (tab["vrai"] <= tab["ic_haut"])
    pred600 = predictive_t(p600, xs)
    # cohérence séquentielle : 2023-2024 puis 2025-2026
    anc = v[v["date_vente"] < "2025-01-01"]; rec = v[v["date_vente"] >= "2025-01-01"]
    p_anc = ajuster_nig(anc, colonnes=cols)
    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"}
    p_seq = ajuster_nig(rec, prior_fn=lambda c: prior_seq, colonnes=cols)
    # petit échantillon : 40 ventes autour du sujet ; a priori emprunté aux autres secteurs
    p_vois = ajuster_nig(se["voisin"])
    micro = se["micro"]; Xm = cas1.matrice(micro); cols_m = colonnes_actives(Xm)
    petits = {}
    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)]:
        p = ajuster_nig(micro, prior_fn=pf, colonnes=cols_m)
        pr = predictive_t(p, x_sujet(s, cols_m))
        petits[nom] = {"post": p, "pred": pr, "largeur95_pct": (pr["haut95"] - pr["bas95"]) / pr["mediane"] * 100}
    # Hull seul (186), a priori faible
    p_hull = ajuster_nig(se["hull"]); pred_hull = predictive_t(p_hull, x_sujet(s, p_hull["colonnes"]))
    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])
    pred_hull["largeur95_pct"] = (pred_hull["haut95"] - pred_hull["bas95"]) / pred_hull["mediane"] * 100
    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),
            "n_rec": len(rec), "p_vois": p_vois, "micro": micro, "cols_micro": cols_m, "petits": petits, "p_hull": p_hull,
            "pred_hull": pred_hull, "n_micro": len(micro)}

def exemple_hierarchique(v: pd.DataFrame, s: dict, ref: dict, n_iter: int = 3000, brulage: int = 1000) -> dict:
    """Étapes 3-4 : Gibbs sur les 600 ventes ; effets relatifs à Hull vs effets fixes (ch. 18) ; prédictive du sujet ;
    expérience Masson-Angers réduit à 5, 10, 20 ventes ; sensibilité à l'a priori sur τ."""
    Xc, moy, z = matrice_hier(v); y = np.log(v["prix_vente"].astype(float)).values; cols = list(Xc.columns)
    g = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage)                 # Gibbs par blocs (résultats du chapitre)
    g_naif = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, bloc=False)  # un paramètre à la fois (diagnostics)
    res = resume_gibbs(g, cols); rel = effets_relatifs(g); res_naif = resume_gibbs(g_naif, cols)
    diag = pd.DataFrame({"rhat_bloc": res["rhat"], "ess_bloc": res["ess"], "rhat_naif": res_naif["rhat"], "ess_naif": res_naif["ess"]})
    fixes = ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], ["coef", "se_hc1"]]; fixes.index = SECTEURS[1:]
    rel["fixe"] = fixes["coef"]; rel["fixe_se"] = fixes["se_hc1"]; rel["vrai"] = [VRAIS[f"sect_{x}"] for x in rel.index]
    # prédictive du sujet (Hull) : ln P* = (x_s − x̄)'β + α_Hull + ε
    xs = (cas1.matrice(cas1.sujet_df(s))[cols].iloc[0] - moy).values
    beta = g["beta"].reshape(-1, len(cols)); alpha = g["alpha"].reshape(-1, 5); s2 = g["sigma2"].reshape(-1)
    rng = np.random.default_rng(GRAINE)
    lnp = beta @ xs + alpha[:, 0] + rng.standard_normal(len(s2)) * np.sqrt(s2)
    P = np.exp(lnp)
    pred = {"tirages": P, "mediane": float(np.median(P)), "moyenne": float(P.mean()), "moyenne_log": float(lnp.mean()),
            "bas90": float(np.quantile(P, 0.05)), "haut90": float(np.quantile(P, 0.95)),
            "bas95": float(np.quantile(P, 0.025)), "haut95": float(np.quantile(P, 0.975)),
            "sd_log": float(lnp.std()), "sd_moyenne_log": float((beta @ xs + alpha[:, 0]).std())}
    # rétrécissement a posteriori : facteur B_s à la médiane de σ², τ²
    s2m = float(np.median(g["sigma2"])); t2m = float(np.median(g["tau2"]))
    B = {sct: facteur_retrecissement(g["meta"]["ns"][j], s2m, t2m) for j, sct in enumerate(SECTEURS)}
    # expérience Masson-Angers : 5, 10, 20, 47 ventes (sous-ensembles emboîtés)
    ma = v[v["secteur"] == "Masson-Angers"]; autres = v[v["secteur"] != "Masson-Angers"]
    perm = np.random.default_rng(GRAINE).permutation(len(ma))
    dist_ma = float(ma["dist_centre_km"].median())
    exp_rows = []
    for n_ma in (5, 10, 20, len(ma)):
        vv = pd.concat([autres, ma.iloc[perm[:n_ma]]]).reset_index(drop=True)
        # effets fixes (MCO, HC1) et prédictive conjuguée (plat) pour un sujet hypothétique à Masson-Angers
        pf = ajuster_nig(vv); i = pf["colonnes"].index("sect_Masson-Angers")
        xma = x_sujet(s, pf["colonnes"], secteur="Masson-Angers", dist_centre=dist_ma)
        pr_f = predictive_t(pf, xma)
        # hiérarchique
        Xv, moy_v, zv = matrice_hier(vv); yv = np.log(vv["prix_vente"].astype(float)).values
        gv = gibbs_hierarchique(yv, Xv.values, zv, n_iter=n_iter, brulage=brulage)
        rr = effets_relatifs(gv).loc["Masson-Angers"]
        xma_c = (cas1.matrice(cas1.sujet_df({**s, "secteur": "Masson-Angers", "dist_centre_km": dist_ma}))[cols].iloc[0] - moy_v).values
        bv = gv["beta"].reshape(-1, len(cols)); av = gv["alpha"].reshape(-1, 5); sv = gv["sigma2"].reshape(-1)
        lnv = bv @ xma_c + av[:, 4] + np.random.default_rng(GRAINE + n_ma).standard_normal(len(sv)) * np.sqrt(sv)
        Pv = np.exp(lnv)
        s2v = float(np.median(gv["sigma2"])); t2v = float(np.median(gv["tau2"]))
        # moyenne brute des résidus du secteur (sans rétrécissement) à la postérieure : r̄ − α0 relatif
        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],
                         "hier": rr["moyenne"], "hier_sd": rr["sd"], "hier_bas": rr["q025"], "hier_haut": rr["q975"], "hier_rhat": rr["rhat"],
                         "tau_med": math.sqrt(t2v), "B": facteur_retrecissement(n_ma, s2v, t2v),
                         "pred_f_med": pr_f["mediane"], "pred_f_bas": pr_f["bas95"], "pred_f_haut": pr_f["haut95"],
                         "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)),
                         "vrai": VRAIS["sect_Masson-Angers"]})
    experience = pd.DataFrame(exp_rows).set_index("n_ma")
    experience["largeur_f_pct"] = (experience["pred_f_haut"] - experience["pred_f_bas"]) / experience["pred_f_med"] * 100
    experience["largeur_h_pct"] = (experience["pred_h_haut"] - experience["pred_h_bas"]) / experience["pred_h_med"] * 100
    # sensibilité à l'a priori sur τ (600 ventes et n_MA = 5)
    v5 = pd.concat([autres, ma.iloc[perm[:5]]]).reset_index(drop=True)
    X5, _, z5 = matrice_hier(v5); y5 = np.log(v5["prix_vente"].astype(float)).values
    tau_eb = float(np.std(ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], "coef"].tolist() + [0.0], ddof=1))
    sens_rows = []
    for nom, kw in [("uniforme sur τ", {"prior_tau": "uniforme"}), ("IG(0,001 ; 0,001) sur τ²", {"prior_tau": "ig"}),
                    ("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})]:
        g600 = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, **kw)
        g5 = gibbs_hierarchique(y5, X5.values, z5, n_iter=n_iter, brulage=brulage, **kw)
        r600 = effets_relatifs(g600).loc["Masson-Angers"]; r5 = effets_relatifs(g5).loc["Masson-Angers"]
        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)),
                          "tau_q95_600": float(np.quantile(np.sqrt(g600["tau2"]), 0.95)),
                          "ma_600": r600["moyenne"], "ma_600_sd": r600["sd"], "tau_med_5": float(np.median(np.sqrt(g5["tau2"]))),
                          "ma_5": r5["moyenne"], "ma_5_sd": r5["sd"], "ma_5_bas": r5["q025"], "ma_5_haut": r5["q975"]})
    sensibilite = pd.DataFrame(sens_rows).set_index("prior")
    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,
            "sigma2_med": s2m, "tau2_med": t2m, "experience": experience, "dist_ma": dist_ma, "sensibilite": sensibilite, "tau_eb": tau_eb,
            "perm_ma": perm}

def criteres(v: pd.DataFrame, hier: dict, ex: dict) -> pd.DataFrame:
    """Comparaison de trois modèles : pooling complet (sans secteur), effets fixes (indicatrices), hiérarchique."""
    y = hier["y"]
    # pooling complet : matrice de référence sans indicatrices de secteur
    X = cas1.matrice(v); cols_pool = [c for c in X.columns if not c.startswith("sect_")]
    p_pool = nig_posterieure(X[cols_pool].values, y, prior_reference(cols_pool)); p_pool["colonnes"] = cols_pool
    p_fixe = ex["p600"]
    lignes = []
    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"])),
                         ("effets fixes de secteur (ch. 18)", log_vrais_pointwise_nig(p_fixe), loo_exact_nig(p_fixe["X"], y, p_fixe["prior"])),
                         ("effets aléatoires de secteur (hiérarchique)", log_vrais_pointwise_hier(hier["g"], y, hier["Xc"].values, hier["z"]), None)]:
        w = waic(lp); lt = loo_tronque(lp)
        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"],
                       "elpd_loo_tronque": lt["elpd_loo"], "elpd_loo_exact": loo["elpd_loo"] if loo else np.nan})
    t = pd.DataFrame(lignes).set_index("modele")
    t["d_elpd_waic"] = t["elpd_waic"] - t["elpd_waic"].max()
    return t

# =============================================================================================================================
# 5. Exercices (chiffres des solutions)
# =============================================================================================================================

def exercices(ex: dict, hier: dict) -> dict:
    r = {}
    # ex 1 : mise à jour normale-normale, variance connue (prix au m² d'un marché mince)
    m0, s0, sig, n, xbar = 3300.0, 200.0, 300.0, 4, 3520.0
    prec = 1 / s0 ** 2 + n / sig ** 2; mn = (m0 / s0 ** 2 + n * xbar / sig ** 2) / prec; sn = math.sqrt(1 / prec)
    w = (n / sig ** 2) / prec
    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),
                "n_equivalent_prior": sig ** 2 / s0 ** 2, "valeur_165": mn * 165, "sd_165": sn * 165}
    prec4 = 1 / 400.0 ** 2 + n / sig ** 2
    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,
                          "n_equivalent_prior": sig ** 2 / 400.0 ** 2}
    # ex 2 : lecture d'intervalles — probabilité a posteriori que la valeur dépasse 650 000 $ sous une postérieure log-normale
    mu, sd = ex["pred600"]["loc"], ex["pred600"]["echelle"]
    r["ex2"] = {"p_sup_650": float(1 - stats.t.cdf((math.log(650_000) - mu) / sd, ex["pred600"]["nu"])),
                "p_inf_550": float(stats.t.cdf((math.log(550_000) - mu) / sd, ex["pred600"]["nu"])),
                "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),
                "z_650": (math.log(650_000) - mu) / sd, "z_550": (math.log(550_000) - mu) / sd}
    # ex 3 : facteur de rétrécissement
    s2, t2 = 0.083 ** 2, 0.15 ** 2
    r["ex3"] = {n_: facteur_retrecissement(n_, s2, t2) for n_ in (5, 10, 20, 47)}
    r["ex3"]["alpha_5"] = facteur_retrecissement(5, s2, t2) * (-0.30) + (1 - facteur_retrecissement(5, s2, t2)) * (-0.08)
    r["ex3"]["tau_005"] = {n_: facteur_retrecissement(n_, s2, 0.05 ** 2) for n_ in (5, 10, 20, 47)}
    r["ex3"]["sd_cond_5"] = math.sqrt(1 / (5 / s2 + 1 / t2)); r["ex3"]["sd_fixe_5"] = math.sqrt(s2 / 5)
    r["ex3"]["n_demi"] = {"tau_015": s2 / t2, "tau_005": s2 / 0.05 ** 2}
    # ex 4 : prédictive Student-t
    p = ex["p600"]; xs = ex["xs"]; q = float(xs @ p["V_n"] @ xs)
    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"],
                "t975": float(stats.t.ppf(0.975, p["nu"])), "z975": 1.959964, "bas95": ex["pred600"]["bas95"], "haut95": ex["pred600"]["haut95"],
                "nu_petit": 20, "t975_20": float(stats.t.ppf(0.975, 20)), "mediane": ex["pred600"]["mediane"],
                "bas95_normal": math.exp(ex["pred600"]["loc"] - 1.959964 * ex["pred600"]["echelle"]), "haut95_normal": math.exp(ex["pred600"]["loc"] + 1.959964 * ex["pred600"]["echelle"]),
                "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"]),
                "elargissement_pct": (float(stats.t.ppf(0.975, 20)) / float(stats.t.ppf(0.975, p["nu"])) - 1) * 100}
    # ex 5 : R̂ sur des moyennes/variances de chaînes données
    moy = np.array([-0.262, -0.258, -0.271, -0.249]); var = np.array([0.0049, 0.0052, 0.0047, 0.0055]); n_ = 1000
    B = n_ * moy.var(ddof=1); W = var.mean(); vp = (n_ - 1) / n_ * W + B / n_
    r["ex5"] = {"B": B, "W": W, "var_plus": vp, "rhat": math.sqrt(vp / W)}
    moy2 = np.array([-0.262, -0.258, -0.371, -0.249]); B2 = n_ * moy2.var(ddof=1); vp2 = (n_ - 1) / n_ * W + B2 / n_
    r["ex5"]["rhat_divergent"] = math.sqrt(vp2 / W); r["ex5"]["B2"] = B2
    r["ex5"]["ess_rho08"] = 4000 / (1 + 2 * 0.8 / (1 - 0.8))
    # ex 6 : script (κ) — valeurs du tableau des petits échantillons
    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")],
                    "mediane": v_["pred"]["mediane"], "largeur95_pct": v_["largeur95_pct"]} for k, v_ in ex["petits"].items()}
    # ex 7 : conditionnelle de α0 avec a priori N(m, s²) — valeurs numériques à τ et α médians
    g = hier["g"]; alpha_bar = float(g["alpha"].reshape(-1, 5).mean()); tau2 = hier["tau2_med"]
    for (m_, s_) in [(13.0, 1.0), (13.0, 0.05)]:
        prec = 5 / tau2 + 1 / s_ ** 2
        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)}
    return r

# =============================================================================================================================
# 6. Main
# =============================================================================================================================

def main(verbose: bool = True, n_iter: int = 3000, brulage: int = 1000) -> dict:
    d = cas1.charger(); v = d["ventes"]; s = d["sujet"]
    ref = cas1.hedonique(v)
    ex = exemple_conjugue(v, s, ref)
    hier = exemple_hierarchique(v, s, ref, n_iter=n_iter, brulage=brulage)
    crit = criteres(v, hier, ex)
    exo = exercices(ex, hier)
    out = {"ventes": v, "sujet": s, "modele": ref, "conjugue": ex, "hier": hier, "criteres": crit, "exercices": exo}
    if verbose: imprimer(out)
    return out

def _f(x, d=4): return f"{x:.{d}f}"
def _d(x): return f"{x:,.0f} $".replace(",", " ")

def imprimer(r: dict) -> None:
    ex = r["conjugue"]; hi = r["hier"]
    print("=== CH39 — modèle conjugué, 600 ventes, a priori faiblement informatif ===")
    t = ex["tab600"]
    print(t.round(4).to_string())
    print(f"couverture des vrais paramètres (IC 95 %) : {int(t['couvre'].sum())}/{int(t['vrai'].notna().sum())}")
    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}")
    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'])}]")
    pr = ex["pred600"]
    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} %")
    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'])}])")
    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}")
    print(f"\n=== petit échantillon : {ex['n_micro']} ventes à 3 km / 18 mois ; colonnes {ex['cols_micro']}")
    for nom, pe in ex["petits"].items():
        po = pe["post"]; i = ex["cols_micro"].index("ln_S"); iq = ex["cols_micro"].index("qualite"); ig = ex["cols_micro"].index("garage")
        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}) ; "
              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} %")
    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} %")
    print("\n=== modèle hiérarchique (Gibbs par blocs) : résumé ===")
    print(hi["resume"].round(4).to_string())
    print("\ndiagnostics : Gibbs par blocs vs un paramètre à la fois"); print(hi["diag"].round(3).to_string())
    print("\neffets relatifs à Hull (hiérarchique vs effets fixes ch. 18 vs vrai) :"); print(hi["relatifs"].round(4).to_string())
    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})")
    print("facteurs de rétrécissement B_s :", {k: round(v_, 4) for k, v_ in hi["B"].items()})
    pd_ = hi["pred"]
    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})")
    print(f"\n=== expérience Masson-Angers (dist. médiane {hi['dist_ma']:.1f} km) ===")
    print(hi["experience"].round(4).to_string())
    print("\n=== sensibilité à l'a priori sur τ ==="); print(hi["sensibilite"].round(4).to_string())
    print("\n=== critères ==="); print(r["criteres"].round(2).to_string())
    print("\n=== exercices ===")
    for k, v_ in r["exercices"].items(): print(k, v_)

if __name__ == "__main__":
    main()