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 clé d'unité, même promesse, même prix) : inscription la plus récente conservée")168 df = df[~dup]169 # R1b — doublon par attributs (la clé peut être absente ou mal saisie sur l'une des deux inscriptions)170 attr = ["secteur", "type", "superficie_habitable_m2", "superficie_terrain_m2", "annee_construction", "date_promesse", "prix_vente"]171 dup2 = df.duplicated(subset=attr, keep="last")172 _journal(j, df, df.index[dup2], "R1", "exclusion", "doublon d'inscription repéré par les attributs (clé d'unité absente ou discordante) : inscription la plus récente conservée")173 df = df[~dup2]; etapes.append(("R1 doublons", len(df)))174 # R3 — meubles : le prix de l'immeuble est le prix déclaré moins la contrepartie des meubles (correction)175 m = df["montant_meubles"].fillna(0) > 0; ids = df.index[m]176 _journal(j, df, ids, "R3", "correction", "contrepartie des meubles retranchée du prix", "prix_vente", df.loc[ids, "prix_vente"].tolist(),177 (df.loc[ids, "prix_vente"] - df.loc[ids, "montant_meubles"]).tolist())178 df["prix_vente"] = df["prix_vente"] - df["montant_meubles"].fillna(0)179 # R4 — prix symbolique180 m = df["prix_vente"] < SEUIL_PRIX_MIN181 _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())182 df = df[~m]; etapes.append(("R4 prix symbolique", len(df)))183 # R5 — unité de superficie : valeur implausible dont la conversion pi² → m² est plausible184 s = df["superficie_habitable_m2"]; conv = s / PI2_PAR_M2185 m = (s > SEUIL_S_MAX_M2) & conv.between(40, 400); ids = df.index[m]186 _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",187 s[m].tolist(), conv[m].round(0).tolist())188 df.loc[m, "superficie_habitable_m2"] = conv[m].round(0)189 m2 = (df["superficie_habitable_m2"] > SEUIL_S_MAX_M2) # implausible et non convertible → exclusion (aucun cas attendu)190 _journal(j, df, df.index[m2], "R5", "exclusion", "superficie implausible non convertible", "superficie_habitable_m2", df.loc[m2, "superficie_habitable_m2"].tolist())191 df = df[~m2]192 # R6 — coordonnées nulles ou hors emprise → invalidées (géocodage de repli à l'étape G1)193 m = (df["lat"].fillna(0) == 0) | (df["lon"].fillna(0) == 0) | ~df["lat"].between(*BBOX["lat"]) | ~df["lon"].between(*BBOX["lon"])194 ids = df.index[m]195 _journal(j, df, ids, "R6", "correction", "coordonnées nulles ou hors de l'emprise du marché : invalidées, à géocoder", "lat;lon",196 [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])197 df["coord_valide"] = ~m; df.loc[m, ["lat", "lon"]] = np.nan198 # R7 — conditions de vente non représentatives déclarées199 m = df["type_transaction"].fillna("ordinaire").str.lower() != "ordinaire"200 _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())201 df = df[~m]; etapes.append(("R7 non représentatives", len(df)))202 df["date_vente"] = df["date_promesse"]203 df["age"] = 2026 - df["annee_construction"]204 return {"df": df, "journal": j, "etapes": etapes}205206def geocodage_repli(df: pd.DataFrame, j: list) -> pd.DataFrame:207 """G1 : sans adresse réelle à géocoder, les coordonnées invalides reçoivent le centre de leur secteur (précision « secteur »)."""208 df = df.copy(); m = ~df["coord_valide"]; ids = df.index[m]209 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"]]210 df["precision_geocodage"] = np.where(m, "centre de secteur", "point d'origine")211 _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,212 [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])213 df["dist_centre_km"] = _dist_centre(df["lat"], df["lon"]).round(2)214 return df215216def imputer(df: pd.DataFrame, j: list | None = None, methode: str = "mediane") -> pd.DataFrame:217 """I1 : valeurs manquantes. methode ∈ {mediane (conditionnelle secteur × type), regression, indicateur, suppression}."""218 df = df.copy()219 champs = {"superficie_terrain_m2": "terrain", "qualite_1a5": "qualite", "annee_construction": "annee"}220 if methode == "suppression":221 m = df[list(champs)].isna().any(axis=1)222 if j is not None: _journal(j, df, df.index[m], "I1", "exclusion", "valeur manquante (suppression par liste)")223 df = df[~m].copy(); df["age"] = 2026 - df["annee_construction"]; return df224 for c in champs:225 m = df[c].isna()226 if not m.any(): continue227 if methode in ("mediane", "indicateur"):228 med = df.groupby(["secteur", "type"])[c].transform("median")229 val = med[m].fillna(df[c].median())230 elif methode == "regression":231 comp = df[~m]232 def X_of(d):233 X = pd.DataFrame({"const": 1.0, "ln_S": np.log(d["superficie_habitable_m2"]), "etages": (d["type"] == "à étages").astype(float)}, index=d.index)234 for s in cas1.SECTEURS[1:]: X[f"s_{s}"] = (d["secteur"] == s).astype(float)235 if c != "annee_construction": X["annee"] = d["annee_construction"].fillna(d["annee_construction"].median())236 if c != "superficie_terrain_m2": X["ln_T"] = np.log(d["superficie_terrain_m2"].fillna(d["superficie_terrain_m2"].median()))237 return X238 y = np.log(comp[c]) if c == "superficie_terrain_m2" else comp[c]239 r = h.mco(X_of(comp).values, y.values); pred = X_of(df[m]).values @ r["beta"]240 val = pd.Series(np.exp(pred) if c == "superficie_terrain_m2" else pred, index=df.index[m])241 if c == "qualite_1a5": val = val.round().clip(1, 5)242 if c == "annee_construction": val = val.round()243 if j is not None:244 _journal(j, df, df.index[m], "I1", "imputation", f"{champs[c]} manquant : imputation ({methode})", c, None, val.round(1).tolist())245 df.loc[m, c] = val.values246 if methode == "indicateur": df[f"manq_{champs[c]}"] = m.astype(float)247 df["age"] = 2026 - df["annee_construction"]248 return df249250def z_robuste(x: np.ndarray) -> np.ndarray:251 x = np.asarray(x, float); med = np.nanmedian(x); mad = np.nanmedian(np.abs(x - med))252 return (x - med) / (C_MAD * mad)253254def scores_robustes(df: pd.DataFrame, j: list | None = None, seuil: float = SEUIL_MAD) -> dict:255 """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."""256 df = df.copy()257 df["ln_prix_m2"] = np.log(df["prix_vente"] / df["superficie_habitable_m2"])258 df["z_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: z_robuste(s.values))259 df["z_superficie"] = z_robuste(np.log(df["superficie_habitable_m2"].values))260 df["z_terrain"] = df.groupby("secteur")["superficie_terrain_m2"].transform(lambda s: z_robuste(np.log(s.values)))261 df["z_classique_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: (s - s.mean()) / s.std(ddof=1))262 m = (df["z_prix_m2"].abs() > seuil) | (df["z_superficie"].abs() > seuil) | (df["z_terrain"].abs() > seuil)263 ids = df.index[m]264 if j is not None:265 _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,266 [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"])])267 return {"df": df, "exclus": df[m], "garde": df[~m]}268269def diagnostics(X: np.ndarray, y: np.ndarray) -> dict:270 """Levier h_ii, résidus studentisés internes r_i et externes t_i, distance de Cook D_i, DFFITS ; vérifie Σh = k."""271 X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape272 r = h.mco(X, y); e = r["residus"]; s2 = r["sigma2"]273 XtX_inv = np.linalg.pinv(X.T @ X)274 hii = np.einsum("ij,jk,ik->i", X, XtX_inv, X)275 r_int = e / np.sqrt(s2 * (1 - hii))276 s2_i = ((n - k) * s2 - e ** 2 / (1 - hii)) / (n - k - 1)277 t_ext = e / np.sqrt(s2_i * (1 - hii))278 cook = r_int ** 2 / k * hii / (1 - hii)279 dffits = t_ext * np.sqrt(hii / (1 - hii))280 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}281282def cook_direct(X: np.ndarray, y: np.ndarray, i: int) -> float:283 """Distance de Cook par réestimation sans l'observation i (contrôle de l'identité algébrique)."""284 X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape285 r = h.mco(X, y); m = np.ones(n, bool); m[i] = False286 b_i = h.mco(X[m], y[m])["beta"]; d = r["beta"] - b_i287 return float(d @ (X.T @ X) @ d / (k * r["sigma2"]))288289def exclusions_regression(df: pd.DataFrame, j: list | None = None, seuil_t: float = SEUIL_T, max_iter: int = 2) -> dict:290 """E1 : ajuste le modèle de référence, exclut |t_i| > seuil, réitère (au plus max_iter fois)."""291 df = df.copy(); exclus = []; iterations = []292 for it in range(1, max_iter + 1):293 X = cas1.matrice(df); y = np.log(df["prix_vente"].astype(float))294 dg = diagnostics(X.values, y.values)295 df["levier"] = dg["h"]; df["t_ext"] = dg["t_ext"]; df["cook"] = dg["cook"]296 m = np.abs(dg["t_ext"]) > seuil_t297 iterations.append({"iteration": it, "n": len(df), "exclus": int(m.sum()), "sigma": float(np.sqrt(dg["mco"]["sigma2"])), "cook_max": float(dg["cook"].max()),298 "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())})299 if not m.any(): break300 ids = df.index[m]301 if j is not None:302 _journal(j, df, ids, "E1", "exclusion", f"résidu studentisé externe |t| > {seuil_t} (itération {it})", "t_ext;cook;levier", None,303 [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"])])304 exclus.append(df[m]); df = df[~m]305 return {"df": df, "exclus": pd.concat(exclus) if exclus else df.iloc[0:0], "iterations": iterations, "diag": dg}306307def pipeline(b: pd.DataFrame | None = None, seuil_mad: float = SEUIL_MAD, seuil_t: float = SEUIL_T, imputation: str = "mediane", ecrire: bool = False) -> dict:308 """Chaîne complète : règles → géocodage de repli → imputation → scores robustes → régression → jeu nettoyé + journal."""309 if b is None: b = charger_brut()["brut"]310 R = regles(b); j = R["journal"]; etapes = R["etapes"]311 df = geocodage_repli(R["df"], j)312 df = imputer(df, j, imputation); etapes.append(("I1 imputation", len(df)))313 S = scores_robustes(df, j, seuil_mad); df = S["garde"]; etapes.append(("S1 scores robustes", len(df)))314 E = exclusions_regression(df, j, seuil_t); df = E["df"]; etapes.append(("E1 régression", len(df)))315 jour = pd.DataFrame(j, columns=["etape", "id_inscription", "action", "motif", "champ", "valeur_avant", "valeur_apres"])316 jour.insert(0, "version_pipeline", "ch14_nettoyage.py 1.0 (2026-10-31)")317 if ecrire: jour.to_csv(D / "cas1_gatineau_journal_exclusions.csv", index=False)318 return {"propre": df, "journal": jour, "etapes": etapes, "regles": R, "scores": S, "regression": E, "apres_regles": R["df"]}319320# ----------------------------------------------------------------------------- 3. effet sur le modèle321322def _modele(df: pd.DataFrame) -> dict:323 d = df.copy()324 if "date_vente" not in d: d["date_vente"] = pd.to_datetime(d["date_promesse"])325 if "age" not in d or d["age"].isna().any(): d["age"] = 2026 - d["annee_construction"]326 if "dist_centre_km" not in d: d["dist_centre_km"] = _dist_centre(d["lat"].fillna(45.428), d["lon"].fillna(-75.714))327 m = cas1.hedonique(d)328 return m329330def brut_coerce(b: pd.DataFrame) -> pd.DataFrame:331 """Le jeu brut tel qu'un analyste pressé le passerait au modèle : lignes complètes seulement, aucune autre correction."""332 d = b.copy(); d["date_vente"] = pd.to_datetime(d["date_promesse"]); d["age"] = 2026 - d["annee_construction"]333 d["dist_centre_km"] = _dist_centre(d["lat"], d["lon"])334 d = d.dropna(subset=["superficie_terrain_m2", "qualite_1a5", "annee_construction"])335 return d336337def comparer_modeles(P: dict, b: pd.DataFrame) -> dict:338 propre_ref = cas1.charger()["ventes"]339 jeux = {"brut (lignes complètes)": brut_coerce(b), "après règles R1–R7 (+G1, I1)": imputer(geocodage_repli(P["regles"]["df"], []), None, "mediane"),340 "après scores robustes S1": P["scores"]["garde"], "après régression E1 (nettoyé)": P["propre"], "jeu propre de référence": propre_ref}341 rows = []; modeles = {}342 for nom, d in jeux.items():343 m = _modele(d); modeles[nom] = m; c = m["coef"]344 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"],345 "age": c.loc["age", "coef"], "qualite": c.loc["qualite", "coef"], "garage": c.loc["garage", "coef"], "mois_pct": 100 * m["croissance_mensuelle"],346 "sect_Aylmer": c.loc["sect_Aylmer", "coef"], "const": c.loc["const", "coef"]})347 tab = pd.DataFrame(rows).set_index("jeu")348 # prédiction pour le sujet349 s = cas1.charger()["sujet"]; sd = cas1.sujet_df(s)350 tab["prediction_sujet"] = [float(cas1.predire(modeles[n], sd)[0]) for n in tab.index]351 # écart nettoyé vs référence en erreurs-types352 ref = modeles["jeu propre de référence"]["coef"]; net = modeles["après régression E1 (nettoyé)"]["coef"]353 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"]354 return {"tableau": tab, "modeles": modeles, "ecarts": ecarts}355356# ----------------------------------------------------------------------------- 4. géocodage (simulation)357358TECHNIQUES = {"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édagogiques359360def simulation_geocodage(seed: int = GRAINE, k_voisins: int = 8) -> dict:361 """Erreur de position gaussienne isotrope d'écart-type σ (par technique) appliquée aux coordonnées propres (projetées en km) :362 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,363 variation du coefficient de distance au centre du modèle hédonique."""364 rng = np.random.default_rng(seed)365 d = cas1.charger(); v = d["ventes"].copy(); s = d["sujet"]366 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]367 dist_s = np.sqrt(((xy - xs) ** 2).sum(1)); dans_rayon = dist_s <= 3.0368 W0 = h.matrice_poids_knn(xy, k_voisins, standardiser=False)369 m_ref = cas1.hedonique(v); beta_ref = m_ref["coef"].loc["dist_centre", "coef"]370 rows = []371 for nom, sig in TECHNIQUES.items():372 err = rng.normal(0, sig, size=xy.shape); xy2 = xy + err373 e = np.sqrt((err ** 2).sum(1))374 d2 = np.sqrt(((xy2 - xs) ** 2).sum(1)); change = (d2 <= 3.0) != dans_rayon375 W1 = h.matrice_poids_knn(xy2, k_voisins, standardiser=False)376 conserves = (W0 * W1).sum(1) / k_voisins377 v2 = v.copy(); v2["dist_centre_km"] = np.sqrt((xy2 ** 2).sum(1))378 m2 = cas1.hedonique(v2)379 rows.append({"technique": nom, "sigma_km": sig, "erreur_mediane_m": 1000 * np.median(e), "erreur_p95_m": 1000 * np.quantile(e, 0.95),380 "pct_rayon_change": 100 * change.mean(), "pct_voisins_conserves": 100 * conserves.mean(), "beta_dist": m2["coef"].loc["dist_centre", "coef"],381 "sigma_modele": m2["sigma"], "erreurs_m": 1000 * e})382 tab = pd.DataFrame(rows).set_index("technique")383 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())}384385# ----------------------------------------------------------------------------- 5. appariement rôle – transactions386387def normaliser_matricule(x) -> str:388 if pd.isna(x): return ""389 return "".join(ch for ch in str(x) if ch.isdigit())390391def appariement_role(b: pd.DataFrame, role: pd.DataFrame, verite: pd.DataFrame | None = None) -> dict:392 """Trois paliers : (1) clé exacte, (2) clé normalisée (chiffres seulement), (3) blocage par secteur + année + terrain ± 2 % avec score393 de concordance (un seul candidat exigé). Taux par palier et, si la vérité est connue (id_vente_origine ↔ matricule canonique), faux appariements."""394 t = b.copy(); r = role.copy()395 r["cle"] = r["matricule"].map(normaliser_matricule)396 t["cle"] = t["matricule"].map(normaliser_matricule)397 t["palier"] = ""; t["matricule_apparie"] = pd.Series([None] * len(t), index=t.index, dtype=object)398 # palier 1 : exact399 exact = t["matricule"].isin(r["matricule"]); t.loc[exact, "palier"] = "1 exact"; t.loc[exact, "matricule_apparie"] = t.loc[exact, "matricule"]400 # palier 2 : normalisé401 m2 = (~exact) & t["cle"].isin(r["cle"]) & (t["cle"] != "")402 cle_to_mat = dict(zip(r["cle"], r["matricule"]))403 t.loc[m2, "palier"] = "2 normalisé"; t.loc[m2, "matricule_apparie"] = t.loc[m2, "cle"].map(cle_to_mat)404 # palier 3 : blocage + concordance (secteur, année, terrain ±2 %) — un seul candidat405 reste = t.index[t["palier"] == ""]406 n3 = 0; ambigus = 0407 for i in reste:408 row = t.loc[i]409 cand = r[(r["secteur"] == row["secteur"]) & (r["role_annee_construction"] == row["annee_construction"])]410 if pd.notna(row["superficie_terrain_m2"]):411 cand = cand[(cand["role_superficie_terrain_m2"] - row["superficie_terrain_m2"]).abs() <= 0.02 * row["superficie_terrain_m2"]]412 else:413 cand = cand.iloc[0:0]414 if len(cand) == 1:415 t.at[i, "palier"] = "3 attributs"; t.at[i, "matricule_apparie"] = cand["matricule"].iloc[0]; n3 += 1416 elif len(cand) > 1:417 ambigus += 1418 t.loc[t["palier"] == "", "palier"] = "non apparié"419 taux = t["palier"].value_counts().reindex(["1 exact", "2 normalisé", "3 attributs", "non apparié"]).fillna(0).astype(int)420 out = {"table": t, "taux": taux, "pct": (100 * taux / len(t)).round(1), "ambigus": ambigus}421 if verite is not None:422 # vérité : matricule canonique de l'unité d'origine423 vrai = dict(zip(role["matricule"], role["matricule"]))424 orig = cas1.charger()["ventes"][["id_vente"]].copy()425 orig["matricule_vrai"] = role["matricule"].values # même ordre que le jeu propre (generer_brut)426 t = t.merge(orig.rename(columns={"id_vente": "id_vente_origine"}), on="id_vente_origine", how="left")427 app = t[t["palier"] != "non apparié"]428 out["faux_appariements"] = int((app["matricule_apparie"] != app["matricule_vrai"]).sum())429 out["faux_par_palier"] = app.assign(faux=app["matricule_apparie"] != app["matricule_vrai"]).groupby("palier")["faux"].sum().to_dict()430 out["table"] = t431 return out432433# ----------------------------------------------------------------------------- 5b. effet de masquage (illustration)434435def illustration_masquage(b: pd.DataFrame, seuil: float = SEUIL_MAD) -> dict:436 """Sur le jeu après R1-R3 seulement (prix symboliques encore présents) : combien d'observations dépassent le seuil selon437 le score classique (moyenne, écart-type) et selon le score robuste (médiane, MAD) du ln du prix au m² par secteur ?"""438 df = b.copy(); df["date_promesse"] = pd.to_datetime(df["date_promesse"]); df["date_acte"] = pd.to_datetime(df["date_acte"])439 df["prix_net"] = df["prix_vente"] - df["montant_meubles"].fillna(0)440 df["lpm"] = np.log(df["prix_net"].clip(lower=1) / df["superficie_habitable_m2"])441 out = []442 for sct, g in df.groupby("secteur"):443 x = g["lpm"].values; zc = (x - x.mean()) / x.std(ddof=1); zr = z_robuste(x)444 out.append({"secteur": sct, "n": len(x), "moyenne": x.mean(), "ecart_type": x.std(ddof=1), "mediane": np.median(x), "mad_x1.4826": C_MAD * np.median(np.abs(x - np.median(x))),445 "n_classique": int((np.abs(zc) > seuil).sum()), "n_robuste": int((np.abs(zr) > seuil).sum()), "n_prix_1": int((g["prix_net"] < SEUIL_PRIX_MIN).sum()),446 "z_classique_min": zc.min(), "z_robuste_min": zr.min()})447 return {"tableau": pd.DataFrame(out).set_index("secteur"), "df": df}448449# ----------------------------------------------------------------------------- 6. imputation comparée et sensibilité450451def _modele_indicateurs(df: pd.DataFrame) -> dict:452 """Modèle de référence augmenté des indicatrices de valeur manquante (méthode de l'indicateur) ; prédiction du sujet avec indicatrices nulles."""453 X = cas1.matrice(df); extra = [c for c in df.columns if c.startswith("manq_")]454 for c in extra: X[c] = df[c].values455 y = np.log(df["prix_vente"].astype(float)); r = h.mco(X.values, y.values)456 coef = pd.DataFrame({"coef": r["beta"], "se_hc1": r["se_hc1"]}, index=X.columns)457 s = cas1.charger()["sujet"]; Xs = cas1.matrice(cas1.sujet_df(s)); 458 for c in extra: Xs[c] = 0.0459 pred = float(np.exp(Xs[X.columns].values @ r["beta"])[0] * np.mean(np.exp(r["residus"])))460 return {"coef": coef, "sigma": float(np.sqrt(r["sigma2"])), "n": r["n"], "r2": r["r2"], "prediction": pred}461462def comparer_imputations(P: dict) -> pd.DataFrame:463 """Quatre traitements des 45 valeurs manquantes, à pipeline égal par ailleurs (S1 puis E1) : suppression par liste, médiane conditionnelle464 (secteur × type), régression auxiliaire, médiane + indicatrices de manquant dans le modèle."""465 base = geocodage_repli(P["regles"]["df"], [])466 rows = []467 for meth in ["suppression", "mediane", "regression", "indicateur"]:468 d = imputer(base, None, meth)469 S = scores_robustes(d, None); E = exclusions_regression(S["garde"], None)470 if meth == "indicateur":471 m = _modele_indicateurs(E["df"]); c = m["coef"]; pred = m["prediction"]472 ind = {k: c.loc[k, "coef"] for k in c.index if k.startswith("manq_")}473 else:474 m = _modele(E["df"]); c = m["coef"]; pred = float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0]); ind = {}475 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"],476 "age": c.loc["age", "coef"], "prediction_sujet": pred, **{f"coef_{k}": v for k, v in ind.items()}})477 return pd.DataFrame(rows).set_index("methode")478479def sensibilite_seuils(b: pd.DataFrame, verite: pd.DataFrame) -> pd.DataFrame:480 """Grille seuil MAD × seuil t : exclusions, σ, β_ln S, et détection des anomalies « à exclure » (ventes liées non déclarées)."""481 a_detecter = set(verite.loc[verite["anomalie"] == "vente_liee_non_declaree", "id_inscription"])482 propres = set(b["id_inscription"]) - set(verite["id_inscription"])483 rows = []484 for zm in [2.5, 3.0, 3.5, 4.0]:485 for zt in [2.5, 3.0, 3.5]:486 P = pipeline(b, zm, zt)487 excl = set(P["journal"].loc[P["journal"]["etape"].isin(["S1", "E1"]), "id_inscription"])488 m = _modele(P["propre"])489 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()),490 "liees_detectees": len(excl & a_detecter), "faux_positifs": len(excl & propres), "sigma": m["sigma"], "ln_S": m["coef"].loc["ln_S", "coef"],491 "prediction_sujet": float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0])})492 return pd.DataFrame(rows)493494# ----------------------------------------------------------------------------- main495496def main(ecrire: bool = True, verbose: bool = True) -> dict:497 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)498 P0 = lambda *a: print(*a) if verbose else None499 g = ecrire_brut() if ecrire else generer_brut()500 b, verite, role = g["brut"], g["verite"], g["role"]501 P0("== Jeu brut :", b.shape, "| anomalies injectées :", len(verite)); P0(verite["anomalie"].value_counts().to_string())502 P = pipeline(b, ecrire=ecrire)503 P0("== Entonnoir :", P["etapes"])504 P0("== Journal par étape/action :\n", P["journal"].groupby(["etape", "action"]).size().to_string())505 P0("== Itérations E1 :", P["regression"]["iterations"])506 # détection : quelles anomalies le pipeline a-t-il traitées ?507 jv = P["journal"].merge(verite, on="id_inscription", how="left")508 P0("== Étape × anomalie :\n", pd.crosstab(jv["etape"], jv["anomalie"].fillna("(ligne propre)")).to_string())509 C = comparer_modeles(P, b)510 P0("== Modèles :\n", C["tableau"].round(4).to_string())511 P0("== Écarts nettoyé − référence (en erreurs-types de la référence) :\n", C["ecarts"].round(3).to_string())512 S = P["scores"]["df"]513 P0("== Scores robustes : max |z_prix_m2| =", round(S["z_prix_m2"].abs().max(), 2), "; exclus S1 :", len(P["scores"]["exclus"]))514 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))515 E = P["regression"]["exclus"]516 P0("== Exclus E1 :\n", E[["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "t_ext", "cook", "levier"]].round(3).to_string(index=False))517 M = illustration_masquage(b); P0("== Masquage (après R1-R3, prix symboliques présents) :\n", M["tableau"].round(3).to_string())518 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"])519 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"))520 I = comparer_imputations(P); P0("== Imputations :\n", I.round(4).to_string())521 T = sensibilite_seuils(b, verite); P0("== Sensibilité aux seuils :\n", T.round(4).to_string(index=False))522 return {"brut": b, "verite": verite, "role": role, "pipeline": P, "modeles": C, "geocodage": G, "appariement": A, "imputations": I, "sensibilite": T, "masquage": M}523524if __name__ == "__main__":525 main(ecrire="--sans-ecriture" not in sys.argv)526