From 91e21feb74c9bece2866b9673255b76462f1732a Mon Sep 17 00:00:00 2001
From: Jared Knowles <jared@civilytics.com>
Date: Tue, 11 Aug 2026 08:55:06 -0400
Subject: [PATCH] feat: add empirical KDE/quantile utilities for real posterior
 draws

---
 package.json          |  1 +
 src/utils/kde.js      | 73 +++++++++++++++++++++++++++++++++++++++++++
 src/utils/kde.test.js | 47 ++++++++++++++++++++++++++++
 3 files changed, 121 insertions(+)
 create mode 100644 src/utils/kde.js
 create mode 100644 src/utils/kde.test.js

diff --git a/package.json b/package.json
index e2e4c70..b61930d 100644
--- a/package.json
+++ b/package.json
@@ -6,6 +6,7 @@
     "dev": "vite",
     "build": "vite build",
     "preview": "vite preview",
+    "test": "node --test 'src/**/*.test.js'",
     "lint": "eslint src/ --ext .js,.jsx,.ts,.tsx",
     "format": "prettier --write \"src/**/*.{js,jsx,css}\""
   },
diff --git a/src/utils/kde.js b/src/utils/kde.js
new file mode 100644
index 0000000..32b28e2
--- /dev/null
+++ b/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
+}
diff --git a/src/utils/kde.test.js b/src/utils/kde.test.js
new file mode 100644
index 0000000..c3c7f5e
--- /dev/null
+++ b/src/utils/kde.test.js
@@ -0,0 +1,47 @@
+import { test } from 'node:test'
+import assert from 'node:assert/strict'
+import { quantile, silvermanBandwidth, kdeCurve } from './kde.js'
+
+test('quantile: median of an odd-length array', () => {
+  assert.equal(quantile([3, 1, 2], 0.5), 2)
+})
+
+test('quantile: linear interpolation between two ranks', () => {
+  // sorted: [10, 20, 30, 40] — p=0.25 -> index 0.75 -> interpolate 10..20
+  assert.equal(quantile([40, 10, 30, 20], 0.25), 17.5)
+})
+
+test('quantile: does not mutate its input array', () => {
+  const input = [5, 3, 4, 1, 2]
+  quantile(input, 0.5)
+  assert.deepEqual(input, [5, 3, 4, 1, 2])
+})
+
+test('silvermanBandwidth: positive, finite floor for identical draws', () => {
+  const h = silvermanBandwidth([7, 7, 7, 7, 7])
+  assert.ok(h > 0 && Number.isFinite(h))
+})
+
+test('silvermanBandwidth: positive, finite floor for a single draw', () => {
+  const h = silvermanBandwidth([7])
+  assert.ok(h > 0 && Number.isFinite(h))
+})
+
+test('kdeCurve: returns n points spanning [min, max]', () => {
+  const draws = [1, 2, 2, 3, 4, 5, 5, 5, 6, 8]
+  const curve = kdeCurve(draws, { min: 0, max: 10, n: 60 })
+  assert.equal(curve.length, 60)
+  assert.equal(curve[0].x, 0)
+  assert.ok(Math.abs(curve[curve.length - 1].x - 10) < 1e-9)
+})
+
+test('kdeCurve: density integrates to ~1 over a wide domain (trapezoidal check)', () => {
+  const draws = [1, 2, 2, 3, 4, 5, 5, 5, 6, 8]
+  const curve = kdeCurve(draws, { min: -20, max: 30, n: 2000 })
+  let area = 0
+  for (let i = 1; i < curve.length; i++) {
+    const dx = curve[i].x - curve[i - 1].x
+    area += (dx * (curve[i].y + curve[i - 1].y)) / 2
+  }
+  assert.ok(Math.abs(area - 1) < 0.01, `expected area ~1, got ${area}`)
+})