spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 14 — Nettoyage, géocodage et qualité des données (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES).341. `generer_brut()` injecte des anomalies déclarées dans le jeu propre `cas1_gatineau_ventes.csv` (graine 20261031) et écrit5 `data/synthetic/cas1_gatineau_ventes_brut.csv` (618 inscriptions), `cas1_gatineau_anomalies_injectees.csv` (table de vérité)6 et `cas1_gatineau_role_extrait.csv` (extrait de rôle fictif pour l'appariement).72. `pipeline()` applique les règles déterministes (R1-R7), le géocodage de repli (G1), l'imputation (I1), les scores robustes (S1),8 les diagnostics de régression (levier, résidus studentisés, distance de Cook ; exclusion E1) et écrit le journal d'exclusions9 `data/synthetic/cas1_gatineau_journal_exclusions.csv`.103. `comparer_modeles()` estime `traite.cas1.hedonique` sur le jeu brut (coercé), après chaque étape et sur le jeu propre.114. `simulation_geocodage()`, `appariement_role()`, `comparer_imputations()`, `sensibilite_seuils()` : sections 14.3 à 14.6.1213Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch14_nettoyage.py14Toutes les lois d'injection sont documentées dans data/synthetic/README-brut.md."""15from __future__ import annotations16import sys, json, pathlib, datetime17import numpy as np, pandas as pd1819sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))20from traite import cas1, hedonique as h2122ROOT = pathlib.Path(__file__).resolve().parents[3]23D = ROOT / "data/synthetic"24GRAINE = 2026103125PI2_PAR_M2 = 10.763926BBOX = {"lat": (45.20, 45.80), "lon": (-76.10, -75.20)} # emprise plausible du jeu (Gatineau et périphérie)27CENTRES = {"Hull": (45.428, -75.714), "Gatineau": (45.477, -75.650), "Aylmer": (45.395, -75.843),28 "Buckingham": (45.588, -75.418), "Masson-Angers": (45.545, -75.420)} # centres de secteur (synth.py)29SEUIL_PRIX_MIN = 5000.0 # ventes < 5 000 $ écartées (convention empruntée au Règlement sur la proportion médiane)30SEUIL_S_MAX_M2 = 450.0 # superficie habitable implausible pour une unifamiliale du jeu → test d'unité31SEUIL_MAD = 3.5 # |z robuste| > 3,532SEUIL_T = 3.0 # |résidu studentisé externe| > 333C_MAD = 1.4826 # 1/Φ^{-1}(0,75)3435# ----------------------------------------------------------------------------- 1. jeu brut3637def _rues(rng, n):38 noms = ["rue des Érables", "rue Laurier", "boulevard Saint-Joseph", "rue Principale", "chemin de la Montagne", "rue Notre-Dame",39 "rue Jacques-Cartier", "avenue Gatineau", "rue de Lorraine", "chemin d'Aylmer", "rue Georges", "rue Maclaren"]40 return [f"{int(rng.integers(10, 1999))}, {noms[int(rng.integers(0, len(noms)))]} (fictive)" for _ in range(n)]4142def generer_brut(seed: int = GRAINE) -> dict:43 """Jeu brut = jeu propre + anomalies déclarées. Retourne brut, table de vérité, extrait de rôle (DataFrames)."""44 rng = np.random.default_rng(seed)45 v = cas1.charger()["ventes"].copy()46 n = len(v)47 # --- identifiants : matricule fictif canonique « SS-FF-NNNN-C » et adresse fictive (jamais réelles)48 code_sect = {"Hull": "05", "Gatineau": "07", "Aylmer": "03", "Buckingham": "11", "Masson-Angers": "12"}49 feuillet = rng.integers(10, 99, n); numero = rng.choice(np.arange(1000, 9999), n, replace=False); verif = rng.integers(0, 9, n)50 v["matricule"] = [f"{code_sect[s]}-{f:02d}-{m:04d}-{c}" for s, f, m, c in zip(v["secteur"], feuillet, numero, verif)]51 v["adresse_fictive"] = _rues(rng, n)52 # --- extrait de rôle fictif (une ligne par unité d'évaluation, format canonique)53 role = v[["matricule", "adresse_fictive", "secteur", "superficie_terrain_m2", "annee_construction", "valeur_role_2025"]].copy()54 role = role.rename(columns={"superficie_terrain_m2": "role_superficie_terrain_m2", "annee_construction": "role_annee_construction"})55 # --- transactions : dates promesse / acte56 delai = rng.integers(21, 91, n)57 b = pd.DataFrame({58 "id_inscription": [f"B{i+1:04d}" for i in range(n)], "id_vente_origine": v["id_vente"], "matricule": v["matricule"],59 "adresse_fictive": v["adresse_fictive"], "date_promesse": v["date_vente"], "date_acte": v["date_vente"] + pd.to_timedelta(delai, unit="D"),60 "secteur": v["secteur"], "type": v["type"], "superficie_habitable_m2": v["superficie_habitable_m2"].astype(float),61 "superficie_terrain_m2": v["superficie_terrain_m2"].astype(float), "annee_construction": v["annee_construction"].astype(float),62 "qualite_1a5": v["qualite_1a5"].astype(float), "garage_places": v["garage_places"], "sous_sol_fini": v["sous_sol_fini"], "piscine": v["piscine"],63 "chambres": v["chambres"], "salles_de_bain": v["salles_de_bain"], "lat": v["lat"].astype(float), "lon": v["lon"].astype(float),64 "valeur_role_2025": v["valeur_role_2025"], "montant_meubles": 0.0, "type_transaction": "ordinaire", "prix_vente": v["prix_vente"].astype(float)})65 # --- tirage des lignes touchées (disjointes) : 109 lignes + 18 doublons ajoutés66 plan = [("prix_symbolique", 5), ("meubles_inclus", 8), ("superficie_pi2", 12), ("dates_inversees", 10), ("coord_nulles", 9),67 ("coord_hors_zone", 6), ("vente_liee_declaree", 8), ("vente_liee_non_declaree", 6), ("terrain_manquant", 20),68 ("qualite_manquante", 15), ("annee_manquante", 10)]69 idx = rng.permutation(n); pos = 0; verite = []70 def prendre(k):71 nonlocal pos; out = idx[pos:pos + k]; pos += k; return out72 for code, k in plan:73 rows = prendre(k)74 for i in rows:75 detail = ""76 if code == "prix_symbolique":77 detail = f"prix {b.at[i,'prix_vente']:.0f} → 1"; b.at[i, "prix_vente"] = 1.078 elif code == "meubles_inclus":79 m = float(np.round(rng.uniform(8000, 35000) / 500) * 500); b.at[i, "montant_meubles"] = m; b.at[i, "prix_vente"] += m80 detail = f"meubles {m:.0f} ajoutés au prix"81 elif code == "superficie_pi2":82 s = b.at[i, "superficie_habitable_m2"]; b.at[i, "superficie_habitable_m2"] = float(np.round(s * PI2_PAR_M2)); detail = f"{s:.0f} m² saisis en pi²"83 elif code == "dates_inversees":84 b.at[i, "date_promesse"], b.at[i, "date_acte"] = b.at[i, "date_acte"], b.at[i, "date_promesse"]; detail = "promesse ↔ acte"85 elif code == "coord_nulles":86 b.at[i, "lat"] = 0.0; b.at[i, "lon"] = 0.0; detail = "lat/lon = 0"87 elif code == "coord_hors_zone":88 if rng.random() < 0.5:89 b.at[i, "lon"] = -b.at[i, "lon"]; detail = "signe de la longitude inversé"90 else:91 b.at[i, "lat"], b.at[i, "lon"] = b.at[i, "lon"], b.at[i, "lat"]; detail = "lat ↔ lon"92 elif code in ("vente_liee_declaree", "vente_liee_non_declaree"):93 f = float(rng.uniform(0.40, 0.60)); p = b.at[i, "prix_vente"]; b.at[i, "prix_vente"] = float(np.round(p * f / 500) * 500)94 if code == "vente_liee_declaree":95 b.at[i, "type_transaction"] = "entre personnes liées" if rng.random() < 0.65 else "succession"96 detail = f"prix × {f:.2f}" + ("" if code == "vente_liee_declaree" else " (non déclarée)")97 elif code == "terrain_manquant":98 b.at[i, "superficie_terrain_m2"] = np.nan; detail = "terrain NA"99 elif code == "qualite_manquante":100 b.at[i, "qualite_1a5"] = np.nan; detail = "qualité NA"101 elif code == "annee_manquante":102 b.at[i, "annee_construction"] = np.nan; detail = "année NA"103 verite.append({"id_inscription": b.at[i, "id_inscription"], "id_vente_origine": b.at[i, "id_vente_origine"], "anomalie": code, "detail": detail})104 # --- doublons d'inscription : 18 copies (12 identiques, 6 avec acte décalé de 1 à 3 jours), ajoutées en fin de fichier105 dup_src = rng.choice(np.setdiff1d(np.arange(n), idx[:pos]), 18, replace=False)106 dups = b.loc[dup_src].copy(); dups["id_inscription"] = [f"B{n + j + 1:04d}" for j in range(18)]107 for j, (k, i) in enumerate(zip(dups.index, dup_src)):108 if j >= 12:109 dups.at[k, "date_acte"] = dups.at[k, "date_acte"] + pd.Timedelta(days=int(rng.integers(1, 4)))110 verite.append({"id_inscription": dups.at[k, "id_inscription"], "id_vente_origine": b.at[i, "id_vente_origine"], "anomalie": "doublon",111 "detail": "copie de " + b.at[i, "id_inscription"] + ("" if j < 12 else " (acte décalé)")})112 b = pd.concat([b, dups], ignore_index=True)113 # --- bruit de format sur la clé matricule des transactions (appariement, §14.4) : la vérité reste dans id_vente_origine114 u = rng.random(len(b)); mat = b["matricule"].tolist()115 for i, (x, m) in enumerate(zip(u, mat)):116 if x < 0.04: mat[i] = np.nan # clé absente117 elif x < 0.09: d = list(m.replace("-", "")); p = int(rng.integers(0, len(d))); d[p] = str((int(d[p]) + 1) % 10); mat[i] = "-".join(["".join(d[0:2]), "".join(d[2:4]), "".join(d[4:8]), d[8]]) # faute de frappe118 elif x < 0.19: mat[i] = m.replace("-", "") # sans séparateurs119 elif x < 0.44: mat[i] = m.replace("-", " ") # espaces120 b["matricule"] = mat121 b = b.sort_values(["date_acte", "id_inscription"]).reset_index(drop=True)122 for c in ("date_promesse", "date_acte"):123 b[c] = pd.to_datetime(b[c]).dt.strftime("%Y-%m-%d")124 verite = pd.DataFrame(verite)125 return {"brut": b, "verite": verite, "role": role.reset_index(drop=True)}126127def ecrire_brut(seed: int = GRAINE) -> dict:128 g = generer_brut(seed)129 g["brut"].to_csv(D / "cas1_gatineau_ventes_brut.csv", index=False)130 g["verite"].to_csv(D / "cas1_gatineau_anomalies_injectees.csv", index=False)131 g["role"].to_csv(D / "cas1_gatineau_role_extrait.csv", index=False)132 return g133134def charger_brut() -> dict:135 b = pd.read_csv(D / "cas1_gatineau_ventes_brut.csv", parse_dates=["date_promesse", "date_acte"])136 v = pd.read_csv(D / "cas1_gatineau_anomalies_injectees.csv")137 r = pd.read_csv(D / "cas1_gatineau_role_extrait.csv")138 return {"brut": b, "verite": v, "role": r}139140# ----------------------------------------------------------------------------- 2. pipeline141142def _dist_centre(lat, lon):143 return np.sqrt(((lat - 45.428) * 111) ** 2 + ((lon + 75.714) * 78) ** 2)144145def _journal(j, df, ids, etape, action, motif, champ="", avant=None, apres=None):146 for k, i in enumerate(ids):147 j.append({"etape": etape, "id_inscription": df.at[i, "id_inscription"], "action": action, "motif": motif, "champ": champ,148 "valeur_avant": None if avant is None else avant[k], "valeur_apres": None if apres is None else apres[k]})149150def regles(b: pd.DataFrame) -> dict:151 """R1 doublons, R2 dates inversées, R3 meubles, R4 prix symbolique, R5 unité de superficie, R6 coordonnées, R7 transaction non représentative."""152 df = b.copy(); j = []; etapes = [("0 lecture", len(df))]153 df["date_promesse"] = pd.to_datetime(df["date_promesse"]); df["date_acte"] = pd.to_datetime(df["date_acte"])154 # R2 — dates : une promesse postérieure à l'acte est impossible → permutation (correction)155 m = df["date_acte"] < df["date_promesse"]156 ids = df.index[m]157 _journal(j, df, ids, "R2", "correction", "date de promesse postérieure à l'acte : champs permutés", "date_promesse",158 df.loc[ids, "date_promesse"].dt.strftime("%Y-%m-%d").tolist(), df.loc[ids, "date_acte"].dt.strftime("%Y-%m-%d").tolist())159 df.loc[m, ["date_promesse", "date_acte"]] = df.loc[m, ["date_acte", "date_promesse"]].values160 df["delai_acte_jours"] = (df["date_acte"] - df["date_promesse"]).dt.days161 # R1 — doublons : même unité (clé normalisée ou, à défaut, attributs), même promesse, même prix → on garde l'inscription la plus récente162 df["cle_unite"] = df["matricule"].astype(str).str.replace(r"\D", "", regex=True)163 sans = df["cle_unite"].isin(["", "nan"])164 df.loc[sans, "cle_unite"] = ("ATTR" + df.loc[sans, "secteur"] + df.loc[sans, "superficie_habitable_m2"].astype(str) + df.loc[sans, "superficie_terrain_m2"].astype(str) + df.loc[sans, "annee_construction"].astype(str))165 df = df.sort_values(["date_acte", "id_inscription"])166 dup = df.duplicated(subset=["cle_unite", "date_promesse", "prix_vente"], keep="last")167 _journal(j, df, df.index[dup], "R1", "exclusion", "doublon d'inscription (même unité, même promesse, même prix) : inscription la plus récente conservée")168 df = df[~dup]; etapes.append(("R1 doublons", len(df)))169 # R3 — meubles : le prix de l'immeuble est le prix déclaré moins la contrepartie des meubles (correction)170 m = df["montant_meubles"].fillna(0) > 0; ids = df.index[m]171 _journal(j, df, ids, "R3", "correction", "contrepartie des meubles retranchée du prix", "prix_vente", df.loc[ids, "prix_vente"].tolist(),172 (df.loc[ids, "prix_vente"] - df.loc[ids, "montant_meubles"]).tolist())173 df["prix_vente"] = df["prix_vente"] - df["montant_meubles"].fillna(0)174 # R4 — prix symbolique175 m = df["prix_vente"] < SEUIL_PRIX_MIN176 _journal(j, df, df.index[m], "R4", "exclusion", f"prix inférieur à {SEUIL_PRIX_MIN:.0f} $ (prix symbolique ou donation)", "prix_vente", df.loc[m, "prix_vente"].tolist())177 df = df[~m]; etapes.append(("R4 prix symbolique", len(df)))178 # R5 — unité de superficie : valeur implausible dont la conversion pi² → m² est plausible179 s = df["superficie_habitable_m2"]; conv = s / PI2_PAR_M2180 m = (s > SEUIL_S_MAX_M2) & conv.between(40, 400); ids = df.index[m]181 _journal(j, df, ids, "R5", "correction", "superficie implausible (> 450 m²) et plausible une fois divisée par 10,7639 : saisie en pi²", "superficie_habitable_m2",182 s[m].tolist(), conv[m].round(0).tolist())183 df.loc[m, "superficie_habitable_m2"] = conv[m].round(0)184 m2 = (df["superficie_habitable_m2"] > SEUIL_S_MAX_M2) # implausible et non convertible → exclusion (aucun cas attendu)185 _journal(j, df, df.index[m2], "R5", "exclusion", "superficie implausible non convertible", "superficie_habitable_m2", df.loc[m2, "superficie_habitable_m2"].tolist())186 df = df[~m2]187 # R6 — coordonnées nulles ou hors emprise → invalidées (géocodage de repli à l'étape G1)188 m = (df["lat"].fillna(0) == 0) | (df["lon"].fillna(0) == 0) | ~df["lat"].between(*BBOX["lat"]) | ~df["lon"].between(*BBOX["lon"])189 ids = df.index[m]190 _journal(j, df, ids, "R6", "correction", "coordonnées nulles ou hors de l'emprise du marché : invalidées, à géocoder", "lat;lon",191 [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])192 df["coord_valide"] = ~m; df.loc[m, ["lat", "lon"]] = np.nan193 # R7 — conditions de vente non représentatives déclarées194 m = df["type_transaction"].fillna("ordinaire").str.lower() != "ordinaire"195 _journal(j, df, df.index[m], "R7", "exclusion", "transaction déclarée non représentative (personnes liées, succession, reprise)", "type_transaction", df.loc[m, "type_transaction"].tolist())196 df = df[~m]; etapes.append(("R7 non représentatives", len(df)))197 df["date_vente"] = df["date_promesse"]198 df["age"] = 2026 - df["annee_construction"]199 return {"df": df, "journal": j, "etapes": etapes}200201def geocodage_repli(df: pd.DataFrame, j: list) -> pd.DataFrame:202 """G1 : sans adresse réelle à géocoder, les coordonnées invalides reçoivent le centre de leur secteur (précision « secteur »)."""203 df = df.copy(); m = ~df["coord_valide"]; ids = df.index[m]204 df.loc[m, "lat"] = [CENTRES[s][0] for s in df.loc[m, "secteur"]]; df.loc[m, "lon"] = [CENTRES[s][1] for s in df.loc[m, "secteur"]]205 df["precision_geocodage"] = np.where(m, "centre de secteur", "point d'origine")206 _journal(j, df, ids, "G1", "imputation", "géocodage de repli au centre du secteur (précision kilométrique) ; variables de distance à utiliser avec prudence", "lat;lon", None,207 [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])208 df["dist_centre_km"] = _dist_centre(df["lat"], df["lon"]).round(2)209 return df210211def imputer(df: pd.DataFrame, j: list | None = None, methode: str = "mediane") -> pd.DataFrame:212 """I1 : valeurs manquantes. methode ∈ {mediane (conditionnelle secteur × type), regression, indicateur, suppression}."""213 df = df.copy()214 champs = {"superficie_terrain_m2": "terrain", "qualite_1a5": "qualite", "annee_construction": "annee"}215 if methode == "suppression":216 m = df[list(champs)].isna().any(axis=1)217 if j is not None: _journal(j, df, df.index[m], "I1", "exclusion", "valeur manquante (suppression par liste)")218 df = df[~m].copy(); df["age"] = 2026 - df["annee_construction"]; return df219 for c in champs:220 m = df[c].isna()221 if not m.any(): continue222 if methode in ("mediane", "indicateur"):223 med = df.groupby(["secteur", "type"])[c].transform("median")224 val = med[m].fillna(df[c].median())225 elif methode == "regression":226 comp = df[~m]227 def X_of(d):228 X = pd.DataFrame({"const": 1.0, "ln_S": np.log(d["superficie_habitable_m2"]), "etages": (d["type"] == "à étages").astype(float)}, index=d.index)229 for s in cas1.SECTEURS[1:]: X[f"s_{s}"] = (d["secteur"] == s).astype(float)230 if c != "annee_construction": X["annee"] = d["annee_construction"].fillna(d["annee_construction"].median())231 if c != "superficie_terrain_m2": X["ln_T"] = np.log(d["superficie_terrain_m2"].fillna(d["superficie_terrain_m2"].median()))232 return X233 y = np.log(comp[c]) if c == "superficie_terrain_m2" else comp[c]234 r = h.mco(X_of(comp).values, y.values); pred = X_of(df[m]).values @ r["beta"]235 val = pd.Series(np.exp(pred) if c == "superficie_terrain_m2" else pred, index=df.index[m])236 if c == "qualite_1a5": val = val.round().clip(1, 5)237 if c == "annee_construction": val = val.round()238 if j is not None:239 _journal(j, df, df.index[m], "I1", "imputation", f"{champs[c]} manquant : imputation ({methode})", c, None, val.round(1).tolist())240 df.loc[m, c] = val.values241 if methode == "indicateur": df[f"manq_{champs[c]}"] = m.astype(float)242 df["age"] = 2026 - df["annee_construction"]243 return df244245def z_robuste(x: np.ndarray) -> np.ndarray:246 x = np.asarray(x, float); med = np.nanmedian(x); mad = np.nanmedian(np.abs(x - med))247 return (x - med) / (C_MAD * mad)248249def scores_robustes(df: pd.DataFrame, j: list | None = None, seuil: float = SEUIL_MAD) -> dict:250 """S1 : scores robustes (médiane / MAD) du ln du prix au m² par secteur, du ln de la superficie habitable et du ln du terrain."""251 df = df.copy()252 df["ln_prix_m2"] = np.log(df["prix_vente"] / df["superficie_habitable_m2"])253 df["z_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: z_robuste(s.values))254 df["z_superficie"] = z_robuste(np.log(df["superficie_habitable_m2"].values))255 df["z_terrain"] = df.groupby("secteur")["superficie_terrain_m2"].transform(lambda s: z_robuste(np.log(s.values)))256 df["z_classique_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: (s - s.mean()) / s.std(ddof=1))257 m = (df["z_prix_m2"].abs() > seuil) | (df["z_superficie"].abs() > seuil) | (df["z_terrain"].abs() > seuil)258 ids = df.index[m]259 if j is not None:260 _journal(j, df, ids, "S1", "exclusion", f"score robuste |z| > {seuil} (prix au m² par secteur, superficie ou terrain)", "z_prix_m2;z_superficie;z_terrain", None,261 [f"{a:.2f};{b_:.2f};{c:.2f}" for a, b_, c in zip(df.loc[ids, "z_prix_m2"], df.loc[ids, "z_superficie"], df.loc[ids, "z_terrain"])])262 return {"df": df, "exclus": df[m], "garde": df[~m]}263264def diagnostics(X: np.ndarray, y: np.ndarray) -> dict:265 """Levier h_ii, résidus studentisés internes r_i et externes t_i, distance de Cook D_i, DFFITS ; vérifie Σh = k."""266 X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape267 r = h.mco(X, y); e = r["residus"]; s2 = r["sigma2"]268 XtX_inv = np.linalg.pinv(X.T @ X)269 hii = np.einsum("ij,jk,ik->i", X, XtX_inv, X)270 r_int = e / np.sqrt(s2 * (1 - hii))271 s2_i = ((n - k) * s2 - e ** 2 / (1 - hii)) / (n - k - 1)272 t_ext = e / np.sqrt(s2_i * (1 - hii))273 cook = r_int ** 2 / k * hii / (1 - hii)274 dffits = t_ext * np.sqrt(hii / (1 - hii))275 return {"h": hii, "r_int": r_int, "t_ext": t_ext, "cook": cook, "dffits": dffits, "n": n, "k": k, "somme_h": float(hii.sum()), "mco": r}276277def cook_direct(X: np.ndarray, y: np.ndarray, i: int) -> float:278 """Distance de Cook par réestimation sans l'observation i (contrôle de l'identité algébrique)."""279 X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape280 r = h.mco(X, y); m = np.ones(n, bool); m[i] = False281 b_i = h.mco(X[m], y[m])["beta"]; d = r["beta"] - b_i282 return float(d @ (X.T @ X) @ d / (k * r["sigma2"]))283284def exclusions_regression(df: pd.DataFrame, j: list | None = None, seuil_t: float = SEUIL_T, max_iter: int = 2) -> dict:285 """E1 : ajuste le modèle de référence, exclut |t_i| > seuil, réitère (au plus max_iter fois)."""286 df = df.copy(); exclus = []; iterations = []287 for it in range(1, max_iter + 1):288 X = cas1.matrice(df); y = np.log(df["prix_vente"].astype(float))289 dg = diagnostics(X.values, y.values)290 df["levier"] = dg["h"]; df["t_ext"] = dg["t_ext"]; df["cook"] = dg["cook"]291 m = np.abs(dg["t_ext"]) > seuil_t292 iterations.append({"iteration": it, "n": len(df), "exclus": int(m.sum()), "sigma": float(np.sqrt(dg["mco"]["sigma2"])), "cook_max": float(dg["cook"].max()),293 "levier_max": float(dg["h"].max()), "n_cook_4n": int((dg["cook"] > 4 / len(df)).sum()), "n_levier_2kn": int((dg["h"] > 2 * dg["k"] / len(df)).sum())})294 if not m.any(): break295 ids = df.index[m]296 if j is not None:297 _journal(j, df, ids, "E1", "exclusion", f"résidu studentisé externe |t| > {seuil_t} (itération {it})", "t_ext;cook;levier", None,298 [f"{a:.2f};{c:.3f};{l:.3f}" for a, c, l in zip(df.loc[ids, "t_ext"], df.loc[ids, "cook"], df.loc[ids, "levier"])])299 exclus.append(df[m]); df = df[~m]300 return {"df": df, "exclus": pd.concat(exclus) if exclus else df.iloc[0:0], "iterations": iterations, "diag": dg}301302def pipeline(b: pd.DataFrame | None = None, seuil_mad: float = SEUIL_MAD, seuil_t: float = SEUIL_T, imputation: str = "mediane", ecrire: bool = False) -> dict:303 """Chaîne complète : règles → géocodage de repli → imputation → scores robustes → régression → jeu nettoyé + journal."""304 if b is None: b = charger_brut()["brut"]305 R = regles(b); j = R["journal"]; etapes = R["etapes"]306 df = geocodage_repli(R["df"], j)307 df = imputer(df, j, imputation); etapes.append(("I1 imputation", len(df)))308 S = scores_robustes(df, j, seuil_mad); df = S["garde"]; etapes.append(("S1 scores robustes", len(df)))309 E = exclusions_regression(df, j, seuil_t); df = E["df"]; etapes.append(("E1 régression", len(df)))310 jour = pd.DataFrame(j, columns=["etape", "id_inscription", "action", "motif", "champ", "valeur_avant", "valeur_apres"])311 jour.insert(0, "version_pipeline", "ch14_nettoyage.py 1.0 (2026-10-31)")312 if ecrire: jour.to_csv(D / "cas1_gatineau_journal_exclusions.csv", index=False)313 return {"propre": df, "journal": jour, "etapes": etapes, "regles": R, "scores": S, "regression": E, "apres_regles": R["df"]}314315# ----------------------------------------------------------------------------- 3. effet sur le modèle316317def _modele(df: pd.DataFrame) -> dict:318 d = df.copy()319 if "date_vente" not in d: d["date_vente"] = pd.to_datetime(d["date_promesse"])320 if "age" not in d or d["age"].isna().any(): d["age"] = 2026 - d["annee_construction"]321 if "dist_centre_km" not in d: d["dist_centre_km"] = _dist_centre(d["lat"].fillna(45.428), d["lon"].fillna(-75.714))322 m = cas1.hedonique(d)323 return m324325def brut_coerce(b: pd.DataFrame) -> pd.DataFrame:326 """Le jeu brut tel qu'un analyste pressé le passerait au modèle : lignes complètes seulement, aucune autre correction."""327 d = b.copy(); d["date_vente"] = pd.to_datetime(d["date_promesse"]); d["age"] = 2026 - d["annee_construction"]328 d["dist_centre_km"] = _dist_centre(d["lat"], d["lon"])329 d = d.dropna(subset=["superficie_terrain_m2", "qualite_1a5", "annee_construction"])330 return d331332def comparer_modeles(P: dict, b: pd.DataFrame) -> dict:333 propre_ref = cas1.charger()["ventes"]334 jeux = {"brut (lignes complètes)": brut_coerce(b), "après règles R1–R7 (+G1, I1)": imputer(geocodage_repli(P["regles"]["df"], []), None, "mediane"),335 "après scores robustes S1": P["scores"]["garde"], "après régression E1 (nettoyé)": P["propre"], "jeu propre de référence": propre_ref}336 rows = []; modeles = {}337 for nom, d in jeux.items():338 m = _modele(d); modeles[nom] = m; c = m["coef"]339 rows.append({"jeu": nom, "n": m["n"], "R2": m["r2"], "sigma": m["sigma"], "ln_S": c.loc["ln_S", "coef"], "se_ln_S": c.loc["ln_S", "se_hc1"], "ln_T": c.loc["ln_T", "coef"],340 "age": c.loc["age", "coef"], "qualite": c.loc["qualite", "coef"], "garage": c.loc["garage", "coef"], "mois_pct": 100 * m["croissance_mensuelle"],341 "sect_Aylmer": c.loc["sect_Aylmer", "coef"], "const": c.loc["const", "coef"]})342 tab = pd.DataFrame(rows).set_index("jeu")343 # prédiction pour le sujet344 s = cas1.charger()["sujet"]; sd = cas1.sujet_df(s)345 tab["prediction_sujet"] = [float(cas1.predire(modeles[n], sd)[0]) for n in tab.index]346 # écart nettoyé vs référence en erreurs-types347 ref = modeles["jeu propre de référence"]["coef"]; net = modeles["après régression E1 (nettoyé)"]["coef"]348 ecarts = pd.DataFrame({"ref": ref["coef"], "nettoye": net["coef"], "se_ref": ref["se_hc1"]}); ecarts["ecart_en_se"] = (ecarts["nettoye"] - ecarts["ref"]) / ecarts["se_ref"]349 return {"tableau": tab, "modeles": modeles, "ecarts": ecarts}350351# ----------------------------------------------------------------------------- 4. géocodage (simulation)352353TECHNIQUES = {"point d'adresse": 0.005, "centroïde de parcelle": 0.015, "interpolation le long de la rue": 0.040, "centroïde de code postal (urbain)": 0.150, "centroïde de code postal (périurbain)": 1.000} # écart-type radial, km — HYPOTHÈSES pédagogiques354355def simulation_geocodage(seed: int = GRAINE, k_voisins: int = 8) -> dict:356 """Erreur de position gaussienne isotrope d'écart-type σ (par technique) appliquée aux coordonnées propres (projetées en km) :357 erreur médiane, 95e centile, part des ventes dont l'appartenance au rayon de 3 km du sujet change, part des k voisins conservés,358 variation du coefficient de distance au centre du modèle hédonique."""359 rng = np.random.default_rng(seed)360 d = cas1.charger(); v = d["ventes"].copy(); s = d["sujet"]361 xy = h.coordonnees_km(v["lat"], v["lon"], 45.428, -75.714); xs = h.coordonnees_km([s["lat"]], [s["lon"]], 45.428, -75.714)[0]362 dist_s = np.sqrt(((xy - xs) ** 2).sum(1)); dans_rayon = dist_s <= 3.0363 W0 = h.matrice_poids_knn(xy, k_voisins, standardiser=False)364 m_ref = cas1.hedonique(v); beta_ref = m_ref["coef"].loc["dist_centre", "coef"]365 rows = []366 for nom, sig in TECHNIQUES.items():367 err = rng.normal(0, sig, size=xy.shape); xy2 = xy + err368 e = np.sqrt((err ** 2).sum(1))369 d2 = np.sqrt(((xy2 - xs) ** 2).sum(1)); change = (d2 <= 3.0) != dans_rayon370 W1 = h.matrice_poids_knn(xy2, k_voisins, standardiser=False)371 conserves = (W0 * W1).sum(1) / k_voisins372 v2 = v.copy(); v2["dist_centre_km"] = np.sqrt((xy2 ** 2).sum(1))373 m2 = cas1.hedonique(v2)374 rows.append({"technique": nom, "sigma_km": sig, "erreur_mediane_m": 1000 * np.median(e), "erreur_p95_m": 1000 * np.quantile(e, 0.95),375 "pct_rayon_change": 100 * change.mean(), "pct_voisins_conserves": 100 * conserves.mean(), "beta_dist": m2["coef"].loc["dist_centre", "coef"],376 "sigma_modele": m2["sigma"], "erreurs_m": 1000 * e})377 tab = pd.DataFrame(rows).set_index("technique")378 return {"tableau": tab.drop(columns="erreurs_m"), "erreurs": {r["technique"]: r["erreurs_m"] for r in rows}, "beta_dist_ref": beta_ref, "n_dans_rayon": int(dans_rayon.sum())}379380# ----------------------------------------------------------------------------- 5. appariement rôle – transactions381382def normaliser_matricule(x) -> str:383 if pd.isna(x): return ""384 return "".join(ch for ch in str(x) if ch.isdigit())385386def appariement_role(b: pd.DataFrame, role: pd.DataFrame, verite: pd.DataFrame | None = None) -> dict:387 """Trois paliers : (1) clé exacte, (2) clé normalisée (chiffres seulement), (3) blocage par secteur + année + terrain ± 2 % avec score388 de concordance (un seul candidat exigé). Taux par palier et, si la vérité est connue (id_vente_origine ↔ matricule canonique), faux appariements."""389 t = b.copy(); r = role.copy()390 r["cle"] = r["matricule"].map(normaliser_matricule)391 t["cle"] = t["matricule"].map(normaliser_matricule)392 t["palier"] = ""; t["matricule_apparie"] = np.nan393 # palier 1 : exact394 exact = t["matricule"].isin(r["matricule"]); t.loc[exact, "palier"] = "1 exact"; t.loc[exact, "matricule_apparie"] = t.loc[exact, "matricule"]395 # palier 2 : normalisé396 m2 = (~exact) & t["cle"].isin(r["cle"]) & (t["cle"] != "")397 cle_to_mat = dict(zip(r["cle"], r["matricule"]))398 t.loc[m2, "palier"] = "2 normalisé"; t.loc[m2, "matricule_apparie"] = t.loc[m2, "cle"].map(cle_to_mat)399 # palier 3 : blocage + concordance (secteur, année, terrain ±2 %) — un seul candidat400 reste = t.index[t["palier"] == ""]401 n3 = 0; ambigus = 0402 for i in reste:403 row = t.loc[i]404 cand = r[(r["secteur"] == row["secteur"]) & (r["role_annee_construction"] == row["annee_construction"])]405 if pd.notna(row["superficie_terrain_m2"]):406 cand = cand[(cand["role_superficie_terrain_m2"] - row["superficie_terrain_m2"]).abs() <= 0.02 * row["superficie_terrain_m2"]]407 else:408 cand = cand.iloc[0:0]409 if len(cand) == 1:410 t.at[i, "palier"] = "3 attributs"; t.at[i, "matricule_apparie"] = cand["matricule"].iloc[0]; n3 += 1411 elif len(cand) > 1:412 ambigus += 1413 t.loc[t["palier"] == "", "palier"] = "non apparié"414 taux = t["palier"].value_counts().reindex(["1 exact", "2 normalisé", "3 attributs", "non apparié"]).fillna(0).astype(int)415 out = {"table": t, "taux": taux, "pct": (100 * taux / len(t)).round(1), "ambigus": ambigus}416 if verite is not None:417 # vérité : matricule canonique de l'unité d'origine418 vrai = dict(zip(role["matricule"], role["matricule"]))419 orig = cas1.charger()["ventes"][["id_vente"]].copy()420 orig["matricule_vrai"] = role["matricule"].values # même ordre que le jeu propre (generer_brut)421 t = t.merge(orig.rename(columns={"id_vente": "id_vente_origine"}), on="id_vente_origine", how="left")422 app = t[t["palier"] != "non apparié"]423 out["faux_appariements"] = int((app["matricule_apparie"] != app["matricule_vrai"]).sum())424 out["faux_par_palier"] = app.assign(faux=app["matricule_apparie"] != app["matricule_vrai"]).groupby("palier")["faux"].sum().to_dict()425 out["table"] = t426 return out427428# ----------------------------------------------------------------------------- 6. imputation comparée et sensibilité429430def comparer_imputations(P: dict) -> pd.DataFrame:431 base = geocodage_repli(P["regles"]["df"], [])432 rows = []433 for meth in ["suppression", "mediane", "regression", "indicateur"]:434 d = imputer(base, None, meth)435 S = scores_robustes(d, None); E = exclusions_regression(S["garde"], None)436 m = _modele(E["df"]); c = m["coef"]437 rows.append({"methode": meth, "n": m["n"], "sigma": m["sigma"], "ln_T": c.loc["ln_T", "coef"], "se_ln_T": c.loc["ln_T", "se_hc1"], "qualite": c.loc["qualite", "coef"],438 "age": c.loc["age", "coef"], "prediction_sujet": float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0])})439 return pd.DataFrame(rows).set_index("methode")440441def sensibilite_seuils(b: pd.DataFrame, verite: pd.DataFrame) -> pd.DataFrame:442 """Grille seuil MAD × seuil t : exclusions, σ, β_ln S, et détection des anomalies « à exclure » (ventes liées non déclarées)."""443 a_detecter = set(verite.loc[verite["anomalie"] == "vente_liee_non_declaree", "id_inscription"])444 propres = set(b["id_inscription"]) - set(verite["id_inscription"])445 rows = []446 for zm in [2.5, 3.0, 3.5, 4.0]:447 for zt in [2.5, 3.0, 3.5]:448 P = pipeline(b, zm, zt)449 excl = set(P["journal"].loc[P["journal"]["etape"].isin(["S1", "E1"]), "id_inscription"])450 m = _modele(P["propre"])451 rows.append({"seuil_mad": zm, "seuil_t": zt, "n_final": m["n"], "exclus_S1": int((P["journal"]["etape"] == "S1").sum()), "exclus_E1": int((P["journal"]["etape"] == "E1").sum()),452 "liees_detectees": len(excl & a_detecter), "faux_positifs": len(excl & propres), "sigma": m["sigma"], "ln_S": m["coef"].loc["ln_S", "coef"],453 "prediction_sujet": float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0])})454 return pd.DataFrame(rows)455456# ----------------------------------------------------------------------------- main457458def main(ecrire: bool = True, verbose: bool = True) -> dict:459 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)460 P0 = lambda *a: print(*a) if verbose else None461 g = ecrire_brut() if ecrire else generer_brut()462 b, verite, role = g["brut"], g["verite"], g["role"]463 P0("== Jeu brut :", b.shape, "| anomalies injectées :", len(verite)); P0(verite["anomalie"].value_counts().to_string())464 P = pipeline(b, ecrire=ecrire)465 P0("== Entonnoir :", P["etapes"])466 P0("== Journal par étape/action :\n", P["journal"].groupby(["etape", "action"]).size().to_string())467 P0("== Itérations E1 :", P["regression"]["iterations"])468 # détection : quelles anomalies le pipeline a-t-il traitées ?469 jv = P["journal"].merge(verite, on="id_inscription", how="left")470 P0("== Étape × anomalie :\n", pd.crosstab(jv["etape"], jv["anomalie"].fillna("(ligne propre)")).to_string())471 C = comparer_modeles(P, b)472 P0("== Modèles :\n", C["tableau"].round(4).to_string())473 P0("== Écarts nettoyé − référence (en erreurs-types de la référence) :\n", C["ecarts"].round(3).to_string())474 S = P["scores"]["df"]475 P0("== Scores robustes : max |z_prix_m2| =", round(S["z_prix_m2"].abs().max(), 2), "; exclus S1 :", len(P["scores"]["exclus"]))476 P0(P["scores"]["exclus"][["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "z_prix_m2", "z_superficie", "z_terrain"]].round(2).to_string(index=False))477 E = P["regression"]["exclus"]478 P0("== Exclus E1 :\n", E[["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "t_ext", "cook", "levier"]].round(3).to_string(index=False))479 G = simulation_geocodage(); P0("== Géocodage (simulation) :\n", G["tableau"].round(3).to_string(), "\nβ_dist référence :", round(G["beta_dist_ref"], 4), "; ventes dans le rayon de 3 km :", G["n_dans_rayon"])480 A = appariement_role(b, role, verite); P0("== Appariement :\n", A["taux"].to_string(), "\n", A["pct"].to_string(), "\nambigus :", A["ambigus"], "; faux appariements :", A.get("faux_appariements"), A.get("faux_par_palier"))481 I = comparer_imputations(P); P0("== Imputations :\n", I.round(4).to_string())482 T = sensibilite_seuils(b, verite); P0("== Sensibilité aux seuils :\n", T.round(4).to_string(index=False))483 return {"brut": b, "verite": verite, "role": role, "pipeline": P, "modeles": C, "geocodage": G, "appariement": A, "imputations": I, "sensibilite": T}484485if __name__ == "__main__":486 main(ecrire="--sans-ecriture" not in sys.argv)487