The gamma and beta functions

@ac-kit/math-analysisPortable

import { gamma, logGamma, beta, logBeta } from "@ac-kit/math-analysis";

Where you meet it. Almost every continuous probability distribution has a normalizing constant — the number you divide by so the total probability comes to 1 — and for the chi-squared, gamma, Student’s t and F distributions, that constant is written with gamma functions. If you are computing a p-value, this is running underneath.

Γ(x)

gamma() is the function itself. Γ(n) = (n − 1)! at every positive integer — note the off-by-one, which is a historical accident everyone has to live with — and it interpolates smoothly between them.

At every non-positive integer it has a pole: the curve shoots off to infinity, and it goes to +∞ on one side and −∞ on the other. There is no sensible single answer there, so gamma returns NaN rather than picking one.

The gamma function and its logarithm

Γ extends the factorial to the reals, with Γ(n) = (n − 1)! at every positive integer. It has a pole at every non-positive integer, where the one-sided limits disagree in sign, so gamma answers NaN there rather than picking one. Drag the bounds across an integer to watch the two branches diverge.

const xs = linspace(lower, upper, samples);

const series = [];
if (scale !== "ln |Γ(x)|") {
  series.push({ name: "Γ(x)", values: xs.map((x) => clip(gamma(x))) });
}
if (scale !== "Γ(x)") {
  series.push({
    name: "ln |Γ(x)|",
    values: xs.map((x) => clip(logGamma(x))),
  });
}

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

When Γ stops fitting

Gamma grows faster than anything else in ordinary use. Γ(171) is about 7 × 10³⁰⁷; Γ(172) is past the top of a JavaScript number and comes back as Infinity.

logGamma() computes ln Γ(x) directly instead of computing Γ(x) and then taking its logarithm. That distinction is the whole point: it never forms the enormous intermediate value, so it keeps working for every argument you can represent.

EvidencelogGamma stays finite past x ≈ 171.6, where gamma overflows binary64 and Math.log(gamma(x)) becomes Infinity.

Where Γ stops fitting and ln Γ carries on

Both curves are ln Γ(x): one taken directly from logGamma, one recovered as Math.log(gamma(x)). They agree exactly until gamma(x) overflows to Infinity, at which point the recovered curve becomes Infinity too and stops being plottable. The direct one keeps going for every argument a caller can represent.

const xs = linspace(2, upper, samples);

const direct = xs.map((x) => logGamma(x));
// `Infinity` breaks the line, which is the point: past the overflow there is
// nothing left to draw.
const viaGamma = xs.map((x) => {
  const value = Math.log(gamma(x));
  return Number.isFinite(value) ? value : Number.NaN;
});

Working in logarithms is not a workaround — it is how these formulas are actually written in practice. A likelihood is a product of many small probabilities, and the log turns every one of those products into a sum that stays comfortably in range.

B(a, b)

beta() computes B(a, b) = Γ(a)Γ(b) / Γ(a + b).

Writing that ratio out by hand overflows for shape parameters far smaller than the ones whose ratio is still perfectly representable — gamma(200) is Infinity, but B(200, 200) is an ordinary small number. That is why it is a function rather than an expression.

The beta function is the normalizing constant of the Beta distribution, which describes a proportion: a conversion rate, a defect rate, a probability you are uncertain about. Through that distribution it also reaches the binomial, Student’s t and the F distribution.

B(a, b) underflows where ln B(a, b) does not

The beta function decays so fast that binary64 cannot hold it: B(a, b) reaches the underflow limit and becomes exactly zero, taking every digit with it. logBeta is not a convenience wrapper around the logarithm of beta(a, b) — it is computed through logGamma directly, so it keeps full relative precision across the whole range. Both curves are plotted as ln B so they are comparable, and the second one falls off a cliff.

const as = linspace(1, upper, samples);

const frame = curveFrame(
  as,
  [
    { name: "logBeta(a, b)", values: as.map((a) => logBeta(a, b)) },
    {
      name: "Math.log(beta(a, b))",
      values: as.map((a) => Math.log(beta(a, b))),
    },
  ],
  { x: "a", y: "logValue", series: "method" },
);

B decays fast in the other direction — it underflows to exactly zero from around a + b ≈ 1500. logBeta() is what to use past that point, and any time the result was going to be logged anyway.

Two naming rules

Every name in this package spells out regularized where it applies, because libraries disagree about whether the bare name means the normalized ratio in [0, 1] or the unnormalised integral. Guessing wrong is silent and off by a factor of Γ(a).

Arguments go parameters first, variable last — regularizedIncompleteGamma(a, x), regularizedIncompleteBeta(a, b, x) — so the gamma and beta families read the same way round.