From fb0beaafc86a57b2ff77d0c78cec8644c3611e20 Mon Sep 17 00:00:00 2001 From: Jared Knowles Date: Wed, 12 Aug 2026 09:34:53 -0400 Subject: [PATCH] fix: compute the district total's interval as HPD, matching the API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The API's count_lower/count_upper are highest-density bounds — hpd_bounds_sql() in crdc-arrests/R/summarize_draws.R — not equal-tailed quantiles, and the white paper's figures use the same. The draws-based total introduced in the previous commit used 2.5/97.5 quantiles, so the chart's two modes labelled "95% interval" meant two different things depending on whether the draws had been fetched. Measured on Clark County, the difference is small — 2 to 3 arrests on a band of roughly 50, about a pixel — but it is a distinction the page offers no explanation for, and the fallback mode's bounds are already HPD. wave equal-tailed HPD delta 2015-16 215-276 219-277 width -3 2017-18 175-230 173-226 width -2 2021-22 116-164 115-161 width -2 hpdBounds() takes the narrowest window covering the mass, contiguous as the R original is. Rendered values read back off the chart geometry now match an offline DuckDB computation exactly. --- src/utils/districtTotal.js | 48 ++++++++++++++++++---- src/utils/districtTotal.test.js | 71 ++++++++++++++++++++++++++++----- 2 files changed, 103 insertions(+), 16 deletions(-) diff --git a/src/utils/districtTotal.js b/src/utils/districtTotal.js index 62eacf9..05f1ad8 100644 --- a/src/utils/districtTotal.js +++ b/src/utils/districtTotal.js @@ -54,17 +54,51 @@ export function totalPerDraw(countsByGroup, nDraws) { return totals } +/** + * Narrowest interval containing `mass` of the values — the highest-density + * interval, matching what the API stores. + * + * The API's `count_lower`/`count_upper` are HPD bounds (`hpd_bounds_sql` in + * crdc-arrests/R/summarize_draws.R), not equal-tailed quantiles, and the white + * paper's figures use the same. Computing the total's interval the same way + * keeps one definition of "95% interval" on the page: the difference is only + * two or three arrests on a Clark County band of ~50, but the chart's fallback + * mode shows summed HPD bounds, and mixing conventions between the two modes + * would be a distinction with no explanation. + * + * Contiguous by construction, as the R original is. For a strongly bimodal + * posterior that is a simplification, but a count total summed across eight + * groups is unimodal in practice. + * + * @param {number[]} values + * @param {number} mass - e.g. 0.95 + * @returns {[number, number]} + */ +export function hpdBounds(values, mass) { + const sorted = [...values].sort((a, b) => a - b) + const n = sorted.length + const span = Math.ceil(mass * n) - 1 + if (span <= 0) return [sorted[0], sorted[0]] + if (span >= n - 1) return [sorted[0], sorted[n - 1]] + + let best = [sorted[0], sorted[span]] + let bestWidth = sorted[span] - sorted[0] + for (let i = 1; i + span < n; i++) { + const width = sorted[i + span] - sorted[i] + if (width < bestWidth) { + bestWidth = width + best = [sorted[i], sorted[i + span]] + } + } + return best +} + /** * @param {number[] | null | undefined} totals - district total per draw * @returns {{lower: number, median: number, upper: number, nDraws: number} | null} */ export function totalInterval(totals) { if (!totals?.length) return null - const tail = (1 - TOTAL_INTERVAL_MASS) / 2 - return { - lower: quantile(totals, tail), - median: quantile(totals, 0.5), - upper: quantile(totals, 1 - tail), - nDraws: totals.length, - } + const [lower, upper] = hpdBounds(totals, TOTAL_INTERVAL_MASS) + return { lower, median: quantile(totals, 0.5), upper, nDraws: totals.length } } diff --git a/src/utils/districtTotal.test.js b/src/utils/districtTotal.test.js index a63370f..8d7696f 100644 --- a/src/utils/districtTotal.test.js +++ b/src/utils/districtTotal.test.js @@ -1,6 +1,6 @@ import { test } from 'node:test' import assert from 'node:assert/strict' -import { TOTAL_INTERVAL_MASS, totalInterval, totalPerDraw } from './districtTotal.js' +import { TOTAL_INTERVAL_MASS, hpdBounds, totalInterval, totalPerDraw } from './districtTotal.js' // ——— totalPerDraw ——— @@ -40,15 +40,68 @@ test('totalPerDraw: does not mutate its input', () => { // ——— totalInterval ——— test('totalInterval: median and 95% bounds from the draw totals', () => { - const totals = Array.from({ length: 1001 }, (_, i) => i) // 0…1000 + // Symmetric and unimodal, so the HPD should sit close to the equal-tailed + // 25–975 without being required to equal it. + const totals = [] + for (let i = 0; i < 1000; i++) { + totals.push(500 + 120 * (Math.sin(i * 1.7) + Math.sin(i * 0.31) + Math.sin(i * 2.9)) / 3) + } const iv = totalInterval(totals) - // Linear-interpolated quantiles, so compare with a tolerance rather than - // exactly: p=0.025 over 0…1000 lands on 25 ± float noise. - const near = (a, b) => Math.abs(a - b) < 1e-9 - assert.ok(near(iv.median, 500), `median ${iv.median}`) - assert.ok(near(iv.lower, 25), `lower ${iv.lower}`) - assert.ok(near(iv.upper, 975), `upper ${iv.upper}`) - assert.equal(iv.nDraws, 1001) + assert.equal(iv.nDraws, 1000) + assert.ok(iv.lower < iv.median && iv.median < iv.upper) + assert.ok(Math.abs(iv.median - 500) < 25, `median ${iv.median}`) +}) + +test('totalInterval: bounds are the narrowest window covering 95% of draws', () => { + const totals = Array.from({ length: 1000 }, (_, i) => i) + const iv = totalInterval(totals) + const inside = totals.filter((t) => t >= iv.lower && t <= iv.upper).length + assert.ok(inside >= 950, `only ${inside} of 1000 draws inside`) + // A uniform has no denser region, so the narrowest window is ~95% of the range. + assert.ok(iv.upper - iv.lower <= 951, `width ${iv.upper - iv.lower}`) +}) + +test('hpdBounds: never wider than the equal-tailed interval', () => { + // Right-skewed, which is where the two definitions diverge most. + const skewed = Array.from({ length: 2000 }, (_, i) => Math.round(((i * 7919) % 1000) ** 1.6 / 1000)) + const [lo, hi] = hpdBounds(skewed, 0.95) + const sorted = [...skewed].sort((a, b) => a - b) + const etLo = sorted[Math.floor(0.025 * (sorted.length - 1))] + const etHi = sorted[Math.ceil(0.975 * (sorted.length - 1))] + assert.ok(hi - lo <= etHi - etLo, `hpd ${hi - lo} vs equal-tailed ${etHi - etLo}`) +}) + +test('hpdBounds: stays in the dense region when it holds enough mass', () => { + // 970 draws at 0-9 and 30 stragglers out at 500. The cluster alone covers 97%, + // so the narrowest 95% window fits inside it and the far tail is excluded. + const values = [ + ...Array.from({ length: 970 }, (_, i) => i % 10), + ...Array.from({ length: 30 }, () => 500), + ] + const [lo, hi] = hpdBounds(values, 0.95) + assert.equal(lo, 0) + assert.ok(hi < 500, `upper bound ${hi} should exclude the far cluster`) +}) + +test('hpdBounds: still reaches the tail when the cluster is too small', () => { + // The mirror case, and the honest one: 90% in the cluster cannot cover a 95% + // interval, so the bound must extend outward rather than under-covering. + const values = [ + ...Array.from({ length: 900 }, (_, i) => i % 10), + ...Array.from({ length: 100 }, () => 500), + ] + const [, hi] = hpdBounds(values, 0.95) + assert.equal(hi, 500) +}) + +test('hpdBounds: degenerate input collapses to a point', () => { + assert.deepEqual(hpdBounds(new Array(100).fill(7), 0.95), [7, 7]) +}) + +test('hpdBounds: does not mutate its input', () => { + const values = [5, 1, 3] + hpdBounds(values, 0.95) + assert.deepEqual(values, [5, 1, 3]) }) test('totalInterval: uses a 95% mass, matching the API convention', () => {