#!/usr/bin/env python3 # Ré-estimation du modèle hédonique DIRECTEMENT SUR LE NŒUD. # # Le pipeline batch historique (laptop, hedonic.py/LightGBM) produit est_2021.. # est_2026 + p10/p90 par unité. Depuis que le connecteur ingest-jdm.mjs ajoute # les ventes récentes en continu dans vraiprix.db, on peut ré-entraîner ici même # (sklearn HistGradientBoostingRegressor ≈ LightGBM ; 28 cœurs / 96 Go). # # Deux modes : # --eval Bench ÉQUITABLE : le modèle laptop a vu les ventes ≤ 2026-07-27 # (build 2026-08-08). On entraîne le challenger sur ce même # périmètre, puis on compare les deux sur les ventes > 2026-07-27 # (jamais vues par aucun des deux). Rien n'est écrit. # --apply Ré-entraîne sur TOUTES les ventes, prédit les 3,7 M d'unités # (médiane + P10/P90 par modèles quantiles), sauvegarde les anciennes # colonnes dans `units_est_prev`, puis met à jour units.est_2026/p10/p90. # # Cible : log(prix). Enrichissements vs les features brutes du rôle : # - te_ratio : médiane locale (code_mun × unite_voisinage, replis municipalité # puis global) de log(prix/valeur_role) — « ratio d'évaluation » local ; # - ppm2_loc : médiane locale du prix/m² ; # - t : date de vente en mois depuis 2021-01 (tendance de marché) ; # - pondération de récence (demi-vie 36 mois) pour coller au marché actuel. # # Usage : python3 scripts/hedonic-retrain.py --eval # python3 scripts/hedonic-retrain.py --apply import argparse import sqlite3 import sys import time import numpy as np import pandas as pd from sklearn.ensemble import HistGradientBoostingRegressor DB = "data/vraiprix.db" BASELINE_CUTOFF = "2026-07-27" # dernière vente vue par le build laptop 2026-08-08 T0 = pd.Timestamp("2021-01-01") AMOUNT_MIN, AMOUNT_MAX = 50_000, 10_000_000 NUM_FEATS = [ "lat", "lng", "annee_construction", "aire_etages_m2", "superficie_terrain_m2", "front_terrain_m", "nb_etages", "nb_logements", "nb_locaux_non_resid", "nb_chambres_locatives", "n_adresses", "valeur_terrain", "valeur_batiment", "valeur_role", "valeur_anterieure", ] CAT_FEATS = ["type_prop", "lien_physique", "genre_construction"] # naive_role / naive_ppm2 : « estimés naïfs » locaux donnés explicitement au GBM # (les arbres capturent mal les relations multiplicatives rôle × ratio local). ENG_FEATS = ["t", "log_role", "te_ratio", "ppm2_loc", "naive_role", "naive_ppm2", "age", "log_aire", "log_terr"] def log(msg): print(f"[hedonic] {msg}", flush=True) def load_sales(con): q = f""" SELECT t.date, t.amount, u.id_provinc, u.code_mun, u.unite_voisinage, u.municipalite, u.lat, u.lng, u.annee_construction, u.aire_etages_m2, u.superficie_terrain_m2, u.front_terrain_m, u.nb_etages, u.nb_logements, u.nb_locaux_non_resid, u.nb_chambres_locatives, u.n_adresses, u.valeur_terrain, u.valeur_batiment, u.valeur_role, u.valeur_anterieure, u.type_prop, u.lien_physique, u.genre_construction FROM transactions t JOIN units u ON u.id_provinc = t.id_provinc WHERE t.amount BETWEEN {AMOUNT_MIN} AND {AMOUNT_MAX} """ df = pd.read_sql_query(q, con) df["t"] = (pd.to_datetime(df["date"]) - T0).dt.days / 30.44 return df # --- encodages locaux (calculés sur l'ENTRAÎNEMENT seulement, replis lissés) --- class LocalEncoder: K = 8 # lissage : poids du repli def fit(self, df): d = df[(df["valeur_role"] > 0)].copy() d["ratio"] = np.log(d["amount"] / d["valeur_role"]) self.g_ratio = d["ratio"].median() m = d.groupby("municipalite")["ratio"].agg(["median", "size"]) self.mun_ratio = ((m["median"] * m["size"] + self.g_ratio * self.K) / (m["size"] + self.K)).to_dict() d["vk"] = d["code_mun"].astype(str) + "|" + d["unite_voisinage"].astype(str) v = d.groupby(["vk", "municipalite"])["ratio"].agg(["median", "size"]).reset_index() v["fb"] = v["municipalite"].map(self.mun_ratio).fillna(self.g_ratio) v["val"] = (v["median"] * v["size"] + v["fb"] * self.K) / (v["size"] + self.K) self.vois_ratio = dict(zip(v["vk"], v["val"])) p = df[(df["aire_etages_m2"] > 20)].copy() p["ppm2"] = p["amount"] / p["aire_etages_m2"] self.g_ppm2 = p["ppm2"].median() m2 = p.groupby("municipalite")["ppm2"].agg(["median", "size"]) self.mun_ppm2 = ((m2["median"] * m2["size"] + self.g_ppm2 * self.K) / (m2["size"] + self.K)).to_dict() p["vk"] = p["code_mun"].astype(str) + "|" + p["unite_voisinage"].astype(str) v2 = p.groupby(["vk", "municipalite"])["ppm2"].agg(["median", "size"]).reset_index() v2["fb"] = v2["municipalite"].map(self.mun_ppm2).fillna(self.g_ppm2) v2["val"] = (v2["median"] * v2["size"] + v2["fb"] * self.K) / (v2["size"] + self.K) self.vois_ppm2 = dict(zip(v2["vk"], v2["val"])) return self def transform(self, df): vk = df["code_mun"].astype(str) + "|" + df["unite_voisinage"].astype(str) mun_r = df["municipalite"].map(self.mun_ratio).fillna(self.g_ratio) df["te_ratio"] = vk.map(self.vois_ratio).fillna(mun_r) mun_p = df["municipalite"].map(self.mun_ppm2).fillna(self.g_ppm2) df["ppm2_loc"] = vk.map(self.vois_ppm2).fillna(mun_p) return df def featurize(df, enc, cat_maps=None): df = enc.transform(df.copy()) df["log_role"] = np.log1p(df["valeur_role"].clip(lower=0)) df["naive_role"] = np.log1p((df["valeur_role"].clip(lower=0) * np.exp(df["te_ratio"])).clip(lower=0)) df["naive_ppm2"] = np.log1p((df["aire_etages_m2"].clip(lower=0) * df["ppm2_loc"]).clip(lower=0)) df["age"] = (2021 + df["t"] / 12.0) - df["annee_construction"] df["log_aire"] = np.log1p(df["aire_etages_m2"].clip(lower=0)) df["log_terr"] = np.log1p(df["superficie_terrain_m2"].clip(lower=0)) if cat_maps is None: cat_maps = {c: {v: i for i, v in enumerate(df[c].astype(str).fillna("NA").unique())} for c in CAT_FEATS} for c in CAT_FEATS: df[c] = df[c].astype(str).fillna("NA").map(cat_maps[c]).fillna(-1).astype(int) cols = NUM_FEATS + ENG_FEATS + CAT_FEATS X = df[cols].astype(float).to_numpy() return X, cols, cat_maps def make_model(loss="absolute_error", quantile=None, cols=None, light=False): # Config volontairement SIMPLE (≈ 8 min d'entraînement) : la variante lourde # (2500 it / 255 feuilles, ~40 min) ne gagnait que ~1-2 pts de MdAPE sur les # segments remplacés — pas rentable. cat_mask = [c in CAT_FEATS for c in cols] kw = dict( max_iter=600 if light else 900, learning_rate=0.06, max_leaf_nodes=127, min_samples_leaf=40, l2_regularization=0.1, max_bins=255, early_stopping=True, validation_fraction=0.05, n_iter_no_change=40, random_state=42, categorical_features=cat_mask, ) if quantile is not None: return HistGradientBoostingRegressor(loss="quantile", quantile=quantile, **kw) return HistGradientBoostingRegressor(loss=loss, **kw) def recency_weights(t, t_max, half_life=24.0): return np.power(0.5, (t_max - t) / half_life) def metrics(y_true, y_pred, label): ape = np.abs(y_pred - y_true) / y_true return { "label": label, "n": len(y_true), "MdAPE": float(np.median(ape) * 100), "MAPE": float(np.mean(ape) * 100), "±10%": float(np.mean(ape <= 0.10) * 100), "±20%": float(np.mean(ape <= 0.20) * 100), } def print_metrics(rows): hdr = f"{'modèle':<26}{'n':>7}{'MdAPE':>8}{'MAPE':>8}{'±10%':>7}{'±20%':>7}" print(hdr); print("-" * len(hdr)) for r in rows: print(f"{r['label']:<26}{r['n']:>7}{r['MdAPE']:>7.2f}%{r['MAPE']:>7.2f}%" f"{r['±10%']:>6.1f}%{r['±20%']:>6.1f}%") def run_eval(con): df = load_sales(con) log(f"ventes jointes exploitables : {len(df)}") train = df[df["date"] <= BASELINE_CUTOFF].copy() test = df[df["date"] > BASELINE_CUTOFF].copy() log(f"train (≤ {BASELINE_CUTOFF}) : {len(train)} | test (>) : {len(test)}") enc = LocalEncoder().fit(train) Xtr, cols, cat_maps = featurize(train, enc) ytr = np.log(train["amount"].to_numpy()) w = recency_weights(train["t"].to_numpy(), train["t"].max()) t0 = time.time() model = make_model(cols=cols) model.fit(Xtr, ytr, sample_weight=w) log(f"challenger entraîné ({model.n_iter_} itérations, {time.time()-t0:.0f}s)") Xte, _, _ = featurize(test, enc, cat_maps) pred = np.exp(model.predict(Xte)) y = test["amount"].to_numpy() # baseline : est_2026 de l'unité jointe base = pd.read_sql_query( "SELECT id_provinc, est_2026, p10, p90 FROM units", con) test = test.merge(base, on="id_provinc", how="left") mask = test["est_2026"].notna().to_numpy() # mélange géométrique laptop × challenger (moyenne sur l'échelle log) base_est = test["est_2026"].to_numpy() blend = np.where(mask, np.exp(0.5 * np.log(np.where(mask, base_est, 1)) + 0.5 * np.log(pred)), pred) rows = [ metrics(y[mask], base_est[mask], "laptop est_2026"), metrics(y[mask], pred[mask], "challenger (nœud)"), metrics(y[mask], blend[mask], "blend 50/50 (géo)"), ] print(); print(f"=== Ventes jamais vues (> {BASELINE_CUTOFF}) — {mask.sum()} obs ===") print_metrics(rows) # par type print("\n--- par type de propriété (MdAPE %) ---") tt = test[mask].copy(); tt["pred"] = pred[mask]; tt["y"] = y[mask] tt["blend"] = blend[mask] for tp, g in tt.groupby(test[mask]["type_prop"]): b = np.median(np.abs(g["est_2026"] - g["y"]) / g["y"]) * 100 c = np.median(np.abs(g["pred"] - g["y"]) / g["y"]) * 100 bl = np.median(np.abs(g["blend"] - g["y"]) / g["y"]) * 100 print(f" {tp:<16} n={len(g):>5} laptop {b:6.2f}% challenger {c:6.2f}% blend {bl:6.2f}%") # couverture des quantiles baseline cov10 = float(np.mean(y[mask] < test["p10"].to_numpy()[mask]) * 100) cov90 = float(np.mean(y[mask] > test["p90"].to_numpy()[mask]) * 100) print(f"\ncouverture P10/P90 laptop sur test : {cov10:.1f}% sous P10 (cible 10) | " f"{cov90:.1f}% au-dessus de P90 (cible 10)") return rows # Types dont l'estimation laptop est remplacée par le challenger (gain validé # sur ventes jamais vues : terrain MdAPE 62→36 %, autre 79→52 %). REPLACE_TYPES = ("terrain", "autre") def quantile_calibration(con, df): """Facteurs d'élargissement des P10/P90 laptop, calibrés sur les ventes jamais vues (> BASELINE_CUTOFF), hors types remplacés. Couverture observée 18,5 %/16,8 % hors bornes → cible 10 %/10 % de chaque côté.""" test = df[df["date"] > BASELINE_CUTOFF].copy() base = pd.read_sql_query( "SELECT id_provinc, est_2026, p10, p90, type_prop tp FROM units", con) test = test.merge(base, on="id_provinc", how="inner") test = test[test["est_2026"].notna() & (test["est_2026"] > 0) & (test["p10"] > 0) & (test["p90"] > 0) & ~test["tp"].isin(REPLACE_TYPES)] z = np.log(test["amount"] / test["est_2026"]) r_lo = np.log(test["p10"] / test["est_2026"]) # < 0 r_hi = np.log(test["p90"] / test["est_2026"]) # > 0 ok = (r_lo < -1e-6) & (r_hi > 1e-6) s_lo = (z[ok] / r_lo[ok]) # >1 ⇒ sous P10 s_hi = (z[ok] / r_hi[ok]) # >1 ⇒ au-dessus de P90 a_lo = float(np.quantile(s_lo, 0.90)) # P(s_lo > a_lo) = 10 % a_hi = float(np.quantile(s_hi, 0.90)) a_lo, a_hi = max(1.0, a_lo), max(1.0, a_hi) log(f"calibration quantiles laptop (n={ok.sum()}) : alpha_lo={a_lo:.3f}, " f"alpha_hi={a_hi:.3f}") return a_lo, a_hi def run_apply(con): df = load_sales(con) a_lo, a_hi = quantile_calibration(con, df) log(f"ré-entraînement final sur {len(df)} ventes (tout l'historique)") enc = LocalEncoder().fit(df) X, cols, cat_maps = featurize(df, enc) yl = np.log(df["amount"].to_numpy()) w = recency_weights(df["t"].to_numpy(), df["t"].max()) t0 = time.time() med = make_model(cols=cols); med.fit(X, yl, sample_weight=w) log(f"modèle médian : {med.n_iter_} it, {time.time()-t0:.0f}s") t0 = time.time() q10 = make_model(cols=cols, quantile=0.10, light=True); q10.fit(X, yl, sample_weight=w) q90 = make_model(cols=cols, quantile=0.90, light=True); q90.fit(X, yl, sample_weight=w) log(f"modèles quantiles P10/P90 : {time.time()-t0:.0f}s") # --- unités des types remplacés, prédites à la date « maintenant » --- t_now = float(df["t"].max()) ph = ",".join("?" * len(REPLACE_TYPES)) units = pd.read_sql_query(f""" SELECT id_provinc, code_mun, unite_voisinage, municipalite, {', '.join(NUM_FEATS)}, {', '.join(CAT_FEATS)} FROM units WHERE type_prop IN ({ph})""", con, params=REPLACE_TYPES) log(f"unités à ré-estimer ({' + '.join(REPLACE_TYPES)}) : {len(units)}") units["t"] = t_now Xu, _, _ = featurize(units, enc, cat_maps) log("prédiction challenger (médiane + P10/P90)…") est = np.exp(med.predict(Xu)) lo = np.exp(q10.predict(Xu)) hi = np.exp(q90.predict(Xu)) lo2 = np.minimum.reduce([lo, est, hi]); hi2 = np.maximum.reduce([lo, est, hi]) out = pd.DataFrame({ "id_provinc": units["id_provinc"], "est": np.clip(est, 1000, None).round(-2), "p10": np.clip(lo2, 1000, None).round(-2), "p90": np.clip(hi2, 1000, None).round(-2), }) log("écriture en base (sauvegarde units_est_prev puis UPDATE)…") cur = con.cursor() cur.execute("PRAGMA busy_timeout=30000") cur.execute("DROP TABLE IF EXISTS units_est_prev") cur.execute("""CREATE TABLE units_est_prev AS SELECT id_provinc, est_2026, p10, p90 FROM units""") con.commit() # 1) types remplacés → challenger cur.execute("DROP TABLE IF EXISTS units_est_new") cur.execute("""CREATE TABLE units_est_new (id_provinc TEXT PRIMARY KEY, est REAL, p10 REAL, p90 REAL)""") cur.executemany("INSERT OR REPLACE INTO units_est_new VALUES (?,?,?,?)", out.itertuples(index=False, name=None)) con.commit() cur.execute("""UPDATE units SET est_2026 = (SELECT est FROM units_est_new n WHERE n.id_provinc = units.id_provinc), p10 = (SELECT p10 FROM units_est_new n WHERE n.id_provinc = units.id_provinc), p90 = (SELECT p90 FROM units_est_new n WHERE n.id_provinc = units.id_provinc) WHERE id_provinc IN (SELECT id_provinc FROM units_est_new)""") con.commit() cur.execute("DROP TABLE units_est_new") con.commit() # 2) autres types → P10/P90 laptop élargis (p' = est × (p/est)^alpha) ph2 = ",".join("?" * len(REPLACE_TYPES)) cur.execute(f"""UPDATE units SET p10 = ROUND(est_2026 * POW(p10 / est_2026, ?), -2), p90 = ROUND(est_2026 * POW(p90 / est_2026, ?), -2) WHERE type_prop NOT IN ({ph2}) AND est_2026 > 0 AND p10 > 0 AND p90 > 0""", (a_lo, a_hi, *REPLACE_TYPES)) con.commit() cur.execute("PRAGMA wal_checkpoint(TRUNCATE)") n = cur.execute( f"SELECT COUNT(*) FROM units WHERE type_prop IN ({ph2}) AND est_2026 IS NOT NULL", REPLACE_TYPES).fetchone()[0] log(f"terminé : {n} unités ré-estimées (challenger), P10/P90 recalibrés ailleurs " f"(alpha {a_lo:.3f}/{a_hi:.3f}). Ancienne version dans units_est_prev.") if __name__ == "__main__": ap = argparse.ArgumentParser() ap.add_argument("--eval", action="store_true") ap.add_argument("--apply", action="store_true") ap.add_argument("--db", default=DB) a = ap.parse_args() con = sqlite3.connect(a.db) try: if a.eval: run_eval(con) elif a.apply: run_apply(con) else: ap.print_help(); sys.exit(1) finally: con.close()