feat: add empirical KDE/quantile utilities for real posterior draws

src/utils/kde.js

@@ -0,0 +1,73 @@
+/**
+ * Empirical density utilities for posterior draw arrays — Gaussian KDE with
+ * Silverman's rule-of-thumb bandwidth, plus a linear-interpolated quantile.
+ * Used in place of distributionApprox.js's analytic fitSkewedInterval/
+ * densityCurve approximation whenever real posterior draws are available.
+ */
+
+const SQRT_2PI = Math.sqrt(2 * Math.PI)
+const MIN_BANDWIDTH = 1e-3
+
+/**
+ * Linear-interpolated quantile (R type-7). Does not mutate `draws`.
+ * @param {number[]} draws
+ * @param {number} p - probability in [0, 1]
+ * @returns {number}
+ */
+export function quantile(draws, p) {
+  const sorted = [...draws].sort((a, b) => a - b)
+  const idx = p * (sorted.length - 1)
+  const lo = Math.floor(idx)
+  const hi = Math.ceil(idx)
+  if (lo === hi) return sorted[lo]
+  const frac = idx - lo
+  return sorted[lo] * (1 - frac) + sorted[hi] * frac
+}
+
+function standardDeviation(draws) {
+  const n = draws.length
+  const mean = draws.reduce((sum, d) => sum + d, 0) / n
+  const variance = draws.reduce((sum, d) => sum + (d - mean) ** 2, 0) / (n - 1)
+  return Math.sqrt(variance)
+}
+
+/**
+ * Silverman's rule-of-thumb bandwidth (robust variant using min(sd, IQR/1.34)),
+ * floored so a near-degenerate draw set (e.g. almost all zeros) never
+ * collapses the kernel to a spike.
+ * @param {number[]} draws
+ * @returns {number}
+ */
+export function silvermanBandwidth(draws) {
+  const n = draws.length
+  if (n < 2) return MIN_BANDWIDTH
+  const sd = standardDeviation(draws)
+  const iqr = quantile(draws, 0.75) - quantile(draws, 0.25)
+  const spread = Math.min(sd, iqr / 1.34) || sd || MIN_BANDWIDTH
+  return Math.max(0.9 * spread * Math.pow(n, -0.2), MIN_BANDWIDTH)
+}
+
+/**
+ * n evenly spaced {x, y} points of a Gaussian KDE over `draws` — same shape
+ * contract as distributionApprox.js's densityCurve, so chart code can switch
+ * between the two without changing its rendering path.
+ * @param {number[]} draws
+ * @param {{min?: number, max?: number, n?: number}} [options]
+ * @returns {Array<{x: number, y: number}>}
+ */
+export function kdeCurve(draws, { min = 0, max, n = 60 } = {}) {
+  const hi = max ?? Math.max(...draws) * 1.1
+  const h = silvermanBandwidth(draws)
+  const step = (hi - min) / (n - 1)
+  const points = []
+  for (let i = 0; i < n; i++) {
+    const x = min + step * i
+    let sum = 0
+    for (const d of draws) {
+      const z = (x - d) / h
+      sum += Math.exp(-0.5 * z * z) / SQRT_2PI
+    }
+    points.push({ x, y: sum / (draws.length * h) })
+  }
+  return points
+}