Function discreteQuantile

src/compute-engine/numerics/distributions.ts:20–56

(
  cdf: (k: number) => number,
  prob: number,
  mean: number,
  stddev: number,
  kMin: number,
  kMax: number,
  deadline?: number
)

 * search seeded from the normal approximation. `kMax` may be `+∞` (Poisson).
 */
export function discreteQuantile(
  cdf: (k: number) => number,
  prob: number,
  mean: number,
  stddev: number,
  kMin: number,
  kMax: number,
  deadline?: number
): number {
  if (prob <= 0) return kMin;
  if (prob >= 1) return kMax;

  // A tiny tolerance makes the `CDF(k) ≥ p` boundary inclusive against
  // round-off: when `p` is itself a CDF value (the `Quantile(CDF(k)) = k`
  // identity) it may be recomputed a few ulps above the machine CDF used in
  // this search, which would otherwise overshoot by one.
  const target = prob - 1e-12;

  // Seed from the normal-approximation inverse-CDF, then correct.
  let k = Math.round(mean + Math.SQRT2 * erfInv(2 * prob - 1) * stddev);
  if (!Number.isFinite(k)) k = Math.round(mean);
  k = Math.max(k, kMin);
  if (Number.isFinite(kMax)) k = Math.min(k, kMax);

  let guard = 0;
  // Step up while the CDF is still below the target probability.
  while (k < kMax && cdf(k) < target) {
    k++;
    if ((++guard & 0x3ff) === 0) checkDeadline(deadline);
  }
  // Step down while the previous integer already reaches the target.
  while (k > kMin && cdf(k - 1) >= target) {
    k--;
    if ((++guard & 0x3ff) === 0) checkDeadline(deadline);
  }
  return k;
}

/** Quantile of Binomial(n, p) at probability `prob` (an integer in [0, n]). */
export function binomialQuantile(