spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 39 — Modèles bayésiens et intervalles (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).3Produit dans assets/figures/ch39/ :4  ch39-fig-posterieures.pdf      densités marginales a posteriori (modèle conjugué, 600 ventes) vs MCO ± 1,96 e.-t. et paramètre vrai5  ch39-fig-prior-informatif.pdf  petit échantillon (40 ventes) : postérieures sous a priori de référence et a priori « marché voisin »6  ch39-fig-traces.pdf            traces de quatre chaînes, Gibbs un paramètre à la fois vs par blocs, avec R̂ et ESS7  ch39-fig-retrecissement.pdf    effet de Masson-Angers selon le nombre de ventes (effets fixes vs hiérarchique) et facteur B(n)8  ch39-fig-predictive.pdf        distribution prédictive du prix du sujet (hiérarchique et conjuguée) vs IC fréquentiste du ch. 189Usage (racine) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch39_figures.py"""10from __future__ import annotations11import sys, pathlib, math12ROOT = pathlib.Path(__file__).resolve().parents[3]13sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))14import numpy as np15import matplotlib.pyplot as plt16from matplotlib.ticker import FuncFormatter17from scipy import stats18from traite.style import setup, palette, savefig, LARGEUR19import ch39_bayes as c392021def _k(x, pos=None): return f"{x/1000:,.0f}".replace(",", " ")22ETIQ = {"ln_S": "Élasticité superficie habitable", "qualite": "Échelon de qualité", "garage": "Place de garage",23        "sect_Masson-Angers": "Secteur Masson-Angers (vs Hull)", "dist_centre": "Distance au centre (> 3 km, par km)", "sigma": "Écart-type résiduel σ"}2425def densite_t(post, j, grille):26    return stats.t.pdf(grille, post["nu"], loc=post["m_n"][j], scale=post["echelle"][j])2728def densite_sigma(post, grille):29    # σ² ~ IG(a_n, b_n) → densité de σ : f(σ) = 2σ · IG(σ²)30    return 2 * grille * stats.invgamma.pdf(grille ** 2, post["a_n"], scale=post["b_n"])3132def fig_posterieures(r):33    ex = r["conjugue"]; p = ex["p600"]; cols = ex["cols"]; coef = r["modele"]["coef"]34    fig, axes = plt.subplots(2, 3, figsize=(LARGEUR, 3.9))35    for ax, nom in zip(axes.ravel(), ["ln_S", "qualite", "garage", "sect_Masson-Angers", "dist_centre", "sigma"]):36        if nom == "sigma":37            g = np.linspace(0.074, 0.093, 300); d = densite_sigma(p, g); vrai = c39.VRAIS["sigma"]38            mco, se = r["modele"]["sigma"], None39        else:40            j = cols.index(nom); m, s = p["m_n"][j], p["echelle"][j]41            g = np.linspace(m - 4.2 * s, m + 4.2 * s, 300); d = densite_t(p, j, g); vrai = c39.VRAIS[nom]42            mco, se = coef.loc[nom, "coef"], coef.loc[nom, "se_hc1"]43        ax.fill_between(g, d, color=palette.cat[0], alpha=0.18, lw=0); ax.plot(g, d, color=palette.cat[0], lw=1.6, label="Postérieure (a priori de référence)")44        ax.axvline(vrai, color=palette.cat[2], ls="--", lw=1.3, label="Paramètre de génération")45        ymax = d.max()46        if se is not None:47            ax.errorbar([mco], [-0.06 * ymax], xerr=[[1.96 * se], [1.96 * se]], fmt="o", color=palette.ink, ms=3.5, capsize=2, lw=1.1, label="MCO ± 1,96 e.-t. HC1")48        else:49            ax.plot([mco], [-0.06 * ymax], "o", color=palette.ink, ms=3.5)50        ax.set_ylim(-0.14 * ymax, 1.12 * ymax); ax.set_yticks([]); ax.set_title(ETIQ[nom], fontsize=8.5)51        ax.grid(False); ax.spines["left"].set_visible(False)52    h, l = axes[0, 0].get_legend_handles_labels()53    fig.legend(h, l, loc="lower center", ncol=3, bbox_to_anchor=(0.5, -0.04), frameon=False)54    fig.tight_layout(); savefig(fig, "ch39/ch39-fig-posterieures.pdf")5556def fig_prior_informatif(r):57    ex = r["conjugue"]; cols_m = ex["cols_micro"]; cols = ex["cols"]; p600 = ex["p600"]58    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.0))59    styles = {"faible": ("Référence (plat), 40 ventes", palette.cat[1], "-"), "kappa_1": ("Marché voisin, κ = 1", palette.cat[0], "-"),60              "kappa_4": ("Marché voisin, κ = 4", palette.cat[0], "--"), "kappa_16": ("Marché voisin, κ = 16", palette.cat[0], ":")}61    for ax, nom in zip(axes, ["ln_S", "qualite"]):62        j6 = cols.index(nom); m6, s6 = p600["m_n"][j6], p600["echelle"][j6]63        jm = cols_m.index(nom); pf = ex["petits"]["faible"]["post"]64        lo = min(m6 - 3 * s6, pf["m_n"][jm] - 3 * pf["echelle"][jm]); hi = max(m6 + 3 * s6, pf["m_n"][jm] + 3 * pf["echelle"][jm])65        g = np.linspace(lo, hi, 400)66        ax.plot(g, densite_t(p600, j6, g), color=palette.ink2, lw=1.2, alpha=0.8, label="600 ventes (référence)")67        for k_, (lab, col, ls) in styles.items():68            ax.plot(g, densite_t(ex["petits"][k_]["post"], jm, g), color=col, ls=ls, lw=1.6, label=lab)69        ax.axvline(c39.VRAIS[nom], color=palette.cat[2], ls="--", lw=1.2, label="Paramètre de génération")70        ax.set_title(ETIQ[nom], fontsize=8.5); ax.set_yticks([]); ax.grid(False); ax.spines["left"].set_visible(False)71    h, l = axes[0].get_legend_handles_labels()72    fig.legend(h, l, loc="lower center", ncol=3, bbox_to_anchor=(0.5, -0.1), frameon=False)73    fig.tight_layout(); savefig(fig, "ch39/ch39-fig-prior-informatif.pdf")7475def fig_traces(r):76    hi = r["hier"]; cols = hi["cols"]; j = cols.index("dist_centre"); jma = c39.SECTEURS.index("Masson-Angers")77    fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 3.9), sharex=True)78    T = 100079    for row, (g, nom) in enumerate([(hi["g_naif"], "Gibbs, un paramètre à la fois"), (hi["g"], "Gibbs par blocs (β, α ensemble)")]):80        for col, (x, lab) in enumerate([(g["alpha"][:, :, jma] - g["alpha"][:, :, 0], "Effet de Masson-Angers (α_MA − α_Hull)"), (g["beta"][:, :, j], "Distance au centre (par km)")]):81            ax = axes[row, col]82            for c in range(x.shape[0]): ax.plot(np.arange(T), x[c, :T], lw=0.5, alpha=0.8, color=palette.cat[c % 6])83            rh = c39.rhat_scinde(x); es = c39.ess(x)84            ax.text(0.98, 0.95, f"R̂ = {rh:.3f}   ESS = {es:,.0f}".replace(",", " "), transform=ax.transAxes, ha="right", va="top", fontsize=7.5,85                    bbox=dict(facecolor=palette.surface, edgecolor="none", pad=1.5))86            if row == 0: ax.set_title(lab, fontsize=8.5)87            if col == 0: ax.set_ylabel(nom, fontsize=7.5)88    for ax in axes[1]: ax.set_xlabel("Itération (après brûlage)")89    fig.tight_layout(); savefig(fig, "ch39/ch39-fig-traces.pdf")9091def fig_retrecissement(r):92    hi = r["hier"]; e = hi["experience"]; n = e.index.values.astype(float); pos = np.arange(len(n), dtype=float)93    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.3), gridspec_kw={"width_ratios": [1.1, 1.0]})94    ax = axes[0]95    ax.errorbar(pos - 0.08, e["fixe"], yerr=[e["fixe"] - e["fixe_bas"], e["fixe_haut"] - e["fixe"]], fmt="s", color=palette.cat[1], capsize=2.5, lw=1.2, ms=4, label="Effet fixe (MCO), IC 95 %")96    ax.errorbar(pos + 0.08, e["hier"], yerr=[e["hier"] - e["hier_bas"], e["hier_haut"] - e["hier"]], fmt="o", color=palette.cat[0], capsize=2.5, lw=1.2, ms=4, label="Effet aléatoire (Gibbs), crédibilité 95 %")97    ax.axhline(c39.VRAIS["sect_Masson-Angers"], color=palette.cat[2], ls="--", lw=1.2, label="Paramètre de génération")98    a0 = float(hi["g"]["alpha0"].mean() - hi["g"]["alpha"][:, :, 0].mean())99    ax.axhline(a0, color=palette.ink2, ls=":", lw=1.1, label="Moyenne des secteurs α₀ − α_Hull")100    ax.set_xticks(pos); ax.set_xticklabels([f"{int(x)}" for x in n]); ax.set_xlim(-0.5, len(n) - 0.5)101    ax.set_xlabel("Nombre de ventes conservées à Masson-Angers"); ax.set_ylabel("Effet relatif à Hull (log)")102    ax.set_ylim(-0.56, 0.22); ax.legend(fontsize=6.3, loc="upper left", ncol=1, handlelength=1.6)103    ax = axes[1]104    grille = np.linspace(1, 60, 300); s2 = hi["sigma2_med"]105    for i, (tau, lab) in enumerate([(0.05, "τ = 0,05"), (0.10, "τ = 0,10"), (math.sqrt(hi["tau2_med"]), f"τ = {math.sqrt(hi['tau2_med']):.2f} (médiane a posteriori)".replace(".", ","))]):106        ax.plot(grille, [c39.facteur_retrecissement(x, s2, tau ** 2) for x in grille], color=palette.seq[2 + 2 * i] if i < 2 else palette.cat[0], lw=1.6, label=lab)107    ax.plot(n, e["B"], "o", color=palette.cat[0], ms=4, label="Expérience (B à τ médian)")108    ax.set_ylim(0.3, 1.02); ax.set_xlim(0, 60); ax.set_xlabel("Ventes dans le secteur $n_s$"); ax.set_ylabel("Facteur de rétrécissement $B_s$")109    ax.legend(fontsize=6.3, loc="lower right", handlelength=1.6)110    fig.tight_layout(); savefig(fig, "ch39/ch39-fig-retrecissement.pdf")111112def fig_predictive(r):113    hi = r["hier"]; pr = hi["pred"]; ex = r["conjugue"]; pt = ex["pred600"]114    P = pr["tirages"]115    fig, ax = plt.subplots(figsize=(LARGEUR, 3.2))116    g = np.linspace(430_000, 900_000, 600)117    kde = stats.gaussian_kde(P, bw_method=0.12)118    ax.fill_between(g, kde(g), color=palette.cat[0], alpha=0.15, lw=0)119    ax.plot(g, kde(g), color=palette.cat[0], lw=1.8, label="Prédictive hiérarchique (tirages de Gibbs)")120    # prédictive conjuguée : log-t → densité en dollars121    dens_t = stats.t.pdf((np.log(g) - pt["loc"]) / pt["echelle"], pt["nu"]) / (pt["echelle"] * g)122    ax.plot(g, dens_t, color=palette.cat[1], ls="--", lw=1.5, label="Prédictive conjuguée (Student-t, 584 ddl)")123    for lo, hi_, lab, al in [(pr["bas95"], pr["haut95"], "Crédibilité 95 %", 0.10), (pr["bas90"], pr["haut90"], "Crédibilité 90 %", 0.16)]:124        m = (g >= lo) & (g <= hi_); ax.fill_between(g[m], kde(g[m]), color=palette.cat[0], alpha=al, lw=0)125    ymax = kde(g).max()126    for x, lab, col, ls in [(c39.IC18["bas95"], "IC fréquentiste 95 % (ch. 18)", palette.ink, "-."), (c39.IC18["haut95"], None, palette.ink, "-."),127                            (c39.IC18["point"], "Prédiction de Duan 637 953 $", palette.cat[2], ":")]:128        ax.axvline(x, color=col, ls=ls, lw=1.1, label=lab)129    ax.annotate(f"médiane a posteriori {_k(pr['mediane'])} k$", (pr["mediane"] - 6000, ymax * 1.06), ha="right", fontsize=7.5)130    ax.annotate("90 %", ((pr["bas90"] + pr["haut90"]) / 2, ymax * 0.35), ha="center", fontsize=7.5, color=palette.ink)131    ax.annotate(f"{_k(pr['bas90'])}", (pr["bas90"], ymax * 0.05), ha="right", fontsize=7, color=palette.ink2)132    ax.annotate(f"{_k(pr['haut90'])}", (pr["haut90"], ymax * 0.05), ha="left", fontsize=7, color=palette.ink2)133    ax.set_ylim(0, ymax * 1.15); ax.set_yticks([]); ax.grid(False); ax.spines["left"].set_visible(False)134    ax.xaxis.set_major_formatter(FuncFormatter(_k)); ax.set_xlabel("Prix du sujet (milliers de dollars)")135    h, l = ax.get_legend_handles_labels()136    fig.legend(h, l, loc="lower center", ncol=2, fontsize=7.5, bbox_to_anchor=(0.5, -0.1), frameon=False)137    fig.tight_layout(); savefig(fig, "ch39/ch39-fig-predictive.pdf")138139def main():140    setup()141    r = c39.main(verbose=False)142    fig_posterieures(r); fig_prior_informatif(r); fig_traces(r); fig_retrecissement(r); fig_predictive(r)143    print("figures ch39 écrites dans assets/figures/ch39/")144    return r145146if __name__ == "__main__":147    main()148