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 variante_mnar_terrain(seed: int = GRAINE) -> pd.DataFrame:480 """Exercice 14.7 : les 20 terrains manquants ne sont plus tirés au hasard (MCAR) mais sont les 20 plus petits terrains du jeu481 parmi les lignes non altérées autrement (MNAR). Retourne le tableau de `comparer_imputations` pour cette variante (+ e.-t. des indicatrices)."""482 g = generer_brut(seed); b = g["brut"].copy(); ver = g["verite"]483 v = cas1.charger()["ventes"].set_index("id_vente")484 ids = ver.loc[ver["anomalie"] == "terrain_manquant", "id_vente_origine"]485 m = b["id_vente_origine"].isin(ids); b.loc[m, "superficie_terrain_m2"] = v.loc[b.loc[m, "id_vente_origine"], "superficie_terrain_m2"].values486 cand = b[~b["id_inscription"].isin(set(ver["id_inscription"]))].nsmallest(20, "superficie_terrain_m2")487 b.loc[cand.index, "superficie_terrain_m2"] = np.nan488 P = pipeline(b); tab = comparer_imputations(P)489 # erreur-type de l'indicatrice de terrain manquant490 base = geocodage_repli(P["regles"]["df"], []); d = imputer(base, None, "indicateur")491 S = scores_robustes(d, None); E = exclusions_regression(S["garde"], None); mi = _modele_indicateurs(E["df"])492 tab.loc["indicateur", "se_manq_terrain"] = mi["coef"].loc["manq_terrain", "se_hc1"]493 return tab494495def sensibilite_seuils(b: pd.DataFrame, verite: pd.DataFrame) -> pd.DataFrame:496 """Grille seuil MAD × seuil t : exclusions, σ, β_ln S, et détection des anomalies « à exclure » (ventes liées non déclarées)."""497 a_detecter = set(verite.loc[verite["anomalie"] == "vente_liee_non_declaree", "id_inscription"])498 propres = set(b["id_inscription"]) - set(verite["id_inscription"])499 rows = []500 for zm in [2.5, 3.0, 3.5, 4.0]:501 for zt in [2.5, 3.0, 3.5]:502 P = pipeline(b, zm, zt)503 excl = set(P["journal"].loc[P["journal"]["etape"].isin(["S1", "E1"]), "id_inscription"])504 m = _modele(P["propre"])505 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()),506 "liees_detectees": len(excl & a_detecter), "faux_positifs": len(excl & propres), "sigma": m["sigma"], "ln_S": m["coef"].loc["ln_S", "coef"],507 "prediction_sujet": float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0])})508 return pd.DataFrame(rows)509510# ----------------------------------------------------------------------------- main511512def main(ecrire: bool = True, verbose: bool = True) -> dict:513 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)514 P0 = lambda *a: print(*a) if verbose else None515 g = ecrire_brut() if ecrire else generer_brut()516 b, verite, role = g["brut"], g["verite"], g["role"]517 P0("== Jeu brut :", b.shape, "| anomalies injectées :", len(verite)); P0(verite["anomalie"].value_counts().to_string())518 P = pipeline(b, ecrire=ecrire)519 P0("== Entonnoir :", P["etapes"])520 P0("== Journal par étape/action :\n", P["journal"].groupby(["etape", "action"]).size().to_string())521 P0("== Itérations E1 :", P["regression"]["iterations"])522 # détection : quelles anomalies le pipeline a-t-il traitées ?523 jv = P["journal"].merge(verite, on="id_inscription", how="left")524 P0("== Étape × anomalie :\n", pd.crosstab(jv["etape"], jv["anomalie"].fillna("(ligne propre)")).to_string())525 C = comparer_modeles(P, b)526 P0("== Modèles :\n", C["tableau"].round(4).to_string())527 P0("== Écarts nettoyé − référence (en erreurs-types de la référence) :\n", C["ecarts"].round(3).to_string())528 S = P["scores"]["df"]529 P0("== Scores robustes : max |z_prix_m2| =", round(S["z_prix_m2"].abs().max(), 2), "; exclus S1 :", len(P["scores"]["exclus"]))530 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))531 E = P["regression"]["exclus"]532 P0("== Exclus E1 :\n", E[["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "t_ext", "cook", "levier"]].round(3).to_string(index=False))533 M = illustration_masquage(b); P0("== Masquage (après R1-R3, prix symboliques présents) :\n", M["tableau"].round(3).to_string())534 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"])535 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"))536 I = comparer_imputations(P); P0("== Imputations :\n", I.round(4).to_string())537 T = sensibilite_seuils(b, verite); P0("== Sensibilité aux seuils :\n", T.round(4).to_string(index=False))538 V = variante_mnar_terrain(); P0("== Variante MNAR (ex. 14.7) :\n", V.round(4).to_string())539 return {"brut": b, "verite": verite, "role": role, "pipeline": P, "modeles": C, "geocodage": G, "appariement": A, "imputations": I, "sensibilite": T, "masquage": M, "mnar": V}540541if __name__ == "__main__":542 main(ecrire="--sans-ecriture" not in sys.argv)543