spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 19 — modèles spatiaux et temporels sur le cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).34Reproduit tous les chiffres cités dans manuscript/chapters/ch19.tex :5 1. matrices de poids (k plus proches voisins, distance inverse, spatio-temporelle de Dubé–Legros) ;6 2. I de Moran des résidus du modèle hédonique de référence (traite.cas1.hedonique) avec et sans indicatrices de secteur ;7 3. corrélogramme (I de Moran par bande de distance) ;8 4. tests du multiplicateur de Lagrange (LM-erreur, LM-décalage et versions robustes) ;9 5. MCO / SAR / SEM / SDM par maximum de vraisemblance concentrée, effets directs et indirects ;10 6. régression géographiquement pondérée (noyau gaussien fixe, bande passante par validation croisée) ;11 7. effets fixes de quartier (secteurs, grilles) et de temps (tendance, trimestres) ;12 8. variogramme empirique des résidus, modèle exponentiel et krigeage ordinaire ;13 9. décalage spatio-temporel (ventes antérieures seulement) ;14 10. validation croisée et indications pour le sujet ; sensibilité au choix de la matrice de poids.1516Contrôle croisé (non requis pour l'exécution) : spreg (ML_Lag, ML_Error, LMtests), esda (Moran) et mgwr (GWR) donnent les17mêmes estimations à 1e-6 près ; voir --controle.18Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch19_spatial.py [--controle]"""19from __future__ import annotations20import sys, pathlib, warnings21import numpy as np, pandas as pd22sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))23from traite import cas124from traite import hedonique as h2526ROOT = pathlib.Path(__file__).resolve().parents[3]27K_VOISINS = 8 # matrice de référence : 8 plus proches voisins, standardisée en ligne28PORTEE_KM = 2.0 # matrice de distance inverse : portée 2 km29FENETRE_MOIS = 12 # matrice spatio-temporelle : ventes antérieures de moins de 12 mois30BANDES_CORRELOGRAMME = [0, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 15]31BORNES_VARIOGRAMME = np.arange(0, 16.1, 1.0)32VOISINS_KRIGEAGE = 4033# primes de secteur du générateur (data/synthetic/README.md : log du rapport des médianes APCIQ T3 2026 à Hull ; Masson-Angers − 0,02)34PRIMES_DGP = {"Hull": 0.0, "Gatineau": float(np.log(464900 / 475000)), "Aylmer": float(np.log(530000 / 475000)),35 "Buckingham": float(np.log(395000 / 475000)), "Masson-Angers": float(np.log(395000 / 475000)) - 0.02}3637_CACHE: dict = {}3839# ----------------------------------------------------------------------------- base4041def charger_base() -> dict:42 """Modèle de référence (traite.cas1.hedonique), spécification « sans secteur », coordonnées en km et sujet."""43 d = cas1.charger(); v = d["ventes"]; s = d["sujet"]44 m = cas1.hedonique(v)45 Xc = m["X"]; y = m["y"]46 Xr = Xc.drop(columns=[c for c in Xc.columns if c.startswith("sect_")]) # sans indicatrices de secteur47 Xn = Xr.drop(columns=["dist_centre"]) # ni secteur ni gradient de distance48 lat0, lon0 = float(v["lat"].mean()), float(v["lon"].mean())49 xy = h.coordonnees_km(v["lat"].values, v["lon"].values, lat0, lon0)50 cs = h.coordonnees_km([s["lat"]], [s["lon"]], lat0, lon0)51 xs_r = cas1.matrice(cas1.sujet_df(s))[Xr.columns].values[0]52 mco_r = h.mco(Xr.values, y); mco_n = h.mco(Xn.values, y)53 return {"ventes": v, "sujet": s, "modele": m, "Xc": Xc, "Xr": Xr, "Xn": Xn, "y": y, "xy": xy, "cs": cs, "xs_r": xs_r,54 "mois": Xc["mois"].values, "mco_r": mco_r, "mco_n": mco_n, "lat0": lat0, "lon0": lon0}5556def matrices(b: dict) -> dict:57 xy, mois = b["xy"], b["mois"]58 return {"knn": h.matrice_poids_knn(xy, K_VOISINS), "dist": h.matrice_poids_distance_inverse(xy, PORTEE_KM),59 "st": h.matrice_poids_spatio_temporelle(xy, mois, PORTEE_KM, FENETRE_MOIS)}6061def _loglik_mco(res: np.ndarray) -> float:62 n = len(res); return float(-n / 2 * (np.log(2 * np.pi * res @ res / n) + 1))6364def _aic_mco(res: np.ndarray, k: int) -> float:65 return 2 * (k + 1) - 2 * _loglik_mco(res)6667def _aicc(res: np.ndarray, k_eff: float) -> float:68 n = len(res); rss = float(res @ res)69 return n * np.log(rss / n) + n * np.log(2 * np.pi) + n * (n + k_eff) / (n - 2 - k_eff)7071# ----------------------------------------------------------------------------- 2-3. Moran et corrélogramme7273def tableau_moran(b: dict, Ws: dict) -> pd.DataFrame:74 series = {"ln P (sans modèle)": b["y"], "Résidus, sans secteur ni distance": b["mco_n"]["residus"],75 "Résidus, sans indicatrices de secteur": b["mco_r"]["residus"], "Résidus, modèle de référence": b["modele"]["residus"]}76 rows = []77 for nom, x in series.items():78 row = {"serie": nom}79 for wn, W in Ws.items():80 r = h.moran_i(x, W); row[f"I_{wn}"] = r["I"]; row[f"z_{wn}"] = r["z"]; row[f"p_{wn}"] = r["p"]81 rows.append(row)82 return pd.DataFrame(rows)8384def correlogramme(b: dict) -> dict:85 return {"sans_secteur": h.moran_par_distance(b["mco_r"]["residus"], b["xy"], BANDES_CORRELOGRAMME),86 "reference": h.moran_par_distance(b["modele"]["residus"], b["xy"], BANDES_CORRELOGRAMME)}8788# ----------------------------------------------------------------------------- 4-5. LM, SAR/SEM/SDM, effets8990def tests_lm(b: dict, W: np.ndarray) -> dict:91 return {"reference": h.tests_lm_spatial(b["y"], b["Xc"].values, W), "sans_secteur": h.tests_lm_spatial(b["y"], b["Xr"].values, W)}9293def modeles_spatiaux(b: dict, W: np.ndarray) -> dict:94 y, Xr, Xc = b["y"], b["Xr"].values, b["Xc"].values95 out = {"mco_r": b["mco_r"], "sar_r": h.sar_mv(y, Xr, W), "sem_r": h.sem_mv(y, Xr, W), "sdm_r": h.sar_mv(y, Xr, W, durbin=True),96 "sar_c": h.sar_mv(y, Xc, W), "sem_c": h.sem_mv(y, Xc, W)}97 for cle in ["sar_r", "sem_r", "sdm_r", "sar_c", "sem_c"]:98 res = out[cle]["innovations"] if "innovations" in out[cle] else out[cle]["residus"]99 out[cle]["moran_residus"] = h.moran_i(res, W)["I"]100 out["mco_r"]["loglik"] = _loglik_mco(b["mco_r"]["residus"]); out["mco_r"]["aic"] = _aic_mco(b["mco_r"]["residus"], Xr.shape[1])101 out["mco_r"]["moran_residus"] = h.moran_i(b["mco_r"]["residus"], W)["I"]102 out["mco_c"] = {"beta": b["modele"]["beta"], "se": b["modele"]["coef"]["se_hc1"].values, "residus": b["modele"]["residus"],103 "loglik": _loglik_mco(b["modele"]["residus"]), "aic": _aic_mco(b["modele"]["residus"], Xc.shape[1]),104 "sigma2": float(b["modele"]["sigma"] ** 2), "moran_residus": h.moran_i(b["modele"]["residus"], W)["I"]}105 return out106107def effets(mod: dict, W: np.ndarray, colonnes) -> pd.DataFrame:108 k = len(colonnes)109 e_sar = h.effets_sar(mod["sar_r"]["rho"], mod["sar_r"]["beta"], W)110 theta = np.concatenate([[0.0], mod["sdm_r"]["beta"][k:]])111 e_sdm = h.effets_sar(mod["sdm_r"]["rho"], mod["sdm_r"]["beta"][:k], W, theta)112 df = pd.DataFrame({"variable": colonnes, "beta_mco": mod["mco_r"]["beta"], "beta_sar": mod["sar_r"]["beta"],113 "sar_direct": e_sar["direct"], "sar_indirect": e_sar["indirect"], "sar_total": e_sar["total"],114 "beta_sdm": mod["sdm_r"]["beta"][:k], "theta_sdm": theta,115 "sdm_direct": e_sdm["direct"], "sdm_indirect": e_sdm["indirect"], "sdm_total": e_sdm["total"]})116 return df.set_index("variable")117118# ----------------------------------------------------------------------------- 6. GWR119120def gwr_resultats(b: dict) -> dict:121 y, X, xy = b["y"], b["Xr"].values, b["xy"]122 sel = h.gwr_bande_cv(y, X, xy)123 g = h.gwr(y, X, xy, sel["bande"])124 gs = h.gwr(y, X, xy, sel["bande"], cibles=b["cs"])125 cols = list(b["Xr"].columns)126 beta = pd.DataFrame(g["beta"], columns=cols)127 quart = beta.describe(percentiles=[0.25, 0.5, 0.75]).T[["min", "25%", "50%", "75%", "max"]]128 inter = pd.Series(g["beta"][:, 0]).groupby(b["ventes"]["secteur"].values).median()129 inter_rel = inter - inter["Hull"]130 dummies = {"Hull": 0.0, **{s: float(b["modele"]["coef"].loc[f"sect_{s}", "coef"]) for s in cas1.SECTEURS[1:]}}131 comp = pd.DataFrame({"gwr_intercept_rel_hull": inter_rel, "indicatrice_reference": pd.Series(dummies), "prime_generateur": pd.Series(PRIMES_DGP)})132 aicc_mco_r = _aicc(b["mco_r"]["residus"], X.shape[1]); aicc_mco_c = _aicc(b["modele"]["residus"], b["Xc"].shape[1])133 return {"selection": sel, "gwr": g, "beta": beta, "quartiles": quart, "comparaison_secteurs": comp, "sujet_beta": gs["beta"][0],134 "aicc_mco_sans_secteur": aicc_mco_r, "aicc_mco_reference": aicc_mco_c}135136# ----------------------------------------------------------------------------- 7. effets fixes137138def effets_fixes(b: dict, W: np.ndarray) -> pd.DataFrame:139 y, xy, v = b["y"], b["xy"], b["ventes"]; n = len(y)140 def ligne(nom, Xm, k_spatial, k_temps):141 r = h.mco(Xm, y); k = Xm.shape[1]142 return {"specification": nom, "parametres": k, "effets_spatiaux": k_spatial, "effets_temps": k_temps, "r2": r["r2"],143 "r2_ajuste": r["r2_ajuste"], "sigma": float(np.sqrt(r["sigma2"])), "aic": _aic_mco(r["residus"], k), "moran": h.moran_i(r["residus"], W)["I"]}144 rows = [ligne("Aucun effet spatial ; tendance linéaire", b["Xr"].values, 0, 1),145 ligne("Gradient de distance seulement ; tendance linéaire", b["Xr"].values, 1, 1) if False else None,146 ligne("Secteurs (référence) ; tendance linéaire", b["Xc"].values, 4, 1)]147 rows = [r for r in rows if r]148 for cell in [5.0, 3.0, 2.0]:149 gx = np.floor(xy[:, 0] / cell); gy = np.floor(xy[:, 1] / cell); cid = pd.factorize(gx * 1000 + gy)[0]150 D = pd.get_dummies(cid, drop_first=True).astype(float).values151 rows.append(ligne(f"Grille de {cell:g} km ({cid.max() + 1} cellules) ; tendance linéaire", np.column_stack([b["Xr"].values, D]), D.shape[1], 1))152 q = pd.PeriodIndex(v["date_vente"], freq="Q").astype(str)153 Dq = pd.get_dummies(q, drop_first=True).astype(float).values154 rows.append(ligne(f"Secteurs ; indicatrices trimestrielles ({Dq.shape[1]})", np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), 4, Dq.shape[1]))155 rt = h.mco(np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), y)156 df = pd.DataFrame(rows)157 df.attrs["indice_trimestriel"] = pd.Series(np.concatenate([[0.0], rt["beta"][-Dq.shape[1]:]]), index=sorted(q.unique()))158 df.attrs["tendance_mensuelle"] = float(b["modele"]["coef"].loc["mois", "coef"])159 return df160161# ----------------------------------------------------------------------------- 8. variogramme et krigeage162163def variogramme_krigeage(b: dict) -> dict:164 e = b["mco_r"]["residus"]; xy = b["xy"]165 vg = h.variogramme_empirique(xy, e, BORNES_VARIOGRAMME)166 vg_ref = h.variogramme_empirique(xy, b["modele"]["residus"], BORNES_VARIOGRAMME)167 p = h.ajuster_variogramme(vg["h"], vg["gamma"], vg["n"], portee_max=5.0)168 loo = h.krigeage_validation_croisee(xy, e, p, VOISINS_KRIGEAGE)169 sujet = h.krigeage_ordinaire(xy, e, b["cs"], p, VOISINS_KRIGEAGE)170 # surface pour la figure : grille régulière (km)171 gx = np.linspace(xy[:, 0].min() - 1, xy[:, 0].max() + 1, 90); gy = np.linspace(xy[:, 1].min() - 1, xy[:, 1].max() + 1, 70)172 GX, GY = np.meshgrid(gx, gy); grille = np.column_stack([GX.ravel(), GY.ravel()])173 kg = h.krigeage_ordinaire(xy, e, grille, p, VOISINS_KRIGEAGE)174 return {"variogramme": vg, "variogramme_reference": vg_ref, "parametres": p, "loo": loo, "sujet": sujet,175 "grille": {"x": gx, "y": gy, "z": kg["prediction"].reshape(GX.shape), "var": kg["variance"].reshape(GX.shape)},176 "variance_residus": float(e.var()), "sigma2_reference": float(b["modele"]["sigma"] ** 2)}177178# ----------------------------------------------------------------------------- 9. spatio-temporel179180def spatio_temporel(b: dict, Ws: dict) -> dict:181 y, Xr, Wst = b["y"], b["Xr"].values, Ws["st"]182 e = b["mco_r"]["residus"]183 ok = Wst.sum(1) > 0 # ventes ayant au moins une vente antérieure voisine184 Wy = Wst @ y185 # décalage spatio-temporel : prédéterminé (matrice triangulaire après tri par date) → MV = MCO sur [W_st y, X]186 Z = np.column_stack([Wy[ok], Xr[ok]]); r = h.mco(Z, y[ok])187 sem_st = h.sem_mv(y, Xr, Wst)188 return {"n_retenues": int(ok.sum()), "n_sans_voisin": int((~ok).sum()), "voisins_moyen": float((Wst[ok] > 0).sum(1).mean()),189 "rho_st": float(r["beta"][0]), "se_rho_st": float(r["se_hc1"][0]), "beta_st": r["beta"][1:], "r2_st": r["r2"],190 "moran_e_st": h.moran_i(e, Wst), "moran_e_knn": h.moran_i(e, Ws["knn"]), "moran_e_dist": h.moran_i(e, Ws["dist"]),191 "sem_st": sem_st, "max_abs_valeur_propre": float(np.abs(np.linalg.eigvals(Wst)).max()),192 "moran_e_apres_lag": h.moran_i(r["residus"], Wst[np.ix_(ok, ok)])["I"]}193194# ----------------------------------------------------------------------------- 10. validation croisée et sujet195196def validation_croisee(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:197 y = b["y"]; n = len(y)198 def loo_mco(Xm):199 r = h.mco(Xm, y); H = Xm @ np.linalg.pinv(Xm.T @ Xm) @ Xm.T; return r["residus"] / (1 - np.diag(H))200 sem = mod["sem_r"]; u = sem["residus"]; B = np.eye(n) - sem["lambda"] * W; Q = B.T @ B201 Eu = np.array([-(Q[i] @ u - Q[i, i] * u[i]) / Q[i, i] for i in range(n)])202 series = {"MCO, modèle de référence (secteurs)": loo_mco(b["Xc"].values),203 "MCO, sans indicatrices de secteur": loo_mco(b["Xr"].values),204 "MCO sans secteur + krigeage des résidus": b["mco_r"]["residus"] - vk["loo"],205 "SEM sans secteur (espérance conditionnelle)": u - Eu,206 "SAR sans secteur (prédiction conditionnelle)": y - (mod["sar_r"]["rho"] * (W @ y) + b["Xr"].values @ mod["sar_r"]["beta"]),207 "GWR sans secteur (validation croisée)": None}208 rows = []209 for nom, r in series.items():210 if r is None:211 rmse = float(np.sqrt(g["gwr"]["cv"] / n)); rows.append({"modele": nom, "rmse_log": rmse, "erreur_abs_moy_pct": np.nan, "moran_erreurs": h.moran_i(g["gwr"]["residus"], W)["I"]})212 else:213 rows.append({"modele": nom, "rmse_log": float(np.sqrt(np.mean(r ** 2))), "erreur_abs_moy_pct": float(np.mean(np.abs(np.exp(r) - 1))), "moran_erreurs": h.moran_i(r, W)["I"]})214 return pd.DataFrame(rows)215216def predictions_sujet(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:217 s, xs = b["sujet"], b["xs_r"]; y = b["y"]; n = len(y)218 duan = lambda res: float(np.mean(np.exp(res)))219 ref = float(cas1.predire(b["modele"], cas1.sujet_df(s))[0])220 mco_r = float(np.exp(xs @ b["mco_r"]["beta"]) * duan(b["mco_r"]["residus"]))221 krig = float(np.exp(xs @ b["mco_r"]["beta"] + vk["sujet"]["prediction"][0]) * duan(b["mco_r"]["residus"] - vk["loo"]))222 # SEM : espérance conditionnelle de u au sujet (matrice augmentée : le sujet n'est voisin d'aucune vente)223 xy_aug = np.vstack([b["xy"], b["cs"]]); Wa = np.zeros((n + 1, n + 1)); Wa[:n, :n] = W224 Wa[n, :n] = h.matrice_poids_knn(xy_aug, K_VOISINS)[n, :n]; Wa[n] /= Wa[n].sum()225 sem = mod["sem_r"]; u_s = float(h.sem_esperance_conditionnelle(sem["lambda"], Wa, sem["residus"], [n])[0])226 sem_pred = float(np.exp(xs @ sem["beta"] + u_s) * duan(sem["innovations"]))227 # SAR : prédiction conditionnelle avec le décalage observé (moyenne des ln P des k ventes voisines)228 idx = np.nonzero(Wa[n, :n])[0]; Wy_s = float(Wa[n, :n] @ y)229 sar = mod["sar_r"]; sar_pred = float(np.exp(sar["rho"] * Wy_s + xs @ sar["beta"]) * duan(sar["residus"]))230 gwr_pred = float(np.exp(xs @ g["sujet_beta"]) * duan(g["gwr"]["residus"]))231 df = pd.DataFrame([232 {"modele": "MCO, modèle de référence (indicatrices de secteur)", "ln_correction": 0.0, "indication": ref},233 {"modele": "MCO sans indicatrices de secteur", "ln_correction": 0.0, "indication": mco_r},234 {"modele": "MCO sans secteur + krigeage du résidu au sujet", "ln_correction": float(vk["sujet"]["prediction"][0]), "indication": krig},235 {"modele": "SEM sans secteur (espérance conditionnelle du résidu)", "ln_correction": u_s, "indication": sem_pred},236 {"modele": "SAR sans secteur (décalage des 8 ventes voisines)", "ln_correction": float(sar["rho"] * Wy_s), "indication": sar_pred},237 {"modele": "GWR sans secteur (coefficients locaux au sujet)", "ln_correction": float(g["sujet_beta"][0] - b["mco_r"]["beta"][0]), "indication": gwr_pred}])238 df.attrs["voisins_sujet"] = list(b["ventes"]["id_vente"].values[idx]); df.attrs["secteurs_voisins"] = list(b["ventes"]["secteur"].values[idx])239 df.attrs["variance_krigeage_sujet"] = float(vk["sujet"]["variance"][0])240 return df241242# ----------------------------------------------------------------------------- 11. sensibilité243244def sensibilite(b: dict) -> pd.DataFrame:245 y, Xr, xy, e = b["y"], b["Xr"].values, b["xy"], b["mco_r"]["residus"]246 specs = [("k plus proches voisins", k, h.matrice_poids_knn(xy, k)) for k in [4, 6, 8, 12, 16]]247 specs += [("distance inverse", d, h.matrice_poids_distance_inverse(xy, d)) for d in [1.0, 2.0, 4.0]]248 rows = []249 for fam, par, W in specs:250 mi = h.moran_i(e, W); sar = h.sar_mv(y, Xr, W); sem = h.sem_mv(y, Xr, W)251 rows.append({"famille": fam, "parametre": par, "voisins_moyen": float((W > 0).sum(1).mean()), "I": mi["I"], "z": mi["z"],252 "rho": sar["rho"], "se_rho": sar["se_rho"], "aic_sar": sar["aic"], "lambda": sem["lambda"], "se_lambda": sem["se_lambda"], "aic_sem": sem["aic"]})253 return pd.DataFrame(rows)254255# ----------------------------------------------------------------------------- contrôle croisé facultatif256257def controle_croise(b: dict, W: np.ndarray, g: dict) -> dict:258 """Compare les estimations maison à spreg/esda/mgwr (installés sur la machine de rédaction)."""259 import spreg260 from libpysal.weights import W as Wcls261 from esda.moran import Moran262 from mgwr.sel_bw import Sel_BW263 from mgwr.gwr import GWR264 y, X = b["y"], b["Xr"]265 neigh = {i: list(np.nonzero(W[i])[0]) for i in range(len(W))}; wts = {i: [W[i, j] for j in neigh[i]] for i in neigh}266 w = Wcls(neigh, wts); w.transform = "r"; Xs = X.drop(columns=["const"]).values; yy = y.reshape(-1, 1)267 with warnings.catch_warnings():268 warnings.simplefilter("ignore")269 ml = spreg.ML_Lag(yy, Xs, w, method="full"); me = spreg.ML_Error(yy, Xs, w, method="full")270 ols = spreg.OLS(yy, Xs, w=w, spat_diag=True, moran=True)271 mo = Moran(b["mco_r"]["residus"], w, permutations=0)272 bw = Sel_BW(b["xy"], yy, Xs, fixed=True, kernel="gaussian").search(criterion="CV")273 gm = GWR(b["xy"], yy, Xs, bw, fixed=True, kernel="gaussian").fit()274 mine_sar = h.sar_mv(y, X.values, W); mine_sem = h.sem_mv(y, X.values, W); lm = h.tests_lm_spatial(y, X.values, W)275 return {"rho_spreg": float(ml.rho), "rho_maison": mine_sar["rho"], "lambda_spreg": float(me.lam), "lambda_maison": mine_sem["lambda"],276 "lm_err_spreg": float(ols.lm_error[0]), "lm_err_maison": lm["LM_erreur"][0], "rlm_lag_spreg": float(ols.rlm_lag[0]), "rlm_lag_maison": lm["RLM_decalage"][0],277 "moran_esda_z": float(mo.z_norm), "moran_maison_z": h.moran_i(b["mco_r"]["residus"], W)["z"],278 "bande_mgwr": float(bw), "bande_maison": g["selection"]["bande"], "aicc_mgwr": float(gm.aicc), "aicc_maison_meme_bande": h.gwr(y, X.values, b["xy"], float(bw))["aicc"]}279280# ----------------------------------------------------------------------------- main281282def main(verbose: bool = True, controle: bool = False) -> dict:283 if "res" in _CACHE and not controle:284 return _CACHE["res"]285 b = charger_base(); Ws = matrices(b); W = Ws["knn"]286 res = {"base": b, "W": Ws, "moran": tableau_moran(b, Ws), "correlogramme": correlogramme(b), "lm": tests_lm(b, W)}287 res["modeles"] = modeles_spatiaux(b, W)288 res["effets"] = effets(res["modeles"], W, list(b["Xr"].columns))289 res["gwr"] = gwr_resultats(b)290 res["effets_fixes"] = effets_fixes(b, W)291 res["variogramme"] = variogramme_krigeage(b)292 res["spatio_temporel"] = spatio_temporel(b, Ws)293 res["validation"] = validation_croisee(b, res["modeles"], res["gwr"], res["variogramme"], W)294 res["sujet"] = predictions_sujet(b, res["modeles"], res["gwr"], res["variogramme"], W)295 res["sensibilite"] = sensibilite(b)296 if controle:297 res["controle"] = controle_croise(b, W, res["gwr"])298 _CACHE["res"] = res299 if verbose:300 imprimer(res)301 return res302303def imprimer(res: dict):304 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)305 b = res["base"]; m = b["modele"]306 print("== Modèle de référence : R² %.3f, σ %.4f, n %d, k %d" % (m["r2"], m["sigma"], m["n"], m["k"]))307 print("== Sans secteur : R² %.3f, σ %.4f ; sans secteur ni distance : R² %.3f" % (b["mco_r"]["r2"], np.sqrt(b["mco_r"]["sigma2"]), b["mco_n"]["r2"]))308 print("== Étendue des coordonnées (km) :", b["xy"].min(0).round(1), b["xy"].max(0).round(1), "; sujet :", b["cs"].round(2))309 for wn, W in res["W"].items():310 print(f" matrice {wn} : voisins moyens {(W > 0).sum(1).mean():.1f}, lignes vides {(W.sum(1) == 0).sum()}")311 print("== I de Moran :\n", res["moran"].round(4).to_string(index=False))312 print("== Corrélogramme (sans secteur) :", [(c["d_min"], c["d_max"], round(c["I"], 3), round(c["z"], 1), c["paires"]) for c in res["correlogramme"]["sans_secteur"]])313 print("== Corrélogramme (référence) :", [(c["d_min"], c["d_max"], round(c["I"], 3), round(c["z"], 1)) for c in res["correlogramme"]["reference"]])314 for k, v in res["lm"].items():315 print(f"== Tests LM ({k}) :", {a: (round(s, 2), round(p, 4)) for a, (s, p) in v.items()})316 mod = res["modeles"]317 print("== SAR sans secteur : ρ %.4f (%.4f), logL %.2f, AIC %.1f, I résidus %.3f" % (mod["sar_r"]["rho"], mod["sar_r"]["se_rho"], mod["sar_r"]["loglik"], mod["sar_r"]["aic"], mod["sar_r"]["moran_residus"]))318 print("== SEM sans secteur : λ %.4f (%.4f), logL %.2f, AIC %.1f, I innovations %.3f" % (mod["sem_r"]["lambda"], mod["sem_r"]["se_lambda"], mod["sem_r"]["loglik"], mod["sem_r"]["aic"], mod["sem_r"]["moran_residus"]))319 print("== SDM sans secteur : ρ %.4f (%.4f), logL %.2f, AIC %.1f, I résidus %.3f" % (mod["sdm_r"]["rho"], mod["sdm_r"]["se_rho"], mod["sdm_r"]["loglik"], mod["sdm_r"]["aic"], mod["sdm_r"]["moran_residus"]))320 print("== MCO sans secteur : logL %.2f, AIC %.1f ; MCO référence : logL %.2f, AIC %.1f" % (mod["mco_r"]["loglik"], mod["mco_r"]["aic"], mod["mco_c"]["loglik"], mod["mco_c"]["aic"]))321 print("== SAR référence : ρ %.4f (%.4f) ; SEM référence : λ %.4f (%.4f)" % (mod["sar_c"]["rho"], mod["sar_c"]["se_rho"], mod["sem_c"]["lambda"], mod["sem_c"]["se_lambda"]))322 cols = list(b["Xr"].columns)323 tab = pd.DataFrame({"MCO réf.": pd.Series(m["beta"], index=b["Xc"].columns).reindex(cols), "MCO": mod["mco_r"]["beta"], "SAR": mod["sar_r"]["beta"], "SEM": mod["sem_r"]["beta"], "SDM β": mod["sdm_r"]["beta"][:len(cols)],324 "SDM θ": np.concatenate([[np.nan], mod["sdm_r"]["beta"][len(cols):]]),325 "se MCO": mod["mco_r"]["se_hc1"], "se SAR": mod["sar_r"]["se"], "se SEM": mod["sem_r"]["se"], "se SDM": mod["sdm_r"]["se"][:len(cols)]}, index=cols)326 print("== Coefficients :\n", tab.round(4).to_string())327 print("== Effets directs/indirects/totaux :\n", res["effets"].round(4).to_string())328 g = res["gwr"]329 print("== GWR : bande %.3f km (CV %.4f), tr(S) %.1f, AICc %.1f (MCO sans secteur %.1f ; référence %.1f), R² %.3f, σ %.4f" % (g["selection"]["bande"], g["selection"]["cv"], g["gwr"]["trace_S"], g["gwr"]["aicc"], g["aicc_mco_sans_secteur"], g["aicc_mco_reference"], g["gwr"]["r2"], g["gwr"]["sigma"]))330 print(" quartiles des coefficients locaux :\n", g["quartiles"].round(4).to_string())331 print(" constante locale médiane par secteur (écart à Hull) vs indicatrice vs générateur :\n", g["comparaison_secteurs"].round(4).to_string())332 print(" coefficients locaux au sujet :", dict(zip(cols, g["sujet_beta"].round(4))))333 print("== Effets fixes :\n", res["effets_fixes"].round(4).to_string(index=False))334 print(" indice trimestriel (log) :", res["effets_fixes"].attrs["indice_trimestriel"].round(3).to_dict())335 vk = res["variogramme"]336 print("== Variogramme (sans secteur) : h", vk["variogramme"]["h"].round(2), "\n γ", vk["variogramme"]["gamma"].round(5), "\n N", vk["variogramme"]["n"])337 print(" référence γ", vk["variogramme_reference"]["gamma"].round(5))338 print(" paramètres :", {k: round(v, 5) for k, v in vk["parametres"].items()}, "; var(e) %.5f ; σ² réf %.5f" % (vk["variance_residus"], vk["sigma2_reference"]))339 print(" krigeage au sujet : %.4f (variance %.5f)" % (vk["sujet"]["prediction"][0], vk["sujet"]["variance"][0]))340 st = res["spatio_temporel"]341 print("== Spatio-temporel : %d ventes retenues (%d sans voisin antérieur), %.1f voisins en moyenne, max |valeur propre| %.1e" % (st["n_retenues"], st["n_sans_voisin"], st["voisins_moyen"], st["max_abs_valeur_propre"]))342 print(" I(e) : W_st %.3f (z %.1f) ; knn %.3f ; distance %.3f" % (st["moran_e_st"]["I"], st["moran_e_st"]["z"], st["moran_e_knn"]["I"], st["moran_e_dist"]["I"]))343 print(" ρ_st (MCO sur ventes antérieures) %.4f (%.4f), R² %.3f, I résidus après décalage %.3f ; SEM-st λ %.4f (%.4f), AIC %.1f" % (st["rho_st"], st["se_rho_st"], st["r2_st"], st["moran_e_apres_lag"], st["sem_st"]["lambda"], st["sem_st"]["se_lambda"], st["sem_st"]["aic"]))344 print("== Validation croisée :\n", res["validation"].round(4).to_string(index=False))345 print("== Sujet :\n", res["sujet"].round(4).to_string(index=False)); print(" voisins du sujet :", res["sujet"].attrs["voisins_sujet"], res["sujet"].attrs["secteurs_voisins"])346 print("== Sensibilité :\n", res["sensibilite"].round(4).to_string(index=False))347 if "controle" in res:348 print("== Contrôle croisé spreg/esda/mgwr :", {k: round(v, 5) for k, v in res["controle"].items()})349350if __name__ == "__main__":351 main(controle="--controle" in sys.argv)352