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 # effet de localisation local = constante locale + coefficient local du gradient × distance (au-delà de 3 km), par secteur,129 # comparé à l'effet de localisation du modèle de référence (indicatrice + gradient) et à celui du générateur (prime − 0,006·max(d − 3, 0))130 sect = b["ventes"]["secteur"].values; dist = b["Xr"]["dist_centre"].values; j = cols.index("dist_centre")131 loc_gwr = g["beta"][:, 0] + g["beta"][:, j] * dist132 coef = b["modele"]["coef"]["coef"]133 loc_ref = np.array([coef.get(f"sect_{sc}", 0.0) for sc in sect]) + coef["dist_centre"] * dist134 loc_dgp = np.array([PRIMES_DGP[sc] for sc in sect]) - 0.006 * dist135 comp = pd.DataFrame({"gwr": loc_gwr, "reference": loc_ref, "generateur": loc_dgp, "secteur": sect}).groupby("secteur").median()136 comp["gwr"] -= float(np.median(g["beta"][sect == "Hull", 0])) # écart à la constante locale médiane de Hull137 comp = comp.rename(columns={"gwr": "gwr_localisation_rel_hull", "reference": "reference_localisation", "generateur": "generateur_localisation"})138 comp.loc["Hull", "gwr_localisation_rel_hull"] = 0.0139 aicc_mco_r = _aicc(b["mco_r"]["residus"], X.shape[1]); aicc_mco_c = _aicc(b["modele"]["residus"], b["Xc"].shape[1])140 return {"selection": sel, "gwr": g, "beta": beta, "quartiles": quart, "comparaison_secteurs": comp, "sujet_beta": gs["beta"][0],141 "aicc_mco_sans_secteur": aicc_mco_r, "aicc_mco_reference": aicc_mco_c}142143# ----------------------------------------------------------------------------- 7. effets fixes144145def effets_fixes(b: dict, W: np.ndarray) -> pd.DataFrame:146 y, xy, v = b["y"], b["xy"], b["ventes"]; n = len(y)147 def ligne(nom, Xm, k_spatial, k_temps):148 r = h.mco(Xm, y); k = Xm.shape[1]149 return {"specification": nom, "parametres": k, "effets_spatiaux": k_spatial, "effets_temps": k_temps, "r2": r["r2"],150 "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"]}151 rows = [ligne("Gradient de distance seul (sans secteur) ; tendance linéaire", b["Xr"].values, 1, 1),152 ligne("Secteurs et gradient (référence) ; tendance linéaire", b["Xc"].values, 5, 1)]153 for cell in [5.0, 3.0, 2.0]:154 gx = np.floor(xy[:, 0] / cell); gy = np.floor(xy[:, 1] / cell); cid = pd.factorize(gx * 1000 + gy)[0]155 D = pd.get_dummies(cid, drop_first=True).astype(float).values156 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))157 q = pd.PeriodIndex(v["date_vente"], freq="Q").astype(str)158 Dq = pd.get_dummies(q, drop_first=True).astype(float).values159 rows.append(ligne(f"Secteurs ; indicatrices trimestrielles ({Dq.shape[1]})", np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), 4, Dq.shape[1]))160 rt = h.mco(np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), y)161 df = pd.DataFrame(rows)162 df.attrs["indice_trimestriel"] = pd.Series(np.concatenate([[0.0], rt["beta"][-Dq.shape[1]:]]), index=sorted(q.unique()))163 df.attrs["tendance_mensuelle"] = float(b["modele"]["coef"].loc["mois", "coef"])164 return df165166# ----------------------------------------------------------------------------- 8. variogramme et krigeage167168def variogramme_krigeage(b: dict) -> dict:169 e = b["mco_r"]["residus"]; xy = b["xy"]170 vg = h.variogramme_empirique(xy, e, BORNES_VARIOGRAMME)171 vg_ref = h.variogramme_empirique(xy, b["modele"]["residus"], BORNES_VARIOGRAMME)172 p = h.ajuster_variogramme(vg["h"], vg["gamma"], vg["n"], portee_max=5.0)173 loo = h.krigeage_validation_croisee(xy, e, p, VOISINS_KRIGEAGE)174 sujet = h.krigeage_ordinaire(xy, e, b["cs"], p, VOISINS_KRIGEAGE)175 # surface pour la figure : grille régulière (km)176 gx = np.linspace(xy[:, 0].min() - 1, xy[:, 0].max() + 1, 90); gy = np.linspace(xy[:, 1].min() - 1, xy[:, 1].max() + 1, 70)177 GX, GY = np.meshgrid(gx, gy); grille = np.column_stack([GX.ravel(), GY.ravel()])178 kg = h.krigeage_ordinaire(xy, e, grille, p, VOISINS_KRIGEAGE)179 return {"variogramme": vg, "variogramme_reference": vg_ref, "parametres": p, "loo": loo, "sujet": sujet,180 "grille": {"x": gx, "y": gy, "z": kg["prediction"].reshape(GX.shape), "var": kg["variance"].reshape(GX.shape)},181 "variance_residus": float(e.var()), "sigma2_reference": float(b["modele"]["sigma"] ** 2)}182183# ----------------------------------------------------------------------------- 9. spatio-temporel184185def spatio_temporel(b: dict, Ws: dict) -> dict:186 y, Xr, Wst = b["y"], b["Xr"].values, Ws["st"]187 e = b["mco_r"]["residus"]188 ok = Wst.sum(1) > 0 # ventes ayant au moins une vente antérieure voisine189 Wy = Wst @ y190 # décalage spatio-temporel : prédéterminé (matrice triangulaire après tri par date) → MV = MCO sur [W_st y, X]191 Z = np.column_stack([Wy[ok], Xr[ok]]); r = h.mco(Z, y[ok])192 sem_st = h.sem_mv(y, Xr, Wst)193 return {"n_retenues": int(ok.sum()), "n_sans_voisin": int((~ok).sum()), "voisins_moyen": float((Wst[ok] > 0).sum(1).mean()),194 "rho_st": float(r["beta"][0]), "se_rho_st": float(r["se_hc1"][0]), "beta_st": r["beta"][1:], "r2_st": r["r2"],195 "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"]),196 "sem_st": sem_st, "max_abs_valeur_propre": float(np.abs(np.linalg.eigvals(Wst)).max()),197 "moran_e_apres_lag": h.moran_i(r["residus"], Wst[np.ix_(ok, ok)])["I"]}198199# ----------------------------------------------------------------------------- 10. validation croisée et sujet200201def validation_croisee(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:202 y = b["y"]; n = len(y)203 def loo_mco(Xm):204 r = h.mco(Xm, y); H = Xm @ np.linalg.pinv(Xm.T @ Xm) @ Xm.T; return r["residus"] / (1 - np.diag(H))205 sem = mod["sem_r"]; u = sem["residus"]; B = np.eye(n) - sem["lambda"] * W; Q = B.T @ B206 Eu = np.array([-(Q[i] @ u - Q[i, i] * u[i]) / Q[i, i] for i in range(n)])207 series = {"MCO, modèle de référence (secteurs)": loo_mco(b["Xc"].values),208 "MCO, sans indicatrices de secteur": loo_mco(b["Xr"].values),209 "MCO sans secteur + krigeage des résidus": b["mco_r"]["residus"] - vk["loo"],210 "SEM sans secteur (espérance conditionnelle)": u - Eu,211 "SAR sans secteur (prédiction conditionnelle)": y - (mod["sar_r"]["rho"] * (W @ y) + b["Xr"].values @ mod["sar_r"]["beta"]),212 "GWR sans secteur (validation croisée)": None}213 rows = []214 for nom, r in series.items():215 if r is None:216 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"]})217 else:218 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"]})219 return pd.DataFrame(rows)220221def predictions_sujet(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:222 s, xs = b["sujet"], b["xs_r"]; y = b["y"]; n = len(y)223 duan = lambda res: float(np.mean(np.exp(res)))224 ref = float(cas1.predire(b["modele"], cas1.sujet_df(s))[0])225 mco_r = float(np.exp(xs @ b["mco_r"]["beta"]) * duan(b["mco_r"]["residus"]))226 krig = float(np.exp(xs @ b["mco_r"]["beta"] + vk["sujet"]["prediction"][0]) * duan(b["mco_r"]["residus"] - vk["loo"]))227 # SEM : espérance conditionnelle de u au sujet (matrice augmentée : le sujet n'est voisin d'aucune vente)228 xy_aug = np.vstack([b["xy"], b["cs"]]); Wa = np.zeros((n + 1, n + 1)); Wa[:n, :n] = W229 Wa[n, :n] = h.matrice_poids_knn(xy_aug, K_VOISINS)[n, :n]; Wa[n] /= Wa[n].sum()230 sem = mod["sem_r"]; u_s = float(h.sem_esperance_conditionnelle(sem["lambda"], Wa, sem["residus"], [n])[0])231 sem_pred = float(np.exp(xs @ sem["beta"] + u_s) * duan(sem["innovations"]))232 # SAR : prédiction conditionnelle avec le décalage observé (moyenne des ln P des k ventes voisines)233 idx = np.nonzero(Wa[n, :n])[0]; Wy_s = float(Wa[n, :n] @ y)234 sar = mod["sar_r"]; sar_pred = float(np.exp(sar["rho"] * Wy_s + xs @ sar["beta"]) * duan(sar["residus"]))235 gwr_pred = float(np.exp(xs @ g["sujet_beta"]) * duan(g["gwr"]["residus"]))236 df = pd.DataFrame([237 {"modele": "MCO, modèle de référence (indicatrices de secteur)", "ln_correction": 0.0, "indication": ref},238 {"modele": "MCO sans indicatrices de secteur", "ln_correction": 0.0, "indication": mco_r},239 {"modele": "MCO sans secteur + krigeage du résidu au sujet", "ln_correction": float(vk["sujet"]["prediction"][0]), "indication": krig},240 {"modele": "SEM sans secteur (espérance conditionnelle du résidu)", "ln_correction": u_s, "indication": sem_pred},241 {"modele": "SAR sans secteur (décalage des 8 ventes voisines)", "ln_correction": float(sar["rho"] * Wy_s), "indication": sar_pred},242 {"modele": "GWR sans secteur (coefficients locaux au sujet)", "ln_correction": float(g["sujet_beta"][0] - b["mco_r"]["beta"][0]), "indication": gwr_pred}])243 df["ecart_log_vs_mco_sans_secteur"] = np.log(df["indication"] / mco_r)244 df.attrs["voisins_sujet"] = list(b["ventes"]["id_vente"].values[idx]); df.attrs["secteurs_voisins"] = list(b["ventes"]["secteur"].values[idx])245 df.attrs["variance_krigeage_sujet"] = float(vk["sujet"]["variance"][0])246 return df247248# ----------------------------------------------------------------------------- 11. sensibilité249250def sensibilite(b: dict) -> pd.DataFrame:251 y, Xr, xy, e = b["y"], b["Xr"].values, b["xy"], b["mco_r"]["residus"]252 specs = [("k plus proches voisins", k, h.matrice_poids_knn(xy, k)) for k in [4, 6, 8, 12, 16]]253 specs += [("distance inverse", d, h.matrice_poids_distance_inverse(xy, d)) for d in [2.0, 3.0, 4.0]]254 rows = []255 for fam, par, W in specs:256 mi = h.moran_i(e, W); sar = h.sar_mv(y, Xr, W); sem = h.sem_mv(y, Xr, W)257 rows.append({"famille": fam, "parametre": par, "voisins_moyen": float((W > 0).sum(1).mean()), "lignes_vides": int((W.sum(1) == 0).sum()), "I": mi["I"], "z": mi["z"],258 "rho": sar["rho"], "se_rho": sar["se_rho"], "aic_sar": sar["aic"], "lambda": sem["lambda"], "se_lambda": sem["se_lambda"], "aic_sem": sem["aic"]})259 return pd.DataFrame(rows)260261262# ----------------------------------------------------------------------------- exercice 19.7263264def exercice_7(b: dict | None = None) -> dict:265 """SEM sans secteur avec 4, 8 et 16 voisins : λ, espérance conditionnelle du résidu au sujet et indication ;266 matrice spatio-temporelle avec fenêtre de 24 mois : ventes sans voisine antérieure et ρ_st."""267 b = b or charger_base(); y, Xr, xy, xs, n = b["y"], b["Xr"].values, b["xy"], b["xs_r"], len(b["y"])268 duan = lambda res: float(np.mean(np.exp(res)))269 rows = []270 for k in [4, 8, 16]:271 W = h.matrice_poids_knn(xy, k); sem = h.sem_mv(y, Xr, W)272 xy_aug = np.vstack([xy, b["cs"]]); Wa = np.zeros((n + 1, n + 1)); Wa[:n, :n] = W273 Wa[n, :n] = h.matrice_poids_knn(xy_aug, k)[n, :n]; Wa[n] /= Wa[n].sum()274 u_s = float(h.sem_esperance_conditionnelle(sem["lambda"], Wa, sem["residus"], [n])[0])275 rows.append({"k": k, "lambda": sem["lambda"], "u_sujet": u_s, "indication": float(np.exp(xs @ sem["beta"] + u_s) * duan(sem["innovations"]))})276 sem_k = pd.DataFrame(rows)277 out = {"sem_par_k": sem_k}278 for fen in [12, 24]:279 Wst = h.matrice_poids_spatio_temporelle(xy, b["mois"], PORTEE_KM, fen); ok = Wst.sum(1) > 0280 Z = np.column_stack([(Wst @ y)[ok], Xr[ok]]); r = h.mco(Z, y[ok])281 out[f"st_{fen}"] = {"sans_voisin": int((~ok).sum()), "voisins_moyen": float((Wst[ok] > 0).sum(1).mean()), "rho_st": float(r["beta"][0]), "se": float(r["se_hc1"][0])}282 return out283284# ----------------------------------------------------------------------------- contrôle croisé facultatif285286def controle_croise(b: dict, W: np.ndarray, g: dict) -> dict:287 """Compare les estimations maison à spreg/esda/mgwr (installés sur la machine de rédaction)."""288 import spreg289 from libpysal.weights import W as Wcls290 from esda.moran import Moran291 from mgwr.sel_bw import Sel_BW292 from mgwr.gwr import GWR293 y, X = b["y"], b["Xr"]294 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}295 w = Wcls(neigh, wts); w.transform = "r"; Xs = X.drop(columns=["const"]).values; yy = y.reshape(-1, 1)296 with warnings.catch_warnings():297 warnings.simplefilter("ignore")298 ml = spreg.ML_Lag(yy, Xs, w, method="full"); me = spreg.ML_Error(yy, Xs, w, method="full")299 ols = spreg.OLS(yy, Xs, w=w, spat_diag=True, moran=True)300 mo = Moran(b["mco_r"]["residus"], w, permutations=0)301 bw = Sel_BW(b["xy"], yy, Xs, fixed=True, kernel="gaussian").search(criterion="CV")302 gm = GWR(b["xy"], yy, Xs, bw, fixed=True, kernel="gaussian").fit()303 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)304 return {"rho_spreg": float(ml.rho), "rho_maison": mine_sar["rho"], "lambda_spreg": float(me.lam), "lambda_maison": mine_sem["lambda"],305 "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],306 "moran_esda_z": float(mo.z_norm), "moran_maison_z": h.moran_i(b["mco_r"]["residus"], W)["z"],307 "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"]}308309# ----------------------------------------------------------------------------- main310311def main(verbose: bool = True, controle: bool = False) -> dict:312 if "res" in _CACHE and not controle:313 return _CACHE["res"]314 b = charger_base(); Ws = matrices(b); W = Ws["knn"]315 res = {"base": b, "W": Ws, "moran": tableau_moran(b, Ws), "correlogramme": correlogramme(b), "lm": tests_lm(b, W)}316 res["modeles"] = modeles_spatiaux(b, W)317 res["effets"] = effets(res["modeles"], W, list(b["Xr"].columns))318 res["gwr"] = gwr_resultats(b)319 res["effets_fixes"] = effets_fixes(b, W)320 res["variogramme"] = variogramme_krigeage(b)321 res["spatio_temporel"] = spatio_temporel(b, Ws)322 res["validation"] = validation_croisee(b, res["modeles"], res["gwr"], res["variogramme"], W)323 res["sujet"] = predictions_sujet(b, res["modeles"], res["gwr"], res["variogramme"], W)324 res["sensibilite"] = sensibilite(b)325 if controle:326 res["controle"] = controle_croise(b, W, res["gwr"])327 _CACHE["res"] = res328 if verbose:329 imprimer(res)330 return res331332def imprimer(res: dict):333 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)334 b = res["base"]; m = b["modele"]335 print("== Modèle de référence : R² %.3f, σ %.4f, n %d, k %d" % (m["r2"], m["sigma"], m["n"], m["k"]))336 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"]))337 print("== Étendue des coordonnées (km) :", b["xy"].min(0).round(1), b["xy"].max(0).round(1), "; sujet :", b["cs"].round(2))338 for wn, W in res["W"].items():339 print(f" matrice {wn} : voisins moyens {(W > 0).sum(1).mean():.1f}, lignes vides {(W.sum(1) == 0).sum()}")340 print("== I de Moran :\n", res["moran"].round(4).to_string(index=False))341 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"]])342 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"]])343 for k, v in res["lm"].items():344 print(f"== Tests LM ({k}) :", {a: (round(s, 2), round(p, 4)) for a, (s, p) in v.items()})345 mod = res["modeles"]346 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"]))347 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"]))348 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"]))349 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"]))350 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"]))351 cols = list(b["Xr"].columns)352 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)],353 "SDM θ": np.concatenate([[np.nan], mod["sdm_r"]["beta"][len(cols):]]),354 "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)355 print("== Coefficients :\n", tab.round(4).to_string())356 print("== Effets directs/indirects/totaux :\n", res["effets"].round(4).to_string())357 g = res["gwr"]358 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"]))359 print(" quartiles des coefficients locaux :\n", g["quartiles"].round(4).to_string())360 print(" constante locale médiane par secteur (écart à Hull) vs indicatrice vs générateur :\n", g["comparaison_secteurs"].round(4).to_string())361 print(" coefficients locaux au sujet :", dict(zip(cols, g["sujet_beta"].round(4))))362 print("== Effets fixes :\n", res["effets_fixes"].round(4).to_string(index=False))363 print(" indice trimestriel (log) :", res["effets_fixes"].attrs["indice_trimestriel"].round(3).to_dict())364 vk = res["variogramme"]365 print("== Variogramme (sans secteur) : h", vk["variogramme"]["h"].round(2), "\n γ", vk["variogramme"]["gamma"].round(5), "\n N", vk["variogramme"]["n"])366 print(" référence γ", vk["variogramme_reference"]["gamma"].round(5))367 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"]))368 print(" krigeage au sujet : %.4f (variance %.5f)" % (vk["sujet"]["prediction"][0], vk["sujet"]["variance"][0]))369 st = res["spatio_temporel"]370 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"]))371 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"]))372 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"]))373 print("== Validation croisée :\n", res["validation"].round(4).to_string(index=False))374 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"])375 print("== Sensibilité :\n", res["sensibilite"].round(4).to_string(index=False))376 if "controle" in res:377 print("== Contrôle croisé spreg/esda/mgwr :", {k: round(v, 5) for k, v in res["controle"].items()})378379if __name__ == "__main__":380 main(controle="--controle" in sys.argv)381