spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 37 — Apprentissage automatique et explicabilité : exemple intégral sur les 600 ventes du cas 13(Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).45Le modèle hédonique de référence est `traite.cas1.hedonique` (semi-log, MCO, HC1) ; il n'est PAS redéfini ici. Ce script lui6oppose un arbre de régression, une forêt aléatoire (Breiman 2001), un gradient boosting (Friedman 2001 ; scikit-learn) et XGBoost7(Chen et Guestrin 2016), sous trois schémas de validation : plis aléatoires (les MÊMES dix plis que `ch18_hedonique.validation_croisee`,8graine 20261031), blocs spatiaux stratifiés par secteur (k-moyennes sur les coordonnées), et partition temporelle (ventes de 20269prédites à partir de 2023-2025). Il ajoute la validation « un secteur exclu » (extrapolation hors support), un scénario où la loi de10génération est perturbée par une interaction et une rupture déclarées (le boosting y gagne), les valeurs SHAP (TreeSHAP) du boosting,11la dépendance partielle, l'importance par permutation, la prédiction du sujet par chaque modèle et les chiffres des exercices.12Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch37_ml.py"""13from __future__ import annotations14import sys, math, pathlib, warnings15import numpy as np, pandas as pd16sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))17from traite import cas118from traite import hedonique as h1920warnings.filterwarnings("ignore", category=FutureWarning); warnings.filterwarnings("ignore", category=UserWarning)21ROOT = pathlib.Path(__file__).resolve().parents[3]22GRAINE = 2026103123N_ARBRES_FORET = 5002425# ----------------------------------------------------------------------------- paramètres de génération (README synthétique)26_g = (1 + 9.06 / 100) ** 0.5 - 127VRAIS = {"ln_S": 0.55, "ln_T": 0.10, "age": -0.0055, "age2": 0.00003, "qualite": 0.06, "garage": 0.035, "sous_sol": 0.03, "piscine": 0.02,28         "etages": 0.02, "dist_centre": -0.006, "sigma": 0.085, "mois": math.log1p(_g) / 12,29         "sect": {"Hull": 0.0, "Gatineau": math.log(464_900 / 475_000), "Aylmer": math.log(530_000 / 475_000),30                  "Buckingham": math.log(395_000 / 475_000), "Masson-Angers": math.log(395_000 / 475_000) - 0.02}}31# scénario non linéaire déclaré (section 37.6) : prime « patrimoine rénové » et « mur de la mi-vie »32SCENARIO = {"prime_qualite_age": 0.12, "qualite_min": 4, "age_min_prime": 40, "mur": -0.10, "age_mur": (25, 45)}3334# ----------------------------------------------------------------------------- régresseurs35VARIABLES_ML = ["superficie", "terrain", "age", "qualite", "garage", "sous_sol", "piscine", "etages", "dist_centre", "lat", "lon", "mois"] + \36               [f"sect_{s}" for s in cas1.SECTEURS]37ETIQUETTES = {"superficie": "Superficie habitable (m²)", "terrain": "Terrain (m²)", "age": "Âge (années)", "qualite": "Qualité (1-5)",38              "garage": "Garage (places)", "sous_sol": "Sous-sol fini", "piscine": "Piscine", "etages": "À étages", "dist_centre": "Distance au centre (km)",39              "lat": "Latitude", "lon": "Longitude", "mois": "Mois depuis 2024-07", "sect_Hull": "Secteur Hull", "sect_Gatineau": "Secteur Gatineau",40              "sect_Aylmer": "Secteur Aylmer", "sect_Buckingham": "Secteur Buckingham", "sect_Masson-Angers": "Secteur Masson-Angers"}4142def matrice_ml(df: pd.DataFrame, s: dict | None = None) -> pd.DataFrame:43    """Régresseurs bruts pour les méthodes à base d'arbres : aucune transformation (les arbres sont invariants aux transformations44    monotones), secteur en indicatrices complètes, coordonnées brutes, temps en mois depuis la date de référence du jeu."""45    X = pd.DataFrame(index=df.index)46    X["superficie"] = df["superficie_habitable_m2"].astype(float); X["terrain"] = df["superficie_terrain_m2"].astype(float)47    X["age"] = df["age"].astype(float); X["qualite"] = df["qualite_1a5"].astype(float); X["garage"] = df["garage_places"].astype(float)48    X["sous_sol"] = df["sous_sol_fini"].astype(float); X["piscine"] = df["piscine"].astype(float)49    X["etages"] = (df["type"] == "à étages").astype(float); X["dist_centre"] = df["dist_centre_km"].astype(float)50    X["lat"] = df["lat"].astype(float); X["lon"] = df["lon"].astype(float)51    d = df["date_vente"] if "date_vente" in df else pd.Series([pd.Timestamp(s["date_evaluation"])] * len(df), index=df.index)52    X["mois"] = ((pd.to_datetime(d) - cas1.REF_DATE).dt.days / 30.4375).astype(float)53    for sct in cas1.SECTEURS:54        X[f"sect_{sct}"] = (df["secteur"] == sct).astype(float)55    return X[VARIABLES_ML]5657# ----------------------------------------------------------------------------- modèles58class Hedonique:59    """Enveloppe du modèle hédonique de référence (traite.hedonique.mco sur traite.cas1.matrice) avec l'interface fit/predict."""60    def __init__(self, colonnes_exclues=()): self.excl = list(colonnes_exclues)61    def fit(self, X: pd.DataFrame, y):62        self.cols = [c for c in X.columns if c not in self.excl]63        self.r = h.mco(X[self.cols].values, np.asarray(y, float)); self.beta = self.r["beta"]; return self64    def predict(self, X: pd.DataFrame): return X[self.cols].values @ self.beta6566def _foret(**k):67    from sklearn.ensemble import RandomForestRegressor68    return RandomForestRegressor(**{**dict(n_estimators=N_ARBRES_FORET, min_samples_leaf=2, max_features=0.5, random_state=GRAINE, n_jobs=-1), **k})6970def _boosting(**k):71    from sklearn.ensemble import GradientBoostingRegressor72    return GradientBoostingRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, min_samples_leaf=5, random_state=GRAINE), **k})7374def _xgboost(**k):75    import xgboost as xgb76    return xgb.XGBRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, colsample_bytree=0.8, min_child_weight=5,77                                      reg_lambda=1.0, random_state=GRAINE, n_jobs=4, verbosity=0), **k})7879def _arbre(**k):80    from sklearn.tree import DecisionTreeRegressor81    return DecisionTreeRegressor(**{**dict(max_depth=3, min_samples_leaf=5, random_state=GRAINE), **k})8283HYPERPARAMETRES = {"Arbre (profondeur 3)": "profondeur 3, feuilles ≥ 5", "Forêt aléatoire": f"{N_ARBRES_FORET} arbres, feuilles ≥ 2, 50 % des variables par division, bootstrap",84                   "Gradient boosting": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, feuilles ≥ 5",85                   "XGBoost": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, 80 % des colonnes, poids min. 5, λ = 1"}86MODELES = {"Hédonique (référence)": ("hedo", Hedonique), "Arbre (profondeur 3)": ("ml", _arbre), "Forêt aléatoire": ("ml", _foret),87           "Gradient boosting": ("ml", _boosting), "XGBoost": ("ml", _xgboost)}8889# ----------------------------------------------------------------------------- plis90def plis_aleatoires(n: int, nfolds: int = 10, graine: int = GRAINE) -> list[np.ndarray]:91    """Les mêmes dix plis que ch18_hedonique.validation_croisee (même graine, même permutation, même découpage)."""92    rng = np.random.default_rng(graine); idx = rng.permutation(n); return list(np.array_split(idx, nfolds))9394def plis_spatiaux(v: pd.DataFrame, par_secteur: int = 2, graine: int = GRAINE) -> tuple[list[np.ndarray], np.ndarray]:95    """Blocs spatiaux stratifiés : dans chaque secteur, k-moyennes (k = 2) sur les coordonnées en km → 10 blocs contigus.96    Chaque secteur reste représenté à l'estimation (par l'autre moitié), mais le voisinage immédiat d'une vente testée ne l'est jamais."""97    from sklearn.cluster import KMeans98    km = h.coordonnees_km(v["lat"].to_numpy(), v["lon"].to_numpy()); bloc = np.zeros(len(v), int); b = 099    for sct in cas1.SECTEURS:100        ii = np.where(v["secteur"].to_numpy() == sct)[0]101        lab = KMeans(par_secteur, n_init=10, random_state=graine).fit_predict(km[ii]); bloc[ii] = b + lab; b += par_secteur102    return [np.where(bloc == j)[0] for j in range(b)], bloc103104def plis_secteurs(v: pd.DataFrame) -> list[np.ndarray]:105    """Un secteur exclu à la fois (extrapolation hors support spatial)."""106    return [np.where(v["secteur"].to_numpy() == sct)[0] for sct in cas1.SECTEURS]107108def partition_temporelle(v: pd.DataFrame, date: str = "2026-01-01") -> tuple[np.ndarray, np.ndarray]:109    test = np.where(v["date_vente"] >= pd.Timestamp(date))[0]; return np.setdiff1d(np.arange(len(v)), test), test110111# ----------------------------------------------------------------------------- métriques112def metriques(pred_log: np.ndarray, y: np.ndarray, facteur: float = 1.0) -> dict:113    """RMSE du logarithme, R² hors échantillon, biais moyen (log), RMSE en dollars, MdAPE, statistiques de ratio IAAO114    (traite.hedonique.statistiques_ratio). `facteur` : correction de retransformation (Duan) appliquée aux prédictions en niveau."""115    m = ~np.isnan(pred_log); e = pred_log[m] - y[m]; P = np.exp(y[m]); V = np.exp(pred_log[m]) * facteur116    r = h.statistiques_ratio(V, P)117    return {"n": int(m.sum()), "rmse_log": float(np.sqrt(np.mean(e ** 2))), "biais_log": float(e.mean()),118            "r2_hors": float(1 - np.sum(e ** 2) / np.sum((y[m] - y[m].mean()) ** 2)), "rmse_dollars": float(np.sqrt(np.mean((V - P) ** 2))),119            "mdape": float(np.median(np.abs(V / P - 1)) * 100), "COD": r["COD"], "PRD": r["PRD"], "PRB": r["PRB"], "PRB_se": r["PRB_se"],120            "ratio_median": r["ratio_median"]}121122def valider(constructeur, genre: str, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> np.ndarray:123    """Prédictions hors échantillon (chaque pli prédit par un modèle estimé sans lui)."""124    X = Xh if genre == "hedo" else Xm; n = len(y); pred = np.full(n, np.nan); tous = np.arange(n)125    for f in plis:126        tr = np.setdiff1d(tous, f); m = constructeur(); m.fit(X.iloc[tr], y[tr]); pred[f] = np.asarray(m.predict(X.iloc[f]), float)127    return pred128129def facteur_duan(pred_log: np.ndarray, y: np.ndarray) -> float:130    """Facteur de retransformation estimé sur les résidus hors échantillon (honnête pour les modèles qui surajustent en échantillon)."""131    m = ~np.isnan(pred_log); return float(np.mean(np.exp(y[m] - pred_log[m])))132133def comparaison(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:134    """Trois schémas de validation pour les cinq modèles. Retourne les prédictions hors échantillon et les métriques."""135    n = len(y); pa = plis_aleatoires(n); ps, bloc = plis_spatiaux(v); tr_t, te_t = partition_temporelle(v)136    out = {"plis": {"aleatoire": pa, "spatial": ps, "bloc_spatial": bloc, "temporel": (tr_t, te_t)}, "pred": {}, "metriques": {}, "duan": {}}137    for nom, (genre, ctor) in MODELES.items():138        p_a = valider(ctor, genre, Xh, Xm, y, pa); F = facteur_duan(p_a, y)139        p_s = valider(ctor, genre, Xh, Xm, y, ps)140        p_t = valider(ctor, genre, Xh, Xm, y, [te_t])            # un seul pli de test : les ventes de 2026141        out["pred"][nom] = {"aleatoire": p_a, "spatial": p_s, "temporel": p_t}; out["duan"][nom] = F142        out["metriques"][nom] = {"aleatoire": metriques(p_a, y, F), "spatial": metriques(p_s, y, F), "temporel": metriques(p_t, y, F)}143    # biais temporel par mois de 2026 pour le boosting et l'hédonique144    mois = Xm["mois"].to_numpy()145    out["biais_mensuel_2026"] = {nom: pd.Series(out["pred"][nom]["temporel"][te_t] - y[te_t]).groupby(np.floor(mois[te_t])).mean().to_dict()146                                 for nom in ("Hédonique (référence)", "Gradient boosting", "Forêt aléatoire")}147    return out148149def boosting_detendance(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:150    """Remède à l'extrapolation temporelle : boosting estimé sur ln P − β̂_mois·mois (tendance de l'hédonique estimée sur l'entraînement),151    sans la variable mois, puis tendance rajoutée. Variante « feuilles ≥ 1 » pour montrer l'effet de la feuille extrême."""152    tr, te = partition_temporelle(v); n = len(y); cols = list(Xh.columns)153    hd = Hedonique().fit(Xh.iloc[tr], y[tr]); bm = hd.beta[cols.index("mois")]; mois = Xm["mois"].to_numpy()154    Xs = Xm.drop(columns=["mois"])155    g = _boosting().fit(Xs.iloc[tr], y[tr] - bm * mois[tr]); p = np.full(n, np.nan); p[te] = g.predict(Xs.iloc[te]) + bm * mois[te]156    g1 = _boosting(min_samples_leaf=1).fit(Xm.iloc[tr], y[tr]); p1 = np.full(n, np.nan); p1[te] = g1.predict(Xm.iloc[te])157    return {"beta_mois_entrainement": float(bm), "pred_detendance": p, "metriques_detendance": metriques(p, y), "pred_feuille1": p1,158            "metriques_feuille1": metriques(p1, y), "n_entrainement_apres_mois_17": int((mois[tr] > 17).sum())}159160def validation_secteurs(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:161    """Un secteur exclu à la fois. L'hédonique n'a aucune information sur un secteur absent : son indicatrice est retirée et le secteur162    est prédit au niveau du secteur de référence (Hull) ; si Hull est exclu, Gatineau devient la référence. Ce choix est un JUGEMENT163    de l'analyste, explicite ; les arbres font implicitement le leur à partir des coordonnées et de la distance au centre."""164    n = len(y); plis = plis_secteurs(v); out = {"pred": {}, "metriques": {}, "biais_secteur": {}}165    for nom, (genre, ctor) in MODELES.items():166        pred = np.full(n, np.nan)167        for sct, f in zip(cas1.SECTEURS, plis):168            tr = np.setdiff1d(np.arange(n), f)169            if genre == "hedo":170                drop = f"sect_{sct}" if sct != "Hull" else "sect_Gatineau"; m = Hedonique([drop]).fit(Xh.iloc[tr], y[tr]); pred[f] = m.predict(Xh.iloc[f])171            else:172                m = ctor(); m.fit(Xm.iloc[tr], y[tr]); pred[f] = m.predict(Xm.iloc[f])173        out["pred"][nom] = pred; out["metriques"][nom] = metriques(pred, y)174        out["biais_secteur"][nom] = {sct: float(np.mean(pred[f] - y[f])) for sct, f in zip(cas1.SECTEURS, plis)}175    return out176177# ----------------------------------------------------------------------------- scénario non linéaire déclaré178def perturbation_scenario(v: pd.DataFrame) -> np.ndarray:179    """ln P† = ln P + 0,12·1[qualité ≥ 4 et âge ≥ 40] − 0,10·1[25 ≤ âge ≤ 45] (hypothèse pédagogique déclarée)."""180    q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO181    return S["prime_qualite_age"] * ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])) + S["mur"] * ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1]))182183def scenario_nonlineaire(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:184    """Mêmes plis aléatoires ; hédonique de référence (mal spécifié désormais), forêt, boosting, XGBoost, et hédonique respécifié185    avec les deux indicatrices (ce que l'analyste ferait après avoir lu la dépendance partielle de l'âge et les interactions SHAP)."""186    delta = perturbation_scenario(v); y2 = y + delta; pa = plis_aleatoires(len(y))187    q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO188    Xh2 = Xh.copy(); Xh2["prime"] = ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])).astype(float); Xh2["mur"] = ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1])).astype(float)189    out = {"delta": delta, "y": y2, "n_prime": int(Xh2["prime"].sum()), "n_mur": int(Xh2["mur"].sum()), "pred": {}, "metriques": {}, "Xh2": Xh2}190    for nom, (genre, ctor) in MODELES.items():191        if nom == "Arbre (profondeur 3)": continue192        p = valider(ctor, genre, Xh, Xm, y2, pa); out["pred"][nom] = p; out["metriques"][nom] = metriques(p, y2, facteur_duan(p, y2))193    p = valider(Hedonique, "hedo", Xh2, Xm, y2, pa); out["pred"]["Hédonique respécifié"] = p; out["metriques"]["Hédonique respécifié"] = metriques(p, y2, facteur_duan(p, y2))194    r = Hedonique().fit(Xh2, y2); cols = list(Xh2.columns)195    out["coef_respecifie"] = {"prime": float(r.beta[cols.index("prime")]), "mur": float(r.beta[cols.index("mur")]), "se_prime": float(r.r["se_hc1"][cols.index("prime")]),196                              "se_mur": float(r.r["se_hc1"][cols.index("mur")]), "sigma": float(np.sqrt(r.r["sigma2"])), "r2": float(r.r["r2"])}197    out["gbm_scenario"] = _boosting().fit(Xm, y2)198    out["hedo_scenario"] = Hedonique().fit(Xh, y2)199    return out200201# ----------------------------------------------------------------------------- explicabilité202def shap_boosting(gbm, Xm: pd.DataFrame, xs: pd.DataFrame) -> dict:203    """TreeSHAP (Lundberg et al. 2020) du gradient boosting : valeurs SHAP des 600 ventes, du sujet, valeur de base, importance204    moyenne |φ|, interactions SHAP moyennes. Vérifie l'axiome d'efficacité (somme des φ + base = prédiction)."""205    import shap206    ex = shap.TreeExplainer(gbm); sv = np.asarray(ex.shap_values(Xm)); base = float(np.ravel(ex.expected_value)[0])207    svs = np.asarray(ex.shap_values(xs))[0]208    si = np.asarray(ex.shap_interaction_values(Xm)); inter = np.abs(si).mean(0).copy(); np.fill_diagonal(inter, 0.0)209    inter_df = pd.DataFrame(inter, index=Xm.columns, columns=Xm.columns)210    paires = inter_df.where(np.triu(np.ones_like(inter, bool), 1)).stack().sort_values(ascending=False)211    return {"valeurs": sv, "base": base, "sujet": pd.Series(svs, index=Xm.columns), "importance": pd.Series(np.abs(sv).mean(0), index=Xm.columns).sort_values(ascending=False),212            "interactions": paires, "efficacite_max_err": float(np.abs(sv.sum(1) + base - gbm.predict(Xm)).max()),213            "efficacite_sujet_err": float(abs(svs.sum() + base - gbm.predict(xs)[0]))}214215def importance_permutation(Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], repetitions: int = 10, graine: int = GRAINE) -> pd.DataFrame:216    """Importance par permutation (Breiman 2001) mesurée HORS échantillon : pour chaque pli, la variable j est permutée dans les217    ventes testées et l'accroissement de l'erreur quadratique (log) est moyenné. Boosting et forêt."""218    rng = np.random.default_rng(graine); n = len(y); res = {}219    for nom, ctor in (("Gradient boosting", _boosting), ("Forêt aléatoire", _foret)):220        gain = np.zeros(Xm.shape[1]); base_err = 0.0221        for f in plis:222            tr = np.setdiff1d(np.arange(n), f); m = ctor(); m.fit(Xm.iloc[tr].to_numpy(), y[tr]); Xt = Xm.iloc[f].to_numpy(); e0 = np.mean((m.predict(Xt) - y[f]) ** 2); base_err += e0 * len(f) / n223            for j in range(Xm.shape[1]):224                acc = 0.0225                for _ in range(repetitions):226                    Xp = Xt.copy(); Xp[:, j] = rng.permutation(Xp[:, j]); acc += np.mean((m.predict(Xp) - y[f]) ** 2)227                gain[j] += (acc / repetitions - e0) * len(f) / n228        res[nom] = pd.Series(gain, index=Xm.columns); res[nom + " (EQM de base)"] = base_err229    return res230231def dependance_partielle(modele, X: pd.DataFrame, variable: str, grille: np.ndarray) -> np.ndarray:232    """PDP (Friedman 2001) : moyenne des prédictions lorsque la variable est fixée à chaque valeur de la grille, les autres variables233    gardant leurs valeurs observées. Retourne la courbe en logarithme du prix."""234    out = np.empty(len(grille))235    for i, g in enumerate(grille):236        Xg = X.copy(); Xg[variable] = g; out[i] = float(np.mean(modele.predict(Xg)))237    return out238239def pdp_vraie(v: pd.DataFrame, variable: str, grille: np.ndarray, scenario: bool = False) -> np.ndarray:240    """Dépendance partielle de la LOI DE GÉNÉRATION (connue) : utile pour juger les modèles sur données synthétiques."""241    a0 = v["age"].to_numpy(float); S0 = v["superficie_habitable_m2"].to_numpy(float); q = v["qualite_1a5"].to_numpy(float); out = np.empty(len(grille))242    for i, g in enumerate(grille):243        a = np.full_like(a0, g) if variable == "age" else a0; S = np.full_like(S0, g) if variable == "superficie" else S0244        f = VRAIS["ln_S"] * np.log(S / 140.0) + VRAIS["age"] * a + VRAIS["age2"] * a ** 2245        if scenario:246            f = f + SCENARIO["prime_qualite_age"] * ((q >= SCENARIO["qualite_min"]) & (a >= SCENARIO["age_min_prime"])) + SCENARIO["mur"] * ((a >= SCENARIO["age_mur"][0]) & (a <= SCENARIO["age_mur"][1]))247        out[i] = float(np.mean(f))248    return out249250def courbes_apprentissage(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], profondeurs=(1, 2, 4)) -> dict:251    """Erreur quadratique (log) d'entraînement et de validation croisée selon le nombre d'arbres du boosting (staged_predict), pour252    plusieurs profondeurs : surapprentissage et arrêt précoce."""253    n = len(y); out = {}254    for d in profondeurs:255        tr_err = np.zeros(1200); te_err = np.zeros(1200)256        for f in plis:257            tr = np.setdiff1d(np.arange(n), f); m = _boosting(max_depth=d).fit(Xm.iloc[tr], y[tr])258            for k, (pt, pv) in enumerate(zip(m.staged_predict(Xm.iloc[tr]), m.staged_predict(Xm.iloc[f]))):259                tr_err[k] += np.mean((pt - y[tr]) ** 2) / len(plis); te_err[k] += np.mean((pv - y[f]) ** 2) * len(f) / n260        out[d] = {"entrainement": np.sqrt(tr_err), "validation": np.sqrt(te_err), "optimum": int(np.argmin(te_err)) + 1, "rmse_min": float(np.sqrt(te_err.min()))}261    return out262263def sensibilite_hyperparametres(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> pd.DataFrame:264    """Grille d'hyperparamètres du boosting (profondeur, taux d'apprentissage, nombre d'arbres, taille minimale des feuilles) et de la265    forêt (feuilles, variables par division) : RMSE log et COD en validation croisée aléatoire."""266    lignes = []267    grille_b = [dict(max_depth=1, n_estimators=1200), dict(max_depth=2, n_estimators=300), dict(max_depth=2, n_estimators=1200), dict(max_depth=2, n_estimators=3000),268                dict(max_depth=3, n_estimators=1200), dict(max_depth=4, n_estimators=1200), dict(max_depth=2, n_estimators=1200, learning_rate=0.1),269                dict(max_depth=2, n_estimators=1200, min_samples_leaf=1), dict(max_depth=2, n_estimators=1200, subsample=1.0)]270    for k in grille_b:271        p = valider(lambda: _boosting(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))272        lignes.append({"modele": "Gradient boosting", "profondeur": k.get("max_depth", 2), "arbres": k.get("n_estimators", 1200), "nu": k.get("learning_rate", 0.03),273                       "feuille_min": k.get("min_samples_leaf", 5), "sous_ech": k.get("subsample", 0.8), "variables": "toutes", **m})274    for k in [dict(min_samples_leaf=1, max_features=1.0), dict(min_samples_leaf=2, max_features=0.5), dict(min_samples_leaf=5, max_features=0.33), dict(min_samples_leaf=2, max_features=0.5, n_estimators=100)]:275        p = valider(lambda: _foret(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))276        lignes.append({"modele": "Forêt aléatoire", "profondeur": "max", "arbres": k.get("n_estimators", N_ARBRES_FORET), "nu": np.nan, "feuille_min": k["min_samples_leaf"],277                       "sous_ech": "bootstrap", "variables": f"{int(round(k['max_features'] * 100))} %", **m})278    return pd.DataFrame(lignes)279280# ----------------------------------------------------------------------------- sujet281def prediction_sujet(modeles_ajustes: dict, duan: dict, Xh_s: pd.DataFrame, Xm_s: pd.DataFrame, ref: dict, s: dict) -> pd.DataFrame:282    """Prédiction du sujet par chaque modèle estimé sur les 600 ventes : logarithme, niveau avec correction de retransformation283    (facteur hors échantillon ; pour l'hédonique, traite.cas1.predire donne la valeur de référence 637 953 $)."""284    rows = []285    for nom, m in modeles_ajustes.items():286        lp = float(m.predict(Xh_s if nom.startswith("Hédonique") else Xm_s)[0]); F = duan[nom]287        rows.append({"modele": nom, "ln_pred": lp, "naif": math.exp(lp), "facteur": F, "prediction": math.exp(lp) * F})288    df = pd.DataFrame(rows).set_index("modele")289    df.loc["Hédonique (référence)", "prediction_cas1"] = float(cas1.predire(ref, cas1.sujet_df(s))[0])290    return df291292# ----------------------------------------------------------------------------- exercices293def exercices(v: pd.DataFrame, comp: dict, secteurs: dict, scen: dict, sh: dict) -> dict:294    e = {}295    # 37.1 : division d'un nœud sur huit ventes (ln P, superficie) — réduction de la somme des carrés296    sub = v.sort_values("id_vente").iloc[:8]; lnP = np.log(sub["prix_vente"].to_numpy(float)); S = sub["superficie_habitable_m2"].to_numpy(float)297    ordre = np.argsort(S); S, lnP = S[ordre], lnP[ordre]; sse0 = float(np.sum((lnP - lnP.mean()) ** 2)); best = None; essais = []298    for k in range(1, 8):299        g, d = lnP[:k], lnP[k:]; sse = float(np.sum((g - g.mean()) ** 2) + np.sum((d - d.mean()) ** 2)); seuil = float((S[k - 1] + S[k]) / 2)300        essais.append({"seuil": seuil, "sse": sse, "gain": sse0 - sse})301        if best is None or sse < best["sse"]: best = {"seuil": seuil, "sse": sse, "gain": sse0 - sse, "k": k, "c_gauche": float(g.mean()), "c_droite": float(d.mean())}302    e["1"] = {"ids": list(sub.sort_values("superficie_habitable_m2")["id_vente"]), "S": S.tolist(), "lnP": lnP.tolist(), "prix": sub.sort_values("superficie_habitable_m2")["prix_vente"].tolist(),303              "sse0": sse0, "essais": essais, "best": best}304    # 37.2 : variance d'une moyenne d'arbres corrélés305    s2, rho = 0.12 ** 2, 0.30; e["2"] = {"sigma": 0.12, "rho": rho, **{f"B{B}": rho * s2 + (1 - rho) * s2 / B for B in (1, 10, 100, 500)}, "limite": rho * s2, "ecart_type_limite": math.sqrt(rho * s2),306                                        "B_pour_1pct": (1 - rho) / 0.01 / rho if rho else np.inf}307    # 37.3 : boosting à la main308    nu = 0.1; F = 12.80; yv = 13.00; traj = [F]309    for _ in range(3): F = F + nu * (yv - F); traj.append(F)310    e["3"] = {"nu": nu, "F": traj, "residus": [yv - t for t in traj], "m_1pct": math.log(0.01) / math.log(1 - nu), "F_nu1": yv}311    # 37.4 : valeurs de Shapley d'un jeu à trois joueurs (S, Q, L) ; v(S) en milliers de dollars au-dessus de la base312    vS = {frozenset(): 0.0, frozenset("S"): 60.0, frozenset("Q"): 30.0, frozenset("L"): 20.0, frozenset("SQ"): 100.0, frozenset("SL"): 85.0, frozenset("QL"): 45.0, frozenset("SQL"): 130.0}313    import itertools314    joueurs = ["S", "Q", "L"]; phi = {}315    for j in joueurs:316        tot = 0.0317        for perm in itertools.permutations(joueurs):318            avant = frozenset(perm[:perm.index(j)]); tot += vS[avant | {j}] - vS[avant]319        phi[j] = tot / 6320    e["4"] = {"v": {"+".join(sorted(k)) if k else "∅": val for k, val in vS.items()}, "phi": phi, "somme": sum(phi.values()), "total": vS[frozenset("SQL")]}321    # 37.5 : dix prédictions d'un AVM fictif322    prix = np.array([312000, 358500, 401000, 445000, 489500, 532000, 598000, 655000, 742000, 895000], float)323    pred = np.array([335000, 349000, 388000, 468000, 471000, 561000, 575000, 702000, 718000, 839000], float)324    r = h.statistiques_ratio(pred, prix); el = np.log(pred / prix)325    e["5"] = {"prix": prix.tolist(), "pred": pred.tolist(), "ratios": (pred / prix).tolist(), "rmse_log": float(np.sqrt(np.mean(el ** 2))), "mdape": float(np.median(np.abs(pred / prix - 1)) * 100),326              "mape": float(np.mean(np.abs(pred / prix - 1)) * 100), **{k: float(val) for k, val in r.items()}}327    # 37.6 : validation par secteur vs aléatoire (chiffres du script)328    e["6"] = {nom: {"aleatoire": comp["metriques"][nom]["aleatoire"]["COD"], "spatial": comp["metriques"][nom]["spatial"]["COD"], "secteur": secteurs["metriques"][nom]["COD"]}329              for nom in ("Hédonique (référence)", "Forêt aléatoire", "Gradient boosting")}330    e["6"]["biais_secteur"] = secteurs["biais_secteur"]331    # 37.7 : scénario — respécification332    e["7"] = {nom: scen["metriques"][nom]["COD"] for nom in scen["metriques"]}; e["7"]["coef"] = scen["coef_respecifie"]; e["7"]["interactions_top"] = sh["scenario_interactions"].head(3).to_dict()333    return e334335# ----------------------------------------------------------------------------- orchestration336def main(verbose: bool = True, complet: bool = True) -> dict:337    d = cas1.charger(); v = d["ventes"].reset_index(drop=True); s = d["sujet"]338    ref = cas1.hedonique(v); Xh = cas1.matrice(v); Xm = matrice_ml(v); y = np.log(v["prix_vente"].astype(float).to_numpy())339    sd = cas1.sujet_df(s); Xh_s = cas1.matrice(sd)[Xh.columns]; Xm_s = matrice_ml(sd, s)340    comp = comparaison(v, Xh, Xm, y)341    det = boosting_detendance(v, Xh, Xm, y)342    sect = validation_secteurs(v, Xh, Xm, y)343    scen = scenario_nonlineaire(v, Xh, Xm, y)344    ajustes = {nom: (ctor().fit(Xh if genre == "hedo" else Xm, y)) for nom, (genre, ctor) in MODELES.items()}345    sujet = prediction_sujet(ajustes, comp["duan"], Xh_s, Xm_s, ref, s)346    sh = shap_boosting(ajustes["Gradient boosting"], Xm, Xm_s)347    sh_sc = shap_boosting(scen["gbm_scenario"], Xm, Xm_s); sh["scenario_interactions"] = sh_sc["interactions"]; sh["scenario_importance"] = sh_sc["importance"]348    pa = comp["plis"]["aleatoire"]349    perm = importance_permutation(Xh, Xm, y, pa) if complet else None350    grille_S = np.linspace(80, 260, 37); grille_A = np.linspace(1, 95, 48)351    pdp = {"superficie": {"grille": grille_S, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "superficie", grille_S),352                          "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "superficie", grille_S),353                          "Hédonique (référence)": dependance_partielle(ajustes["Hédonique (référence)"], Xh, "ln_S", np.log(grille_S / 140.0)), "vraie": pdp_vraie(v, "superficie", grille_S)},354           "age": {"grille": grille_A, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "age", grille_A),355                   "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "age", grille_A), "vraie": pdp_vraie(v, "age", grille_A)}}356    # hédonique : l'âge entre par deux colonnes (age, age2_100) → PDP par substitution conjointe357    def pdp_hedo_age(modele, Xh_, grille):358        out = np.empty(len(grille))359        for i, g in enumerate(grille):360            Xg = Xh_.copy(); Xg["age"] = g; Xg["age2_100"] = g ** 2 / 100.0; out[i] = float(np.mean(modele.predict(Xg)))361        return out362    pdp["age"]["Hédonique (référence)"] = pdp_hedo_age(ajustes["Hédonique (référence)"], Xh, grille_A)363    pdp["age_scenario"] = {"grille": grille_A, "Gradient boosting": dependance_partielle(scen["gbm_scenario"], Xm, "age", grille_A),364                           "Hédonique (référence)": pdp_hedo_age(scen["hedo_scenario"], Xh, grille_A), "vraie": pdp_vraie(v, "age", grille_A, scenario=True)}365    courbes = courbes_apprentissage(Xm, y, pa) if complet else None366    sens = sensibilite_hyperparametres(Xm, y, pa) if complet else None367    ex = exercices(v, comp, sect, scen, sh)368    out = {"ventes": v, "sujet": s, "modele_ref": ref, "Xh": Xh, "Xm": Xm, "y": y, "Xh_s": Xh_s, "Xm_s": Xm_s, "comparaison": comp, "detendance": det, "secteurs": sect,369           "scenario": scen, "ajustes": ajustes, "sujet_pred": sujet, "shap": sh, "permutation": perm, "pdp": pdp, "courbes": courbes, "sensibilite": sens, "exercices": ex}370    if verbose:371        pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)372        print("== Modèle hédonique de référence : R² %.3f, σ %.4f, n %d, k %d" % (ref["r2"], ref["sigma"], ref["n"], ref["k"]))373        for schema in ("aleatoire", "spatial", "temporel"):374            print(f"== Validation {schema} :")375            for nom, m in comp["metriques"].items():376                mm = m[schema]; print(f"   {nom:24s} n={mm['n']:3d} RMSElog={mm['rmse_log']:.4f} R²hors={mm['r2_hors']:.3f} biais={mm['biais_log']:+.4f} RMSE$={mm['rmse_dollars']:,.0f} MdAPE={mm['mdape']:.2f} COD={mm['COD']:.2f} PRD={mm['PRD']:.3f} PRB={mm['PRB']:.3f} méd={mm['ratio_median']:.3f} Duan={comp['duan'][nom]:.4f}")377        print("   biais mensuel 2026 :", {k: {int(kk): round(vv, 3) for kk, vv in vals.items()} for k, vals in comp["biais_mensuel_2026"].items()})378        md = det["metriques_detendance"]; m1 = det["metriques_feuille1"]379        print(f"== Boosting détendancé (test 2026) : biais={md['biais_log']:+.4f} RMSElog={md['rmse_log']:.4f} COD={md['COD']:.2f} MdAPE={md['mdape']:.2f} ; β_mois entraînement {det['beta_mois_entrainement']:.5f} ; n(mois>17)={det['n_entrainement_apres_mois_17']}")380        print(f"   Boosting feuilles ≥ 1 (test 2026) : biais={m1['biais_log']:+.4f} RMSElog={m1['rmse_log']:.4f} COD={m1['COD']:.2f}")381        print("== Un secteur exclu :")382        for nom, m in sect["metriques"].items():383            print(f"   {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} biais/secteur=", {k: round(vv, 3) for k, vv in sect["biais_secteur"][nom].items()})384        print(f"== Scénario non linéaire (n prime={scen['n_prime']}, n mur={scen['n_mur']}) :")385        for nom, m in scen["metriques"].items():386            print(f"   {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} PRB={m['PRB']:.3f}")387        print("   coefficients respécifiés :", {k: round(vv, 4) for k, vv in scen["coef_respecifie"].items()})388        print("== Sujet :"); print(sujet.round(4).to_string())389        print(f"== SHAP (boosting) : base {sh['base']:.4f} (exp {math.exp(sh['base']):,.0f} $) ; efficacité max err {sh['efficacite_max_err']:.2e}")390        print("   importance moyenne |φ| :", sh["importance"].round(4).to_dict())391        print("   sujet φ :", sh["sujet"].sort_values().round(4).to_dict())392        print("   interactions (base) :", sh["interactions"].head(4).round(4).to_dict()); print("   interactions (scénario) :", sh["scenario_interactions"].head(4).round(4).to_dict())393        if perm is not None:394            print("== Importance par permutation (ΔEQM log, hors échantillon) :")395            for nom in ("Gradient boosting", "Forêt aléatoire"): print(f"   {nom} (EQM base {perm[nom + ' (EQM de base)']:.5f}) :", perm[nom].sort_values(ascending=False).round(5).head(8).to_dict())396        if courbes is not None:397            print("== Courbes d'apprentissage :", {d_: (c["optimum"], round(c["rmse_min"], 4), round(float(c["entrainement"][-1]), 4)) for d_, c in courbes.items()})398        if sens is not None:399            print("== Sensibilité aux hyperparamètres :"); print(sens[["modele", "profondeur", "arbres", "nu", "feuille_min", "sous_ech", "variables", "rmse_log", "mdape", "COD", "PRB"]].round(4).to_string(index=False))400        print("== Exercices :")401        print("   37.1", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["1"]["best"].items()}, "SSE0", round(ex["1"]["sse0"], 4)); print("        ventes", ex["1"]["ids"], ex["1"]["S"], [round(x, 4) for x in ex["1"]["lnP"]])402        print("   37.2", {k: (round(vv, 6) if isinstance(vv, float) else vv) for k, vv in ex["2"].items()})403        print("   37.3", {k: (np.round(vv, 4).tolist() if isinstance(vv, list) else round(vv, 4)) for k, vv in ex["3"].items()})404        print("   37.4", {k: (round(vv, 3) if isinstance(vv, float) else vv) for k, vv in ex["4"].items()})405        print("   37.5", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["5"].items() if k not in ("prix", "pred", "ratios")}, [round(r_, 3) for r_ in ex["5"]["ratios"]])406        print("   37.6", ex["6"]); print("   37.7", ex["7"])407    return out408409if __name__ == "__main__":410    main()411