fix: clamp KDE bandwidth to the plotting domain, not absolute units
Small districts produce rare-event count posteriors that are heavily zero-inflated, and the absolute 1e-3 floor / raw-sd fallback broke on both extremes of that data: - All-identical draws (e.g. 500 zeros, the norm for a group with a handful of students) collapsed to h = 1e-3, a near-delta spike of density ~399. Because SexRidgeColumn shares one maxPdf per column, that single spike flattened every other ridge in the column to sub-pixel height. Measured on a real district (0400315 AZ, unified_m4_mod): the three informative ridges rendered at 0.07-0.13px of a 47.56px row. - When IQR is 0 -- the normal case when most draws are 0 -- min(sd, iqr/1.34) was falsy and the rule fell all the way back to raw sd, which a few extreme draws inflate until the ridge is a flat line claiming maximal uncertainty. Bandwidth is now clamped to [domainWidth/50, domainWidth/6] and the robust rule degrades by picking the smallest *positive* spread estimate instead of discarding robustness entirely. The floor is slightly wider than one render step at n = 60, so a degenerate draw set resolves as a narrow bump; it does not bind on an ordinary posterior (spread wider than ~8% of the domain keeps its own Silverman bandwidth). On the district above the informative ridges now render at 5.66-5.70px, a ~60x improvement. Five new tests cover both failure modes; all five fail against the old formula.
This commit is contained in:
+45
-8
@@ -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++) {
|
||||
|
||||
+62
-2
@@ -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]', () => {
|
||||
|
||||
Reference in New Issue
Block a user