Ré-estimation hédonique sur le nœud : terrain/autre + recalibrage P10/P90

Commit

scripts/hedonic-retrain.py (sklearn HGB ≈ LightGBM, log-prix, encodages
locaux par voisinage, pondération de récence) — bench équitable sur les
ventes jamais vues par le build laptop (> 2026-07-27, n=4849) :

- unifamilial : laptop meilleur (10,3 % vs 13 % MdAPE) → conservé ;
- terrain 62→36 %, autre 79→52 % → challenger appliqué sur ces 926k unités ;
- P10/P90 laptop recalibrés (couverture 18,5/16,8 % → 10,1/9,8 %) ;
- clamp p10 ≤ est ≤ p90 (38 826 violations du build laptop → 0).

Sauvegarde intégrale dans units_est_prev (rollback documenté §10).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Simon-Pierre Boucher committed 2 mo ago · Aug 23, 2026 parent cf84e27

2 changed files +404 −1 Browse files

modified docs/PIPELINE-DONNEES.md +44 −1

+   sessions neuves ; ~1 essai suffit en général).
+ - N'alimente QUE `transactions`. Les estimations hédoniques (`est_*`, `p10/p90`)
+   et `market_index` d'une vente toute neuve restent celles du dernier build
++  complet tant que le pipeline batch n'a pas été rejoué (voir toutefois §10 :
++  ré-estimation partielle possible sur le nœud).
+ - Historique de premier comblement : gap 2026-07-27 → 2026-08-12 rempli
+   province-wide le 2026-08-22 (+4189 ventes, appariement units 100 %).
++
++## 10. Ré-estimation hédonique sur le nœud (`scripts/hedonic-retrain.py`) — 2026-08-22
++
++Grâce aux ventes fraîches du connecteur (§9), on dispose d'un jeu de test
++« jamais vu » par le modèle laptop (ventes > 2026-07-27, build du 2026-08-08).
++`scripts/hedonic-retrain.py` (sklearn HistGradientBoostingRegressor ≈ LightGBM,
++cible log-prix, features rôle + encodages locaux par unité de voisinage +
++tendance temporelle, pondération de récence) permet d'évaluer et de corriger le
++modèle en place sans le pipeline laptop :
++
++- `--eval` : bench équitable laptop vs challenger sur les ventes jamais vues
++  (n=4849, aucune écriture).
++- `--apply` : ré-entraîne sur tout l'historique et applique la stratégie
++  composite validée, avec sauvegarde intégrale dans `units_est_prev`.
++
++### Verdict du bench 2026-08-22 (MdAPE sur ventes jamais vues)
++| Segment | laptop | challenger nœud |
++|---|---|---|
++| unifamilial (n=4420) | **10,3 %** | 12,8-13,0 % |
++| terrain (n=139) | 62,2 % | **36-37 %** |
++| autre (n=150) | 79,3 % | **52-54 %** |
++
++→ Le modèle laptop reste MEILLEUR sur l'unifamilial (conservé) ; le challenger
++est appliqué UNIQUEMENT sur `terrain` + `autre` (926 432 unités, 25 % du parc).
++
++### Autres corrections appliquées le 2026-08-22
++- **Recalibrage P10/P90** du modèle laptop (couverture observée 18,5 %/16,8 %
++  hors bornes vs cible 10/10) : élargissement `p' = est×(p/est)^α`,
++  α_lo=1,754 / α_hi=1,255 → couverture après : **10,1 % / 9,8 %**. Les
++  fourchettes affichées sont plus larges mais honnêtes.
++- **Clamp p10 ≤ est ≤ p90** partout (38 826 violations dans le build laptop,
++  0 après).
++
++### Rollback
++`units_est_prev` contient les (est_2026, p10, p90) d'avant intervention :
++```sql
++UPDATE units SET est_2026=(SELECT est_2026 FROM units_est_prev p WHERE p.id_provinc=units.id_provinc),
++  p10=(SELECT p10 FROM units_est_prev p WHERE p.id_provinc=units.id_provinc),
++  p90=(SELECT p90 FROM units_est_prev p WHERE p.id_provinc=units.id_provinc);
++```
++NB : le prochain build complet laptop écrasera tout — reporter au besoin ces
++améliorations dans `hedonic.py` (recalibrage quantiles + segments terrain/autre).

added scripts/hedonic-retrain.py +360 −0

++#!/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()

spb vrai-prix

Ré-estimation hédonique sur le nœud : terrain/autre + recalibrage P10/P90

Theme