spb/uqo-eval

1// UQO Éval — modélisation hédonique : voisinage spatial (k plus proches voisins), matrice W, Moran.2import { normalCdf } from "./stats";34/** Index spatial par grille pour des requêtes kNN rapides (coordonnées → km locaux). */5export class KnnIndex {6  private xs: Float64Array;7  private ys: Float64Array;8  private cell: number;9  private buckets = new Map<string, number[]>();10  readonly n: number;1112  constructor(lat: ArrayLike<number>, lng: ArrayLike<number>) {13    this.n = lat.length;14    let mlat = 0;15    for (let i = 0; i < this.n; i++) mlat += lat[i];16    mlat /= Math.max(1, this.n);17    const cos = Math.cos((mlat * Math.PI) / 180);18    this.cosLat = cos;19    this.xs = new Float64Array(this.n);20    this.ys = new Float64Array(this.n);21    let minx = Infinity, maxx = -Infinity, miny = Infinity, maxy = -Infinity;22    for (let i = 0; i < this.n; i++) {23      this.xs[i] = lng[i] * 111 * cos;24      this.ys[i] = lat[i] * 111;25      minx = Math.min(minx, this.xs[i]); maxx = Math.max(maxx, this.xs[i]);26      miny = Math.min(miny, this.ys[i]); maxy = Math.max(maxy, this.ys[i]);27    }28    // ~25 points par cellule29    const area = Math.max(1e-6, (maxx - minx) * (maxy - miny));30    this.cell = Math.max(0.05, Math.sqrt((area / Math.max(1, this.n)) * 25));31    for (let i = 0; i < this.n; i++) {32      const key = this.key(this.xs[i], this.ys[i]);33      const b = this.buckets.get(key);34      if (b) b.push(i);35      else this.buckets.set(key, [i]);36    }37  }38  private key(x: number, y: number) {39    return `${Math.floor(x / this.cell)}|${Math.floor(y / this.cell)}`;40  }41  private cosLat = 1;42  /** k voisins les plus proches d'un point (km locaux), en excluant éventuellement un index. */43  queryXY(x: number, y: number, k: number, exclude = -1): { idx: Int32Array; dist: Float64Array } {44    const cx = Math.floor(x / this.cell);45    const cy = Math.floor(y / this.cell);46    const best: { i: number; d: number }[] = [];47    let ring = 0;48    const maxRing = 400;49    while (ring <= maxRing) {50      // les cellules de l'anneau `ring`51      for (let dx = -ring; dx <= ring; dx++) {52        for (let dy = -ring; dy <= ring; dy++) {53          if (Math.max(Math.abs(dx), Math.abs(dy)) !== ring) continue;54          const b = this.buckets.get(`${cx + dx}|${cy + dy}`);55          if (!b) continue;56          for (const i of b) {57            if (i === exclude) continue;58            const d = Math.hypot(this.xs[i] - x, this.ys[i] - y);59            if (best.length < k) {60              best.push({ i, d });61              if (best.length === k) best.sort((a, b2) => a.d - b2.d);62            } else if (d < best[k - 1].d) {63              best[k - 1] = { i, d };64              best.sort((a, b2) => a.d - b2.d);65            }66          }67        }68      }69      // arrêt : k trouvés et l'anneau suivant ne peut pas contenir plus proche70      if (best.length >= k && ring * this.cell >= best[k - 1].d) break;71      if (best.length >= k && ring > 60) break;72      ring++;73    }74    best.sort((a, b2) => a.d - b2.d);75    const m = Math.min(k, best.length);76    const idx = new Int32Array(m);77    const dist = new Float64Array(m);78    for (let j = 0; j < m; j++) {79      idx[j] = best[j].i;80      dist[j] = best[j].d;81    }82    return { idx, dist };83  }84  queryLatLng(lat: number, lng: number, k: number, exclude = -1) {85    return this.queryXY(lng * 111 * this.cosLat, lat * 111, k, exclude);86  }87  /** Voisinage de chaque point de l'index (sans lui-même). */88  selfNeighbours(k: number): Int32Array[] {89    const out: Int32Array[] = new Array(this.n);90    for (let i = 0; i < this.n; i++) out[i] = this.queryXY(this.xs[i], this.ys[i], k, i).idx;91    return out;92  }93}9495/** W·v avec W = kNN ligne-standardisée (moyenne des voisins). */96export function lag(nb: Int32Array[], v: ArrayLike<number>): Float64Array {97  const out = new Float64Array(nb.length);98  for (let i = 0; i < nb.length; i++) {99    const ids = nb[i];100    if (!ids.length) continue;101    let s = 0;102    for (let j = 0; j < ids.length; j++) s += v[ids[j]];103    out[i] = s / ids.length;104  }105  return out;106}107108export interface Moran {109  I: number;110  expected: number;111  z: number;112  p: number;113  k: number;114}115116/** I de Moran des résidus sous W kNN ligne-standardisée (test de normalité asymptotique). */117export function moranI(nb: Int32Array[], e: Float64Array): Moran {118  const n = e.length;119  const k = nb[0]?.length ?? 0;120  if (n < 10 || !k) return { I: NaN, expected: NaN, z: NaN, p: NaN, k };121  let m = 0;122  for (let i = 0; i < n; i++) m += e[i];123  m /= n;124  const We = lag(nb, e);125  let num = 0;126  let den = 0;127  for (let i = 0; i < n; i++) {128    num += (e[i] - m) * (We[i] - m); // Σ_i (e_i−m)·Σ_j w_ij (e_j−m) (lag conserve la moyenne)129    den += (e[i] - m) ** 2;130  }131  const S0 = n; // lignes standardisées132  const I = (n / S0) * (num / den);133  // S1, S2134  const sets = nb.map((ids) => new Set(Array.from(ids)));135  const colSum = new Float64Array(n); // c_i = Σ_j w_ji136  for (let i = 0; i < n; i++) for (const j of nb[i]) colSum[j] += 1 / k;137  let S1 = 0;138  for (let i = 0; i < n; i++)139    for (const j of nb[i]) {140      const wji = sets[j].has(i) ? 1 / k : 0;141      S1 += (1 / k + wji) ** 2;142    }143  // chaque paire symétrique a été visitée deux fois (depuis i et depuis j) : S1 = Σ_visites − ½·Σ_sym144  let sym = 0;145  for (let i = 0; i < n; i++) for (const j of nb[i]) if (sets[j].has(i)) sym += (2 / k) ** 2;146  S1 = S1 - sym / 2;147  let S2 = 0;148  for (let i = 0; i < n; i++) S2 += (1 + colSum[i]) ** 2;149  const EI = -1 / (n - 1);150  const varI = (n * n * S1 - n * S2 + 3 * S0 * S0) / ((n * n - 1) * S0 * S0) - EI * EI;151  const z = varI > 0 ? (I - EI) / Math.sqrt(varI) : NaN;152  const p = Number.isFinite(z) ? 2 * (1 - normalCdf(Math.abs(z))) : NaN;153  return { I, expected: EI, z, p, k };154}155