fix: compute the district total's interval as HPD, matching the API

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

@@ -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
+  const [lower, upper] = hpdBounds(totals, TOTAL_INTERVAL_MASS)
-  return {
+  return { lower, median: quantile(totals, 0.5), upper, nDraws: totals.length }
-    lower: quantile(totals, tail),
-    median: quantile(totals, 0.5),
-    upper: quantile(totals, 1 - tail),
-    nDraws: totals.length,
-  }
 }

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
+  assert.equal(iv.nDraws, 1000)
-  // exactly: p=0.025 over 0…1000 lands on 25 ± float noise.
+  assert.ok(iv.lower < iv.median && iv.median < iv.upper)
-  const near = (a, b) => Math.abs(a - b) < 1e-9
+  assert.ok(Math.abs(iv.median - 500) < 25, `median ${iv.median}`)
-  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}`)
+test('totalInterval: bounds are the narrowest window covering 95% of draws', () => {
-  assert.equal(iv.nDraws, 1001)
+  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', () => {