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 (c) : hédonique où les indicatrices de secteur sont remplacées par un polynôme de degré 2 des coordonnées (km)328    Xh = cas1.matrice(v); km = h.coordonnees_km(v["lat"].to_numpy(), v["lon"].to_numpy()); y = np.log(v["prix_vente"].astype(float).to_numpy())329    Xp = Xh.drop(columns=[f"sect_{s_}" for s_ in cas1.SECTEURS[1:]]).copy()330    Xp["x_km"], Xp["y_km"], Xp["x2"], Xp["y2"], Xp["xy"] = km[:, 0], km[:, 1], km[:, 0] ** 2 / 100, km[:, 1] ** 2 / 100, km[:, 0] * km[:, 1] / 100331    n = len(y); pa = plis_aleatoires(n); p_al = valider(Hedonique, "hedo", Xp, None, y, pa)332    p_sec = np.full(n, np.nan)333    for sct, f in zip(cas1.SECTEURS, plis_secteurs(v)):334        tr = np.setdiff1d(np.arange(n), f); m = Hedonique().fit(Xp.iloc[tr], y[tr]); p_sec[f] = m.predict(Xp.iloc[f])335    e["6c"] = {"aleatoire": metriques(p_al, y, facteur_duan(p_al, y)), "secteur": metriques(p_sec, y),336               "biais_secteur": {sct: float(np.mean(p_sec[f] - y[f])) for sct, f in zip(cas1.SECTEURS, plis_secteurs(v))}}337    # 37.7 (c) : respécification avec des seuils décalés ; (d) : validation imbriquée pour la profondeur du boosting338    q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); y2 = scen["y"]; Xm = matrice_ml(v)339    e["7c"] = {}340    for nom_s, (lo, hi) in {"22-48": (22, 48), "28-42": (28, 42), "25-45 (vrai)": (25, 45)}.items():341        X2 = Xh.copy(); X2["prime"] = ((q >= 4) & (a >= 40)).astype(float); X2["mur"] = ((a >= lo) & (a <= hi)).astype(float)342        p = valider(Hedonique, "hedo", X2, None, y2, pa); r = Hedonique().fit(X2, y2)343        e["7c"][nom_s] = {**metriques(p, y2, facteur_duan(p, y2)), "coef_mur": float(r.beta[list(X2.columns).index("mur")]),344                          "n_mal_classees": int((X2["mur"] != (((a >= 25) & (a <= 45)).astype(float))).sum())}345    rng = np.random.default_rng(GRAINE); pred = np.full(n, np.nan); choix = []346    for f in pa:347        tr = np.setdiff1d(np.arange(n), f); inner = np.array_split(rng.permutation(tr), 3); err = {}348        for d_ in (1, 2, 3):349            se = 0.0350            for fi in inner:351                tri = np.setdiff1d(tr, fi); m = _boosting(max_depth=d_).fit(Xm.iloc[tri], y[tri]); se += np.sum((m.predict(Xm.iloc[fi]) - y[fi]) ** 2)352            err[d_] = se353        d_opt = min(err, key=err.get); choix.append(d_opt); m = _boosting(max_depth=d_opt).fit(Xm.iloc[tr], y[tr]); pred[f] = m.predict(Xm.iloc[f])354    e["7d"] = {"profondeurs_choisies": choix, **metriques(pred, y, facteur_duan(pred, y))}355    # 37.6 : validation par secteur vs aléatoire (chiffres du script)356    e["6"] = {nom: {"aleatoire": comp["metriques"][nom]["aleatoire"]["COD"], "spatial": comp["metriques"][nom]["spatial"]["COD"], "secteur": secteurs["metriques"][nom]["COD"]}357              for nom in ("Hédonique (référence)", "Forêt aléatoire", "Gradient boosting")}358    e["6"]["biais_secteur"] = secteurs["biais_secteur"]359    # 37.7 : scénario — respécification360    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()361    return e362363# ----------------------------------------------------------------------------- orchestration364def main(verbose: bool = True, complet: bool = True) -> dict:365    d = cas1.charger(); v = d["ventes"].reset_index(drop=True); s = d["sujet"]366    ref = cas1.hedonique(v); Xh = cas1.matrice(v); Xm = matrice_ml(v); y = np.log(v["prix_vente"].astype(float).to_numpy())367    sd = cas1.sujet_df(s); Xh_s = cas1.matrice(sd)[Xh.columns]; Xm_s = matrice_ml(sd, s)368    comp = comparaison(v, Xh, Xm, y)369    det = boosting_detendance(v, Xh, Xm, y)370    sect = validation_secteurs(v, Xh, Xm, y)371    scen = scenario_nonlineaire(v, Xh, Xm, y)372    ajustes = {nom: (ctor().fit(Xh if genre == "hedo" else Xm, y)) for nom, (genre, ctor) in MODELES.items()}373    sujet = prediction_sujet(ajustes, comp["duan"], Xh_s, Xm_s, ref, s)374    sh = shap_boosting(ajustes["Gradient boosting"], Xm, Xm_s)375    sh_sc = shap_boosting(scen["gbm_scenario"], Xm, Xm_s); sh["scenario_interactions"] = sh_sc["interactions"]; sh["scenario_importance"] = sh_sc["importance"]376    pa = comp["plis"]["aleatoire"]377    perm = importance_permutation(Xh, Xm, y, pa) if complet else None378    grille_S = np.linspace(80, 260, 37); grille_A = np.linspace(1, 95, 48)379    pdp = {"superficie": {"grille": grille_S, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "superficie", grille_S),380                          "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "superficie", grille_S),381                          "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)},382           "age": {"grille": grille_A, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "age", grille_A),383                   "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "age", grille_A), "vraie": pdp_vraie(v, "age", grille_A)}}384    # hédonique : l'âge entre par deux colonnes (age, age2_100) → PDP par substitution conjointe385    def pdp_hedo_age(modele, Xh_, grille):386        out = np.empty(len(grille))387        for i, g in enumerate(grille):388            Xg = Xh_.copy(); Xg["age"] = g; Xg["age2_100"] = g ** 2 / 100.0; out[i] = float(np.mean(modele.predict(Xg)))389        return out390    pdp["age"]["Hédonique (référence)"] = pdp_hedo_age(ajustes["Hédonique (référence)"], Xh, grille_A)391    pdp["age_scenario"] = {"grille": grille_A, "Gradient boosting": dependance_partielle(scen["gbm_scenario"], Xm, "age", grille_A),392                           "Hédonique (référence)": pdp_hedo_age(scen["hedo_scenario"], Xh, grille_A), "vraie": pdp_vraie(v, "age", grille_A, scenario=True)}393    courbes = courbes_apprentissage(Xm, y, pa) if complet else None394    sens = sensibilite_hyperparametres(Xm, y, pa) if complet else None395    ex = exercices(v, comp, sect, scen, sh)396    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,397           "scenario": scen, "ajustes": ajustes, "sujet_pred": sujet, "shap": sh, "permutation": perm, "pdp": pdp, "courbes": courbes, "sensibilite": sens, "exercices": ex}398    if verbose:399        pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)400        print("== Modèle hédonique de référence : R² %.3f, σ %.4f, n %d, k %d" % (ref["r2"], ref["sigma"], ref["n"], ref["k"]))401        for schema in ("aleatoire", "spatial", "temporel"):402            print(f"== Validation {schema} :")403            for nom, m in comp["metriques"].items():404                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}")405        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()})406        md = det["metriques_detendance"]; m1 = det["metriques_feuille1"]407        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']}")408        print(f"   Boosting feuilles ≥ 1 (test 2026) : biais={m1['biais_log']:+.4f} RMSElog={m1['rmse_log']:.4f} COD={m1['COD']:.2f}")409        print("== Un secteur exclu :")410        for nom, m in sect["metriques"].items():411            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()})412        print(f"== Scénario non linéaire (n prime={scen['n_prime']}, n mur={scen['n_mur']}) :")413        for nom, m in scen["metriques"].items():414            print(f"   {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} PRB={m['PRB']:.3f}")415        print("   coefficients respécifiés :", {k: round(vv, 4) for k, vv in scen["coef_respecifie"].items()})416        print("== Sujet :"); print(sujet.round(4).to_string())417        print(f"== SHAP (boosting) : base {sh['base']:.4f} (exp {math.exp(sh['base']):,.0f} $) ; efficacité max err {sh['efficacite_max_err']:.2e}")418        print("   importance moyenne |φ| :", sh["importance"].round(4).to_dict())419        print("   sujet φ :", sh["sujet"].sort_values().round(4).to_dict())420        print("   interactions (base) :", sh["interactions"].head(4).round(4).to_dict()); print("   interactions (scénario) :", sh["scenario_interactions"].head(4).round(4).to_dict())421        if perm is not None:422            print("== Importance par permutation (ΔEQM log, hors échantillon) :")423            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())424        if courbes is not None:425            print("== Courbes d'apprentissage :", {d_: (c["optimum"], round(c["rmse_min"], 4), round(float(c["entrainement"][-1]), 4)) for d_, c in courbes.items()})426        if sens is not None:427            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))428        print("== Exercices :")429        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"]])430        print("   37.2", {k: (round(vv, 6) if isinstance(vv, float) else vv) for k, vv in ex["2"].items()})431        print("   37.3", {k: (np.round(vv, 4).tolist() if isinstance(vv, list) else round(vv, 4)) for k, vv in ex["3"].items()})432        print("   37.4", {k: (round(vv, 3) if isinstance(vv, float) else vv) for k, vv in ex["4"].items()})433        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"]])434        print("   37.6", ex["6"]); print("   37.7", ex["7"])435        print("   37.1 essais", [(round(x["seuil"], 1), round(x["gain"], 3)) for x in ex["1"]["essais"]])436        print("   37.6c", {k: (round(vv["rmse_log"], 4), round(vv["COD"], 2)) for k, vv in ex["6c"].items() if k != "biais_secteur"}, {k: round(vv, 3) for k, vv in ex["6c"]["biais_secteur"].items()})437        print("   37.7c", {k: (round(vv["rmse_log"], 4), round(vv["COD"], 2), round(vv["coef_mur"], 4), vv["n_mal_classees"]) for k, vv in ex["7c"].items()})438        print("   37.7d", ex["7d"]["profondeurs_choisies"], round(ex["7d"]["rmse_log"], 4), round(ex["7d"]["COD"], 2))439    return out440441if __name__ == "__main__":442    main()443