Incomplete gamma and beta

@ac-kit/math-analysisPortable

import { regularizedIncompleteGamma, regularizedIncompleteGammaUpper, regularizedIncompleteBeta, inverseRegularizedIncompleteGamma, inverseRegularizedIncompleteBeta } from "@ac-kit/math-analysis";

Where you meet them. Every time software reports a p-value it evaluates one of these. A chi-squared test goes through the gamma family; a t-test, an F-test and a binomial confidence interval all go through the beta family. If you have ever read p = 0.003 off a screen, this is the code that produced it.

Regularized means the normalized version, in [0, 1]. The word is in the name because libraries disagree about what the bare name means, and the difference is a factor of Γ(a) that no type system will catch for you.

The gamma family

regularizedIncompleteGamma() computes P(a, x) — the proportion of a Gamma(a, 1) distribution at or below x. It is the cumulative function of the gamma, chi-squared and Poisson distributions.

regularizedIncompleteGammaUpper() computes Q(a, x) = 1 − P(a, x), the proportion above x. In practice this is the p-value, and it is a separate function for the same reason erfc is: when the answer is small, subtracting it out of 1 has already thrown it away.

The regularized incomplete gamma family

P(a, x) is the chi-squared and Poisson cumulative distribution in disguise: it rises from 0 to 1 as x sweeps the domain, with the shape parameter a setting where the rise happens. Q = 1 − P is the survival function. Both names spell out 'regularized' because libraries disagree about whether the bare name means this ratio in [0, 1] or the unnormalised integral — and the argument order puts the parameter first and the variable last, so this family and the beta family read the same way.

const xs = linspace(0, upper, samples);
const evaluate =
  branch === "lower P(a, x)"
    ? regularizedIncompleteGamma
    : regularizedIncompleteGammaUpper;

const series = [];
for (let a = 1; a <= shapes; a++) {
  series.push({
    name: `a = ${a}`,
    values: xs.map((x) => evaluate(a, x)),
  });
}

const frame = curveFrame(xs, series, {
  x: "x",
  y: "probability",
  series: "shape",
});

The beta family

regularizedIncompleteBeta() computes Iₓ(a, b), the cumulative function of the Beta distribution — the distribution of a proportion. Through it come Student’s t, the F distribution and the binomial tail, which is why it appears far more often than its name suggests.

The regularized incomplete beta function

I(a, b; x) is the beta distribution's cumulative function, and through it Student's t, the F distribution and the binomial tail. It is evaluated by a continued fraction rather than by its defining integral, with the identity I(a, b; x) = 1 − I(b, a; 1 − x) used to pick whichever side converges quickly — the fraction is slow exactly where the reflected one is fast.

const xs = linspace(0, 1, samples);

const series = [];
for (let multiple = 1; multiple <= spread; multiple++) {
  const shape = b * multiple;
  series.push({
    name: `b = ${shape}`,
    values: xs.map((x) => regularizedIncompleteBeta(a, shape, x)),
  });
}

const frame = curveFrame(xs, series, {
  x: "x",
  y: "probability",
  series: "shape",
});

Running them backwards

A cumulative function answers “which percentile is this value at?”. Very often you want the other direction: “which value sits at the 95th percentile?” That is the quantile, and it is what turns a confidence level into a critical value.

inverseRegularizedIncompleteGamma() and inverseRegularizedIncompleteBeta() do that inversion. Neither has a closed form, so both refine an estimate until it stops moving — see convergence and iteration limits for how that termination is controlled.

An inverse is only as good as the round trip it survives: put a probability in, get a value out, feed that value back through the forward function, and see how far you have drifted.

EvidenceInverting a probability and evaluating it back agrees with the input to better than 1e-12.

Round-tripping the quantile functions

An inverse is only as good as the round trip it survives. Each point takes a probability p, asks the inverse for the x with that cumulative probability, then feeds x back through the forward function and measures how far the answer drifts from p. A few ulps across the whole unit interval is what correct looks like; the residual climbing near the ends is where the forward function's own conditioning gives out, not where the inverse is wrong.

// Endpoints excluded: p = 0 and p = 1 map to the domain's own boundaries.
const ps = linspace(0.001, 0.999, samples);

const residuals = ps.map((p) => {
  const x =
    family === "gamma"
      ? inverseRegularizedIncompleteGamma(a, p)
      : inverseRegularizedIncompleteBeta(a, b, p);
  const back =
    family === "gamma"
      ? regularizedIncompleteGamma(a, x)
      : regularizedIncompleteBeta(a, b, x);
  return Math.max(Math.abs(back - p), Number.MIN_VALUE);
});

const frame = curveFrame(
  ps,
  [{ name: "|P(inverse(p)) − p|", values: residuals }],
  { x: "p", y: "residual", series: "quantity" },
);

Incomplete gamma and beta

The gamma family

The beta family

Running them backwards

Further reading