spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 19 — modèles spatiaux et temporels (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).3Toutes les cartes sont tracées sur les coordonnées synthétiques (km est/nord autour du centroïde des ventes), sans fond de carte.4Sorties (assets/figures/ch19/) :5  ch19-fig-carte-residus.pdf      résidus MCO sans indicatrices de secteur vs modèle de référence6  ch19-fig-correlogramme.pdf      I de Moran par bande de distance7  ch19-fig-gwr-coefficients.pdf   coefficients locaux de la GWR (localisation, superficie, âge)8  ch19-fig-gwr-bande.pdf          score de validation croisée selon la bande passante9  ch19-fig-variogramme.pdf        variogramme empirique et modèle exponentiel10  ch19-fig-krigeage.pdf           surface krigée des résidus11  ch19-fig-sensibilite-W.pdf      I, ρ et λ selon la matrice de poids12Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch19_figures.py"""13import sys, pathlib14ROOT = 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 np17import matplotlib.pyplot as plt18from matplotlib.ticker import FuncFormatter19from traite.style import setup, palette, savefig, LARGEUR, DEMI20from traite import hedonique as h21import ch19_spatial as c192223SECTEURS = ["Hull", "Gatineau", "Aylmer", "Buckingham", "Masson-Angers"]2425def _axes_km(ax, xy):26    ax.set_aspect("equal"); ax.set_xlim(xy[:, 0].min() - 1.5, xy[:, 0].max() + 1.5); ax.set_ylim(xy[:, 1].min() - 1.5, xy[:, 1].max() + 1.5)27    ax.set_xlabel("Est (km, coordonnées synthétiques)"); ax.grid(True, lw=0.4)2829def _etiquettes_secteurs(ax, v, xy, dy=0.0):30    for s in SECTEURS:31        m = (v["secteur"] == s).values32        cx, cy = np.median(xy[m, 0]), np.median(xy[m, 1])33        ax.text(cx, cy + dy, s, fontsize=7, color=palette.ink2, ha="center", va="center",34                bbox=dict(boxstyle="round,pad=0.15", fc=palette.surface, ec="none", alpha=0.75))3536def fig_carte_residus(res):37    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]38    series = [("Sans indicatrices de secteur (I = %.2f)" % res["moran"].loc[2, "I_knn"], b["mco_r"]["residus"]),39              ("Modèle de référence (I = %.2f)" % res["moran"].loc[3, "I_knn"], b["modele"]["residus"])]40    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.6), sharey=True)41    vmax = 0.25; cmap = palette.cmap_div()42    for ax, (titre, e) in zip(axes, series):43        sc = ax.scatter(xy[:, 0], xy[:, 1], c=np.clip(e, -vmax, vmax), cmap=cmap, vmin=-vmax, vmax=vmax, s=9, lw=0)44        ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=110, color=palette.ink, zorder=5)45        _axes_km(ax, xy); ax.set_title(titre, fontsize=8.5); _etiquettes_secteurs(ax, v, xy, dy=4.5)46    axes[0].set_ylabel("Nord (km)")47    cb = fig.colorbar(sc, ax=axes, orientation="horizontal", fraction=0.05, pad=0.16, aspect=40)48    cb.set_label("Résidu du modèle semi-logarithmique (ln P observé − ln P ajusté)", fontsize=8)49    cb.ax.tick_params(labelsize=7.5)50    savefig(fig, "ch19/ch19-fig-carte-residus.pdf")5152def fig_correlogramme(res):53    cg = res["correlogramme"]54    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))55    for i, (cle, nom) in enumerate([("sans_secteur", "Résidus sans indicatrices de secteur"), ("reference", "Résidus du modèle de référence")]):56        pts = cg[cle]; x = np.array([(p["d_min"] + p["d_max"]) / 2 for p in pts]); I = np.array([p["I"] for p in pts])57        sd = np.sqrt(np.array([p["Var_I"] for p in pts]))58        ax.errorbar(x + (i - 0.5) * 0.12, I, yerr=1.96 * sd, fmt="o-", color=palette.cat[i], ms=4.5, lw=1.4, capsize=2, label=nom)59    ax.axhline(cg["reference"][0]["E_I"], color=palette.ink2, lw=0.9, ls="--")60    ax.text(14.6, cg["reference"][0]["E_I"] + 0.012, "E[I] = −1/(n−1)", fontsize=7, color=palette.ink2, ha="right")61    ax.set_xlabel("Distance entre paires de ventes (km, centre de la bande)"); ax.set_ylabel("I de Moran de la bande")62    ax.legend(loc="upper right", fontsize=7.5)63    savefig(fig, "ch19/ch19-fig-correlogramme.pdf")6465def fig_gwr_coefficients(res):66    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]; g = res["gwr"]["gwr"]; cols = list(b["Xr"].columns)67    j = cols.index("dist_centre"); dist = b["Xr"]["dist_centre"].values68    sect = v["secteur"].values69    loc = g["beta"][:, 0] + g["beta"][:, j] * dist - np.median(g["beta"][sect == "Hull", 0])70    panneaux = [("Effet de localisation local\n(écart à Hull, en log)", loc, palette.cmap_div(), -0.32, 0.32),71                ("Coefficient local de ln(superficie)\n(élasticité)", g["beta"][:, cols.index("ln_S")], palette.cmap_seq(), None, None),72                ("Coefficient local de l'âge\n(× 100, % par année)", 100 * g["beta"][:, cols.index("age")], palette.cmap_seq(), None, None)]73    fig, axes = plt.subplots(1, 3, figsize=(LARGEUR, 2.75))74    for ax, (titre, z, cmap, vmin, vmax) in zip(axes, panneaux):75        sc = ax.scatter(xy[:, 0], xy[:, 1], c=z, cmap=cmap, vmin=vmin, vmax=vmax, s=9, lw=0)76        ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=70, color=palette.ink, zorder=5)77        _axes_km(ax, xy); ax.set_title(titre, fontsize=7.5); ax.tick_params(labelsize=6.5); ax.set_xlabel("Est (km)", fontsize=7)78        cb = fig.colorbar(sc, ax=ax, orientation="horizontal", fraction=0.06, pad=0.22, aspect=28); cb.ax.tick_params(labelsize=6.5)79    axes[0].set_ylabel("Nord (km)", fontsize=7)80    savefig(fig, "ch19/ch19-fig-gwr-coefficients.pdf")8182def fig_gwr_bande(res):83    sel = res["gwr"]["selection"]84    fig, ax = plt.subplots(figsize=(DEMI, 2.5))85    ax.plot(sel["grille"], sel["cv_grille"], marker="o", ms=3.5, color=palette.cat[0])86    ax.axvline(sel["bande"], color=palette.ink2, ls="--", lw=0.9)87    ax.text(sel["bande"] * 1.1, sel["cv"] * 1.6, "h* = %.2f km" % sel["bande"], fontsize=7, color=palette.ink2, va="bottom")88    ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("Bande passante h (km, échelle log)"); ax.set_ylabel("Score de validation croisée\n(somme des carrés des erreurs, échelle log)")89    ax.set_xticks([0.3, 0.5, 1, 2, 5, 10]); ax.get_xaxis().set_major_formatter(FuncFormatter(lambda x, p: ("%g" % x)))90    savefig(fig, "ch19/ch19-fig-gwr-bande.pdf")9192def fig_variogramme(res):93    vk = res["variogramme"]; vg, vr, p = vk["variogramme"], vk["variogramme_reference"], vk["parametres"]94    fig, ax = plt.subplots(figsize=(LARGEUR, 3.1))95    taille = 12 + 60 * vg["n"] / vg["n"].max()96    ax.scatter(vg["h"], vg["gamma"], s=taille, color=palette.cat[0], alpha=0.85, lw=0, label="Résidus sans indicatrices de secteur (taille selon le nombre de paires)", zorder=3)97    hh = np.linspace(0.05, 16, 200)98    ax.plot(hh, h.variogramme_exponentiel(hh, p["pepite"], p["palier_partiel"], p["portee"]), color=palette.cat[0], lw=1.6,99            label="Modèle exponentiel ajusté (pépite %.4f, palier %.4f, portée pratique %.0f km)" % (p["pepite"], p["palier"], p["portee_pratique"]))100    ax.scatter(vr["h"], vr["gamma"], s=14, marker="s", color=palette.cat[1], lw=0, label="Résidus du modèle de référence", zorder=3)101    ax.axhline(vk["sigma2_reference"], color=palette.cat[1], lw=0.9, ls=":")102    ax.axhline(vk["variance_residus"], color=palette.cat[0], lw=0.9, ls=":")103    ax.text(15.8, vk["variance_residus"] + 0.0004, "variance des résidus sans secteur", fontsize=7, color=palette.ink2, ha="right")104    ax.text(15.8, vk["sigma2_reference"] - 0.0012, "σ² du modèle de référence", fontsize=7, color=palette.ink2, ha="right")105    ax.set_xlabel("Distance h entre paires de ventes (km)"); ax.set_ylabel("Semi-variance γ(h) (log²)")106    ax.set_ylim(0, 0.02); ax.legend(loc="upper left", fontsize=7)107    savefig(fig, "ch19/ch19-fig-variogramme.pdf")108109def fig_krigeage(res):110    b = res["base"]; v, xy, cs = b["ventes"], b["xy"], b["cs"]; gr = res["variogramme"]["grille"]111    GX, GY = np.meshgrid(gr["x"], gr["y"]); pts = np.column_stack([GX.ravel(), GY.ravel()])112    from scipy.spatial import cKDTree113    dmin = cKDTree(xy).query(pts)[0].reshape(GX.shape)114    Z = np.ma.masked_where(dmin > 2.5, gr["z"])115    fig, ax = plt.subplots(figsize=(LARGEUR, 3.6))116    vmax = 0.25117    cf = ax.contourf(GX, GY, np.clip(Z, -vmax, vmax), levels=np.linspace(-vmax, vmax, 21), cmap=palette.cmap_div(), vmin=-vmax, vmax=vmax)118    ax.scatter(xy[:, 0], xy[:, 1], s=3, color=palette.ink, alpha=0.45, lw=0)119    ax.scatter(cs[:, 0], cs[:, 1], marker="*", s=120, color=palette.ink, zorder=5)120    ax.annotate("sujet : résidu krigé %.3f" % res["variogramme"]["sujet"]["prediction"][0], (cs[0, 0], cs[0, 1]), xytext=(14, -34), textcoords="offset points",121                fontsize=7, color=palette.ink, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.7))122    _axes_km(ax, xy); ax.set_ylabel("Nord (km)"); _etiquettes_secteurs(ax, v, xy, dy=5.0)123    cb = fig.colorbar(cf, ax=ax, fraction=0.035, pad=0.02); cb.set_label("Résidu krigé (log)", fontsize=8); cb.ax.tick_params(labelsize=7.5)124    savefig(fig, "ch19/ch19-fig-krigeage.pdf")125126def fig_sensibilite(res):127    s = res["sensibilite"]; k = s[s["famille"].str.startswith("k")]128    fig, ax = plt.subplots(figsize=(LARGEUR, 2.9))129    ax.errorbar(k["parametre"], k["rho"], yerr=1.96 * k["se_rho"], fmt="o-", color=palette.cat[0], ms=4.5, capsize=2, label="ρ du SAR (± 1,96 é.-t.)")130    ax.errorbar(k["parametre"], k["lambda"], yerr=1.96 * k["se_lambda"], fmt="s-", color=palette.cat[1], ms=4.5, capsize=2, label="λ du SEM (± 1,96 é.-t.)")131    ax.plot(k["parametre"], k["I"], "^-", color=palette.cat[2], ms=4.5, label="I de Moran des résidus MCO")132    ax.set_xlabel("Nombre k de plus proches voisins dans W"); ax.set_ylabel("Coefficient"); ax.set_xticks(k["parametre"].astype(int))133    ax.set_ylim(0, 0.95); ax.legend(loc="upper left", fontsize=7.5)134    savefig(fig, "ch19/ch19-fig-sensibilite-W.pdf")135136def main():137    setup()138    res = c19.main(verbose=False)139    for f in [fig_carte_residus, fig_correlogramme, fig_gwr_coefficients, fig_gwr_bande, fig_variogramme, fig_krigeage, fig_sensibilite]:140        f(res); print("OK", f.__name__)141142if __name__ == "__main__":143    main()144