spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 37 — Apprentissage automatique et explicabilité (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).3Produit dans assets/figures/ch37/ :4  ch37-fig-arbre.pdf          arbre de régression de profondeur 3 estimé sur les 600 ventes (feuilles en $)5  ch37-fig-apprentissage.pdf  erreur d'entraînement et de validation croisée du boosting selon le nombre d'arbres et la profondeur6  ch37-fig-erreurs.pdf        distribution des erreurs hors échantillon (log, en %) par modèle et par schéma de validation7  ch37-fig-predit-observe.pdf prix prédit hors échantillon (plis aléatoires) contre prix observé : hédonique, forêt, boosting8  ch37-fig-shap-beeswarm.pdf  valeurs SHAP des 600 ventes (boosting) : importance et sens des effets9  ch37-fig-pdp.pdf            dépendance partielle : superficie et âge (hédonique, forêt, boosting, loi de génération) ; âge dans le scénario non linéaire10  ch37-fig-waterfall.pdf      cascade SHAP de la prédiction du sujet par le boosting (en dollars)11Usage (racine) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch37_figures.py"""12from __future__ import annotations13import sys, pathlib, math, re14ROOT = pathlib.Path(__file__).resolve().parents[3]15sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))16import numpy as np, pandas as pd17import matplotlib.pyplot as plt18from matplotlib.ticker import FuncFormatter19from traite.style import setup, palette, savefig, LARGEUR, fmt_milliers20import ch37_ml as c372122COURT = {"Hédonique (référence)": "Hédonique", "Arbre (profondeur 3)": "Arbre", "Forêt aléatoire": "Forêt", "Gradient boosting": "Boosting", "XGBoost": "XGBoost",23         "Hédonique respécifié": "Hédonique respécifié"}24COUL = {"Hédonique (référence)": palette.cat[0], "Arbre (profondeur 3)": palette.cat[4], "Forêt aléatoire": palette.cat[1], "Gradient boosting": palette.cat[2], "XGBoost": palette.cat[3]}25def _k(x, pos=None): return f"{x/1000:,.0f}".replace(",", " ")2627def fig_arbre(res):28    from sklearn.tree import plot_tree29    arbre = res["ajustes"]["Arbre (profondeur 3)"]30    noms = [c37.ETIQUETTES[c].replace(" (m²)", " (m²)") for c in res["Xm"].columns]31    fig, ax = plt.subplots(figsize=(LARGEUR, 3.6))32    plot_tree(arbre, feature_names=noms, filled=False, impurity=False, rounded=True, fontsize=6, ax=ax, precision=3)33    for t in ax.texts:34        s = t.get_text()35        s = re.sub(r"samples = (\d+)", lambda m: f"n = {m.group(1)}", s)36        s = re.sub(r"value = (-?\d+\.\d+)", lambda m: f"ln P = {float(m.group(1)):.3f}\n≈ {math.exp(float(m.group(1)))/1000:,.0f} k$".replace(",", " "), s)37        s = s.replace(" <= ", " ≤ ")38        t.set_text(s); t.set_color(palette.ink)39        bb = t.get_bbox_patch()40        if bb is not None: bb.set_edgecolor(palette.ink2); bb.set_facecolor(palette.surface); bb.set_linewidth(0.6)41    ax.set_axis_off(); ax.grid(False)42    savefig(fig, "ch37/ch37-fig-arbre.pdf")4344def fig_apprentissage(res):45    c = res["courbes"]; ref = res["comparaison"]["metriques"]["Hédonique (référence)"]["aleatoire"]["rmse_log"]46    fig, ax = plt.subplots(figsize=(LARGEUR, 3.1))47    it = np.arange(1, 1201)48    for i, (d, cc) in enumerate(sorted(c.items())):49        ax.plot(it, cc["validation"], color=palette.cat[i], lw=1.6, label=f"validation croisée, profondeur {d}")50        ax.plot(it, cc["entrainement"], color=palette.cat[i], lw=1.0, ls="--", alpha=0.85)51        ax.plot([cc["optimum"]], [cc["rmse_min"]], marker="o", ms=4, color=palette.cat[i])52    ax.axhline(ref, color=palette.ink2, lw=0.9, ls=":")53    ax.text(1180, ref + 0.002, f"hédonique de référence (validation croisée) : {ref:.4f}", ha="right", va="bottom", fontsize=7.5, color=palette.ink2)54    ax.plot([], [], color=palette.ink2, lw=1.0, ls="--", label="erreur d'entraînement (tiretés)")55    ax.set_xlabel("Nombre d'arbres (itérations du boosting, ν = 0,03)"); ax.set_ylabel("Racine de l'erreur quadratique moyenne (log du prix)")56    ax.set_ylim(0.0, 0.20); ax.set_xlim(0, 1200); ax.legend(loc="upper right", ncol=1)57    savefig(fig, "ch37/ch37-fig-apprentissage.pdf")5859def fig_erreurs(res):60    comp = res["comparaison"]; y = res["y"]; noms = list(c37.MODELES)61    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 3.0), sharey=True)62    titres = {"aleatoire": "Plis aléatoires (10)", "spatial": "Blocs spatiaux par secteur (10)", "temporel": "Prospectif : ventes de 2026"}63    for ax, sch in zip(axes, ["aleatoire", "spatial", "temporel"]):64        data = []65        for nom in noms:66            p = comp["pred"][nom][sch]; m = ~np.isnan(p); data.append((np.exp(p[m] - y[m]) - 1) * 100)67        bp = ax.boxplot(data, widths=0.6, patch_artist=True, showfliers=True, flierprops=dict(marker=".", ms=2.5, markeredgecolor=palette.ink2, alpha=0.5),68                        medianprops=dict(color=palette.ink, lw=1.2), whiskerprops=dict(color=palette.ink2, lw=0.8), capprops=dict(color=palette.ink2, lw=0.8))69        for patch, nom in zip(bp["boxes"], noms): patch.set_facecolor(COUL[nom]); patch.set_alpha(0.55); patch.set_edgecolor(palette.ink2)70        ax.axhline(0, color=palette.ink2, lw=0.8, ls=":")71        ax.set_xticks(range(1, len(noms) + 1)); ax.set_xticklabels([COURT[n] for n in noms], rotation=35, ha="right")72        ax.set_title(titres[sch], fontsize=8.5)73        for i, nom in enumerate(noms):74            mm = comp["metriques"][nom][sch]; ax.text(i + 1, 43, f"COD\n{mm['COD']:.1f}".replace(".", ","), ha="center", va="top", fontsize=6.8, color=palette.ink)75    axes[0].set_ylabel("Erreur de prédiction hors échantillon (%)"); axes[0].set_ylim(-45, 45)76    savefig(fig, "ch37/ch37-fig-erreurs.pdf")7778def fig_predit_observe(res):79    comp = res["comparaison"]; P = np.exp(res["y"])80    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 2.9), sharex=True, sharey=True)81    for ax, nom in zip(axes, ["Hédonique (référence)", "Forêt aléatoire", "Gradient boosting"]):82        pred = np.exp(comp["pred"][nom]["aleatoire"]) * comp["duan"][nom]83        ax.scatter(P, pred, s=7, color=COUL[nom], alpha=0.55, edgecolors="none")84        lim = [180_000, 900_000]; ax.plot(lim, lim, color=palette.ink2, lw=0.9)85        ax.plot(lim, [l * 1.18 for l in lim], color=palette.ink2, lw=0.6, ls="--"); ax.plot(lim, [l / 1.18 for l in lim], color=palette.ink2, lw=0.6, ls="--")86        m = comp["metriques"][nom]["aleatoire"]87        ax.set_title(COURT[nom], fontsize=8.5)88        ax.text(0.03, 0.97, f"RMSE log {m['rmse_log']:.3f}\nMdAPE {m['mdape']:.1f} %\nCOD {m['COD']:.2f}".replace(".", ","), transform=ax.transAxes, va="top", ha="left", fontsize=7.2)89        ax.set_xlim(lim); ax.set_ylim(lim); ax.set_xlabel("Prix observé (k$)")90        ax.xaxis.set_major_formatter(FuncFormatter(_k)); ax.yaxis.set_major_formatter(FuncFormatter(_k))91    axes[0].set_ylabel("Prix prédit hors échantillon (k$)")92    savefig(fig, "ch37/ch37-fig-predit-observe.pdf")9394def fig_beeswarm(res, top: int = 11):95    sh = res["shap"]; Xm = res["Xm"]; sv = sh["valeurs"]; imp = sh["importance"]96    cols = list(imp.index[:top]); rng = np.random.default_rng(20261031)97    fig, ax = plt.subplots(figsize=(LARGEUR, 3.7)); cmap = palette.cmap_seq()98    for i, c in enumerate(cols[::-1]):99        j = list(Xm.columns).index(c); x = sv[:, j] * 100; val = Xm[c].to_numpy()100        lo, hi = np.nanpercentile(val, [5, 95]); col = np.clip((val - lo) / (hi - lo if hi > lo else 1), 0, 1)101        # essaim : décalage vertical proportionnel à la densité locale102        ordre = np.argsort(x); yj = np.zeros(len(x)); bins = np.digitize(x, np.linspace(x.min(), x.max() + 1e-9, 60))103        for b in np.unique(bins):104            idx = np.where(bins == b)[0]; k = len(idx); yj[idx] = (np.arange(k) - (k - 1) / 2) * min(0.32 / max(k, 1) * 2, 0.06) + rng.normal(0, 0.004, k)105        ax.scatter(x, i + yj, c=cmap(col), s=6, edgecolors="none", alpha=0.85)106        ax.text(ax.get_xlim()[1] if False else 20.5, i, f"{imp[c]*100:.2f}".replace(".", ","), va="center", ha="left", fontsize=7.2, color=palette.ink2)107    ax.axvline(0, color=palette.ink2, lw=0.8)108    ax.set_yticks(range(len(cols))); ax.set_yticklabels([c37.ETIQUETTES[c] for c in cols[::-1]])109    ax.set_xlabel("Valeur SHAP φ (points de pourcentage du logarithme du prix)"); ax.set_xlim(-22, 26)110    ax.text(20.5, len(cols) - 0.35, "moy. |φ|", fontsize=7.2, color=palette.ink2, ha="left", va="bottom")111    sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(0, 1)); sm.set_array([])112    cb = fig.colorbar(sm, ax=ax, pad=0.01, fraction=0.03, ticks=[0, 1]); cb.ax.set_yticklabels(["faible", "élevée"]); cb.set_label("Valeur de la variable", fontsize=7.5)113    cb.outline.set_visible(False)114    savefig(fig, "ch37/ch37-fig-shap-beeswarm.pdf")115116def _centrer(c): return c - np.mean(c)117118def fig_pdp(res):119    pdp = res["pdp"]; v = res["ventes"]120    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 2.9))121    series = [("vraie", "loi de génération", palette.ink, ":"), ("Hédonique (référence)", "hédonique", palette.cat[0], "-"), ("Forêt aléatoire", "forêt aléatoire", palette.cat[1], "-"), ("Gradient boosting", "boosting", palette.cat[2], "-")]122    for ax, (cle, var, xlab) in zip(axes, [("superficie", "superficie_habitable_m2", "Superficie habitable (m²)"), ("age", "age", "Âge (années)"), ("age_scenario", "age", "Âge (années), scénario non linéaire")]):123        d = pdp[cle]; g = d["grille"]124        for k, lab, col, ls in series:125            if k in d: ax.plot(g, _centrer(d[k]) * 100, color=col, ls=ls, lw=1.7 if k == "vraie" else 1.4, label=lab)126        q = np.percentile(v[var], np.arange(5, 100, 5)); ax.plot(q, np.full_like(q, ax.get_ylim()[0] if False else -38), "|", color=palette.ink2, ms=5, alpha=0.8)127        ax.set_xlabel(xlab); ax.set_ylim(-42, 42)128    axes[0].set_ylabel("Dépendance partielle, écart au niveau moyen (%)"); axes[0].legend(loc="upper left", fontsize=7)129    axes[2].text(0.97, 0.04, "prime +12 % (qualité ≥ 4, âge ≥ 40)\n« mur » −10 % (25 à 45 ans)", transform=axes[2].transAxes, ha="right", va="bottom", fontsize=6.8, color=palette.ink2)130    savefig(fig, "ch37/ch37-fig-pdp.pdf")131132def fig_waterfall(res, top: int = 8):133    sh = res["shap"]; phi = sh["sujet"].copy(); base = sh["base"]134    ordre = phi.abs().sort_values(ascending=False); princip = list(ordre.index[:top]); autres = float(phi.drop(princip).sum())135    etapes = [(c37.ETIQUETTES[c], float(phi[c])) for c in princip] + [("Autres variables (somme)", autres)]136    fig, ax = plt.subplots(figsize=(LARGEUR, 3.3))137    cur = base; ys = []138    for i, (lab, ph) in enumerate(etapes):139        a, b = math.exp(cur), math.exp(cur + ph); col = palette.div["pos"] if ph >= 0 else palette.div["neg"]140        ax.barh(i, b - a, left=a, color=col, height=0.62, edgecolor="none")141        ax.text(max(a, b) + 4000, i, f"{ph*100:+.1f} %".replace(".", ","), va="center", ha="left", fontsize=7.5, color=palette.ink)142        if i + 1 < len(etapes): ax.plot([b, b], [i + 0.31, i + 0.69], color=palette.ink2, lw=0.6)143        cur += ph; ys.append(lab)144    fin = math.exp(cur)145    ax.axvline(math.exp(base), color=palette.ink2, lw=0.9, ls=":"); ax.axvline(fin, color=palette.ink, lw=0.9, ls="--")146    ax.text(math.exp(base), -0.9, f"base : {math.exp(base):,.0f} $".replace(",", " "), ha="center", va="top", fontsize=7.5, color=palette.ink2)147    ax.text(fin, len(etapes) - 0.2, f"prédiction : {fin:,.0f} $".replace(",", " "), ha="center", va="bottom", fontsize=7.5, color=palette.ink)148    ax.set_yticks(range(len(etapes))); ax.set_yticklabels(ys); ax.invert_yaxis()149    ax.set_xlabel("Prix prédit (k$) — chaque contribution est un facteur multiplicatif exp(φ)")150    ax.xaxis.set_major_formatter(FuncFormatter(_k)); ax.set_xlim(430_000, 700_000); ax.set_ylim(len(etapes) - 0.3, -1.5)151    savefig(fig, "ch37/ch37-fig-waterfall.pdf")152153def main():154    setup(); res = c37.main(verbose=False)155    fig_arbre(res); fig_apprentissage(res); fig_erreurs(res); fig_predit_observe(res); fig_beeswarm(res); fig_pdp(res); fig_waterfall(res)156    print("figures ch37 écrites dans", ROOT / "assets/figures/ch37")157    return res158159if __name__ == "__main__":160    main()161