crdc-demo/src/utils/kde.js

/**
 * 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)
// Absolute last-resort floor, used only when no plotting domain is known.
const MIN_BANDWIDTH = 1e-3
// Bandwidth is clamped relative to the plotting domain, not to absolute units,
// because these draws are rate-per-1,000 values whose scale varies by orders of
// magnitude between districts. domainWidth/50 is slightly wider than one render
// step at the charts' n = 60 (step = domainWidth/59), so a degenerate draw set
// (e.g. 500 identical zeros, common for small districts) resolves as a narrow
// bump instead of a delta spike that flattens every other ridge sharing the
// column's maxPdf. It is also loose enough not to bind on an ordinary posterior:
// a spread wider than ~8% of the domain keeps its own Silverman bandwidth.
// domainWidth/6 stops a handful of extreme draws from inflating sd until the
// curve is a flat line.
const BANDWIDTH_FLOOR_DIVISOR = 50
const BANDWIDTH_CEILING_DIVISOR = 6

/**
 * 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 the smallest
 * *positive* spread estimate among sd and IQR/1.34), clamped to a fraction of
 * the plotting domain.
 *
 * Both ends of the clamp matter for the zero-inflated count posteriors small
 * districts produce. Without the floor, an all-identical draw set (sd = IQR = 0)
 * collapses to a delta-function spike. Without the ceiling, a group whose IQR is
 * 0 (the normal case when most draws are 0) falls back to raw sd, which a
 * handful of extreme draws inflates until the ridge is a featureless flat line.
 *
 * @param {number[]} draws
 * @param {number} [domainWidth] - width of the x-range the curve will be drawn
 *   over. Omit only when no domain is known; the clamp then degrades to the
 *   absolute MIN_BANDWIDTH floor.
 * @returns {number}
 */
export function silvermanBandwidth(draws, domainWidth = 0) {
  const width = domainWidth > 0 ? domainWidth : 0
  const floor = width ? width / BANDWIDTH_FLOOR_DIVISOR : MIN_BANDWIDTH
  const ceiling = width ? width / BANDWIDTH_CEILING_DIVISOR : Infinity

  const n = draws.length
  if (n < 2) return floor

  const sd = standardDeviation(draws)
  const iqr = quantile(draws, 0.75) - quantile(draws, 0.25)
  // Degrade gracefully: keep the robust rule when IQR is informative, use sd
  // when it isn't, and let the floor handle a fully degenerate draw set —
  // rather than treating a zero spread as "no estimate available".
  const candidates = [sd, iqr / 1.34].filter((v) => v > 0)
  const spread = candidates.length ? Math.min(...candidates) : 0

  const raw = 0.9 * spread * Math.pow(n, -0.2)
  return Math.min(Math.max(raw, floor), ceiling)
}

/**
 * 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
  // Guarded against a degenerate/inverted domain so the bandwidth clamp can
  // never be handed a negative width.
  const domainWidth = Math.max(hi - min, 0)
  const h = silvermanBandwidth(draws, domainWidth)
  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
}