spb/traite-evaluation-qc

1"""Outils économétriques : hédonique semi-log, indices (dummies temporelles, ventes répétées BMN / Case–Shiller),2statistiques de ratio IAAO (COD, PRD, PRB), autocorrélation spatiale (I de Moran), matrices de poids kNN/distance."""3from __future__ import annotations4import numpy as np56# ----------------------------------------------------------------------------- MCO78def mco(X: np.ndarray, y: np.ndarray, robuste: bool = True) -> dict:9    """Moindres carrés ordinaires avec erreurs-types classiques et HC1 (White). X doit contenir la constante."""10    X = np.asarray(X, float); y = np.asarray(y, float)11    n, k = X.shape12    XtX_inv = np.linalg.pinv(X.T @ X)13    beta = XtX_inv @ X.T @ y14    res = y - X @ beta15    s2 = res @ res / (n - k)16    V_cl = s2 * XtX_inv17    meat = (X * (res[:, None] ** 2)).T @ X18    V_hc1 = XtX_inv @ meat @ XtX_inv * n / (n - k)19    tss = ((y - y.mean()) ** 2).sum()20    r2 = 1 - (res @ res) / tss21    return {"beta": beta, "se": np.sqrt(np.diag(V_cl)), "se_hc1": np.sqrt(np.diag(V_hc1)), "residus": res,22            "sigma2": s2, "r2": r2, "r2_ajuste": 1 - (1 - r2) * (n - 1) / (n - k), "n": n, "k": k, "V_hc1": V_hc1}2324def indice_dummies_temps(beta_temps: np.ndarray, base: float = 100.0) -> np.ndarray:25    """Indice hédonique à partir des coefficients des indicatrices de période (période 0 = référence, coefficient 0)."""26    b = np.concatenate([[0.0], np.asarray(beta_temps, float)])27    return base * np.exp(b)2829# ----------------------------------------------------------------------------- ventes répétées3031def ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:32    """Estimateur de Bailey–Muth–Nourse (1963) : ln(P2/P1) = Σ_t D_t β_t + ε, D = −1 à t1, +1 à t2, période 0 omise."""33    p1 = np.asarray(prix_1, float); p2 = np.asarray(prix_2, float)34    t1 = np.asarray(t_1, int); t2 = np.asarray(t_2, int)35    y = np.log(p2 / p1)36    D = np.zeros((len(y), n_periodes))37    for i in range(len(y)):38        D[i, t1[i]] -= 1; D[i, t2[i]] += 139    D = D[:, 1:]40    beta, *_ = np.linalg.lstsq(D, y, rcond=None)41    res = y - D @ beta42    idx = base * np.exp(np.concatenate([[0.0], beta]))43    return {"beta": beta, "indice": idx, "residus": res, "D": D, "y": y}4445def ventes_repetees_case_shiller(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:46    """Case–Shiller (1987, 1989) pondéré en trois étapes : (1) BMN, (2) régression de ε² sur l'intervalle,47    (3) moindres carrés pondérés par 1/√(variance ajustée)."""48    e1 = ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes, base)49    gap = (np.asarray(t_2) - np.asarray(t_1)).astype(float)50    Z = np.column_stack([np.ones_like(gap), gap])51    a, *_ = np.linalg.lstsq(Z, e1["residus"] ** 2, rcond=None)52    var_hat = np.maximum(Z @ a, 1e-8)53    w = 1 / np.sqrt(var_hat)54    beta, *_ = np.linalg.lstsq(e1["D"] * w[:, None], e1["y"] * w, rcond=None)55    idx = base * np.exp(np.concatenate([[0.0], beta]))56    return {"beta": beta, "indice": idx, "coef_variance": a, "poids": w}5758# ----------------------------------------------------------------------------- ratios IAAO5960def statistiques_ratio(valeurs_estimees, prix) -> dict:61    """Statistiques d'étude de ratios (IAAO, Standard on Ratio Studies) :62    ratio médian, COD = 100 × moyenne|r − médiane| / médiane, PRD = moyenne(r) / (Σ V / Σ P),63    PRB = pente de (r − médiane)/médiane sur log2 de la valeur proxy ((V/médiane + P)/2)."""64    V = np.asarray(valeurs_estimees, float); P = np.asarray(prix, float)65    r = V / P66    med = np.median(r)67    cod = 100.0 * np.mean(np.abs(r - med)) / med68    prd = np.mean(r) / (V.sum() / P.sum())69    proxy = 0.5 * (V / med + P)70    x = np.log2(proxy)71    yv = (r - med) / med72    X = np.column_stack([np.ones_like(x), x])73    b, *_ = np.linalg.lstsq(X, yv, rcond=None)74    res = yv - X @ b75    se = np.sqrt((res @ res) / (len(x) - 2) / ((x - x.mean()) ** 2).sum())76    return {"n": len(r), "ratio_median": med, "ratio_moyen": r.mean(), "COD": cod, "PRD": prd, "PRB": b[1], "PRB_se": se,77            "ratio_pondere": V.sum() / P.sum()}7879# ----------------------------------------------------------------------------- spatial8081def matrice_poids_knn(coords: np.ndarray, k: int = 8, standardiser: bool = True) -> np.ndarray:82    """Matrice W (n×n) des k plus proches voisins (distance euclidienne), standardisée en ligne."""83    c = np.asarray(coords, float)84    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))85    np.fill_diagonal(d, np.inf)86    W = np.zeros_like(d)87    idx = np.argsort(d, axis=1)[:, :k]88    for i in range(len(c)):89        W[i, idx[i]] = 1.090    if standardiser:91        W = W / W.sum(1, keepdims=True)92    return W9394def matrice_poids_distance_inverse(coords: np.ndarray, portee: float, puissance: float = 1.0, standardiser: bool = True) -> np.ndarray:95    c = np.asarray(coords, float)96    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))97    with np.errstate(divide="ignore"):98        W = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)99    if standardiser:100        s = W.sum(1, keepdims=True); s[s == 0] = 1101        W = W / s102    return W103104def moran_i(x: np.ndarray, W: np.ndarray) -> dict:105    """I de Moran et test sous hypothèse de normalité : I = (n/S0) · (z'Wz)/(z'z)."""106    z = np.asarray(x, float); z = z - z.mean()107    n = len(z); S0 = W.sum()108    I = (n / S0) * (z @ W @ z) / (z @ z)109    EI = -1.0 / (n - 1)110    S1 = 0.5 * ((W + W.T) ** 2).sum()111    S2 = ((W.sum(0) + W.sum(1)) ** 2).sum()112    VI = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / ((n ** 2 - 1) * S0 ** 2) - EI ** 2113    zI = (I - EI) / np.sqrt(VI)114    from math import erf, sqrt115    p = 2 * (1 - 0.5 * (1 + erf(abs(zI) / sqrt(2))))116    return {"I": I, "E_I": EI, "Var_I": VI, "z": zI, "p": p}117118def decalage_spatial(W: np.ndarray, x: np.ndarray) -> np.ndarray:119    return W @ np.asarray(x, float)120121def sar_lag_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, grille_rho=None) -> dict:122    """Modèle autorégressif spatial y = ρWy + Xβ + ε estimé par maximum de vraisemblance concentrée (recherche sur ρ)."""123    y = np.asarray(y, float); X = np.asarray(X, float); n = len(y)124    Wy = W @ y125    b0, *_ = np.linalg.lstsq(X, y, rcond=None); e0 = y - X @ b0126    bL, *_ = np.linalg.lstsq(X, Wy, rcond=None); eL = Wy - X @ bL127    eig = np.linalg.eigvals(W).real128    rhos = grille_rho if grille_rho is not None else np.linspace(-0.9, 0.99, 380)129    best = None130    for rho in rhos:131        e = e0 - rho * eL132        s2 = e @ e / n133        ll = np.log(np.abs(1 - rho * eig)).sum() - n / 2 * np.log(2 * np.pi * s2) - n / 2134        if best is None or ll > best[0]:135            best = (ll, rho, s2)136    ll, rho, s2 = best137    beta = b0 - rho * bL138    return {"rho": rho, "beta": beta, "sigma2": s2, "loglik": ll}139140# ----------------------------------------------------------------------------- spatial, suite (chapitre 19)141# Fonctions ajoutées le 2026-10-31 pour le chapitre 19 (tests : analysis/python/tests/test_ch19.py).142# Elles n'altèrent pas les fonctions précédentes. Déterminants jacobiens calculés avec les valeurs propres143# complexes de W (|det(I − ρW)| = Π |1 − ρ λ_i|), ce qui est exact pour une matrice standardisée en ligne.144145def coordonnees_km(lat, lon, lat0: float | None = None, lon0: float | None = None) -> np.ndarray:146    """Projection équirectangulaire locale : (lat, lon) → (est, nord) en kilomètres autour de (lat0, lon0)."""147    lat = np.asarray(lat, float); lon = np.asarray(lon, float)148    lat0 = float(np.mean(lat)) if lat0 is None else lat0149    lon0 = float(np.mean(lon)) if lon0 is None else lon0150    x = (lon - lon0) * 111.32 * np.cos(np.radians(lat0)); y = (lat - lat0) * 110.57151    return np.column_stack([x, y])152153def _distances(c: np.ndarray) -> np.ndarray:154    c = np.asarray(c, float)155    return np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))156157def _logdet(eig: np.ndarray, rho: float) -> float:158    return float(np.log(np.abs(1.0 - rho * eig)).sum())159160def _intervalle_rho(eig: np.ndarray) -> tuple[float, float]:161    """Intervalle admissible (1/λ_min, 1/λ_max) sur les valeurs propres réelles de W (λ_max = 1 si standardisée)."""162    re = eig.real[np.abs(eig.imag) < 1e-8]163    lo = 1.0 / re.min() if len(re) and re.min() < -1e-8 else -0.999164    hi = 1.0 / re.max() if len(re) and re.max() > 1e-8 else 0.999165    return max(lo, -0.999) + 1e-6, min(hi, 0.999) - 1e-6166167def _optimiser_concentree(f, bornes, n_grille: int = 60):168    """Maximise une vraisemblance concentrée scalaire : grille puis raffinement (section dorée, scipy)."""169    from scipy.optimize import minimize_scalar170    g = np.linspace(bornes[0], bornes[1], n_grille)171    vals = np.array([f(r) for r in g]); i = int(np.argmax(vals))172    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]173    r = minimize_scalar(lambda t: -f(t), bounds=(lo, hi), method="bounded", options={"xatol": 1e-7})174    return float(r.x), float(-r.fun)175176def sar_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, durbin: bool = False) -> dict:177    """SAR (y = ρWy + Xβ + ε) ou SDM (durbin=True : y = ρWy + Xβ + WXθ + ε) par maximum de vraisemblance concentrée.178    X contient la constante en première colonne (non décalée dans le SDM). Erreurs-types asymptotiques (matrice d'information179    analytique, Anselin 1988, ch. 6). Retourne aussi la log-vraisemblance, l'AIC et les résidus."""180    y = np.asarray(y, float); X0 = np.asarray(X, float); n = len(y)181    Xr = np.column_stack([X0, W @ X0[:, 1:]]) if durbin else X0182    k = Xr.shape[1]183    eig = np.linalg.eigvals(W)184    Wy = W @ y185    XtX_inv = np.linalg.pinv(Xr.T @ Xr)186    b0 = XtX_inv @ Xr.T @ y; e0 = y - Xr @ b0187    bL = XtX_inv @ Xr.T @ Wy; eL = Wy - Xr @ bL188    def ll_c(rho):189        e = e0 - rho * eL; s2 = e @ e / n190        return _logdet(eig, rho) - n / 2 * np.log(2 * np.pi * s2) - n / 2191    rho, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))192    beta = b0 - rho * bL; e = e0 - rho * eL; s2 = float(e @ e / n)193    # matrice d'information (β, ρ, σ²)194    A_inv = np.linalg.inv(np.eye(n) - rho * W); WA = W @ A_inv; WAXb = WA @ (Xr @ beta)195    I = np.zeros((k + 2, k + 2))196    I[:k, :k] = Xr.T @ Xr / s2197    I[:k, k] = I[k, :k] = Xr.T @ WAXb / s2198    I[k, k] = np.trace(WA.T @ WA) + np.trace(WA @ WA) + WAXb @ WAXb / s2199    I[k, k + 1] = I[k + 1, k] = np.trace(WA) / s2200    I[k + 1, k + 1] = n / (2 * s2 ** 2)201    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))202    return {"rho": rho, "beta": beta, "se": se[:k], "se_rho": float(se[k]), "sigma2": s2, "loglik": ll,203            "aic": 2 * (k + 2) - 2 * ll, "residus": e, "n": n, "k": k, "X": Xr, "durbin": durbin, "ajuste": y - e}204205def sem_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:206    """Modèle à erreur spatiale (SEM) y = Xβ + u, u = λWu + ε, par maximum de vraisemblance concentrée sur λ.207    Retourne λ, β (MCG spatiaux), σ², log-vraisemblance, AIC, résidus u (= y − Xβ) et innovations ε (= u − λWu)."""208    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape209    eig = np.linalg.eigvals(W); Wy = W @ y; WX = W @ X210    def ajuster(lam):211        ys = y - lam * Wy; Xs = X - lam * WX212        b = np.linalg.pinv(Xs.T @ Xs) @ Xs.T @ ys; eps = ys - Xs @ b213        return b, eps, float(eps @ eps / n), Xs214    def ll_c(lam):215        _, eps, s2, _ = ajuster(lam)216        return _logdet(eig, lam) - n / 2 * np.log(2 * np.pi * s2) - n / 2217    lam, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))218    beta, eps, s2, Xs = ajuster(lam)219    B_inv = np.linalg.inv(np.eye(n) - lam * W); WB = W @ B_inv220    I = np.zeros((k + 2, k + 2))221    I[:k, :k] = Xs.T @ Xs / s2222    I[k, k] = np.trace(WB.T @ WB) + np.trace(WB @ WB)223    I[k, k + 1] = I[k + 1, k] = np.trace(WB) / s2224    I[k + 1, k + 1] = n / (2 * s2 ** 2)225    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))226    u = y - X @ beta227    return {"lambda": lam, "beta": beta, "se": se[:k], "se_lambda": float(se[k]), "sigma2": s2, "loglik": ll,228            "aic": 2 * (k + 2) - 2 * ll, "residus": u, "innovations": eps, "n": n, "k": k, "ajuste": X @ beta}229230def tests_lm_spatial(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:231    """Tests du multiplicateur de Lagrange sur les résidus MCO (Anselin 1988 ; versions robustes d'Anselin, Bera, Florax et Yoon 1996) :232    LM-erreur, LM-décalage, RLM-erreur, RLM-décalage (χ² à 1 degré de liberté) et SARMA (2 d.l.)."""233    from scipy.stats import chi2234    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape235    XtX_inv = np.linalg.pinv(X.T @ X); b = XtX_inv @ X.T @ y; e = y - X @ b; s2 = e @ e / n236    M = np.eye(n) - X @ XtX_inv @ X.T237    T = np.trace((W.T + W) @ W)238    WXb = W @ (X @ b)239    RJ = (WXb @ M @ WXb) / s2 + T240    d_rho = (e @ W @ y) / s2; d_lam = (e @ W @ e) / s2241    lm_err = d_lam ** 2 / T; lm_lag = d_rho ** 2 / RJ242    rlm_lag = (d_rho - d_lam) ** 2 / (RJ - T)243    rlm_err = (d_lam - T / RJ * d_rho) ** 2 / (T * (1 - T / RJ))244    sarma = lm_err + rlm_lag245    p = lambda s, df=1: float(chi2.sf(s, df))246    return {"LM_erreur": (float(lm_err), p(lm_err)), "LM_decalage": (float(lm_lag), p(lm_lag)),247            "RLM_erreur": (float(rlm_err), p(rlm_err)), "RLM_decalage": (float(rlm_lag), p(rlm_lag)), "SARMA": (float(sarma), p(sarma, 2))}248249def effets_sar(rho: float, beta: np.ndarray, W: np.ndarray, theta: np.ndarray | None = None) -> dict:250    """Effets direct, indirect et total (LeSage et Pace 2009) : S_r = (I − ρW)^{-1}(β_r I + θ_r W) ;251    direct = tr(S_r)/n, total = 1'S_r 1/n, indirect = total − direct. theta = coefficients des WX (SDM), None pour le SAR."""252    n = W.shape[0]; A_inv = np.linalg.inv(np.eye(n) - rho * W)253    beta = np.asarray(beta, float); theta = np.zeros_like(beta) if theta is None else np.asarray(theta, float)254    direct, total = [], []255    for br, tr in zip(beta, theta):256        S = A_inv @ (br * np.eye(n) + tr * W)257        direct.append(np.trace(S) / n); total.append(S.sum() / n)258    direct = np.array(direct); total = np.array(total)259    return {"direct": direct, "indirect": total - direct, "total": total}260261def matrice_poids_spatio_temporelle(coords: np.ndarray, temps: np.ndarray, portee: float, fenetre: float,262                                     puissance: float = 1.0, unidirectionnel: bool = True, standardiser: bool = True) -> np.ndarray:263    """Matrice spatio-temporelle à la Dubé et Legros (2013) : w_ij = s_ij · t_ij, avec s_ij = 1/d_ij^p pour 0 < d_ij ≤ portée264    et t_ij = 1 si la vente j précède la vente i de moins de `fenetre` unités de temps (unidirectionnel : seules les ventes265    antérieures informent). temps : dates en unités numériques (mois). Standardisée en ligne (lignes vides laissées nulles)."""266    d = _distances(coords); t = np.asarray(temps, float)267    dt = t[:, None] - t[None, :]           # > 0 si j est antérieure à i268    with np.errstate(divide="ignore"):269        S = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)270    if unidirectionnel:271        Tm = ((dt > 0) & (dt <= fenetre)).astype(float)272    else:273        Tm = ((dt != 0) & (np.abs(dt) <= fenetre)).astype(float)274    W = S * Tm275    if standardiser:276        s = W.sum(1, keepdims=True); s[s == 0] = 1.0; W = W / s277    return W278279def moran_par_distance(x: np.ndarray, coords: np.ndarray, bornes) -> list[dict]:280    """Corrélogramme : I de Moran calculé avec une matrice binaire des paires dont la distance tombe dans chaque bande [b_k, b_{k+1})."""281    d = _distances(coords); out = []282    for lo, hi in zip(bornes[:-1], bornes[1:]):283        W = ((d >= lo) & (d < hi)).astype(float); np.fill_diagonal(W, 0.0)284        npaires = int(W.sum() / 2)285        if npaires < 10:286            continue287        r = moran_i(x, W); r.update({"d_min": lo, "d_max": hi, "paires": npaires}); out.append(r)288    return out289290def gwr(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bande: float, cibles: np.ndarray | None = None) -> dict:291    """Régression géographiquement pondérée à noyau gaussien fixe : w_ij = exp(−½ (d_ij / h)²).292    Coefficients locaux à chaque cible (par défaut les observations), valeurs ajustées, trace de la matrice chapeau,293    score de validation croisée (somme des carrés des erreurs leave-one-out) et AICc (Fotheringham, Brunsdon et Charlton 2002)."""294    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape295    c = np.asarray(coords, float); cib = c if cibles is None else np.asarray(cibles, float)296    d = np.sqrt(((cib[:, None, :] - c[None, :, :]) ** 2).sum(-1))297    Wk = np.exp(-0.5 * (d / bande) ** 2)298    betas = np.zeros((len(cib), k))299    for i in range(len(cib)):300        Xw = X * Wk[i][:, None]301        betas[i] = np.linalg.solve(X.T @ Xw + 1e-10 * np.eye(k), Xw.T @ y)302    out = {"beta": betas, "bande": bande}303    if cibles is None:304        yhat = (X * betas).sum(1)305        # matrice chapeau et validation croisée leave-one-out306        S_diag = np.zeros(n); cv = 0.0307        for i in range(n):308            Xw = X * Wk[i][:, None]309            Ainv = np.linalg.inv(X.T @ Xw + 1e-10 * np.eye(k))310            S_diag[i] = X[i] @ Ainv @ Xw[i]311            w_loo = Wk[i].copy(); w_loo[i] = 0.0312            b_loo = np.linalg.solve(X.T @ (X * w_loo[:, None]) + 1e-10 * np.eye(k), (X * w_loo[:, None]).T @ y)313            cv += (y[i] - X[i] @ b_loo) ** 2314        res = y - yhat; rss = float(res @ res); tr_S = float(S_diag.sum())315        aicc = n * np.log(rss / n) + n * np.log(2 * np.pi) + n * (n + tr_S) / (n - 2 - tr_S)316        out.update({"ajuste": yhat, "residus": res, "trace_S": tr_S, "cv": float(cv), "aicc": float(aicc),317                    "r2": 1 - rss / float(((y - y.mean()) ** 2).sum()), "sigma": float(np.sqrt(rss / (n - tr_S)))})318    return out319320def gwr_bande_cv(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bornes=(0.3, 15.0), n_grille: int = 25) -> dict:321    """Bande passante (km) minimisant le score de validation croisée du noyau gaussien : grille log puis section dorée."""322    from scipy.optimize import minimize_scalar323    g = np.exp(np.linspace(np.log(bornes[0]), np.log(bornes[1]), n_grille))324    cv = np.array([gwr(y, X, coords, h)["cv"] for h in g]); i = int(np.argmin(cv))325    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]326    r = minimize_scalar(lambda h: gwr(y, X, coords, h)["cv"], bounds=(lo, hi), method="bounded", options={"xatol": 1e-3})327    return {"bande": float(r.x), "cv": float(r.fun), "grille": g, "cv_grille": cv}328329def variogramme_empirique(coords: np.ndarray, z: np.ndarray, bornes) -> dict:330    """Semi-variogramme empirique (Matheron) : γ(h) = Σ_{paires dans la bande} (z_i − z_j)² / (2 N(h))."""331    d = _distances(coords); z = np.asarray(z, float)332    dz2 = (z[:, None] - z[None, :]) ** 2333    iu = np.triu_indices(len(z), 1); dd = d[iu]; zz = dz2[iu]334    h, g, nb = [], [], []335    for lo, hi in zip(bornes[:-1], bornes[1:]):336        m = (dd >= lo) & (dd < hi)337        if m.sum() < 30:338            continue339        h.append(dd[m].mean()); g.append(zz[m].sum() / (2 * m.sum())); nb.append(int(m.sum()))340    return {"h": np.array(h), "gamma": np.array(g), "n": np.array(nb)}341342def variogramme_exponentiel(h, pepite: float, palier_partiel: float, portee: float) -> np.ndarray:343    """γ(h) = c₀ + c₁ (1 − exp(−h/a)) pour h > 0 ; portée pratique ≈ 3a."""344    h = np.asarray(h, float)345    return np.where(h > 0, pepite + palier_partiel * (1 - np.exp(-h / portee)), 0.0)346347def ajuster_variogramme(h, gamma, n_paires, portee_init: float = 2.0, portee_max: float = 10.0) -> dict:348    """Ajustement du modèle exponentiel par moindres carrés pondérés (poids N(h)/γ_modèle², Cressie 1985) ; le paramètre349    de portée a est borné par portee_max (portée pratique 3a) pour éviter les solutions dégénérées quand le palier n'est pas atteint."""350    from scipy.optimize import least_squares351    h = np.asarray(h, float); g = np.asarray(gamma, float); nb = np.asarray(n_paires, float)352    def res(p):353        gm = variogramme_exponentiel(h, *p); return np.sqrt(nb) * (g - gm) / np.maximum(gm, 1e-9)354    p0 = [0.5 * g[0], max(g[-1] - 0.5 * g[0], 1e-6), portee_init]355    r = least_squares(res, p0, bounds=([0, 0, 1e-3], [np.inf, np.inf, portee_max]))356    c0, c1, a = r.x357    return {"pepite": float(c0), "palier_partiel": float(c1), "portee": float(a), "portee_pratique": float(3 * a), "palier": float(c0 + c1)}358359def krigeage_ordinaire(coords_obs: np.ndarray, z: np.ndarray, coords_cibles: np.ndarray, params: dict, voisins: int | None = 40) -> dict:360    """Krigeage ordinaire avec covariance exponentielle C(h) = c₁ exp(−h/a) (+ c₀ en h = 0), voisinage des `voisins` points361    les plus proches (None = tous). Retourne la prédiction et la variance de krigeage à chaque cible."""362    co = np.asarray(coords_obs, float); z = np.asarray(z, float); ct = np.asarray(coords_cibles, float)363    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]364    cov = lambda h: c1 * np.exp(-h / a)365    dct = np.sqrt(((ct[:, None, :] - co[None, :, :]) ** 2).sum(-1))366    pred = np.zeros(len(ct)); var = np.zeros(len(ct))367    for i in range(len(ct)):368        idx = np.argsort(dct[i])[:voisins] if voisins else np.arange(len(co))369        sub = co[idx]; dd = _distances(sub); m = len(idx)370        K = np.zeros((m + 1, m + 1)); K[:m, :m] = cov(dd); K[np.arange(m), np.arange(m)] = c0 + c1371        K[m, :m] = 1.0; K[:m, m] = 1.0372        rhs = np.concatenate([cov(dct[i, idx]), [1.0]])373        sol = np.linalg.solve(K, rhs); w = sol[:m]; mu = sol[m]374        pred[i] = w @ z[idx]; var[i] = (c0 + c1) - w @ rhs[:m] - mu375    return {"prediction": pred, "variance": var}376377def krigeage_validation_croisee(coords: np.ndarray, z: np.ndarray, params: dict, voisins: int = 40) -> np.ndarray:378    """Prédiction leave-one-out de chaque observation par krigeage ordinaire de ses voisins (hors elle-même)."""379    c = np.asarray(coords, float); z = np.asarray(z, float); d = _distances(c); out = np.zeros(len(z))380    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]; cov = lambda h: c1 * np.exp(-h / a)381    for i in range(len(z)):382        idx = np.argsort(d[i])[1:voisins + 1]; m = len(idx); dd = d[np.ix_(idx, idx)]383        K = np.zeros((m + 1, m + 1)); K[:m, :m] = cov(dd); K[np.arange(m), np.arange(m)] = c0 + c1; K[m, :m] = 1; K[:m, m] = 1384        rhs = np.concatenate([cov(d[i, idx]), [1.0]]); sol = np.linalg.solve(K, rhs); out[i] = sol[:m] @ z[idx]385    return out386387def sem_esperance_conditionnelle(lam: float, W_aug: np.ndarray, u_obs: np.ndarray, idx_cibles) -> np.ndarray:388    """Espérance conditionnelle des perturbations u aux cibles sachant les u observés, dans un SEM de matrice augmentée W_aug389    (observations puis cibles) : avec Q = (I − λW)'(I − λW) (précision à σ² près), E[u_c | u_o] = −Q_cc^{-1} Q_co u_o."""390    n = W_aug.shape[0]; B = np.eye(n) - lam * W_aug; Q = B.T @ B391    idx_c = np.asarray(idx_cibles); idx_o = np.setdiff1d(np.arange(n), idx_c)392    Qcc = Q[np.ix_(idx_c, idx_c)]; Qco = Q[np.ix_(idx_c, idx_o)]393    return -np.linalg.solve(Qcc, Qco @ np.asarray(u_obs, float))394