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 re.min() < 0 else -0.999164    return max(lo, -0.999) + 1e-6, 1.0 / re.max() - 1e-6165166def _optimiser_concentree(f, bornes, n_grille: int = 60):167    """Maximise une vraisemblance concentrée scalaire : grille puis raffinement (section dorée, scipy)."""168    from scipy.optimize import minimize_scalar169    g = np.linspace(bornes[0], bornes[1], n_grille)170    vals = np.array([f(r) for r in g]); i = int(np.argmax(vals))171    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]172    r = minimize_scalar(lambda t: -f(t), bounds=(lo, hi), method="bounded", options={"xatol": 1e-7})173    return float(r.x), float(-r.fun)174175def sar_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, durbin: bool = False) -> dict:176    """SAR (y = ρWy + Xβ + ε) ou SDM (durbin=True : y = ρWy + Xβ + WXθ + ε) par maximum de vraisemblance concentrée.177    X contient la constante en première colonne (non décalée dans le SDM). Erreurs-types asymptotiques (matrice d'information178    analytique, Anselin 1988, ch. 6). Retourne aussi la log-vraisemblance, l'AIC et les résidus."""179    y = np.asarray(y, float); X0 = np.asarray(X, float); n = len(y)180    Xr = np.column_stack([X0, W @ X0[:, 1:]]) if durbin else X0181    k = Xr.shape[1]182    eig = np.linalg.eigvals(W)183    Wy = W @ y184    XtX_inv = np.linalg.pinv(Xr.T @ Xr)185    b0 = XtX_inv @ Xr.T @ y; e0 = y - Xr @ b0186    bL = XtX_inv @ Xr.T @ Wy; eL = Wy - Xr @ bL187    def ll_c(rho):188        e = e0 - rho * eL; s2 = e @ e / n189        return _logdet(eig, rho) - n / 2 * np.log(2 * np.pi * s2) - n / 2190    rho, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))191    beta = b0 - rho * bL; e = e0 - rho * eL; s2 = float(e @ e / n)192    # matrice d'information (β, ρ, σ²)193    A_inv = np.linalg.inv(np.eye(n) - rho * W); WA = W @ A_inv; WAXb = WA @ (Xr @ beta)194    I = np.zeros((k + 2, k + 2))195    I[:k, :k] = Xr.T @ Xr / s2196    I[:k, k] = I[k, :k] = Xr.T @ WAXb / s2197    I[k, k] = np.trace(WA.T @ WA) + np.trace(WA @ WA) + WAXb @ WAXb / s2198    I[k, k + 1] = I[k + 1, k] = np.trace(WA) / s2199    I[k + 1, k + 1] = n / (2 * s2 ** 2)200    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))201    return {"rho": rho, "beta": beta, "se": se[:k], "se_rho": float(se[k]), "sigma2": s2, "loglik": ll,202            "aic": 2 * (k + 2) - 2 * ll, "residus": e, "n": n, "k": k, "X": Xr, "durbin": durbin, "ajuste": y - e}203204def sem_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:205    """Modèle à erreur spatiale (SEM) y = Xβ + u, u = λWu + ε, par maximum de vraisemblance concentrée sur λ.206    Retourne λ, β (MCG spatiaux), σ², log-vraisemblance, AIC, résidus u (= y − Xβ) et innovations ε (= u − λWu)."""207    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape208    eig = np.linalg.eigvals(W); Wy = W @ y; WX = W @ X209    def ajuster(lam):210        ys = y - lam * Wy; Xs = X - lam * WX211        b = np.linalg.pinv(Xs.T @ Xs) @ Xs.T @ ys; eps = ys - Xs @ b212        return b, eps, float(eps @ eps / n), Xs213    def ll_c(lam):214        _, eps, s2, _ = ajuster(lam)215        return _logdet(eig, lam) - n / 2 * np.log(2 * np.pi * s2) - n / 2216    lam, ll = _optimiser_concentree(ll_c, _intervalle_rho(eig))217    beta, eps, s2, Xs = ajuster(lam)218    B_inv = np.linalg.inv(np.eye(n) - lam * W); WB = W @ B_inv219    I = np.zeros((k + 2, k + 2))220    I[:k, :k] = Xs.T @ Xs / s2221    I[k, k] = np.trace(WB.T @ WB) + np.trace(WB @ WB)222    I[k, k + 1] = I[k + 1, k] = np.trace(WB) / s2223    I[k + 1, k + 1] = n / (2 * s2 ** 2)224    V = np.linalg.inv(I); se = np.sqrt(np.diag(V))225    u = y - X @ beta226    return {"lambda": lam, "beta": beta, "se": se[:k], "se_lambda": float(se[k]), "sigma2": s2, "loglik": ll,227            "aic": 2 * (k + 2) - 2 * ll, "residus": u, "innovations": eps, "n": n, "k": k, "ajuste": X @ beta}228229def tests_lm_spatial(y: np.ndarray, X: np.ndarray, W: np.ndarray) -> dict:230    """Tests du multiplicateur de Lagrange sur les résidus MCO (Anselin 1988 ; versions robustes d'Anselin, Bera, Florax et Yoon 1996) :231    LM-erreur, LM-décalage, RLM-erreur, RLM-décalage (χ² à 1 degré de liberté) et SARMA (2 d.l.)."""232    from scipy.stats import chi2233    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape234    XtX_inv = np.linalg.pinv(X.T @ X); b = XtX_inv @ X.T @ y; e = y - X @ b; s2 = e @ e / n235    M = np.eye(n) - X @ XtX_inv @ X.T236    T = np.trace((W.T + W) @ W)237    WXb = W @ (X @ b)238    RJ = (WXb @ M @ WXb) / s2 + T239    d_rho = (e @ W @ y) / s2; d_lam = (e @ W @ e) / s2240    lm_err = d_lam ** 2 / T; lm_lag = d_rho ** 2 / RJ241    rlm_lag = (d_rho - d_lam) ** 2 / (RJ - T)242    rlm_err = (d_lam - T / RJ * d_rho) ** 2 / (T * (1 - T / RJ))243    sarma = lm_err + rlm_lag244    p = lambda s, df=1: float(chi2.sf(s, df))245    return {"LM_erreur": (float(lm_err), p(lm_err)), "LM_decalage": (float(lm_lag), p(lm_lag)),246            "RLM_erreur": (float(rlm_err), p(rlm_err)), "RLM_decalage": (float(rlm_lag), p(rlm_lag)), "SARMA": (float(sarma), p(sarma, 2))}247248def effets_sar(rho: float, beta: np.ndarray, W: np.ndarray, theta: np.ndarray | None = None) -> dict:249    """Effets direct, indirect et total (LeSage et Pace 2009) : S_r = (I − ρW)^{-1}(β_r I + θ_r W) ;250    direct = tr(S_r)/n, total = 1'S_r 1/n, indirect = total − direct. theta = coefficients des WX (SDM), None pour le SAR."""251    n = W.shape[0]; A_inv = np.linalg.inv(np.eye(n) - rho * W)252    beta = np.asarray(beta, float); theta = np.zeros_like(beta) if theta is None else np.asarray(theta, float)253    direct, total = [], []254    for br, tr in zip(beta, theta):255        S = A_inv @ (br * np.eye(n) + tr * W)256        direct.append(np.trace(S) / n); total.append(S.sum() / n)257    direct = np.array(direct); total = np.array(total)258    return {"direct": direct, "indirect": total - direct, "total": total}259260def matrice_poids_spatio_temporelle(coords: np.ndarray, temps: np.ndarray, portee: float, fenetre: float,261                                     puissance: float = 1.0, unidirectionnel: bool = True, standardiser: bool = True) -> np.ndarray:262    """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ée263    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 ventes264    antérieures informent). temps : dates en unités numériques (mois). Standardisée en ligne (lignes vides laissées nulles)."""265    d = _distances(coords); t = np.asarray(temps, float)266    dt = t[:, None] - t[None, :]           # > 0 si j est antérieure à i267    with np.errstate(divide="ignore"):268        S = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)269    if unidirectionnel:270        Tm = ((dt > 0) & (dt <= fenetre)).astype(float)271    else:272        Tm = ((dt != 0) & (np.abs(dt) <= fenetre)).astype(float)273    W = S * Tm274    if standardiser:275        s = W.sum(1, keepdims=True); s[s == 0] = 1.0; W = W / s276    return W277278def moran_par_distance(x: np.ndarray, coords: np.ndarray, bornes) -> list[dict]:279    """Corrélogramme : I de Moran calculé avec une matrice binaire des paires dont la distance tombe dans chaque bande [b_k, b_{k+1})."""280    d = _distances(coords); out = []281    for lo, hi in zip(bornes[:-1], bornes[1:]):282        W = ((d >= lo) & (d < hi)).astype(float); np.fill_diagonal(W, 0.0)283        npaires = int(W.sum() / 2)284        if npaires < 10:285            continue286        r = moran_i(x, W); r.update({"d_min": lo, "d_max": hi, "paires": npaires}); out.append(r)287    return out288289def gwr(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bande: float, cibles: np.ndarray | None = None) -> dict:290    """Régression géographiquement pondérée à noyau gaussien fixe : w_ij = exp(−½ (d_ij / h)²).291    Coefficients locaux à chaque cible (par défaut les observations), valeurs ajustées, trace de la matrice chapeau,292    score de validation croisée (somme des carrés des erreurs leave-one-out) et AICc (Fotheringham, Brunsdon et Charlton 2002)."""293    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape294    c = np.asarray(coords, float); cib = c if cibles is None else np.asarray(cibles, float)295    d = np.sqrt(((cib[:, None, :] - c[None, :, :]) ** 2).sum(-1))296    Wk = np.exp(-0.5 * (d / bande) ** 2)297    betas = np.zeros((len(cib), k))298    for i in range(len(cib)):299        Xw = X * Wk[i][:, None]300        betas[i] = np.linalg.solve(X.T @ Xw + 1e-10 * np.eye(k), Xw.T @ y)301    out = {"beta": betas, "bande": bande}302    if cibles is None:303        yhat = (X * betas).sum(1)304        # matrice chapeau et validation croisée leave-one-out305        S_diag = np.zeros(n); cv = 0.0306        for i in range(n):307            Xw = X * Wk[i][:, None]308            Ainv = np.linalg.inv(X.T @ Xw + 1e-10 * np.eye(k))309            S_diag[i] = X[i] @ Ainv @ Xw[i]310            w_loo = Wk[i].copy(); w_loo[i] = 0.0311            b_loo = np.linalg.solve(X.T @ (X * w_loo[:, None]) + 1e-10 * np.eye(k), (X * w_loo[:, None]).T @ y)312            cv += (y[i] - X[i] @ b_loo) ** 2313        res = y - yhat; rss = float(res @ res); tr_S = float(S_diag.sum())314        aicc = n * np.log(rss / n) + n * np.log(2 * np.pi) + n * (n + tr_S) / (n - 2 - tr_S)315        out.update({"ajuste": yhat, "residus": res, "trace_S": tr_S, "cv": float(cv), "aicc": float(aicc),316                    "r2": 1 - rss / float(((y - y.mean()) ** 2).sum()), "sigma": float(np.sqrt(rss / (n - tr_S)))})317    return out318319def gwr_bande_cv(y: np.ndarray, X: np.ndarray, coords: np.ndarray, bornes=(0.3, 15.0), n_grille: int = 25) -> dict:320    """Bande passante (km) minimisant le score de validation croisée du noyau gaussien : grille log puis section dorée."""321    from scipy.optimize import minimize_scalar322    g = np.exp(np.linspace(np.log(bornes[0]), np.log(bornes[1]), n_grille))323    cv = np.array([gwr(y, X, coords, h)["cv"] for h in g]); i = int(np.argmin(cv))324    lo = g[max(i - 1, 0)]; hi = g[min(i + 1, n_grille - 1)]325    r = minimize_scalar(lambda h: gwr(y, X, coords, h)["cv"], bounds=(lo, hi), method="bounded", options={"xatol": 1e-3})326    return {"bande": float(r.x), "cv": float(r.fun), "grille": g, "cv_grille": cv}327328def variogramme_empirique(coords: np.ndarray, z: np.ndarray, bornes) -> dict:329    """Semi-variogramme empirique (Matheron) : γ(h) = Σ_{paires dans la bande} (z_i − z_j)² / (2 N(h))."""330    d = _distances(coords); z = np.asarray(z, float)331    dz2 = (z[:, None] - z[None, :]) ** 2332    iu = np.triu_indices(len(z), 1); dd = d[iu]; zz = dz2[iu]333    h, g, nb = [], [], []334    for lo, hi in zip(bornes[:-1], bornes[1:]):335        m = (dd >= lo) & (dd < hi)336        if m.sum() < 30:337            continue338        h.append(dd[m].mean()); g.append(zz[m].sum() / (2 * m.sum())); nb.append(int(m.sum()))339    return {"h": np.array(h), "gamma": np.array(g), "n": np.array(nb)}340341def variogramme_exponentiel(h, pepite: float, palier_partiel: float, portee: float) -> np.ndarray:342    """γ(h) = c₀ + c₁ (1 − exp(−h/a)) pour h > 0 ; portée pratique ≈ 3a."""343    h = np.asarray(h, float)344    return np.where(h > 0, pepite + palier_partiel * (1 - np.exp(-h / portee)), 0.0)345346def ajuster_variogramme(h, gamma, n_paires, portee_init: float = 2.0) -> dict:347    """Ajustement du modèle exponentiel par moindres carrés pondérés (poids N(h)/γ_modèle², Cressie 1985)."""348    from scipy.optimize import least_squares349    h = np.asarray(h, float); g = np.asarray(gamma, float); nb = np.asarray(n_paires, float)350    def res(p):351        gm = variogramme_exponentiel(h, *p); return np.sqrt(nb) * (g - gm) / np.maximum(gm, 1e-9)352    p0 = [0.5 * g[0], max(g[-1] - 0.5 * g[0], 1e-6), portee_init]353    r = least_squares(res, p0, bounds=([0, 0, 1e-3], [np.inf, np.inf, np.inf]))354    c0, c1, a = r.x355    return {"pepite": float(c0), "palier_partiel": float(c1), "portee": float(a), "portee_pratique": float(3 * a), "palier": float(c0 + c1)}356357def krigeage_ordinaire(coords_obs: np.ndarray, z: np.ndarray, coords_cibles: np.ndarray, params: dict, voisins: int | None = 40) -> dict:358    """Krigeage ordinaire avec covariance exponentielle C(h) = c₁ exp(−h/a) (+ c₀ en h = 0), voisinage des `voisins` points359    les plus proches (None = tous). Retourne la prédiction et la variance de krigeage à chaque cible."""360    co = np.asarray(coords_obs, float); z = np.asarray(z, float); ct = np.asarray(coords_cibles, float)361    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]362    cov = lambda h: c1 * np.exp(-h / a)363    dct = np.sqrt(((ct[:, None, :] - co[None, :, :]) ** 2).sum(-1))364    pred = np.zeros(len(ct)); var = np.zeros(len(ct))365    for i in range(len(ct)):366        idx = np.argsort(dct[i])[:voisins] if voisins else np.arange(len(co))367        sub = co[idx]; dd = _distances(sub); m = len(idx)368        K = np.zeros((m + 1, m + 1)); K[:m, :m] = cov(dd); K[np.arange(m), np.arange(m)] = c0 + c1369        K[m, :m] = 1.0; K[:m, m] = 1.0370        rhs = np.concatenate([cov(dct[i, idx]), [1.0]])371        sol = np.linalg.solve(K, rhs); w = sol[:m]; mu = sol[m]372        pred[i] = w @ z[idx]; var[i] = (c0 + c1) - w @ rhs[:m] - mu373    return {"prediction": pred, "variance": var}374375def krigeage_validation_croisee(coords: np.ndarray, z: np.ndarray, params: dict, voisins: int = 40) -> np.ndarray:376    """Prédiction leave-one-out de chaque observation par krigeage ordinaire de ses voisins (hors elle-même)."""377    c = np.asarray(coords, float); z = np.asarray(z, float); d = _distances(c); out = np.zeros(len(z))378    c0, c1, a = params["pepite"], params["palier_partiel"], params["portee"]; cov = lambda h: c1 * np.exp(-h / a)379    for i in range(len(z)):380        idx = np.argsort(d[i])[1:voisins + 1]; m = len(idx); dd = d[np.ix_(idx, idx)]381        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] = 1382        rhs = np.concatenate([cov(d[i, idx]), [1.0]]); sol = np.linalg.solve(K, rhs); out[i] = sol[:m] @ z[idx]383    return out384385def sem_esperance_conditionnelle(lam: float, W_aug: np.ndarray, u_obs: np.ndarray, idx_cibles) -> np.ndarray:386    """Espérance conditionnelle des perturbations u aux cibles sachant les u observés, dans un SEM de matrice augmentée W_aug387    (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."""388    n = W_aug.shape[0]; B = np.eye(n) - lam * W_aug; Q = B.T @ B389    idx_c = np.asarray(idx_cibles); idx_o = np.setdiff1d(np.arange(n), idx_c)390    Qcc = Q[np.ix_(idx_c, idx_c)]; Qco = Q[np.ix_(idx_c, idx_o)]391    return -np.linalg.solve(Qcc, Qco @ np.asarray(u_obs, float))392