TypeScript 90.1%
JavaScript 3.2%
Python 3.2%
CSS 1.8%
HTML 1.8%
1// UQO Éval — modélisation hédonique : estimateurs (MCO, Huber, LAD, ridge, LASSO, SAR-2SLS) et covariances.2import { cholesky, cholInverse, cholSolve, gram, hcat, matVec, solveSPD, takeCols, zeros, type Mat } from "./linalg";3import { lag, KnnIndex } from "./spatial";4import type { SeType } from "./spec";56export interface Fit {7 beta: Float64Array;8 fitted: Float64Array;9 resid: Float64Array;10 /** (X'WX)^-1 — « pain » du sandwich ; null pour les méthodes sans inférence. */11 bread: Float64Array | null;12 /** Matrice de design effectivement utilisée pour l'inférence (X ou X̂ en 2SLS). */13 Xinf: Mat | null;14 weights: Float64Array | null;15 sigma2: number;16 df: number;17 k: number;18 extra: Record<string, number>;19}2021function residuals(X: Mat, y: Float64Array, beta: Float64Array) {22 const fitted = matVec(X, beta);23 const resid = new Float64Array(X.n);24 for (let i = 0; i < X.n; i++) resid[i] = y[i] - fitted[i];25 return { fitted, resid };26}2728/** MCO (ou MCP si poids). */29export function ols(X: Mat, y: Float64Array, weights?: Float64Array): Fit {30 const { G, c } = gram(X, y, weights);31 const { x: beta, inv } = solveSPD(G, X.k, c);32 const { fitted, resid } = residuals(X, y, beta);33 let ss = 0;34 for (let i = 0; i < X.n; i++) ss += (weights ? weights[i] : 1) * resid[i] * resid[i];35 const df = X.n - X.k;36 return { beta, fitted, resid, bread: inv, Xinf: X, weights: weights ?? null, sigma2: ss / Math.max(1, df), df, k: X.k, extra: {} };37}3839/** Covariance des coefficients : classique, HC1 ou grappes. */40export function covariance(fit: Fit, type: SeType, clusters?: Int32Array, nClusters = 0): Float64Array | null {41 if (!fit.bread || !fit.Xinf) return null;42 const X = fit.Xinf;43 const { n, k, a } = X;44 const B = fit.bread;45 const w = fit.weights;46 if (type === "classic") {47 const V = new Float64Array(k * k);48 for (let i = 0; i < k * k; i++) V[i] = fit.sigma2 * B[i];49 return V;50 }51 // viande : Σ_g s_g s_g' avec s_g = Σ_{i∈g} w_i e_i x_i (HC1 : chaque obs = une grappe)52 const M = new Float64Array(k * k);53 const useCl = type === "cluster" && clusters && nClusters > 1;54 const nG = useCl ? nClusters : n;55 const S = useCl ? new Float64Array(nG * k) : null;56 for (let i = 0; i < n; i++) {57 const s = (w ? w[i] : 1) * fit.resid[i];58 const off = i * k;59 if (useCl && S) {60 const g = clusters![i] * k;61 for (let j = 0; j < k; j++) S[g + j] += s * a[off + j];62 } else {63 for (let j = 0; j < k; j++) {64 const sj = s * a[off + j];65 if (sj === 0) continue;66 for (let l = 0; l < k; l++) M[j * k + l] += sj * s * a[off + l];67 }68 }69 }70 if (useCl && S) for (let g = 0; g < nG; g++) for (let j = 0; j < k; j++) for (let l = 0; l < k; l++) M[j * k + l] += S[g * k + j] * S[g * k + l];71 // V = B M B × correction petits échantillons72 const BM = new Float64Array(k * k);73 for (let i = 0; i < k; i++) for (let j = 0; j < k; j++) {74 let s = 0;75 for (let l = 0; l < k; l++) s += B[i * k + l] * M[l * k + j];76 BM[i * k + j] = s;77 }78 const V = new Float64Array(k * k);79 const corr = useCl ? (nG / (nG - 1)) * ((n - 1) / (n - k)) : n / (n - k);80 for (let i = 0; i < k; i++) for (let j = 0; j < k; j++) {81 let s = 0;82 for (let l = 0; l < k; l++) s += BM[i * k + l] * B[l * k + j];83 V[i * k + j] = s * corr;84 }85 return V;86}8788const mad = (e: Float64Array) => {89 const s = Float64Array.from(e, Math.abs).sort();90 return s[Math.floor(s.length / 2)] / 0.6745;91};9293/** Régression robuste de Huber (IRLS, c = 1,345). */94export function huber(X: Mat, y: Float64Array, iters = 30): Fit {95 let fit = ols(X, y);96 let w = new Float64Array(X.n).fill(1);97 let downweighted = 0;98 for (let it = 0; it < iters; it++) {99 const s = Math.max(1e-9, mad(fit.resid));100 const c = 1.345 * s;101 const w2 = new Float64Array(X.n);102 downweighted = 0;103 for (let i = 0; i < X.n; i++) {104 const ae = Math.abs(fit.resid[i]);105 w2[i] = ae <= c ? 1 : c / ae;106 if (w2[i] < 1) downweighted++;107 }108 const next = ols(X, y, w2);109 let delta = 0;110 for (let j = 0; j < X.k; j++) delta = Math.max(delta, Math.abs(next.beta[j] - fit.beta[j]) / (1e-9 + Math.abs(fit.beta[j])));111 fit = next;112 w = w2;113 if (delta < 1e-6) break;114 }115 fit.extra = { scale: mad(fit.resid), downweighted, downweightedPct: (100 * downweighted) / X.n };116 fit.weights = w;117 return fit;118}119120/** Régression médiane (LAD) par IRLS — coefficients seulement. */121export function lad(X: Mat, y: Float64Array, iters = 50): Fit {122 let fit = ols(X, y);123 let sy = 0;124 for (let i = 0; i < X.n; i++) sy += Math.abs(y[i]);125 const eps = 1e-6 * (sy / X.n);126 for (let it = 0; it < iters; it++) {127 const w = new Float64Array(X.n);128 for (let i = 0; i < X.n; i++) w[i] = 1 / Math.max(Math.abs(fit.resid[i]), eps);129 const next = ols(X, y, w);130 let delta = 0;131 for (let j = 0; j < X.k; j++) delta = Math.max(delta, Math.abs(next.beta[j] - fit.beta[j]) / (1e-9 + Math.abs(fit.beta[j])));132 fit = next;133 if (delta < 1e-7) break;134 }135 let sae = 0;136 for (let i = 0; i < X.n; i++) sae += Math.abs(fit.resid[i]);137 return { ...fit, bread: null, Xinf: null, weights: null, extra: { mae: sae / X.n } };138}139140/* ------------------------------ pénalisation ------------------------------ */141interface Std {142 mean: Float64Array;143 sd: Float64Array;144 ymean: number;145}146function standardize(X: Mat, y: Float64Array): { Z: Mat; yc: Float64Array; std: Std } {147 const { n, k, a } = X;148 const mean = new Float64Array(k);149 const sd = new Float64Array(k);150 for (let j = 1; j < k; j++) {151 let s = 0;152 for (let i = 0; i < n; i++) s += a[i * k + j];153 mean[j] = s / n;154 let v = 0;155 for (let i = 0; i < n; i++) v += (a[i * k + j] - mean[j]) ** 2;156 sd[j] = Math.sqrt(v / n) || 1;157 }158 const Z = zeros(n, k - 1);159 for (let i = 0; i < n; i++) for (let j = 1; j < k; j++) Z.a[i * (k - 1) + j - 1] = (a[i * k + j] - mean[j]) / sd[j];160 let ym = 0;161 for (let i = 0; i < n; i++) ym += y[i];162 ym /= n;163 const yc = new Float64Array(n);164 for (let i = 0; i < n; i++) yc[i] = y[i] - ym;165 return { Z, yc, std: { mean, sd, ymean: ym } };166}167168function destandardize(b: Float64Array, std: Std, k: number): Float64Array {169 const beta = new Float64Array(k);170 let b0 = std.ymean;171 for (let j = 1; j < k; j++) {172 beta[j] = b[j - 1] / std.sd[j];173 b0 -= beta[j] * std.mean[j];174 }175 beta[0] = b0;176 return beta;177}178179/** Ridge sur données standardisées : (G/n + λI)^-1 c/n. */180function ridgeStd(G: Float64Array, c: Float64Array, p: number, n: number, lambda: number): Float64Array {181 const A = new Float64Array(p * p);182 const cc = new Float64Array(p);183 for (let i = 0; i < p * p; i++) A[i] = G[i] / n;184 for (let j = 0; j < p; j++) {185 A[j * p + j] += lambda;186 cc[j] = c[j] / n;187 }188 const L = cholesky(A, p, 1e-10);189 if (!L) throw new Error("ridge : matrice singulière");190 return cholSolve(L, p, cc);191}192193/** LASSO par descente de coordonnées sur la matrice de Gram (données standardisées). */194function lassoStd(G: Float64Array, c: Float64Array, p: number, n: number, lambda: number, b0?: Float64Array, iters = 200): Float64Array {195 const b = b0 ? Float64Array.from(b0) : new Float64Array(p);196 const soft = (z: number, l: number) => (z > l ? z - l : z < -l ? z + l : 0);197 for (let it = 0; it < iters; it++) {198 let maxd = 0;199 for (let j = 0; j < p; j++) {200 let z = c[j] / n;201 const row = j * p;202 for (let l = 0; l < p; l++) if (l !== j && b[l] !== 0) z -= (G[row + l] / n) * b[l];203 const nb = soft(z, lambda) / (G[row + j] / n || 1);204 maxd = Math.max(maxd, Math.abs(nb - b[j]));205 b[j] = nb;206 }207 if (maxd < 1e-7) break;208 }209 return b;210}211212export interface PenalizedResult {213 fit: Fit;214 lambda: number;215 grid: { lambda: number; cvErr: number }[];216 selected: number[]; // indices de colonnes (X) non nulles (LASSO)217 postFit?: Fit; // MCO post-LASSO sur les colonnes retenues218}219220export function penalized(X: Mat, y: Float64Array, method: "ridge" | "lasso", folds = 5): PenalizedResult {221 const { Z, yc, std } = standardize(X, y);222 const p = Z.k;223 const n = Z.n;224 const full = gram(Z, yc);225 // grille λ226 let lmax = 0;227 for (let j = 0; j < p; j++) lmax = Math.max(lmax, Math.abs(full.c[j]) / n);228 const grid: number[] = [];229 const nl = 24;230 if (method === "lasso") for (let i = 0; i < nl; i++) grid.push(lmax * Math.pow(1e-3, i / (nl - 1)));231 else for (let i = 0; i < nl; i++) grid.push(10 * Math.pow(1e-5, i / (nl - 1)));232 // validation croisée par blocs (l'échantillon est déjà mélangé)233 const cvErr = new Float64Array(grid.length);234 const foldOf = (i: number) => i % folds;235 for (let f = 0; f < folds; f++) {236 const tr: number[] = [];237 const te: number[] = [];238 for (let i = 0; i < n; i++) (foldOf(i) === f ? te : tr).push(i);239 const Ztr = zeros(tr.length, p);240 const ytr = new Float64Array(tr.length);241 for (let r = 0; r < tr.length; r++) {242 Ztr.a.set(Z.a.subarray(tr[r] * p, tr[r] * p + p), r * p);243 ytr[r] = yc[tr[r]];244 }245 const g = gram(Ztr, ytr);246 let warm: Float64Array | undefined;247 for (let li = 0; li < grid.length; li++) {248 const b = method === "lasso" ? lassoStd(g.G, g.c, p, tr.length, grid[li], warm) : ridgeStd(g.G, g.c, p, tr.length, grid[li]);249 warm = b;250 let se = 0;251 for (const i of te) {252 let pred = 0;253 for (let j = 0; j < p; j++) pred += Z.a[i * p + j] * b[j];254 se += (yc[i] - pred) ** 2;255 }256 cvErr[li] += se / n;257 }258 }259 let best = 0;260 for (let li = 1; li < grid.length; li++) if (cvErr[li] < cvErr[best]) best = li;261 const lambda = grid[best];262 const bStd = method === "lasso" ? lassoStd(full.G, full.c, p, n, lambda, undefined, 500) : ridgeStd(full.G, full.c, p, n, lambda);263 const beta = destandardize(bStd, std, X.k);264 const { fitted, resid } = residuals(X, y, beta);265 let ss = 0;266 for (let i = 0; i < n; i++) ss += resid[i] * resid[i];267 const selected = [0];268 for (let j = 1; j < X.k; j++) if (Math.abs(beta[j]) > 1e-12) selected.push(j);269 const dfEff = method === "lasso" ? selected.length : X.k;270 const fit: Fit = { beta, fitted, resid, bread: null, Xinf: null, weights: null, sigma2: ss / Math.max(1, n - dfEff), df: n - dfEff, k: X.k, extra: { lambda, nonzero: selected.length - 1 } };271 const res: PenalizedResult = { fit, lambda, grid: grid.map((l, i) => ({ lambda: l, cvErr: cvErr[i] })), selected };272 if (method === "lasso" && selected.length > 1) res.postFit = ols(takeCols(X, selected), y);273 return res;274}275276/* --------------------------------- SAR 2SLS --------------------------------- */277export interface SarResult {278 fit: Fit; // beta = [ρ, β...] ; Xinf = X̂ = [Ŵy, X]279 rho: number;280 Wy: Float64Array;281 nb: Int32Array[];282 index: KnnIndex;283}284285/** y = ρWy + Xβ + ε, instruments Z = [X, WX, W²X] (Kelejian & Prucha, 1998). */286export function sar(X: Mat, y: Float64Array, lat: Float64Array, lng: Float64Array, k = 8): SarResult {287 const index = new KnnIndex(lat, lng);288 const nb = index.selfNeighbours(k);289 const Wy = lag(nb, y);290 const { n, k: p } = X;291 // WX et W²X (hors constante, hors colonnes quasi-constantes)292 const cols: number[] = [];293 for (let j = 1; j < p; j++) cols.push(j);294 const Xnc = takeCols(X, cols);295 const WX = zeros(n, Xnc.k);296 const W2X = zeros(n, Xnc.k);297 for (let j = 0; j < Xnc.k; j++) {298 const col = new Float64Array(n);299 for (let i = 0; i < n; i++) col[i] = Xnc.a[i * Xnc.k + j];300 const w1 = lag(nb, col);301 const w2 = lag(nb, w1);302 for (let i = 0; i < n; i++) {303 WX.a[i * Xnc.k + j] = w1[i];304 W2X.a[i * Xnc.k + j] = w2[i];305 }306 }307 const Zm = hcat(hcat(X, WX), W2X);308 // 1re étape : Ŵy = Z (Z'Z)^-1 Z' Wy (ridge infime pour la colinéarité des instruments)309 const gz = gram(Zm, Wy);310 let tr = 0;311 for (let j = 0; j < Zm.k; j++) tr += gz.G[j * Zm.k + j];312 const { x: gamma } = solveSPD(gz.G, Zm.k, gz.c, (tr / Zm.k) * 1e-8);313 const WyHat = matVec(Zm, gamma);314 // 2e étape : régresser y sur X̂ = [Ŵy, X]315 const WyM: Mat = { n, k: 1, a: WyHat };316 const Xhat = hcat(WyM, X);317 const g2 = gram(Xhat, y);318 const { x: beta, inv } = solveSPD(g2.G, Xhat.k, g2.c);319 // résidus structurels avec le Wy observé320 const Xfull = hcat({ n, k: 1, a: Wy }, X);321 const { fitted, resid } = residuals(Xfull, y, beta);322 let ss = 0;323 for (let i = 0; i < n; i++) ss += resid[i] * resid[i];324 const df = n - Xhat.k;325 const fit: Fit = { beta, fitted, resid, bread: inv, Xinf: Xhat, weights: null, sigma2: ss / Math.max(1, df), df, k: Xhat.k, extra: { rho: beta[0], knn: k } };326 return { fit, rho: beta[0], Wy, nb, index };327}328329/** Prédiction SAR hors échantillon : Wy des k voisins d'entraînement. */330export function sarPredict(res: SarResult, Xte: Mat, latTe: Float64Array, lngTe: Float64Array, yTrain: Float64Array, beta: Float64Array): Float64Array {331 const out = new Float64Array(Xte.n);332 for (let i = 0; i < Xte.n; i++) {333 const q = res.index.queryLatLng(latTe[i], lngTe[i], res.nb[0]?.length ?? 8);334 let wy = 0;335 for (let j = 0; j < q.idx.length; j++) wy += yTrain[q.idx[j]];336 wy = q.idx.length ? wy / q.idx.length : 0;337 let s = beta[0] * wy;338 for (let j = 0; j < Xte.k; j++) s += Xte.a[i * Xte.k + j] * beta[j + 1];339 out[i] = s;340 }341 return out;342}343344export { cholInverse };345