diff --git a/src/utils/kde.js b/src/utils/kde.js index 32b28e2..903a158 100644 --- a/src/utils/kde.js +++ b/src/utils/kde.js @@ -6,7 +6,20 @@ */ 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`. @@ -32,19 +45,40 @@ function standardDeviation(draws) { } /** - * 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. + * 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) { +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 MIN_BANDWIDTH + if (n < 2) return floor + 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) + // 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) } /** @@ -57,7 +91,10 @@ export function silvermanBandwidth(draws) { */ export function kdeCurve(draws, { min = 0, max, n = 60 } = {}) { const hi = max ?? Math.max(...draws) * 1.1 - const h = silvermanBandwidth(draws) + // 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++) { diff --git a/src/utils/kde.test.js b/src/utils/kde.test.js index c3c7f5e..e18f937 100644 --- a/src/utils/kde.test.js +++ b/src/utils/kde.test.js @@ -18,13 +18,73 @@ test('quantile: does not mutate its input array', () => { }) test('silvermanBandwidth: positive, finite floor for identical draws', () => { - const h = silvermanBandwidth([7, 7, 7, 7, 7]) + const h = silvermanBandwidth([7, 7, 7, 7, 7], 10) assert.ok(h > 0 && Number.isFinite(h)) + assert.equal(h, 10 / 50, 'degenerate spread falls back to the domain-relative floor') }) test('silvermanBandwidth: positive, finite floor for a single draw', () => { - const h = silvermanBandwidth([7]) + const h = silvermanBandwidth([7], 10) assert.ok(h > 0 && Number.isFinite(h)) + assert.equal(h, 10 / 50) +}) + +test('silvermanBandwidth: positive, finite floor when no domain is supplied', () => { + assert.ok(silvermanBandwidth([7, 7, 7, 7, 7]) > 0) + assert.ok(silvermanBandwidth([7]) > 0) +}) + +test('silvermanBandwidth: clamped to [domainWidth/50, domainWidth/6]', () => { + const domainWidth = 30 + // sd is inflated by 50 extreme draws; without the ceiling this is ~15.6. + const outlierDraws = [...Array(450).fill(0), ...Array(50).fill(200)] + assert.equal(silvermanBandwidth(outlierDraws, domainWidth), domainWidth / 6) + // Fully degenerate: sd = IQR = 0. + assert.equal(silvermanBandwidth(Array(500).fill(0), domainWidth), domainWidth / 50) +}) + +test('silvermanBandwidth: leaves an ordinary spread untouched by the clamp', () => { + const draws = Array.from({ length: 500 }, (_, i) => 3 + Math.sin(i) * 1.2 + (i % 7) * 0.15) + const h = silvermanBandwidth(draws, 10) + assert.ok(h > 10 / 50 && h < 10 / 6, `expected an unclamped bandwidth, got ${h}`) +}) + +test('kdeCurve: all-identical draws do not produce a delta-function spike', () => { + // 500 identical zeros is the common case for a small district's rare-event + // count posterior. The old absolute 1e-3 floor gave maxY ~399 here (peak + // ~3989x a uniform density over the same domain), which flattened every other + // ridge sharing the column's maxPdf to sub-pixel height. + const domainWidth = 10 + const curve = kdeCurve(Array(500).fill(0), { min: 0, max: domainWidth, n: 60 }) + const maxY = Math.max(...curve.map((p) => p.y)) + assert.ok(Number.isFinite(maxY) && maxY > 0) + assert.ok( + maxY * domainWidth < 25, + `peak density should stay within ~25x a uniform density over the domain, got ${maxY * domainWidth}x`, + ) +}) + +test('kdeCurve: zero-inflated draws keep their shape (not over-smoothed to flat)', () => { + const domainWidth = 10 + const draws = [...Array(450).fill(0), ...Array(50).fill(1)] + const curve = kdeCurve(draws, { min: 0, max: domainWidth, n: 60 }) + const ys = curve.map((p) => p.y) + const maxY = Math.max(...ys) + const meanY = ys.reduce((sum, y) => sum + y, 0) / ys.length + assert.ok(maxY / meanY > 3, `expected a peaked curve, got peak/mean ${maxY / meanY}`) +}) + +test('kdeCurve: outlier-inflated sd does not flatten the curve', () => { + // IQR is 0 here (most draws are 0), so the rule falls back to sd — which these + // 50 extreme draws inflate to ~60. Unclamped that gives h ~15.6 on a domain of + // 30, i.e. peak/mean ~1.6: a near-flat line claiming maximal uncertainty. + const domainWidth = 30 + const draws = [...Array(450).fill(0), ...Array(50).fill(200)] + const curve = kdeCurve(draws, { min: 0, max: domainWidth, n: 60 }) + const ys = curve.map((p) => p.y) + const maxY = Math.max(...ys) + const meanY = ys.reduce((sum, y) => sum + y, 0) / ys.length + assert.ok(maxY / meanY > 3, `expected a peaked curve, got peak/mean ${maxY / meanY}`) }) test('kdeCurve: returns n points spanning [min, max]', () => {