spb/uqo-eval

1// UQO Éval — modélisation hédonique : qualité d'ajustement, tests de spécification, VIF.2import { gram, mean, median, solveSPD, sortedCopy, type Mat } from "./linalg";3import { chi2Sf, normalQuantile } from "./stats";45export interface FitMetrics {6  r2: number;7  adjR2: number;8  rmse: number; // en $9  mae: number; // en $10  mdape: number; // fraction11  within10: number;12  within20: number;13  aic: number;14  bic: number;15  sigma: number; // écart-type des résidus (échelle de la variable dépendante)16  smear: number; // facteur de Duan (log) ou 117  n: number;18  k: number;19}2021/** Prédictions en dollars (retransformation de Duan pour le log). */22export function toDollars(fitted: Float64Array, logModel: boolean, smear: number): Float64Array {23  const out = new Float64Array(fitted.length);24  for (let i = 0; i < fitted.length; i++) out[i] = logModel ? Math.exp(fitted[i]) * smear : fitted[i];25  return out;26}2728export function smearing(resid: Float64Array, logModel: boolean): number {29  if (!logModel) return 1;30  let s = 0;31  for (let i = 0; i < resid.length; i++) s += Math.exp(resid[i]);32  return s / resid.length;33}3435export function fitMetrics(y: Float64Array, fitted: Float64Array, resid: Float64Array, prices: Float64Array, logModel: boolean, k: number, smearOverride?: number): FitMetrics {36  const n = y.length;37  const smear = smearOverride ?? smearing(resid, logModel);38  const pred = toDollars(fitted, logModel, smear);39  const ym = mean(y);40  let sst = 0;41  let sse = 0;42  for (let i = 0; i < n; i++) {43    sst += (y[i] - ym) ** 2;44    sse += resid[i] ** 2;45  }46  const r2 = sst > 0 ? 1 - sse / sst : NaN;47  let se$ = 0;48  let ae$ = 0;49  const ape: number[] = [];50  let w10 = 0;51  let w20 = 0;52  for (let i = 0; i < n; i++) {53    const d = pred[i] - prices[i];54    se$ += d * d;55    ae$ += Math.abs(d);56    const a = Math.abs(d) / prices[i];57    ape.push(a);58    if (a <= 0.1) w10++;59    if (a <= 0.2) w20++;60  }61  return {62    r2,63    adjR2: 1 - ((1 - r2) * (n - 1)) / Math.max(1, n - k),64    rmse: Math.sqrt(se$ / n),65    mae: ae$ / n,66    mdape: median(ape),67    within10: w10 / n,68    within20: w20 / n,69    aic: n * Math.log(sse / n) + 2 * k,70    bic: n * Math.log(sse / n) + k * Math.log(n),71    sigma: Math.sqrt(sse / Math.max(1, n - k)),72    smear,73    n,74    k,75  };76}7778/** Breusch-Pagan (version Koenker, robuste à la non-normalité) : n·R² de e² sur X. */79export function breuschPagan(X: Mat, resid: Float64Array): { stat: number; df: number; p: number } {80  const n = X.n;81  const e2 = new Float64Array(n);82  for (let i = 0; i < n; i++) e2[i] = resid[i] * resid[i];83  const { G, c } = gram(X, e2);84  let tr = 0;85  for (let j = 0; j < X.k; j++) tr += G[j * X.k + j];86  const { x: b } = solveSPD(G, X.k, c, (tr / X.k) * 1e-10);87  const m = mean(e2);88  let sst = 0;89  let sse = 0;90  for (let i = 0; i < n; i++) {91    let f = 0;92    for (let j = 0; j < X.k; j++) f += X.a[i * X.k + j] * b[j];93    sst += (e2[i] - m) ** 2;94    sse += (e2[i] - f) ** 2;95  }96  const r2 = sst > 0 ? 1 - sse / sst : 0;97  const stat = n * r2;98  const df = X.k - 1;99  return { stat, df, p: chi2Sf(stat, df) };100}101102export function jarqueBera(resid: Float64Array): { stat: number; p: number; skew: number; kurt: number } {103  const n = resid.length;104  const m = mean(resid);105  let m2 = 0, m3 = 0, m4 = 0;106  for (let i = 0; i < n; i++) {107    const d = resid[i] - m;108    m2 += d * d;109    m3 += d * d * d;110    m4 += d * d * d * d;111  }112  m2 /= n; m3 /= n; m4 /= n;113  const skew = m3 / Math.pow(m2, 1.5);114  const kurt = m4 / (m2 * m2);115  const stat = (n / 6) * (skew * skew + ((kurt - 3) ** 2) / 4);116  return { stat, p: chi2Sf(stat, 2), skew, kurt };117}118119/** VIF des colonnes demandées : VIF_j = [(X'X)^-1]_jj · Σ(x_ij − x̄_j)². */120export function vif(X: Mat, XtXinv: Float64Array, cols: number[]): number[] {121  const { n, k, a } = X;122  return cols.map((j) => {123    let s = 0;124    for (let i = 0; i < n; i++) s += a[i * k + j];125    const m = s / n;126    let ss = 0;127    for (let i = 0; i < n; i++) ss += (a[i * k + j] - m) ** 2;128    return XtXinv[j * k + j] * ss;129  });130}131132/** Points du diagramme quantile-quantile (résidus standardisés). */133export function qqPoints(resid: Float64Array, sigma: number, maxPts = 300): { q: number; e: number }[] {134  const s = sortedCopy(resid);135  const n = s.length;136  const step = Math.max(1, Math.floor(n / maxPts));137  const out: { q: number; e: number }[] = [];138  for (let i = 0; i < n; i += step) {139    const p = (i + 0.5) / n;140    out.push({ q: normalQuantile(p), e: s[i] / sigma });141  }142  return out;143}144145export function histogram(v: Float64Array, bins = 30): { lo: number; hi: number; n: number }[] {146  const s = sortedCopy(v);147  const lo = s[Math.floor(s.length * 0.005)];148  const hi = s[Math.min(s.length - 1, Math.ceil(s.length * 0.995))];149  const w = (hi - lo) / bins || 1;150  const out = Array.from({ length: bins }, (_, i) => ({ lo: lo + i * w, hi: lo + (i + 1) * w, n: 0 }));151  for (let i = 0; i < v.length; i++) {152    const b = Math.min(bins - 1, Math.max(0, Math.floor((v[i] - lo) / w)));153    out[b].n++;154  }155  return out;156}157