distribution (yearly casualties in car accidents)

percentile.js

const input = [
  [2005, 465],
  [2006, 439],
  [2007, 415],
  [2008, 433],
  [2009, 346],
  [2010, 375],
  [2011, 382],
  [2012, 290],
  [2013, 309],
  [2014, 319],
  [2015, 356],
  [2016, 377],
  [2017, 364],
  [2018, 316],
  [2019, 355],
  [2020, 305],
  [2021, 363],
  [2022, 348],
  [2023, 358],
  [2024, 436],
]

const data = input.map(([y, n]) => n)

// Fit parameters
const mean = data.reduce((a, b) => a + b, 0) / data.length
const variance = data.reduce((sum, x) => sum + Math.pow(x - mean, 2), 0) / data.length
const r = (mean * mean) / (variance - mean)
const p = r / (r + mean)

// Fixed logGamma - the input should be (n) not (n-1) for the standard gamma
function logGamma(z) {
  // For Gamma(z), not Gamma(z+1)
  if (z < 0.5) {
    return Math.log(Math.PI / Math.sin(Math.PI * z)) - logGamma(1 - z)
  }

  z -= 1
  const g = 7
  const c = [
    0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
    12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
  ]

  let x = c[0]
  for (let i = 1; i < g + 2; i++) {
    x += c[i] / (z + i)
  }

  const t = z + g + 0.5
  return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x)
}

// Sanity check
console.log(`logGamma(5) = ${logGamma(5)}, expected = ${Math.log(24)}`) // Gamma(5) = 4! = 24
console.log(`logGamma(10) = ${logGamma(10)}, expected = ${Math.log(362880)}`) // Gamma(10) = 9!

function negBinomPMF(k, r, p) {
  const logP = logGamma(k + r) - logGamma(k + 1) - logGamma(r) + r * Math.log(p) + k * Math.log(1 - p)
  return Math.exp(logP)
}

function negBinomCDF(k, r, p) {
  let sum = 0
  for (let i = 0; i <= k; i++) {
    sum += negBinomPMF(i, r, p)
  }
  return sum
}

function negBinomQuantile(prob, r, p) {
  let k = 0
  let cumulative = 0

  while (k < 10000) {
    cumulative += negBinomPMF(k, r, p)
    if (cumulative >= prob) return k
    k++
  }
  return k
}

function getConfidenceInterval(confidence, r, p) {
  const alpha = 1 - confidence
  return {
    lower: negBinomQuantile(alpha / 2, r, p),
    upper: negBinomQuantile(1 - alpha / 2, r, p),
    median: negBinomQuantile(0.5, r, p),
  }
}

const ci90 = getConfidenceInterval(0.9, r, p)
const ci95 = getConfidenceInterval(0.95, r, p)
const ci99 = getConfidenceInterval(0.99, r, p)

console.log(`Samples: ${data.length}`)
console.log(`Median: ${ci95.median}`)
console.log(`90% CI: [${ci90.lower}, ${ci90.upper}]`)
console.log(`95% CI: [${ci95.lower}, ${ci95.upper}]`)
console.log(`99% CI: [${ci99.lower}, ${ci99.upper}]`)

function getPercentile(k, r, p) {
  return negBinomCDF(k, r, p) * 100
}

const observed = 450
const percentile = getPercentile(observed, r, p)
console.log(`${observed} casualties is at the ${percentile.toFixed(1)}th percentile`)
console.log(`P(X > ${observed}) = ${(100 - percentile).toFixed(2)}%`)