Prelude · Statistics

Distributions

Densities, CDFs, quantiles and samplers for twenty probability distributions, with the special functions (ln Γ, incomplete gamma and beta, erf) beneath them.

open Prelude.DistributionsDistributions.hb →

This module provides the probability distributions used across statistics and simulation: the normal, uniform, Cauchy, exponential, Laplace, logistic, Gumbel, Pareto, Weibull, Rayleigh, log-normal, χ², gamma, beta, Student’s t and F for continuous data, and Bernoulli, geometric, Poisson and binomial for counts. For each you get a density (for a discrete distribution, the probability mass), a cumulative distribution function and, for most, a quantile function. Beneath them sit the special functions they are built from: the log-gamma function, the regularised incomplete gamma and beta functions, and the error function.

Two functions carry most of the weight. The regularised lower incomplete gamma \(P(a, x)\) gives erf, χ², gamma and Poisson; the regularised incomplete beta \(I_x(a, b)\) gives Student’s t, F, beta and binomial. Both are evaluated by a fixed number of iterations — a series where it converges and a continued fraction where it does not. erf is defined here rather than being an intrinsic, so that every executor runs exactly the same arithmetic and returns the same digits.

Conventions. Parameters come first and the point of evaluation last: normalCdf(m, s, x), normalQuantile(m, s, p). Every argument may be a tensor and broadcasts, so a vector of x gives a vector of densities; the statements below are written for real numbers and apply entrywise. A quantile takes a probability p in (0, 1). Nothing here generates random numbers: to sample, pass uniform draws in [0, 1) to a quantile function (normalSample is exactly that), which turns them into draws from the distribution. Parameter names follow the source — for example the gamma distribution takes a shape and a scale (not a rate), while the exponential takes a rate; each entry says which.

Limitations. Choices between regions are made with 0/1 masks, not branches, so both sides of every choice are evaluated (on arguments clamped to their own domain so the discarded side cannot overflow). Iteration counts are fixed, which gives predictable cost and good accuracy for moderate parameters but no convergence test. Upper-tail probabilities obtained as one minus a CDF lose relative accuracy once they fall below about \(10^{-16}\). Log-gamma is valid for positive arguments.

open Prelude.Distributions

Guards

The masks used throughout this module evaluate both alternatives and multiply the unwanted one by zero; since \(0 \cdot \infty\) is NaN, not zero, each alternative is protected from division by zero, logarithms of zero, and arguments outside the region where it converges. These helpers are not meant to be called directly.

Throughout this page, \([\,\cdot\,]\) is the Iverson bracket (\([x \ge a] = 1\) if true, \(0\) otherwise), and \(v_+ = \operatorname{positive}(v) = \max(v, 10^{-300})\) is the floored argument a logarithm or a power is taken of.

Building blocks

• guard(v) — Let \(v \in \mathbb{R}\). Then \(\operatorname{guard}(v) = v + 10^{-300}\,[\lvert v \rvert < 10^{-300}]\). Nudges a value whose magnitude is below 1e-300 away from zero, so it can be a denominator.
• positive(v) — Let \(v \in \mathbb{R}\). Then \(\operatorname{positive}(v) = v_+ = \max(v, 10^{-300})\). Clamps a value to at least 1e-300, so its logarithm is finite.
• clipLow(a, x) — Let \(a, x \in \mathbb{R}\). Then \(\operatorname{clipLow}(a, x) = \min(x, a + 1)\). Clamps x to at most a + 1, the region where the incomplete-gamma series converges.
• clipHigh(a, x) — Let \(a, x \in \mathbb{R}\). Then \(\operatorname{clipHigh}(a, x) = \max(x, a + 1)\). Clamps x to at least a + 1, the region where the incomplete-gamma continued fraction converges.

The log-gamma function

The gamma function \(\Gamma(a)\) grows too fast to compute directly, so its logarithm is used everywhere. It is computed with the Lanczos approximation (\(g = 7\), nine coefficients): a fixed rational sum with no iteration, accurate to about fifteen digits over the positive half-line.

lgamma

lgamma(a) · a argument, positive scalar or tensor → \(\ln \Gamma(a)\) same shape

Let \(a > 0\). Then

\[\operatorname{lgamma}(a) = \tfrac12 \ln(2\pi) + \left(a - \tfrac12\right) \ln(a + 6.5) - (a + 6.5) + \ln \operatorname{lanczos}(a - 1) \approx \ln \Gamma(a),\]

where \(\tfrac12 \ln(2\pi)\) is the constant 0.9189385332046727 and lanczos is the rational sum below.

The natural logarithm of the gamma function. For a whole number, \(\ln \Gamma(n+1) = \ln n!\).

Intuition \(\Gamma(a)\) extends the factorial to real arguments (\(\Gamma(n+1) = n!\)), and it overflows a double once \(a\) passes about 171. Its logarithm grows only like \(a \ln a\), so products and ratios of gamma functions become sums and differences that stay finite.

Use it when you need a factorial, a binomial coefficient or a normalising constant with large arguments: \(\ln\binom{n}{k}\) is lgamma(n + 1) − lgamma(k + 1) − lgamma(n − k + 1), exactly as binomialPdf does it. Exponentiate only at the end. Arguments must be positive: there is no reflection formula, so zero or negative input is not supported.

combinatoricsnumerical-stabilitylikelihood-evaluation

def lgamma(a)    ≔ 0.9189385332046727 + (a − 0.5) · ln(a + 6.5) − (a + 6.5) + ln(lanczos(a − 1))

lbeta

lbeta(a, b) · a first argument, positive · b second argument, positive → \(\ln B(a, b)\) same shape

Let \(a > 0\) and \(b > 0\). Then

\[\operatorname{lbeta}(a, b) = \operatorname{lgamma}(a) + \operatorname{lgamma}(b) - \operatorname{lgamma}(a + b) \approx \ln B(a, b).\]

The logarithm of the beta function, \(\ln B(a, b) = \ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a + b)\).

Intuition \(B(a, b)\) is the area under \(x^{a-1}(1-x)^{b-1}\) on [0, 1]: the constant that turns that curve into a probability density. For whole numbers it is \((a-1)!\,(b-1)!/(a+b-1)!\), a continuous cousin of the reciprocal binomial coefficient.

Use it when you normalise a beta density, compute a beta-binomial likelihood, or do Bayesian updates of a proportion: the marginal likelihood of \(k\) successes in \(n\) trials under a Beta\((a, b)\) prior contains lbeta(a + k, b + n − k) − lbeta(a, b). Stay in logs. betaPdf, fPdf and betaI already call it for you.

bayesian-updatinglikelihood-evaluationnumerical-stability

def lbeta(a, b)  ≔ lgamma(a) + lgamma(b) − lgamma(a + b)

Building blocks

• lanczos(z) — Let \(z > -1\). Then \(\operatorname{lanczos}(z) = c_0 + \sum_{k=1}^{8} \frac{c_k}{z + k}\), where \(c_0 = 0.99999999999980993\), \(c_1 = 676.5203681218851\), \(c_2 = -1259.1392167224028\), \(c_3 = 771.32342877765313\), \(c_4 = -176.61502916214059\), \(c_5 = 12.507343278686905\), \(c_6 = -0.13857109526572012\), \(c_7 = 9.9843695780195716 \times 10^{-6}\) and \(c_8 = 1.5056327351493116 \times 10^{-7}\). The Lanczos rational sum with g = 7; used by lgamma.

The regularised incomplete gamma function

The regularised lower incomplete gamma function is

\[P(a, x) = \frac{1}{\Gamma(a)} \int_0^x t^{a-1} e^{-t}\, dt\]

and \(Q(a, x) = 1 - P(a, x)\) is its upper counterpart. For \(x < a + 1\) it is summed as a power series, each term the previous times \(x/(a+n)\) (201 terms: the first and 200 more); for \(x \ge a + 1\) it is evaluated as a continued fraction by the modified Lentz method (200 turns).

lowerGamma

lowerGamma(a, x) · a shape, positive · x upper limit, non-negative → \(P(a, x)\) in [0, 1] same shape

Let \(a > 0\) and \(x \ge 0\). Then

\[\operatorname{lowerGamma}(a, x) = [x < a + 1]\; \operatorname{gammaSeries}\big(a, \min(x, a + 1)\big) + [x \ge a + 1]\; \Big(1 - \operatorname{gammaFraction}\big(a, \max(x, a + 1)\big)\Big) \approx P(a, x),\]

where, away from the guard thresholds,

\[\operatorname{gammaSeries}(a, x) = \frac{e^{-x} x^{a}}{\Gamma(a)} \sum_{n=0}^{200} \frac{x^{n}}{a(a+1)\cdots(a+n)}, \qquad \operatorname{gammaFraction}(a, x) = \frac{e^{-x} x^{a}}{\Gamma(a)} \cdot \cfrac{1}{x + 1 - a - \cfrac{1 \cdot (1 - a)}{x + 3 - a - \cfrac{2 (2 - a)}{x + 5 - a - \cdots}}}\]

with the fraction cut after 200 turns (see the building blocks for the exact recurrences).

The regularised lower incomplete gamma function, choosing the series or the continued fraction by the mask \(x < a + 1\). It is the CDF of a gamma distribution with shape a and unit scale.

Intuition For whole \(a\), \(P(a, x)\) is the probability that a unit-rate Poisson process has produced at least \(a\) events by time \(x\), i.e. the CDF of the wait for the \(a\)-th event. For any \(a > 0\) it is the CDF of a gamma variable with shape \(a\) and unit scale, and χ², gamma, Poisson and erf are all this function after a change of variable.

Use it when you need a gamma-family CDF the module does not name, or a building block for your own distribution. Prefer gammaCdf, chiCdf, poissonCdf or erf when one fits, since they handle scaling and supports. The iteration counts are fixed (201 series terms, 200 fraction turns) with no convergence test, so check very large shapes against a reference value.

waiting-timestail-probabilitiesprobabilistic-modelling

def lowerGamma(a, x) ≔ (x < a + 1) · gammaSeries(a, clipLow(a, x)) + (x ≥ a + 1) · (1 − gammaFraction(a, clipHigh(a, x)))

upperGamma

upperGamma(a, x) · a shape, positive · x lower limit, non-negative → \(Q(a, x)\) in [0, 1] same shape

Let \(a > 0\) and \(x \ge 0\). Then

\[\operatorname{upperGamma}(a, x) = 1 - \operatorname{lowerGamma}(a, x) \approx Q(a, x).\]

The regularised upper incomplete gamma function, \(1 - P(a, x)\). Because it is formed by subtraction, very small values lose relative accuracy.

Intuition \(Q(a, x)\) asks the complementary question: the probability that the \(a\)-th event has not happened by time \(x\). For whole \(a\) it equals the probability that a Poisson count with mean \(x\) is at most \(a - 1\), which is how poissonCdf uses it.

Use it when you want an upper-tail (survival) probability of a gamma or χ² variable, or a Poisson CDF. It is computed as \(1 - P\), so its absolute error is around \(10^{-16}\): a result of \(10^{-12}\) has only a few correct digits, and smaller tails come out as 0. It is fine for ordinary p-values, not for rare-event probabilities.

tail-probabilitiescount-datahypothesis-testing

def upperGamma(a, x) ≔ 1 − lowerGamma(a, x)

Building blocks

• gsTerm(a, x, n, t) — Let \(a, x, n, t \in \mathbb{R}\). Then \(\operatorname{gsTerm}(a, x, n, t) = t\,x / \operatorname{guard}(a + n)\). The next series term: the previous term t times x / (a + n).
• gsLoop(k, a, x, n, t, s) — Let \(k \in \mathbb{N}\) and \(a, x, n, t, s \in \mathbb{R}\), and write \(t' = \operatorname{gsTerm}(a, x, n + 1, t)\). Then \(\operatorname{gsLoop}(0, a, x, n, t, s) = s\) and \(\operatorname{gsLoop}(k, a, x, n, t, s) = \operatorname{gsLoop}(k - 1, a, x, n + 1, t', s + t')\). Adds k further series terms, carrying the term and the running sum.
• gammaSeries(a, x) — Let \(a > 0\) and \(x \in \mathbb{R}\). Let \(t_0 = 1/\operatorname{guard}(a)\) and \(t_n = t_{n-1}\,x / \operatorname{guard}(a + n)\) for \(n = 1, \dots, 200\). Then \(\operatorname{gammaSeries}(a, x) = \Big(\sum_{n=0}^{200} t_n\Big) \exp\big(-x + a \ln x_+ - \operatorname{lgamma}(a)\big)\). P(a, x) by the power series (201 terms) times its exponential prefactor; used where x < a + 1.
• cfA(a, i) — Let \(a, i \in \mathbb{R}\). Then \(\operatorname{cfA}(a, i) = -i\,(i - a)\). The i-th partial numerator of the continued fraction.
• cfB(b) — Let \(b \in \mathbb{R}\). Then \(\operatorname{cfB}(b) = b + 2\). The next partial denominator: the previous one plus 2.
• cfD(a, i, b, d) — Let \(a, i, b, d \in \mathbb{R}\). Then \(\operatorname{cfD}(a, i, b, d) = 1 / \operatorname{guard}\big(-i(i - a)\,d + b + 2\big)\). The Lentz D update, carried already inverted.
• cfC(a, i, b, c) — Let \(a, i, b, c \in \mathbb{R}\). Then \(\operatorname{cfC}(a, i, b, c) = \operatorname{guard}\big(b + 2 - i(i - a) / \operatorname{guard}(c)\big)\). The Lentz C update.
• cfLoop(k, a, i, b, c, d, h) — Let \(k \in \mathbb{N}\) and \(a, i, b, c, d, h \in \mathbb{R}\), and write \(c' = \operatorname{cfC}(a, i + 1, b, c)\) and \(d' = \operatorname{cfD}(a, i + 1, b, d)\). Then \(\operatorname{cfLoop}(0, a, i, b, c, d, h) = h\) and \(\operatorname{cfLoop}(k, a, i, b, c, d, h) = \operatorname{cfLoop}(k - 1, a, i + 1, b + 2, c', d', h\,d'\,c')\). Runs k Lentz turns, multiplying the running value h by C·D each turn.
• gammaFraction(a, x) — Let \(a > 0\) and \(x \in \mathbb{R}\). Let \(b_0 = x + 1 - a\), \(c_0 = 10^{300}\), \(d_0 = h_0 = 1/\operatorname{guard}(b_0)\), and for \(i = 1, \dots, 200\) let \(b_i = b_{i-1} + 2\), \(d_i = 1/\operatorname{guard}(-i(i - a)\,d_{i-1} + b_i)\), \(c_i = \operatorname{guard}(b_i - i(i - a)/\operatorname{guard}(c_{i-1}))\) and \(h_i = h_{i-1}\,d_i\,c_i\). Then \(\operatorname{gammaFraction}(a, x) = h_{200} \exp\big(-x + a \ln x_+ - \operatorname{lgamma}(a)\big)\). Q(a, x) by the continued fraction (200 turns) times its exponential prefactor; used where x ≥ a + 1.

The error function

The error function is an incomplete gamma function in disguise: \(\operatorname{erf}(x) = \operatorname{sign}(x)\, P(\tfrac12, x^2)\).

erf

erf(x) · x argument scalar or tensor → \(\operatorname{erf}(x)\) in (−1, 1) same shape

Let \(x \in \mathbb{R}\). Then

\[\operatorname{erf}(x) = \operatorname{sign}(x)\; \operatorname{lowerGamma}\!\left(\tfrac12,\; x^2\right) \approx \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2}\, dt .\]

The error function, computed from the lower incomplete gamma function. It is odd: \(\operatorname{erf}(-x) = -\operatorname{erf}(x)\).

Intuition \(\operatorname{erf}(x)\) is the probability that a normal variable with mean 0 and variance \(\tfrac12\) lands in \([-x, x]\). Rescaled, \(\operatorname{erf}(z/\sqrt2)\) is the chance a standard normal lies within \(z\) standard deviations of its mean: about 0.683 at \(z = 1\) and 0.954 at \(z = 2\).

Use it when you need a symmetric “within \(k\) sigma” probability, a Gaussian integral from physics (diffusion or heat-conduction profiles), or the same digits on every executor. For a one-sided probability with a mean and a standard deviation, normalCdf is this function already rescaled.

measurement-errortail-probabilitiesuncertainty-quantification

def erf(x)  ≔ sign(x) · lowerGamma(0.5, x · x)

erfc

erfc(x) · x argument scalar or tensor → \(1 - \operatorname{erf}(x)\) in (0, 2) same shape

Let \(x \in \mathbb{R}\). Then

\[\operatorname{erfc}(x) = 1 - \operatorname{erf}(x).\]

The complementary error function. Because it is formed by subtraction, for large positive x it loses digits steadily and returns exactly 0 beyond about 6.

Intuition The mass outside \([-x, x]\) for the same variance-\(\tfrac12\) normal: \(\operatorname{erfc}(z/\sqrt2)\) is the probability that a standard normal falls more than \(z\) standard deviations from its mean, in either direction.

Use it when you want a two-sided normal tail, such as a two-sided z-test p-value, \(\operatorname{erfc}(\lvert z \rvert/\sqrt2)\), or a bit-error rate in a communications model. Because it is computed as \(1 - \operatorname{erf}\), it is unreliable for very small tails and is exactly 0 beyond about 6, so do not use it for far-tail or rare-event probabilities.

hypothesis-testingtail-probabilitieserror-analysis

def erfc(x) ≔ 1 − erf(x)

The normal distribution

The normal distribution with mean m and standard deviation s has density \(\frac{1}{s\sqrt{2\pi}} e^{-(x-m)^2 / 2s^2}\). Its quantile uses Wichura’s algorithm AS 241 (PPND16): three rational approximations — one for the centre, \(\lvert p - \tfrac12 \rvert \le 0.425\), and two for the tails — accurate to about sixteen digits with no iteration.

normalPdf

normalPdf(m, s, x) · m mean · s standard deviation, positive · x point → density broadcast shape

Let \(m \in \mathbb{R}\) be the mean, \(s > 0\) the standard deviation and \(x \in \mathbb{R}\) the point. Then

\[\operatorname{normalPdf}(m, s, x) = \frac{1}{s\sqrt{2\pi}} \exp\!\left(-\frac{(x - m)^2}{2 s^2}\right),\]

with \(\sqrt{2\pi}\) written as the constant 2.5066282746310002.

The normal probability density at x.

Intuition The bell curve is what you get when many small independent effects add up (the central limit theorem), which makes it the default model for measurement noise and for averages. \(m\) sets the centre and \(s\) the width, with about 68% of the mass within one \(s\) of \(m\). The density is a height, not a probability: it exceeds 1 when \(s\) is small.

Use it when you evaluate the likelihood of data under a Gaussian model, weight points in a Gaussian kernel or mixture, or draw the curve. For the probability of a range use normalCdf. When you combine many observations, add log-densities written directly as \(-\tfrac12\big((x-m)/s\big)^2 - \ln s\) plus a constant, rather than multiplying densities that underflow to 0.

likelihood-evaluationprobabilistic-modellingmeasurement-error

def normalPdf(m, s, x) ≔ exp(0 − (x − m) · (x − m) / (2 · s · s)) / (s · 2.5066282746310002)

normalCdf

normalCdf(m, s, x) · m mean · s standard deviation, positive · x point → probability broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{normalCdf}(m, s, x) = \Phi\!\left(\frac{x - m}{s}\right) = \frac12\left(1 + \operatorname{erf}\frac{x - m}{s\sqrt2}\right).\]

The probability that a normal variable is at most x.

Intuition The area under the bell curve to the left of \(x\): the fraction of the population at or below that value. Any interval probability is a difference of two of these, \(P(a < X \le b) = F(b) - F(a)\).

Use it when you need the share of a normal population inside a limit (a specification limit, a pass mark, a service threshold), a one-sided p-value for a z statistic, or a percentile rank. Upper tails are \(1 - F\) and lose accuracy below about \(10^{-16}\), so extreme-tail risks come out as 0. To go from a probability back to a value, use normalQuantile.

tail-probabilitieshypothesis-testingpercentiles

def normalCdf(m, s, x) ≔ 0.5 · (1 + erf((x − m) / (s · 1.4142135623730951)))

standardQuantile

standardQuantile(p) · p probability in (0, 1) scalar or tensor → standard normal quantile same shape

Let \(p \in [0, 1]\) and \(q = p - \tfrac12\). Then

\[\operatorname{standardQuantile}(p) = [\,\lvert q \rvert \le 0.425\,]\; \operatorname{midQ}(p) + [\,\lvert q \rvert > 0.425\,]\; \operatorname{tailQ}(p) \approx \Phi^{-1}(p),\]

where midQ and tailQ are the central and tail rational approximations of AS 241 given in the building blocks.

The inverse of the standard normal CDF, \(\Phi^{-1}(p)\), by AS 241. At exactly 0 or 1 it returns a large finite value rather than an infinity.

Intuition Given a probability, the z-score that cuts it off from below. It is the ruler that converts “the top 2.5%” into “1.96 standard deviations”, with no mean or scale attached.

Use it when you need a z critical value (\(z_{1-\alpha/2}\) for a two-sided interval or test), normal scores for a Q-Q plot or a rank-based transform, or a starting point for another quantile (the module seeds studentQuantile, chiQuantile and gammaQuantile from it). It is a fixed rational formula, so it is cheap. Keep \(p\) strictly inside (0, 1): at the ends it returns a large finite number, not an infinity.

confidence-intervalshypothesis-testingpercentiles

def standardQuantile(p) ≔ (abs(qOf(p)) ≤ 0.425) · midQ(p) + (abs(qOf(p)) > 0.425) · tailQ(p)

normalQuantile

normalQuantile(m, s, p) · m mean · s standard deviation, positive · p probability in (0, 1) → quantile broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(0 < p < 1\). Then

\[\operatorname{normalQuantile}(m, s, p) = m + s\, \operatorname{standardQuantile}(p) \approx m + s\,\Phi^{-1}(p).\]

The value below which a normal variable falls with probability p.

Intuition The standard z-score placed on your scale: shifted to the mean \(m\) and stretched by \(s\). It answers “below which value does a fraction \(p\) of the population fall?”

Use it when you set a threshold from a target probability: a tolerance or control limit, a reorder point for a service level, or a normal value-at-risk on a loss distribution. The bounds of a central \(1 - \alpha\) interval are the quantiles at \(\alpha/2\) and \(1 - \alpha/2\). To turn a batch of uniform draws into samples, normalSample is the same call under the name that says so.

confidence-intervalsrisk-measurementpercentiles

def normalQuantile(m, s, p) ≔ m + s · standardQuantile(p)

normalSample

normalSample(m, s, u) · m mean · s standard deviation, positive · u uniform draws in [0, 1) tensor → normal draws shape of u

Let \(m \in \mathbb{R}\), \(s > 0\), and let \(u \in [0, 1)^n\) be independent uniform draws. Then, for \(i = 1, \dots, n\),

\[\operatorname{normalSample}(m, s, u)_i = \operatorname{normalQuantile}(m, s, u_i) = m + s\, \operatorname{standardQuantile}(u_i).\]

Turns uniform draws into normal draws by the inverse-CDF method. Supply the uniforms yourself (for example from a seeded generator in C#); the same uniforms always give the same normals.

Intuition Inverse-transform sampling: if \(U\) is uniform on [0, 1), then \(F^{-1}(U)\) has CDF \(F\). Each uniform is read as “a random percentile” and mapped to the value at that percentile, with no rejection loop and no hidden random state.

Use it when you simulate from a normal model (Monte Carlo risk, noise injection, perturbation studies) and want reproducible results: seed the uniform generator in C#, and the same uniforms always give the same normals. Stratified or quasi-random uniforms (Latin hypercube, Sobol points) pass straight through with their structure intact. A uniform of exactly 0 maps to a large negative finite value rather than minus infinity.

simulationrandom-variate-generationuncertainty-quantification

def normalSample(m, s, u) ≔ normalQuantile(m, s, u)

Building blocks

• qOf(p) — Let \(p \in \mathbb{R}\). Then \(\operatorname{qOf}(p) = q = p - \tfrac12\). The centred probability p − 0.5.
• midR(p) — Let \(p \in \mathbb{R}\) and \(q = p - \tfrac12\). Then \(\operatorname{midR}(p) = 0.180625 - q^2\). The variable of the central rational approximation.
• midNum(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{midNum}(r) = \sum_{k=0}^{7} a_k r^k\) (evaluated in Horner form), where \(a_0 = 3.387132872796366608\), \(a_1 = 133.14166789178437745\), \(a_2 = 1971.5909503065514427\), \(a_3 = 13731.693765509461125\), \(a_4 = 45921.953931549871457\), \(a_5 = 67265.770927008700853\), \(a_6 = 33430.575583588128105\) and \(a_7 = 2509.0809287301226727\). Numerator polynomial of the central approximation.
• midDen(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{midDen}(r) = \sum_{k=0}^{7} b_k r^k\) (Horner form), where \(b_0 = 1\), \(b_1 = 42.313330701600911252\), \(b_2 = 687.1870074920579083\), \(b_3 = 5394.1960214247511077\), \(b_4 = 21213.794301586595867\), \(b_5 = 39307.89580009271061\), \(b_6 = 28729.085735721942674\) and \(b_7 = 5226.495278852854561\). Denominator polynomial of the central approximation.
• midQ(p) — Let \(p \in \mathbb{R}\), \(q = p - \tfrac12\) and \(r = 0.180625 - q^2\). Then \(\operatorname{midQ}(p) = q\, \operatorname{midNum}(r) / \operatorname{midDen}(r)\). The standard normal quantile for the central region.
• tailR(p) — Let \(p \in \mathbb{R}\). Then \(\operatorname{tailR}(p) = \rho = \sqrt{-\ln\big(\min(p, 1 - p)\big)_+}\). The tail variable, the square root of minus the log of the smaller tail probability.
• nearR(p) — Let \(p \in \mathbb{R}\) and \(\rho = \operatorname{tailR}(p)\). Then \(\operatorname{nearR}(p) = \rho - 1.6\). The tail variable shifted for the near-tail approximation.
• nearNum(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{nearNum}(r) = \sum_{k=0}^{7} c_k r^k\) (Horner form), where \(c_0 = 1.42343711074968357734\), \(c_1 = 4.6303378461565452959\), \(c_2 = 5.7694972214606914055\), \(c_3 = 3.64784832476320460504\), \(c_4 = 1.27045825245236838258\), \(c_5 = 0.24178072517745061177\), \(c_6 = 0.0227238449892691845833\) and \(c_7 = 7.7454501427834140764 \times 10^{-4}\). Numerator polynomial of the near-tail approximation.
• nearDen(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{nearDen}(r) = \sum_{k=0}^{7} d_k r^k\) (Horner form), where \(d_0 = 1\), \(d_1 = 2.05319162663775882187\), \(d_2 = 1.6763848301838038494\), \(d_3 = 0.689767334985100004550\), \(d_4 = 0.148103976427480074590\), \(d_5 = 0.0151986665636164571966\), \(d_6 = 5.47593808499534494600 \times 10^{-4}\) and \(d_7 = 1.05075007164441684324 \times 10^{-9}\). Denominator polynomial of the near-tail approximation.
• farR(p) — Let \(p \in \mathbb{R}\) and \(\rho = \operatorname{tailR}(p)\). Then \(\operatorname{farR}(p) = \rho - 5\). The tail variable shifted for the far-tail approximation.
• farNum(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{farNum}(r) = \sum_{k=0}^{7} e_k r^k\) (Horner form), where \(e_0 = 6.65790464350110377720\), \(e_1 = 5.46378491116411436990\), \(e_2 = 1.78482653991729133580\), \(e_3 = 0.296560571828504891230\), \(e_4 = 0.0265321895265761230930\), \(e_5 = 0.00124266094738807843860\), \(e_6 = 2.71155556874348757815 \times 10^{-5}\) and \(e_7 = 2.01033439929228813265 \times 10^{-7}\). Numerator polynomial of the far-tail approximation.
• farDen(r) — Let \(r \in \mathbb{R}\). Then \(\operatorname{farDen}(r) = \sum_{k=0}^{7} f_k r^k\) (Horner form), where \(f_0 = 1\), \(f_1 = 0.599832206555887937690\), \(f_2 = 0.136929880922735805310\), \(f_3 = 0.0148753612908506148525\), \(f_4 = 7.86869131145613259100 \times 10^{-4}\), \(f_5 = 1.84631831751005468180 \times 10^{-5}\), \(f_6 = 1.4215117583164458887 \times 10^{-7}\) and \(f_7 = 2.04426310338993978564 \times 10^{-15}\). Denominator polynomial of the far-tail approximation.
• tailQ(p) — Let \(p \in \mathbb{R}\), \(q = p - \tfrac12\) and \(\rho = \operatorname{tailR}(p)\). Then \(\operatorname{tailQ}(p) = \operatorname{sign}(q)\Big([\rho \le 5]\,\frac{\operatorname{nearNum}(\rho - 1.6)}{\operatorname{nearDen}(\rho - 1.6)} + [\rho > 5]\,\frac{\operatorname{farNum}(\rho - 5)}{\operatorname{farDen}(\rho - 5)}\Big)\). The standard normal quantile for both tails, signed by the side of 0.5.

Closed-form distributions

These distributions have densities, CDFs and quantiles in closed form, so they need no series at all. Each family has three functions — …Pdf, …Cdf and …Quantile — with the same leading parameters. Supports are enforced with masks: outside the support the density and CDF are 0 (or the CDF is 1 above it), with no error. Every quantile, applied to uniform draws, samples the distribution. As in the Guards section, \(v_+ = \max(v, 10^{-300})\).

uniformPdf

uniformPdf(a, b, x) · a lower bound · b upper bound, greater than a · x point → density broadcast shape

Let \(a < b\) be the bounds and \(x \in \mathbb{R}\). Then

\[\operatorname{uniformPdf}(a, b, x) = \frac{[x \ge a]\,[x \le b]}{b - a}.\]

\(1/(b - a)\) for \(a \le x \le b\), and 0 elsewhere.

Intuition “Anywhere in \([a, b]\), nothing preferred”: a flat density whose height \(1/(b-a)\) makes the area one. It is the model you reach for when all you know is the bounds.

Use it when you model a quantity known only to lie in a range (rounding error, a tolerance band, an arrival somewhere in a known window), write a flat prior, or evaluate the likelihood of bounded data. The density is 0 outside the bounds, so a log-likelihood there is minus infinity. When values near the middle are likelier, a beta rescaled to the interval (betaPdf) is the usual next step.

probabilistic-modellinglikelihood-evaluationmeasurement-error

def uniformPdf(a, b, x)      ≔ (x ≥ a) · (x ≤ b) / (b − a)

uniformCdf

uniformCdf(a, b, x) · a lower bound · b upper bound · x point → probability broadcast shape

Let \(a < b\) and \(x \in \mathbb{R}\). Then

\[\operatorname{uniformCdf}(a, b, x) = \min\!\left(\max\!\left(\frac{x - a}{b - a},\, 0\right),\, 1\right).\]

\((x - a)/(b - a)\) clamped to [0, 1].

Intuition The fraction of the interval already covered at \(x\): a straight ramp from 0 at \(a\) to 1 at \(b\), flat on either side.

Use it when you need the probability that a uniformly distributed quantity is below a threshold (a timer that fires somewhere in a window, a random offset), or you map a value in \([a, b]\) to a percentile in [0, 1] before passing it to another distribution’s quantile. Values outside \([a, b]\) are clamped to 0 or 1, not rejected.

tail-probabilitiespercentiles

def uniformCdf(a, b, x)      ≔ min(max((x − a) / (b − a), 0), 1)

uniformQuantile

uniformQuantile(a, b, p) · a lower bound · b upper bound · p probability → quantile broadcast shape

Let \(a < b\) and \(p \in [0, 1]\). Then

\[\operatorname{uniformQuantile}(a, b, p) = a + p\,(b - a).\]

It rescales draws in [0, 1) to [a, b).

Intuition A linear rescaling: the point a fraction \(p\) of the way from \(a\) to \(b\). Applied to standard uniform draws, it gives uniform draws on any interval.

Use it when you draw random values in a range (initial weights, jitter, random search over a hyperparameter interval) from the uniforms your C# generator supplies. It does not clamp: \(p\) outside [0, 1] gives a value outside \([a, b]\). For a range searched on a log scale, such as a learning rate, apply it to the logarithms of the bounds and exponentiate the result.

random-variate-generationsimulationhyperparameter-choice

def uniformQuantile(a, b, p) ≔ a + p · (b − a)

cauchyPdf

cauchyPdf(x0, g, x) · x0 location (median) · g scale, positive · x point → density broadcast shape

Let \(x_0 \in \mathbb{R}\) (the argument x0) be the location, \(\gamma > 0\) (the argument g) the scale and \(x \in \mathbb{R}\). Then

\[\operatorname{cauchyPdf}(x_0, \gamma, x) = \frac{\gamma}{\pi\,\big(\gamma^2 + (x - x_0)^2\big)}.\]

The Cauchy distribution has no mean or variance, which makes it a standard test case for robust statistics.

Intuition The ratio of two independent standard normals is Cauchy, and so is the point where a ray fired at a uniformly random angle hits a line. Its tails decay like \(1/x^2\), so slowly that the mean and variance do not exist: averages of Cauchy samples never settle down. \(x_0\) is the median and the peak; \(\gamma\) is the half-width at half the peak height.

Use it when you model resonance line shapes (the Lorentzian profile in spectroscopy), need a heavy-tailed likelihood that shrugs off outliers, or want a stress test for an estimator. Summarise it with the median and quantiles, never the sample mean. For heavy tails with a tunable weight, use studentPdf (one degree of freedom is the standard Cauchy).

heavy-tailsrobust-statisticslikelihood-evaluation

def cauchyPdf(x0, g, x)      ≔ g / (3.141592653589793 · (g · g + (x − x0) · (x − x0)))

cauchyCdf

cauchyCdf(x0, g, x) · x0 location · g scale, positive · x point → probability broadcast shape

Let \(x_0 \in \mathbb{R}\), \(\gamma > 0\) (the argument g) and \(x \in \mathbb{R}\). Then

\[\operatorname{cauchyCdf}(x_0, \gamma, x) = \frac12 + \frac1\pi \arctan\frac{x - x_0}{\gamma}.\]

Intuition An arctangent: half the mass lies on each side of \(x_0\), and the quartiles sit exactly at \(x_0 \pm \gamma\). It creeps towards 0 and 1 so slowly that about 6% of the mass lies more than \(10\gamma\) from the centre.

Use it when you need the probability that a Cauchy-distributed quantity exceeds a threshold, or a percentile rank under a heavy-tailed model. Set beside normalCdf at the same threshold, it shows how much more probability a heavy-tailed model gives to extreme values.

heavy-tailstail-probabilities

def cauchyCdf(x0, g, x)      ≔ 0.5 + atan((x − x0) / g) / 3.141592653589793

cauchyQuantile

cauchyQuantile(x0, g, p) · x0 location · g scale, positive · p probability → quantile broadcast shape

Let \(x_0 \in \mathbb{R}\), \(\gamma > 0\) (the argument g) and \(0 < p < 1\). Then

\[\operatorname{cauchyQuantile}(x_0, \gamma, p) = x_0 + \gamma \tan\!\Big(\pi \big(p - \tfrac12\big)\Big).\]

Intuition A tangent: uniform angles in \((-\pi/2, \pi/2)\) projected onto a line. The tangent blows up near the ends of that range, so quantiles near 0 or 1 are enormous, which is the heavy tail seen from the other side.

Use it when you generate Cauchy draws from uniforms (outlier-heavy test data, robustness experiments), or you need interval bounds under a Cauchy model: the central 50% interval is \(x_0 \pm \gamma\). There is no guard, so at \(p = 0\) or very close to 1 the tangent returns huge values; keep the uniforms strictly inside (0, 1).

random-variate-generationheavy-tailsrobust-statistics

def cauchyQuantile(x0, g, p) ≔ x0 + g · tan(3.141592653589793 · (p − 0.5))

exponentialPdf

exponentialPdf(r, x) · r rate, positive · x point → density broadcast shape

Let \(r > 0\) be the rate and \(x \in \mathbb{R}\). Then

\[\operatorname{exponentialPdf}(r, x) = [x \ge 0]\; r\, e^{-rx}.\]

\(r\,e^{-rx}\) for \(x \ge 0\), and 0 for negative x. The mean is \(1/r\).

Intuition The waiting time until the first event when events happen independently at a constant rate \(r\). It is memoryless: however long you have already waited, the remaining wait has the same distribution. The mean wait is \(1/r\), so the parameter is events per unit time, not a scale.

Use it when you model inter-arrival times (requests, customers, decays), lifetimes with a constant failure rate, or evaluate the likelihood of observed gaps (the maximum-likelihood rate is one over the mean gap). If the failure rate rises or falls with age, use weibullPdf; for the wait until the \(k\)-th event, use gammaPdf.

waiting-timesreliability-analysislikelihood-evaluation

def exponentialPdf(r, x)      ≔ (x ≥ 0) · r · exp(0 − r · x)

exponentialCdf

exponentialCdf(r, x) · r rate, positive · x point → probability broadcast shape

Let \(r > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{exponentialCdf}(r, x) = [x \ge 0]\,\big(1 - e^{-rx}\big).\]

Intuition The probability that the event has happened by time \(x\). Its complement \(e^{-rx}\) is the survival probability, which halves every \(\ln 2 / r\) time units, like radioactive decay.

Use it when you need the chance that a wait ends within a deadline (a queue served within a time limit under a constant rate, a component failing before its warranty ends) or, as \(1 - F\), the probability of surviving past a time. For the time by which a given fraction of waits have ended, use exponentialQuantile.

waiting-timesreliability-analysistail-probabilities

def exponentialCdf(r, x)      ≔ (x ≥ 0) · (1 − exp(0 − r · x))

exponentialQuantile

exponentialQuantile(r, p) · r rate, positive · p probability → quantile broadcast shape

Let \(r > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{exponentialQuantile}(r, p) = -\frac{\ln (1 - p)_+}{r},\]

which is \(-\ln(1 - p)/r\) for \(p < 1\) and a large finite value at \(p = 1\).

Intuition The time by which a fraction \(p\) of waits have ended. At \(p = \tfrac12\) it is the median wait \(\ln 2 / r\), shorter than the mean \(1/r\) because a few long waits pull the mean up.

Use it when you simulate arrivals or lifetimes from uniforms (add successive draws to get the event times of a Poisson process), or set a timeout that a target fraction of waits will beat. At \(p = 1\) the guarded logarithm gives a large finite value instead of infinity.

random-variate-generationwaiting-timessimulation

def exponentialQuantile(r, p) ≔ 0 − ln(positive(1 − p)) / r

laplacePdf

laplacePdf(m, b, x) · m location · b scale, positive · x point → density broadcast shape

Let \(m \in \mathbb{R}\) be the location, \(b > 0\) the scale and \(x \in \mathbb{R}\). Then

\[\operatorname{laplacePdf}(m, b, x) = \frac{1}{2b}\, e^{-\lvert x - m \rvert / b}.\]

The double-exponential density.

Intuition Two exponentials placed back to back around \(m\): a sharp peak, and tails that fall like \(e^{-\lvert x - m \rvert/b}\), heavier than the normal’s. It is also the difference of two independent exponential waits with the same rate. Maximising a Laplace likelihood minimises absolute error, which is why its location estimate is the median.

Use it when you model errors with occasional large deviations (L1 regression, robust fitting), put a sparsity-inducing prior on weights (the Bayesian reading of the lasso), or add noise for differential privacy (the Laplace mechanism uses scale = sensitivity / ε). For rarer, larger outliers use cauchyPdf or studentPdf.

robust-statisticslikelihood-evaluationdifferential-privacy

def laplacePdf(m, b, x)      ≔ exp(0 − abs(x − m) / b) / (2 · b)

laplaceCdf

laplaceCdf(m, b, x) · m location · b scale, positive · x point → probability broadcast shape

Let \(m \in \mathbb{R}\), \(b > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{laplaceCdf}(m, b, x) = \frac12 + \frac12 \operatorname{sign}(x - m)\,\big(1 - e^{-\lvert x - m \rvert / b}\big).\]

Intuition Half the mass lies on each side of \(m\), and each tail beyond a distance \(d\) holds \(\tfrac12 e^{-d/b}\): exponential decay in both directions.

Use it when you need the probability that Laplace noise exceeds a bound, e.g. how often a differentially private count is off by more than \(d\), which is \(e^{-d/b}\) for both sides together. It also gives percentile ranks of residuals under an L1 error model.

tail-probabilitiesdifferential-privacyrobust-statistics

def laplaceCdf(m, b, x)      ≔ 0.5 + 0.5 · sign(x − m) · (1 − exp(0 − abs(x − m) / b))

laplaceQuantile

laplaceQuantile(m, b, p) · m location · b scale, positive · p probability → quantile broadcast shape

Let \(m \in \mathbb{R}\), \(b > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{laplaceQuantile}(m, b, p) = m - b\, \operatorname{sign}\!\big(p - \tfrac12\big) \ln \Big(1 - 2\,\big\lvert p - \tfrac12 \big\rvert\Big)_+ .\]

Intuition The inverse of that two-sided decay: the distance from \(m\) grows with the logarithm of one over the tail probability, so extreme quantiles grow only slowly.

Use it when you draw Laplace noise from uniforms (the Laplace mechanism, simulated robust residuals) or need a symmetric error bound: the interval \(m \pm b \ln(1/\alpha)\) holds \(1 - \alpha\) of the mass. At \(p = 0\) or 1 the guarded logarithm gives large finite values.

random-variate-generationdifferential-privacyconfidence-intervals

def laplaceQuantile(m, b, p) ≔ m − b · sign(p − 0.5) · ln(positive(1 − 2 · abs(p − 0.5)))

logisticPdf

logisticPdf(m, s, x) · m location · s scale, positive · x point → density broadcast shape

Let \(m \in \mathbb{R}\) be the location, \(s > 0\) the scale, \(x \in \mathbb{R}\) and \(z = (x - m)/s\). Then

\[\operatorname{logisticPdf}(m, s, x) = \frac{e^{-z}}{s\,(1 + e^{-z})^2}.\]

Intuition The derivative of the sigmoid: a bell shape like the normal’s, with somewhat heavier exponential tails. It is the noise behind logistic regression: if \(y = 1\) exactly when a score plus logistic noise is positive, then \(P(y = 1)\) is a sigmoid of the score. Its variance is \(\pi^2 s^2/3\), so \(s\) is not the standard deviation.

Use it when you need a normal-like likelihood whose CDF has a closed form, model latent-utility noise, or fit growth and dose-response data. For probabilities use logisticCdf. The formula exponentiates \(-(x - m)/s\), so more than about 709 scales below \(m\) that exponential overflows and the result is NaN.

classificationprobabilistic-modellinglikelihood-evaluation

def logisticPdf(m, s, x)      ≔ exp(0 − (x − m) / s) / (s · (1 + exp(0 − (x − m) / s)) · (1 + exp(0 − (x − m) / s)))

logisticCdf

logisticCdf(m, s, x) · m location · s scale, positive · x point → probability broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{logisticCdf}(m, s, x) = \frac{1}{1 + e^{-(x - m)/s}}.\]

The logistic (sigmoid) function.

Intuition The sigmoid squashes any real score into a probability: it crosses one half at \(m\), and \(s\) sets how gradually it rises. Each step of \(s\) in \(x\) multiplies the odds by \(e\).

Use it when you turn a linear score into a class probability (logistic regression, Elo-style win probabilities), or model an S-shaped curve such as dose-response or adoption over time. It stays finite for large scores in either direction. For the reverse map from probability to score, use logisticQuantile.

classificationbinary-outcomesdose-response

def logisticCdf(m, s, x)      ≔ 1 / (1 + exp(0 − (x − m) / s))

logisticQuantile

logisticQuantile(m, s, p) · m location · s scale, positive · p probability → quantile broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{logisticQuantile}(m, s, p) = m + s \ln \frac{p_+}{(1 - p)_+}.\]

The logit, \(m + s \ln\big(p/(1 - p)\big)\) for \(0 < p < 1\).

Intuition The logit: a probability turned into log-odds, stretched by \(s\) and shifted by \(m\). Probabilities near 0 and 1 are pulled far apart, which is why log-odds are the natural scale for adding up evidence.

Use it when you convert a probability into a score (a base rate turned into a model’s bias term, a calibrated probability back into its linear predictor), or draw logistic noise from uniforms. Probabilities of exactly 0 or 1 are floored at \(10^{-300}\), giving about \(m \pm 690\,s\) rather than an infinity.

classificationrandom-variate-generation

def logisticQuantile(m, s, p) ≔ m + s · ln(positive(p) / positive(1 − p))

gumbelPdf

gumbelPdf(m, b, x) · m location · b scale, positive · x point → density broadcast shape

Let \(m \in \mathbb{R}\) be the location, \(b > 0\) the scale, \(x \in \mathbb{R}\) and \(z = (x - m)/b\). Then

\[\operatorname{gumbelPdf}(m, b, x) = \frac1b\, e^{-(z + e^{-z})}.\]

The Gumbel (maximum extreme value) density.

Intuition The distribution of the largest of many independent values with light (exponential-type) tails: the highest flood of the year, the worst latency of the day, the peak load on a structure. It leans right, because a maximum can overshoot far more easily than undershoot. \(m\) is the mode and \(b\) the spread.

Use it when you fit block maxima (annual peak rainfall, wind speed, daily worst response time) by likelihood, or model random utilities in discrete-choice models. For maxima of heavy-tailed data (insurance claims, market losses) it understates the tail; a Pareto tail (paretoPdf) is safer. For minima, apply it to the negated values.

extreme-valueslikelihood-evaluationrisk-measurement

def gumbelPdf(m, b, x)      ≔ exp(0 − (x − m) / b − exp(0 − (x − m) / b)) / b

gumbelCdf

gumbelCdf(m, b, x) · m location · b scale, positive · x point → probability broadcast shape

Let \(m \in \mathbb{R}\), \(b > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{gumbelCdf}(m, b, x) = \exp\!\big(-e^{-(x - m)/b}\big).\]

Intuition A double exponential, \(\exp(-e^{-z})\). If one period’s maximum has this CDF, the maximum over \(n\) independent periods is Gumbel again with \(m\) moved up by \(b \ln n\), so the family is stable under taking maxima.

Use it when you need the probability that a period’s maximum stays below a design level (a levee height, a capacity limit), or the return period of a level, \(1/(1 - F)\) periods on average between exceedances.

extreme-valuestail-probabilitiesrisk-measurement

def gumbelCdf(m, b, x)      ≔ exp(0 − exp(0 − (x − m) / b))

gumbelQuantile

gumbelQuantile(m, b, p) · m location · b scale, positive · p probability → quantile broadcast shape

Let \(m \in \mathbb{R}\), \(b > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{gumbelQuantile}(m, b, p) = m - b \ln \big(-\ln p_+\big)_+ ,\]

which is \(m - b \ln(-\ln p)\) for \(0 < p < 1\).

Intuition Turns a non-exceedance probability into a level: the value a period’s maximum stays below with probability \(p\). For long return periods the level grows with the logarithm of the period: each tenfold longer period adds about the same amount, \(b \ln 10\).

Use it when you compute the design level for a return period of \(T\) periods, gumbelQuantile(m, b, 1 − 1/T), or draw Gumbel noise from uniforms: with \(m = 0\) and \(b = 1\), the argmax of log-probabilities plus independent Gumbel draws samples a categorical choice (the Gumbel-max trick).

extreme-valuesrisk-measurementrandom-variate-generation

def gumbelQuantile(m, b, p) ≔ m − b · ln(positive(0 − ln(positive(p))))

paretoPdf

paretoPdf(xm, a, x) · xm minimum (scale), positive · a shape, positive · x point → density broadcast shape

Let \(x_m > 0\) (the argument xm) be the minimum, \(a > 0\) the shape and \(x \in \mathbb{R}\). Then

\[\operatorname{paretoPdf}(x_m, a, x) = [x \ge x_m]\; \frac{a\, x_m^{a}}{x_+^{\,a+1}}.\]

\(\dfrac{a\, x_m^a}{x^{a+1}}\) for \(x \ge x_m\), and 0 below it.

Intuition The power law behind “a few account for most”: wealth, city sizes, file sizes, insurance losses above a threshold. The density falls like \(x^{-(a+1)}\), so doubling a value divides its density by the same factor \(2^{a+1}\) at every scale. Smaller \(a\) means heavier tails: the mean is infinite for \(a \le 1\) and the variance for \(a \le 2\).

Use it when you fit the tail of heavy-tailed data above a threshold \(x_m\) (large-loss models, bursty network traffic), or evaluate a power-law likelihood. The density is 0 below \(x_m\), so choose the threshold before fitting. When data are skewed but all moments are finite, lognormalPdf is often the better model.

heavy-tailsrisk-measurementlikelihood-evaluation

def paretoPdf(xm, a, x)      ≔ (x ≥ xm) · a · xm^a / positive(x)^(a + 1)

paretoCdf

paretoCdf(xm, a, x) · xm minimum (scale), positive · a shape, positive · x point → probability broadcast shape

Let \(x_m > 0\) (the argument xm), \(a > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{paretoCdf}(x_m, a, x) = [x \ge x_m]\,\Big(1 - \big(x_m / x_+\big)^{a}\Big).\]

Intuition The survival probability is a pure power, \((x_m/x)^a\): the share of values above \(x\) drops by the factor \(2^a\) each time \(x\) doubles. The 80/20 rule corresponds to \(a \approx 1.16\).

Use it when you estimate exceedance probabilities for heavy-tailed quantities: how likely a claim above a large amount is, or a traffic burst above a link’s capacity. On a plot of \(\ln(1 - F)\) against \(\ln x\), a straight line of slope \(-a\) is the signature of a Pareto tail, which is a quick check before relying on this model.

heavy-tailstail-probabilitiesrisk-measurement

def paretoCdf(xm, a, x)      ≔ (x ≥ xm) · (1 − (xm / positive(x))^a)

paretoQuantile

paretoQuantile(xm, a, p) · xm minimum (scale), positive · a shape, positive · p probability → quantile broadcast shape

Let \(x_m > 0\) (the argument xm), \(a > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{paretoQuantile}(x_m, a, p) = \frac{x_m}{(1 - p)_+^{\,1/a}},\]

which is \(x_m (1 - p)^{-1/a}\) for \(p < 1\).

Intuition The level exceeded with probability \(1 - p\). As \(p \to 1\) it grows like a power of \(1/(1 - p)\), not a logarithm, which is why heavy-tail risk estimates explode at high confidence levels.

Use it when you compute a value-at-risk or a capacity level for heavy-tailed losses or loads, or draw power-law samples from uniforms (synthetic wealth, file sizes, workloads). At \(p = 1\) the result is huge, and for small \(a\) the floored power underflows and the result is infinite, so draw uniforms in [0, 1).

risk-measurementheavy-tailsrandom-variate-generation

def paretoQuantile(xm, a, p) ≔ xm / positive(1 − p)^(1 / a)

weibullPdf

weibullPdf(l, k, x) · l scale λ, positive · k shape, positive · x point → density broadcast shape

Let \(\lambda > 0\) (the argument l) be the scale, \(k > 0\) the shape and \(x \in \mathbb{R}\). Then

\[\operatorname{weibullPdf}(\lambda, k, x) = [x \ge 0]\; \frac{k}{\lambda}\left(\frac{x_+}{\lambda}\right)^{k-1} \exp\!\left(-\left(\frac{x_+}{\lambda}\right)^{k}\right).\]

Intuition A lifetime whose hazard (instantaneous failure rate) changes as a power of age: \(k < 1\) is infant mortality (failures thin out with age), \(k = 1\) is the exponential with rate \(1/\lambda\), and \(k > 1\) is wear-out. It also arises as the weakest link, the minimum of many independent strengths. About 63.2% of lifetimes end before \(\lambda\), whatever \(k\) is.

Use it when you model component lifetimes, material strength or wind speeds, and fit them by likelihood from failure data. For a constant hazard, exponentialPdf is the simpler special case; for data produced by multiplicative growth, compare with lognormalPdf.

reliability-analysiswaiting-timeslikelihood-evaluation

def weibullPdf(l, k, x)      ≔ (x ≥ 0) · (k / l) · (positive(x) / l)^(k − 1) · exp(0 − (positive(x) / l)^k)

weibullCdf

weibullCdf(l, k, x) · l scale λ, positive · k shape, positive · x point → probability broadcast shape

Let \(\lambda > 0\) (the argument l), \(k > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{weibullCdf}(\lambda, k, x) = [x \ge 0]\left(1 - \exp\!\left(-\left(\frac{x_+}{\lambda}\right)^{k}\right)\right).\]

Intuition The fraction of units failed by age \(x\). Because \(\ln(-\ln(1 - F))\) is a straight line in \(\ln x\) with slope \(k\), a Weibull plot of failure data reads the shape straight off the slope.

Use it when you estimate the probability that a part fails within its warranty or maintenance interval, the reliability \(1 - F\) at a mission time, or the share of wind speeds below a turbine’s cut-in speed.

reliability-analysistail-probabilities

def weibullCdf(l, k, x)      ≔ (x ≥ 0) · (1 − exp(0 − (positive(x) / l)^k))

weibullQuantile

weibullQuantile(l, k, p) · l scale λ, positive · k shape, positive · p probability → quantile broadcast shape

Let \(\lambda > 0\) (the argument l), \(k > 0\) and \(p \in [0, 1]\). Then

\[\operatorname{weibullQuantile}(\lambda, k, p) = \lambda\,\big(-\ln (1 - p)_+\big)^{1/k}.\]

Intuition The age by which a fraction \(p\) of units have failed. Reliability engineers call weibullQuantile(l, k, 0.1) the B10 life: the age at which 10% have failed.

Use it when you set maintenance or replacement intervals from a failure-fraction target, quote B-lives, or draw Weibull lifetimes from uniforms for a reliability simulation. With \(k = 1\) it reduces to the exponential quantile with rate \(1/\lambda\).

reliability-analysisrandom-variate-generationsimulation

def weibullQuantile(l, k, p) ≔ l · (0 − ln(positive(1 − p)))^(1 / k)

rayleighPdf

rayleighPdf(s, x) · s scale σ, positive · x point → density broadcast shape

Let \(\sigma > 0\) (the argument s) be the scale and \(x \in \mathbb{R}\). Then

\[\operatorname{rayleighPdf}(\sigma, x) = [x \ge 0]\; \frac{x}{\sigma^2}\, e^{-x^2 / 2\sigma^2}.\]

Intuition The length of a 2-D vector whose two components are independent normals with mean 0 and standard deviation \(\sigma\): a shot’s distance from the bullseye, the amplitude of a noisy signal, a speed built from two orthogonal components. The density is 0 at the centre and peaks at \(\sigma\), because a ring of small radius has little room in it.

Use it when you model radial miss distances, the envelope of narrowband noise or fading amplitudes, or wave heights, and evaluate their likelihood. It is the Weibull with \(k = 2\) and \(\lambda = \sigma\sqrt2\). If the vector has a non-zero mean (a direct signal path), the Rice distribution applies instead, and this module does not provide it.

signal-amplitudemeasurement-errorlikelihood-evaluation

def rayleighPdf(s, x)      ≔ (x ≥ 0) · x · exp(0 − x · x / (2 · s · s)) / (s · s)

rayleighCdf

rayleighCdf(s, x) · s scale σ, positive · x point → probability broadcast shape

Let \(\sigma > 0\) (the argument s) and \(x \in \mathbb{R}\). Then

\[\operatorname{rayleighCdf}(\sigma, x) = [x \ge 0]\,\big(1 - e^{-x^2 / 2\sigma^2}\big).\]

Intuition The probability that the vector’s length is at most \(x\), i.e. that a point lands inside the circle of radius \(x\) around the centre.

Use it when you compute a hit probability within a target radius, the chance that a fading signal’s amplitude drops below a receiver threshold (the outage probability), or, as \(1 - F\), the share of wave heights above a level.

tail-probabilitiessignal-amplitude

def rayleighCdf(s, x)      ≔ (x ≥ 0) · (1 − exp(0 − x · x / (2 · s · s)))

rayleighQuantile

rayleighQuantile(s, p) · s scale σ, positive · p probability → quantile broadcast shape

Let \(\sigma > 0\) (the argument s) and \(p \in [0, 1]\). Then

\[\operatorname{rayleighQuantile}(\sigma, p) = \sigma \sqrt{-2 \ln (1 - p)_+}.\]

Intuition The radius of the circle that contains a fraction \(p\) of the points. At \(p = \tfrac12\) it is the circular error probable, \(\sigma\sqrt{2 \ln 2} \approx 1.177\,\sigma\).

Use it when you size a target or tolerance circle for a positioning error (the CEP of a GPS fix), set a detection threshold on a noise envelope for a target false-alarm rate (\(p = 1 -\) that rate), or draw Rayleigh variates from uniforms.

confidence-intervalssignal-amplituderandom-variate-generation

def rayleighQuantile(s, p) ≔ s · sqrt(0 − 2 · ln(positive(1 − p)))

lognormalPdf

lognormalPdf(m, s, x) · m mean of ln X · s standard deviation of ln X, positive · x point → density broadcast shape

Let \(m \in \mathbb{R}\) be the mean and \(s > 0\) the standard deviation of \(\ln X\), and let \(x \in \mathbb{R}\). Then

\[\operatorname{lognormalPdf}(m, s, x) = [x > 0]\; \frac{\operatorname{normalPdf}(m, s, \ln x_+)}{x_+} = [x > 0]\; \frac{1}{x_+\, s\sqrt{2\pi}} \exp\!\left(-\frac{(\ln x_+ - m)^2}{2s^2}\right).\]

A variable whose logarithm is normal with mean m and standard deviation s. The density is the normal density of \(\ln x\) divided by x, for \(x > 0\).

Intuition What you get when many small independent effects multiply instead of add: incomes, compounded prices, particle sizes, response times. Its logarithm is normal, so it is positive and skewed to the right. \(m\) and \(s\) describe \(\ln X\), not \(X\): the median is \(e^{m}\) and the mean \(e^{m + s^2/2}\).

Use it when you model positive quantities whose spread grows with their size, or evaluate a likelihood for them; the fitted \(m\) and \(s\) are the mean and standard deviation of the logged data. Passing the raw data’s mean and standard deviation is the classic mistake. For power-law tails use paretoPdf.

probabilistic-modellinglikelihood-evaluationskewed-data

def lognormalPdf(m, s, x)      ≔ (x > 0) · normalPdf(m, s, ln(positive(x))) / positive(x)

lognormalCdf

lognormalCdf(m, s, x) · m mean of ln X · s standard deviation of ln X, positive · x point → probability broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{lognormalCdf}(m, s, x) = [x > 0]\; \operatorname{normalCdf}(m, s, \ln x_+).\]

The normal CDF of \(\ln x\), for \(x > 0\).

Intuition The normal CDF read on a log axis: being below \(x\) is the same event as \(\ln X\) being below \(\ln x\), because the logarithm preserves order.

Use it when you need the share of a skewed positive quantity under a threshold: the fraction of requests faster than a latency target, of loans below a size, of particles passing a sieve. As \(1 - F\) it gives the exceedance probability for a latency or loss budget, subject to the usual \(10^{-16}\) limit on upper tails.

tail-probabilitiesskewed-data

def lognormalCdf(m, s, x)      ≔ (x > 0) · normalCdf(m, s, ln(positive(x)))

lognormalQuantile

lognormalQuantile(m, s, p) · m mean of ln X · s standard deviation of ln X, positive · p probability → quantile broadcast shape

Let \(m \in \mathbb{R}\), \(s > 0\) and \(0 < p < 1\). Then

\[\operatorname{lognormalQuantile}(m, s, p) = \exp\big(\operatorname{normalQuantile}(m, s, p)\big) \approx e^{m + s\,\Phi^{-1}(p)}.\]

The exponential of the normal quantile.

Intuition The normal quantile of \(\ln X\), exponentiated. Percentiles carry over directly, but intervals become lopsided: the upper bound sits further from the median than the lower bound.

Use it when you quote percentiles of a log-normal model (p95 or p99 latency, a price level under geometric Brownian motion at a horizon), build bounds for a positive quantity, or draw log-normal samples from uniforms.

percentilesconfidence-intervalsrandom-variate-generation

def lognormalQuantile(m, s, p) ≔ exp(normalQuantile(m, s, p))

bernoulliPdf

bernoulliPdf(q, x) · q probability of success · x outcome, 0 or 1 → probability broadcast shape

Let \(q \in [0, 1]\) be the probability of success and \(x \in \mathbb{R}\). Then

\[\operatorname{bernoulliPdf}(q, x) = [x = 1]\, q + [x = 0]\,(1 - q).\]

The probability mass: q at 1, \(1 - q\) at 0, and 0 for any other value.

Intuition One yes/no trial that succeeds with probability \(q\): a coin flip, a click, a defective part. It is the building block of the geometric and binomial distributions.

Use it when you evaluate the likelihood of binary outcomes under predicted probabilities: the sum of the logs over observations is minus the log-loss (binary cross-entropy) a classifier minimises. Outcomes must be exactly 0 or 1; any other value gets probability 0. If a predicted \(q\) can be exactly 0 or 1, clamp it before taking logs. For the number of successes in \(n\) trials use binomialPdf.

binary-outcomeslikelihood-evaluationclassification

def bernoulliPdf(q, x)      ≔ (x = 1) · q + (x = 0) · (1 − q)

bernoulliCdf

bernoulliCdf(q, x) · q probability of success · x point → probability broadcast shape

Let \(q \in [0, 1]\) and \(x \in \mathbb{R}\). Then

\[\operatorname{bernoulliCdf}(q, x) = [x \ge 0]\,(1 - q) + [x \ge 1]\, q .\]

0 below 0, \(1 - q\) on [0, 1), and 1 from 1 upwards.

Intuition A two-step staircase: probability \(1 - q\) sits at 0 and the rest at 1, so the CDF jumps to \(1 - q\) at 0 and to 1 at 1.

Use it when a generic code path needs a CDF for every family, for example a randomised probability integral transform that checks the calibration of binary predictions. In most binary work you want bernoulliPdf (likelihoods) or bernoulliQuantile (simulation) instead.

binary-outcomesmodel-calibration

def bernoulliCdf(q, x)      ≔ (x ≥ 0) · (1 − q) + (x ≥ 1) · q

bernoulliQuantile

bernoulliQuantile(q, p) · q probability of success · p probability → outcome 0 or 1 broadcast shape

Let \(q \in [0, 1]\) and \(p \in [0, 1]\). Then

\[\operatorname{bernoulliQuantile}(q, p) = [\,p > 1 - q\,].\]

1 when \(p > 1 - q\), otherwise 0. Applied to uniform draws it produces Bernoulli trials.

Intuition Cuts [0, 1) at \(1 - q\): a uniform landing above the cut is a success, which happens with probability \(q\). Raising \(q\) lowers the cut, so with the same uniform a success never turns into a failure.

Use it when you simulate yes/no events (a customer converts, a packet is dropped, a unit is dropped out) from supplied uniforms. Reusing the same uniforms across scenarios gives common random numbers, which makes the differences between scenarios far less noisy than independent runs would.

simulationrandom-variate-generationbinary-outcomes

def bernoulliQuantile(q, p) ≔ (p > 1 − q)

geometricPdf

geometricPdf(q, x) · q probability of success per trial · x number of failures before the first success, a whole number → probability broadcast shape

Let \(q \in (0, 1]\) be the probability of success per trial and \(x \in \mathbb{R}\) (normally \(x \in \mathbb{N}\)). Then

\[\operatorname{geometricPdf}(q, x) = [x \ge 0]\; q\,(1 - q)^{x}.\]

\(q\,(1 - q)^x\) for \(x = 0, 1, 2, \ldots\). The support starts at 0 (failures counted, not trials). A fractional x is not rounded.

Intuition The number of failures before the first success in independent trials that each succeed with probability \(q\): calls before a sale, retries before a packet gets through. It is the discrete waiting time and, like the exponential, memoryless. Counting starts at 0, and the mean is \((1 - q)/q\).

Use it when you model retry counts or run lengths, or evaluate the likelihood of observed failure counts. If your data count attempts including the success, subtract 1 first. A fractional \(x\) is not rounded, so pass whole numbers. For successes in a fixed number of trials use binomialPdf.

waiting-timescount-datalikelihood-evaluation

def geometricPdf(q, x)      ≔ (x ≥ 0) · q · (1 − q)^x

geometricCdf

geometricCdf(q, x) · q probability of success per trial · x point → probability broadcast shape

Let \(q \in (0, 1]\) and \(x \in \mathbb{R}\). Then

\[\operatorname{geometricCdf}(q, x) = [x \ge 0]\,\Big(1 - (1 - q)^{\lfloor x \rfloor + 1}\Big).\]

Intuition The probability that the first success comes within the first \(\lfloor x \rfloor + 1\) trials: one minus the chance that all of them fail.

Use it when you ask how many attempts are enough: the probability that a policy allowing up to 5 retries (6 attempts) eventually succeeds is geometricCdf(q, 5), and the same call gives the chance a user converts within a number of visits. Unlike geometricPdf, a fractional \(x\) is floored here.

waiting-timescount-datatail-probabilities

def geometricCdf(q, x)      ≔ (x ≥ 0) · (1 − (1 − q)^(floor(x) + 1))

geometricQuantile

geometricQuantile(q, p) · q probability of success per trial · p probability → whole number broadcast shape

Let \(q \in (0, 1)\) and \(p \in (0, 1)\). Then

\[\operatorname{geometricQuantile}(q, p) = \left\lceil \frac{\ln (1 - p)_+}{\ln (1 - q)_+} - 1 \right\rceil .\]

The smallest count k with \(P(X \le k) \ge p\). At \(p = 0\) the formula gives \(-1\) rather than 0.

Intuition The smallest failure count \(k\) whose CDF reaches \(p\), solved in closed form: take logarithms of \((1 - q)^{k+1} \le 1 - p\) and round up, with no search.

Use it when you size a retry budget for a target success probability, or draw geometric counts from uniforms. One edge case follows from the formula: at \(p = 0\) (a uniform draw of exactly 0) the ceiling of \(0 - 1\) gives \(-1\), not 0, so clamp the result at 0 when sampling.

random-variate-generationwaiting-timessimulation

def geometricQuantile(q, p) ≔ ceil(ln(positive(1 − p)) / ln(positive(1 − q)) − 1)

The regularised incomplete beta function

The regularised incomplete beta function is

\[I_x(a, b) = \frac{1}{B(a, b)} \int_0^x t^{a-1} (1 - t)^{b-1}\, dt\]

It is evaluated by Lentz’s continued fraction (240 sub-steps). The fraction converges quickly for \(x < (a+1)/(a+b+2)\); on the other side the symmetry \(I_x(a, b) = 1 - I_{1-x}(b, a)\) is used.

betaI

betaI(a, b, x) · a first shape, positive · b second shape, positive · x point in [0, 1] → \(I_x(a, b)\) in [0, 1] broadcast shape

Let \(a > 0\), \(b > 0\), \(x \in [0, 1]\) and \(\xi = (a + 1)/(a + b + 2)\). Then

\[\operatorname{betaI}(a, b, x) = [x < \xi]\; \operatorname{betaNear}\big(a, b, \min(x, \xi)\big) + [x \ge \xi]\; \operatorname{betaFar}\big(a, b, \max(x, \xi)\big) \approx I_x(a, b),\]

where

\[\operatorname{betaNear}(a, b, x) = \frac{x^{a} (1 - x)^{b}}{a\, B(a, b)}\, \operatorname{betacf}(a, b, x), \qquad \operatorname{betaFar}(a, b, x) = 1 - \frac{(1 - x)^{b}\, x^{a}}{b\, B(a, b)}\, \operatorname{betacf}(b, a, 1 - x),\]

away from the guard thresholds, and \(\operatorname{betacf}\) is the continued fraction

\[\operatorname{betacf}(a, b, x) \approx \cfrac{1}{1 + \cfrac{\tau_0}{1 + \cfrac{\tau_1}{1 + \cfrac{\tau_2}{1 + \cdots}}}}, \qquad \tau_{2m-1} = \frac{m(b - m)\,x}{(a + 2m - 1)(a + 2m)}, \quad \tau_{2m} = -\frac{(a + m)(a + b + m)\,x}{(a + 2m)(a + 2m + 1)} \quad (m \ge 1 \text{ and } m \ge 0 \text{ respectively}),\]

cut after 240 sub-steps (see the building blocks for the exact recurrence).

The regularised incomplete beta function, choosing the direct fraction or the symmetric one by the mask \(x < (a+1)/(a+b+2)\). It is the CDF of a beta distribution.

Intuition For whole \(a\) and \(b\), \(I_x(a, b)\) is the probability of at least \(a\) successes in \(a + b - 1\) trials with success probability \(x\); equivalently, the chance that the \(a\)-th smallest of \(a + b - 1\) uniform draws is at most \(x\). That is why beta, binomial, Student’s t and F probabilities all reduce to it.

Use it when you need a probability built on it that the module does not name, such as an order-statistic probability or a binomial upper tail. Prefer betaCdf, binomialCdf, studentCdf or fCdf when one fits. The continued fraction runs a fixed 240 sub-steps with no convergence test, so check very large shapes against a reference value.

proportionstail-probabilitieshypothesis-testing

def betaI(a, b, x)    ≔ (x < betaEdge(a, b)) · betaNear(a, b, min(x, betaEdge(a, b))) + (x ≥ betaEdge(a, b)) · betaFar(a, b, max(x, betaEdge(a, b)))

Building blocks

• bcOdd(j) — Let \(j \in \mathbb{N}\). Then \(\operatorname{bcOdd}(j) = j - 2\lfloor j/2 \rfloor = [\,j \text{ odd}\,]\). 1 when sub-step j is odd, 0 when even; selects the coefficient formula.
• bcFirst(a, b, x, m) — Let \(a, b > 0\), \(x \in [0, 1]\) and \(m \in \mathbb{R}\). Then \(\operatorname{bcFirst}(a, b, x, m) = m(b - m)\,x / \operatorname{guard}\big((a - 1 + 2m)(a + 2m)\big)\). The odd-step coefficient of the continued fraction.
• bcSecond(a, b, x, m) — Let \(a, b > 0\), \(x \in [0, 1]\) and \(m \in \mathbb{R}\). Then \(\operatorname{bcSecond}(a, b, x, m) = -(a + m)(a + b + m)\,x / \operatorname{guard}\big((a + 2m)(a + 1 + 2m)\big)\). The even-step coefficient of the continued fraction.
• bcTerm(a, b, x, j) — Let \(a, b > 0\), \(x \in [0, 1]\), \(j \in \mathbb{N}\) and \(\omega = j - 2\lfloor j/2 \rfloor\). Then \(\operatorname{bcTerm}(a, b, x, j) = \tau_j = \omega\, \operatorname{bcFirst}\big(a, b, x, \tfrac{j+1}{2}\big) + (1 - \omega)\, \operatorname{bcSecond}\big(a, b, x, \tfrac{j}{2}\big)\). The coefficient of sub-step j, chosen by parity.
• bcD(t, d) — Let \(t, d \in \mathbb{R}\). Then \(\operatorname{bcD}(t, d) = 1 / \operatorname{guard}(1 + t\,d)\). The Lentz D update for coefficient t.
• bcC(t, c) — Let \(t, c \in \mathbb{R}\). Then \(\operatorname{bcC}(t, c) = \operatorname{guard}\big(1 + t / \operatorname{guard}(c)\big)\). The Lentz C update for coefficient t.
• bcLoop(k, a, b, x, j, c, d, h, t) — Let \(k \in \mathbb{N}\), \(a, b > 0\), \(x \in [0, 1]\), \(j \in \mathbb{N}\) and \(c, d, h, t \in \mathbb{R}\), and write \(c' = \operatorname{bcC}(t, c)\) and \(d' = \operatorname{bcD}(t, d)\). Then \(\operatorname{bcLoop}(0, a, b, x, j, c, d, h, t) = h\) and \(\operatorname{bcLoop}(k, a, b, x, j, c, d, h, t) = \operatorname{bcLoop}(k - 1, a, b, x, j + 1, c', d', h\,d'\,c', \tau_{j+2})\) with \(\tau_{j+2} = \operatorname{bcTerm}(a, b, x, j + 2)\). Runs k Lentz sub-steps, carrying the current coefficient and the running value h.
• bcSeed(a, b, x) — Let \(a, b > 0\) and \(x \in [0, 1]\). Then \(\operatorname{bcSeed}(a, b, x) = 1 / \operatorname{guard}\big(1 - (a + b)\,x / \operatorname{guard}(a + 1)\big)\). The starting D (and value) of the fraction.
• betacf(a, b, x) — Let \(a, b > 0\) and \(x \in [0, 1]\), and \(\tau_j = \operatorname{bcTerm}(a, b, x, j)\). Let \(c_0 = 1\), \(d_0 = h_0 = \operatorname{bcSeed}(a, b, x)\), and for \(i = 1, \dots, 240\) let \(c_i = \operatorname{guard}(1 + \tau_i / \operatorname{guard}(c_{i-1}))\), \(d_i = 1/\operatorname{guard}(1 + \tau_i\, d_{i-1})\) and \(h_i = h_{i-1}\, d_i\, c_i\). Then \(\operatorname{betacf}(a, b, x) = h_{240}\). The continued fraction itself, 240 sub-steps.
• betaEdge(a, b) — Let \(a, b > 0\). Then \(\operatorname{betaEdge}(a, b) = \xi = (a + 1)/(a + b + 2)\). The switch-over point.
• betaFront(a, b, x) — Let \(a, b > 0\) and \(x \in [0, 1]\). Then \(\operatorname{betaFront}(a, b, x) = \exp\big(a \ln x_+ + b \ln (1 - x)_+ - \operatorname{lbeta}(a, b)\big)\). The prefactor x to the a, times (1 − x) to the b, divided by B(a, b).
• betaNear(a, b, x) — Let \(a, b > 0\) and \(x \in [0, 1]\). Then \(\operatorname{betaNear}(a, b, x) = \operatorname{betaFront}(a, b, x)\, \operatorname{betacf}(a, b, x) / \operatorname{guard}(a)\). The incomplete beta by the direct fraction; used below the switch-over point.
• betaFar(a, b, x) — Let \(a, b > 0\) and \(x \in [0, 1]\). Then \(\operatorname{betaFar}(a, b, x) = 1 - \operatorname{betaFront}(b, a, 1 - x)\, \operatorname{betacf}(b, a, 1 - x) / \operatorname{guard}(b)\). The incomplete beta by the symmetric fraction; used at or above the switch-over point.

Distributions from the incomplete gamma: χ², gamma, Poisson

The χ² distribution with v degrees of freedom is a gamma distribution with shape \(v/2\) and scale 2; the Poisson CDF is an upper incomplete gamma, \(P(X \le k) = Q(k + 1, \lambda)\), which turns a count distribution into a continuous computation.

chiPdf

chiPdf(v, x) · v degrees of freedom, positive · x point → density broadcast shape

Let \(v > 0\) be the degrees of freedom and \(x \in \mathbb{R}\). Then

\[\operatorname{chiPdf}(v, x) = [x > 0]\; \exp\!\left(\left(\frac v2 - 1\right) \ln x_+ - \frac x2 - \frac v2 \ln 2 - \operatorname{lgamma}\!\left(\frac v2\right)\right) \approx [x > 0]\; \frac{x^{v/2-1} e^{-x/2}}{2^{v/2}\,\Gamma(v/2)},\]

with \(\ln 2\) written as the constant 0.6931471805599453.

The χ² density for \(x > 0\), and 0 otherwise.

Intuition The sum of \(v\) squared independent standard normals: the total squared deviation of \(v\) independent standardised errors. Its mean is \(v\), so a statistic near \(v\) is unremarkable. It is the gamma distribution with shape \(v/2\) and scale 2.

Use it when you evaluate the likelihood of a sum of squared standardised residuals or of a scaled sample variance (\((n-1)S^2/\sigma^2\) has \(n - 1\) degrees of freedom), or draw the reference curve of a χ² test. For p-values you need chiCdf, not the density.

hypothesis-testinglikelihood-evaluationgoodness-of-fit

def chiPdf(v, x)   ≔ (x > 0) · exp((v / 2 − 1) · ln(positive(x)) − x / 2 − (v / 2) · 0.6931471805599453 − lgamma(v / 2))

chiCdf

chiCdf(v, x) · v degrees of freedom, positive · x point → probability broadcast shape

Let \(v > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{chiCdf}(v, x) = [x > 0]\; \operatorname{lowerGamma}\!\left(\frac v2,\; \frac{\max(x, 0)}{2}\right) \approx [x > 0]\; P\!\left(\frac v2, \frac x2\right).\]

For the p-value of a χ² statistic use 1 − chiCdf(v, x).

Intuition The probability that a sum of \(v\) squared standard normals is at most \(x\). A statistic far above \(v\) is unlikely under the null model, and the upper tail \(1 - F\) measures how unlikely.

Use it when you turn a χ² statistic into a p-value, 1 − chiCdf(v, x): a contingency-table test (chiSquared and freedom in Categorical), a goodness-of-fit test, or a likelihood-ratio test (\(2\,\Delta \ln L\), with as many degrees of freedom as constrained parameters). p-values below about \(10^{-16}\) come out as 0.

hypothesis-testinggoodness-of-fittail-probabilities

def chiCdf(v, x)   ≔ (x > 0) · lowerGamma(v / 2, max(x, 0) / 2)

gammaPdf

gammaPdf(s, r, x) · s shape, positive · r scale θ, positive · x point → density broadcast shape

Let \(s > 0\) be the shape, \(\theta > 0\) (the argument r) the scale and \(x \in \mathbb{R}\). Then

\[\operatorname{gammaPdf}(s, \theta, x) = [x > 0]\; \exp\!\big((s - 1) \ln x_+ - x/\theta - s \ln \theta - \operatorname{lgamma}(s)\big) \approx [x > 0]\; \frac{x^{s-1} e^{-x/\theta}}{\theta^s\,\Gamma(s)}.\]

The gamma density with shape s and scale θ (mean \(s\theta\)). For a rate β, pass \(1/\beta\).

Intuition The waiting time until the \(s\)-th event of a Poisson process whose mean gap is θ: a sum of \(s\) independent exponential waits. Non-integer \(s\) interpolates smoothly; small \(s\) is sharply skewed and large \(s\) looks normal. The mean is \(s\theta\) and the variance \(s\theta^2\).

Use it when you model positive, right-skewed totals (repair times, rainfall amounts, claim sizes) or use it as the conjugate prior for a Poisson rate. Watch the parameterisation: the argument r is the scale θ, so a rate β must be passed as \(1/\beta\). For a single wait, exponentialPdf takes the rate \(1/\theta\).

waiting-timesbayesian-updatinglikelihood-evaluation

def gammaPdf(s, r, x) ≔ (x > 0) · exp((s − 1) · ln(positive(x)) − x / r − s · ln(r) − lgamma(s))

gammaCdf

gammaCdf(s, r, x) · s shape, positive · r scale θ, positive · x point → probability broadcast shape

Let \(s > 0\), \(\theta > 0\) (the argument r) and \(x \in \mathbb{R}\). Then

\[\operatorname{gammaCdf}(s, \theta, x) = [x > 0]\; \operatorname{lowerGamma}\big(s,\; \max(x, 0)/\theta\big) \approx [x > 0]\; P(s, x/\theta).\]

Intuition The probability that the \(s\)-th event has arrived by time \(x\), or that a gamma-distributed total stays below \(x\).

Use it when you need the chance that \(s\) sequential exponential stages finish within a deadline (a multi-step job, a chain of repairs), a posterior probability from a gamma posterior on a rate (“the rate is below 2 per day”), or, as \(1 - F\), the tail probability of a gamma-modelled claim. As with gammaPdf, pass the scale, not the rate.

waiting-timestail-probabilitiesbayesian-updating

def gammaCdf(s, r, x) ≔ (x > 0) · lowerGamma(s, max(x, 0) / r)

poissonPdf

poissonPdf(l, x) · l mean λ, positive · x count, a whole number → probability broadcast shape

Let \(\lambda > 0\) (the argument l) be the mean and \(x \in \mathbb{R}\) (normally \(x \in \mathbb{N}\)). Then

\[\operatorname{poissonPdf}(\lambda, x) = [x \ge 0]\; \exp\!\big(x \ln \lambda_+ - \lambda - \operatorname{lgamma}(x + 1)\big) \approx [x \ge 0]\; \frac{\lambda^x e^{-\lambda}}{x!}.\]

The Poisson probability mass, computed through lgamma. It is intended for whole-number x; a fractional x is not rounded and gives a smooth interpolation rather than 0.

Intuition The number of events in a fixed window when events occur independently at a constant average rate: calls per hour, defects per metre, arrivals per minute. \(\lambda\) is both the mean and the variance, and the distribution is the limit of a binomial with many trials and a small success probability.

Use it when you evaluate the likelihood of counts (Poisson regression, rare-event rates, anomaly scores on event counts). If counts vary more than their mean (overdispersion), a Poisson fit understates the uncertainty; the usual alternative, the negative binomial, is not in this module. The mass is computed in logs, so large λ and counts are safe; pass whole numbers.

count-datalikelihood-evaluationrare-events

def poissonPdf(l, x) ≔ (x ≥ 0) · exp(x · ln(positive(l)) − l − lgamma(x + 1))

poissonCdf

poissonCdf(l, x) · l mean λ, positive · x point → probability broadcast shape

Let \(\lambda > 0\) (the argument l) and \(x \in \mathbb{R}\). Then

\[\operatorname{poissonCdf}(\lambda, x) = [x \ge 0]\; \operatorname{upperGamma}\big(\lfloor x \rfloor + 1,\; \lambda\big) \approx [x \ge 0]\; Q(\lfloor x \rfloor + 1, \lambda) = P(X \le \lfloor x \rfloor).\]

Intuition The probability of at most \(\lfloor x \rfloor\) events when \(\lambda\) are expected, computed as one upper incomplete gamma rather than a sum of masses.

Use it when you ask whether an observed count is surprisingly high (the p-value for “at least \(k\)” is 1 − poissonCdf(l, k − 1)), plan capacity for peak counts, or set stock against random demand. Note the floor: poissonCdf(l, 2.9) equals poissonCdf(l, 2).

count-datahypothesis-testingcapacity-planning

def poissonCdf(l, x) ≔ (x ≥ 0) · upperGamma(floor(x) + 1, l)

Distributions from the incomplete beta: beta, Student’s t, F, binomial

Student’s t, F and the binomial CDF are all incomplete beta functions after a change of variable, so they share one continued fraction. The binomial identity is the discrete twin of the Poisson one: \(P(X \le k) = I_{1-q}(n - k,\; k + 1)\).

betaPdf

betaPdf(a, b, x) · a first shape, positive · b second shape, positive · x point → density broadcast shape

Let \(a > 0\) and \(b > 0\) be the shapes and \(x \in \mathbb{R}\). Then

\[\operatorname{betaPdf}(a, b, x) = [x > 0]\,[x < 1]\; \exp\!\big((a - 1) \ln x_+ + (b - 1) \ln (1 - x)_+ - \operatorname{lbeta}(a, b)\big) \approx [0 < x < 1]\; \frac{x^{a-1}(1 - x)^{b-1}}{B(a, b)}.\]

The beta density for \(0 < x < 1\), and 0 otherwise.

Intuition A distribution over a probability or proportion in (0, 1). Read \(a - 1\) and \(b - 1\) as successes and failures already seen: for \(a, b > 1\) the peak is at \((a-1)/(a+b-2)\), and the curve narrows as \(a + b\) grows. \(a = b = 1\) is uniform, and \(a, b < 1\) is U-shaped.

Use it when you represent uncertainty about a rate: after \(k\) successes in \(n\) trials, a Beta\((a, b)\) prior becomes Beta\((a + k,\; b + n - k)\), the basis of Bayesian A/B testing and Thompson sampling. It also models bounded fractions such as a budget share. The density is set to 0 at the endpoints, even where the true density diverges.

proportionsbayesian-updatinguncertainty-quantification

def betaPdf(a, b, x)  ≔ (x > 0) · (x < 1) · exp((a − 1) · ln(positive(x)) + (b − 1) · ln(positive(1 − x)) − lbeta(a, b))

betaCdf

betaCdf(a, b, x) · a first shape, positive · b second shape, positive · x point → probability broadcast shape

Let \(a > 0\), \(b > 0\) and \(x \in \mathbb{R}\). Then

\[\operatorname{betaCdf}(a, b, x) = [x > 0]\,\Big([x \ge 1] + [x < 1]\; \operatorname{betaI}\big(a, b, \min(\max(x, 0), 1)\big)\Big).\]

\(I_x(a, b)\) on (0, 1), 0 at or below 0 and 1 at or above 1.

Intuition The probability that an unknown proportion is at most \(x\), for a belief described by Beta\((a, b)\).

Use it when you need a posterior probability such as “the conversion rate is below 5%” or “variant B beats a fixed target”, or the probability that a bounded fraction stays under a limit. Credible-interval bounds go the other way and come from betaQuantile. Arguments outside (0, 1) are safe: 0 at or below 0, 1 at or above 1.

proportionsbayesian-updatingtail-probabilities

def betaCdf(a, b, x)  ≔ (x > 0) · ((x ≥ 1) + (x < 1) · betaI(a, b, min(max(x, 0), 1)))

studentPdf

studentPdf(v, t) · v degrees of freedom, positive · t point → density broadcast shape

Let \(v > 0\) be the degrees of freedom and \(t \in \mathbb{R}\). Then

\[\operatorname{studentPdf}(v, t) = \exp\!\Big(\operatorname{lgamma}\!\big(\tfrac{v+1}{2}\big) - \operatorname{lgamma}\!\big(\tfrac v2\big) - \tfrac12 \ln(v\pi)\Big) \left(1 + \frac{t^2}{v}\right)^{-\frac{v+1}{2}} \approx \frac{\Gamma\!\left(\frac{v+1}{2}\right)}{\sqrt{v\pi}\;\Gamma\!\left(\frac{v}{2}\right)} \left(1 + \frac{t^2}{v}\right)^{-\frac{v+1}{2}}.\]

Student’s t density.

Intuition A standardised mean when the standard deviation is itself estimated from the data: the \(v\) degrees of freedom of that estimate set how much the tails are fattened. With \(v = 1\) it is the standard Cauchy; as \(v\) grows it approaches the standard normal. It takes no location or scale arguments.

Use it when you want a heavy-tailed likelihood that down-weights outliers (robust regression, noise models in Bayesian fits): standardise as \((x - m)/\sigma\) and divide the density by \(\sigma\). It also draws the reference curve of a t test. For probabilities use studentCdf.

heavy-tailsrobust-statisticslikelihood-evaluation

def studentPdf(v, t)  ≔ exp(lgamma((v + 1) / 2) − lgamma(v / 2) − 0.5 · ln(v · 3.141592653589793)) · (1 + t · t / v)^(0 − (v + 1) / 2)

studentCdf

studentCdf(v, t) · v degrees of freedom, positive · t point → probability broadcast shape

Let \(v > 0\) and \(t \in \mathbb{R}\). Then

\[\operatorname{studentCdf}(v, t) = \frac12 + \operatorname{sign}(t)\left(\frac12 - \operatorname{studentTail}(v, t)\right), \qquad \operatorname{studentTail}(v, t) = \frac12\, \operatorname{betaI}\!\left(\frac v2,\; \frac12,\; \frac{v}{v + t^2}\right).\]

The probability that a t variable is at most t, from the one-sided tail \(\tfrac12 I_{v/(v+t^2)}(v/2, \tfrac12)\). A two-sided p-value for a statistic t is \(2\,(1 - \operatorname{studentCdf}(v, \lvert t \rvert))\).

Intuition The probability that a t statistic with \(v\) degrees of freedom is at most \(t\). With few degrees of freedom the tails are heavier, so the same statistic is weaker evidence than it would be under the normal.

Use it when you compute a p-value for a one- or two-sample t test, a regression coefficient, or any mean comparison with a small sample and an estimated variance. The two-sided p-value is \(2\,(1 - F(\lvert t \rvert))\), as in the Examples. As \(v\) grows the answer approaches normalCdf(0, 1, t).

hypothesis-testingcomparing-groupstail-probabilities

def studentCdf(v, t)  ≔ 0.5 + sign(t) · (0.5 − studentTail(v, t))

fPdf

fPdf(d1, d2, x) · d1 numerator degrees of freedom, positive · d2 denominator degrees of freedom, positive · x point → density broadcast shape

Let \(d_1 > 0\) and \(d_2 > 0\) (the arguments d1, d2) be the degrees of freedom and \(x \in \mathbb{R}\). Then

\[\operatorname{fPdf}(d_1, d_2, x) = [x > 0]\; \exp\!\Big(\tfrac{d_1}{2} \ln d_1 + \tfrac{d_2}{2} \ln d_2 + \big(\tfrac{d_1}{2} - 1\big) \ln x_+ - \tfrac{d_1 + d_2}{2} \ln(d_1 x_+ + d_2) - \operatorname{lbeta}\big(\tfrac{d_1}{2}, \tfrac{d_2}{2}\big)\Big),\]

which for \(x > 0\) is

\[\frac{d_1^{d_1/2}\, d_2^{d_2/2}\, x^{d_1/2 - 1}}{(d_1 x + d_2)^{(d_1+d_2)/2}\, B(d_1/2,\, d_2/2)} .\]

The F density, computed in logarithms.

Intuition The ratio of two independent variance estimates, each a χ² divided by its degrees of freedom: \(d_1\) above the line, \(d_2\) below. When the true variances are equal the ratio sits near 1; values far above 1 mean the numerator’s source of variation dominates.

Use it when you evaluate the likelihood or draw the null curve of an ANOVA, regression or variance-ratio statistic. For p-values use fCdf. The density is computed in logs, so large degrees of freedom are safe.

hypothesis-testingcomparing-variances

def fPdf(d1, d2, x)   ≔ (x > 0) · exp((d1 / 2) · ln(d1) + (d2 / 2) · ln(d2) + (d1 / 2 − 1) · ln(positive(x)) − ((d1 + d2) / 2) · ln(d1 · positive(x) + d2) − lbeta(d1 / 2, d2 / 2))

fCdf

fCdf(d1, d2, x) · d1 numerator degrees of freedom · d2 denominator degrees of freedom · x point → probability broadcast shape

Let \(d_1, d_2 > 0\) (the arguments d1, d2) and \(x \in \mathbb{R}\), and write \(x^{*} = \max(x, 0)\). Then

\[\operatorname{fCdf}(d_1, d_2, x) = [x > 0]\; \operatorname{betaI}\!\left(\frac{d_1}{2},\; \frac{d_2}{2},\; \frac{d_1 x^{*}}{d_1 x^{*} + d_2}\right).\]

\(I_{d_1 x/(d_1 x + d_2)}(d_1/2,\; d_2/2)\) for \(x > 0\).

Intuition The probability that a variance ratio is at most \(x\) when the two underlying variances are truly equal.

Use it when you turn an F statistic into a p-value, 1 − fCdf(d1, d2, x): one-way ANOVA (between-group over within-group mean squares), the overall F test of a regression, a nested-model comparison, or a test of equal variances. Order matters: swapping \(d_1\) and \(d_2\) describes \(1/F\), not \(F\).

hypothesis-testingcomparing-variancescomparing-groups

def fCdf(d1, d2, x)   ≔ (x > 0) · betaI(d1 / 2, d2 / 2, d1 · max(x, 0) / (d1 · max(x, 0) + d2))

binomialPdf

binomialPdf(n, q, x) · n number of trials · q probability of success · x number of successes, a whole number → probability broadcast shape

Let \(n \in \mathbb{N}\) be the number of trials, \(q \in [0, 1]\) the probability of success and \(x \in \mathbb{R}\) (normally \(x \in \{0, \dots, n\}\)). Then

\[\operatorname{binomialPdf}(n, q, x) = [x \ge 0]\,[x \le n]\; \exp\!\big(\operatorname{lgamma}(n + 1) - \operatorname{lgamma}(x + 1) - \operatorname{lgamma}(n - x + 1) + x \ln q_+ + (n - x) \ln (1 - q)_+\big),\]

which for whole \(x\) is \(\binom{n}{x} q^x (1 - q)^{n-x}\).

The binomial probability mass for \(0 \le x \le n\), computed through lgamma. Like poissonPdf, it expects whole-number x and does not round.

Intuition The number of successes in \(n\) independent trials that each succeed with probability \(q\): a sum of \(n\) Bernoulli outcomes. The mean is \(nq\); the shape is symmetric at \(q = \tfrac12\) and skewed otherwise.

Use it when you evaluate the likelihood of \(k\) successes out of \(n\) (conversions from visitors, defects in a batch, correct answers on a test), or need exact probabilities where a normal approximation is crude: small \(n\), or \(q\) near 0 or 1. It is computed in logs, so large \(n\) is safe; pass a whole-number \(x\). For a rare event over very many trials, poissonPdf is the simpler limit.

count-dataproportionslikelihood-evaluation

def binomialPdf(n, q, x) ≔ (x ≥ 0) · (x ≤ n) · exp(lgamma(n + 1) − lgamma(x + 1) − lgamma(n − x + 1) + x · ln(positive(q)) + (n − x) · ln(positive(1 − q)))

binomialCdf

binomialCdf(n, q, x) · n number of trials · q probability of success · x point → probability broadcast shape

Let \(n \in \mathbb{N}\), \(q \in [0, 1]\) and \(x \in \mathbb{R}\). Then

\[\operatorname{binomialCdf}(n, q, x) = [x \ge n] + [x \ge 0]\,[x < n]\; \operatorname{betaI}\big(n - \lfloor x \rfloor,\; \lfloor x \rfloor + 1,\; 1 - q\big).\]

\(P(X \le \lfloor x \rfloor)\), via \(I_{1-q}(n - \lfloor x \rfloor,\; \lfloor x \rfloor + 1)\); 1 for \(x \ge n\).

Intuition The probability of at most \(\lfloor x \rfloor\) successes, obtained from one incomplete beta instead of adding up to \(n + 1\) masses.

Use it when you run an exact binomial test (are 9 defects in 100 more than a 5% process allows? The p-value is 1 − binomialCdf(100, 0.05, 8)), check an acceptance-sampling plan, or compute the chance that a quorum or majority is reached. Arguments at or past \(n\) give 1, and negative ones give 0.

hypothesis-testingproportionsquality-control

def binomialCdf(n, q, x) ≔ (x ≥ n) + (x ≥ 0) · (x < n) · betaI(n − floor(x), floor(x) + 1, 1 − q)

Building blocks

• studentTail(v, t) — Let \(v > 0\) and \(t \in \mathbb{R}\). Then \(\operatorname{studentTail}(v, t) = \tfrac12\, \operatorname{betaI}\big(v/2,\; \tfrac12,\; v/(v + t^2)\big)\). The probability beyond the absolute value of t in one tail; used by studentCdf.

Quantiles without a closed form

For Student’s t, χ² and gamma the quantile is found by Newton’s method on the CDF, \(x \leftarrow x - (F(x) - p)/f(x)\), for a fixed 60 steps from a starting point that is already close. The beta quantile uses bisection instead, which cannot fail; the F quantile is obtained from the beta quantile.

studentQuantile

studentQuantile(v, p) · v degrees of freedom, positive · p probability in (0, 1) → quantile broadcast shape

Let \(v > 0\) and \(0 < p < 1\). Let \(x_0 = \operatorname{standardQuantile}(p)\) and

\[x_{k+1} = x_k - \frac{\operatorname{studentCdf}(v, x_k) - p}{\operatorname{guard}\big(\operatorname{studentPdf}(v, x_k)\big)}, \qquad k = 0, 1, \dots, 59 .\]

Then \(\operatorname{studentQuantile}(v, p) = x_{60}\).

The value t with \(P(T \le t) = p\). Newton’s method starts from the standard normal quantile, which t approaches as the degrees of freedom grow. For a two-sided test at level α use studentQuantile(v, 1 − α/2).

Intuition The t critical value: how many estimated standard errors to step away from an estimate to capture probability \(p\). For small \(v\) it exceeds the normal value (2.228 against 1.96 at \(v = 10\) and \(p = 0.975\)); that gap is the price of estimating the variance.

Use it when you build a confidence interval for a mean or a regression coefficient, \(\hat\mu \pm t_{v,\,1-\alpha/2}\,\mathrm{SE}\), or need a critical value for a t test. Each call runs 60 Newton steps, each one an incomplete-beta evaluation, so compute it once per \((v, p)\) rather than once per observation.

confidence-intervalshypothesis-testing

def studentQuantile(v, p) ≔ tLoop(60, v, p, standardQuantile(p))

chiQuantile

chiQuantile(v, p) · v degrees of freedom, positive · p probability in (0, 1) → quantile broadcast shape

Let \(v > 0\) and \(0 < p < 1\). Let \(x_0 = \max\big(\operatorname{wilson}(v, p),\, 10^{-6}\big)\) and

\[x_{k+1} = \max\!\left(x_k - \frac{\operatorname{chiCdf}(v, x_k) - p}{\operatorname{guard}\big(\operatorname{chiPdf}(v, x_k)\big)},\; 10^{-12}\right), \qquad k = 0, 1, \dots, 59 .\]

Then \(\operatorname{chiQuantile}(v, p) = x_{60}\).

The χ² critical value with \(P(X \le x) = p\). Newton’s method starts from the Wilson–Hilferty approximation (the cube root of \(\chi^2/v\) is nearly normal) and keeps every iterate positive.

Intuition The value that a sum of \(v\) squared standard normals stays below with probability \(p\): the rejection threshold of a χ² test, and the squared radius of a Gaussian confidence ellipsoid in \(v\) dimensions.

Use it when you need a critical value (reject independence when chiSquared exceeds chiQuantile(freedom, 0.95)), a confidence interval for a normal variance, \(\big[(n-1)S^2/\chi^2_{1-\alpha/2},\; (n-1)S^2/\chi^2_{\alpha/2}\big]\) on \(n - 1\) degrees of freedom, or the size of a confidence ellipse. Like studentQuantile it runs 60 Newton steps, so compute it once rather than per data point.

confidence-intervalshypothesis-testinggoodness-of-fit

def chiQuantile(v, p)   ≔ chiLoop(60, v, p, max(wilson(v, p), 1e-6))

gammaQuantile

gammaQuantile(s, r, p) · s shape, positive · r scale θ, positive · p probability in (0, 1) → quantile broadcast shape

Let \(s > 0\), \(\theta > 0\) (the argument r) and \(0 < p < 1\). Let \(x_0 = \theta \max\big(\operatorname{wilson}(2s, p)/2,\, 10^{-6}\big)\) and

\[x_{k+1} = \max\!\left(x_k - \frac{\operatorname{gammaCdf}(s, \theta, x_k) - p}{\operatorname{guard}\big(\operatorname{gammaPdf}(s, \theta, x_k)\big)},\; 10^{-12}\right), \qquad k = 0, 1, \dots, 59 .\]

Then \(\operatorname{gammaQuantile}(s, \theta, p) = x_{60}\).

The gamma quantile, by Newton’s method from the Wilson–Hilferty starting point for the matching χ² (2s degrees of freedom), scaled by θ.

Intuition The time by which the \(s\)-th event of a process with mean gap θ has arrived with probability \(p\). It equals θ/2 times the χ² quantile with \(2s\) degrees of freedom, which is where its starting point comes from.

Use it when you need a credible interval for a Poisson rate from a gamma posterior, a duration that a multi-stage process meets with a given probability, or an exact (Garwood) confidence interval for a Poisson mean, whose bounds are gamma quantiles with unit scale. Pass the scale θ, not a rate. At 60 Newton steps per value it is too slow to use as a per-draw sampler in large simulations.

confidence-intervalsbayesian-updatingwaiting-times

def gammaQuantile(s, r, p)   ≔ gammaLoop(60, s, r, p, r · max(wilson(2 · s, p) / 2, 1e-6))

betaQuantile

betaQuantile(a, b, p) · a first shape, positive · b second shape, positive · p probability in [0, 1] → quantile in [0, 1] broadcast shape

Let \(a > 0\), \(b > 0\) and \(p \in [0, 1]\). Let \(\ell_0 = 0\) and \(u_0 = 1\), and for \(k = 0, 1, \dots, 59\) let \(\mu_k = (\ell_k + u_k)/2\), \(\omega_k = [\operatorname{betaCdf}(a, b, \mu_k) < p]\) and

\[\ell_{k+1} = \omega_k\, \mu_k + (1 - \omega_k)\, \ell_k, \qquad u_{k+1} = \omega_k\, u_k + (1 - \omega_k)\, \mu_k .\]

Then \(\operatorname{betaQuantile}(a, b, p) = \mu_{60} = (\ell_{60} + u_{60})/2\).

The beta quantile by bisection on [0, 1]: 60 halvings, which exhausts double precision. Bisection is used because densities such as Beta(½, ½) diverge at the ends, where Newton’s method stalls. The price is several incomplete-beta evaluations per step, so this is the most expensive quantile in the module.

Intuition The proportion below which a Beta\((a, b)\) belief places probability \(p\). It is found by halving [0, 1] sixty times, which cannot fail even for U-shaped densities where Newton’s method stalls.

Use it when you report a Bayesian credible interval for a rate (the posterior quantiles at \(\alpha/2\) and \(1 - \alpha/2\)), compute an exact Clopper–Pearson interval for a binomial proportion (its bounds are beta quantiles), or feed fQuantile. It is the most expensive quantile in the module, so avoid it as a per-draw sampler in large simulations.

confidence-intervalsproportionsbayesian-updating

def betaQuantile(a, b, p)       ≔ betaLoop(60, a, b, p, 0, 1)

fQuantile

fQuantile(d1, d2, p) · d1 numerator degrees of freedom · d2 denominator degrees of freedom · p probability in (0, 1) → quantile broadcast shape

Let \(d_1, d_2 > 0\) (the arguments d1, d2), \(0 < p < 1\) and \(y = \operatorname{betaQuantile}(d_1/2, d_2/2, p)\). Then

\[\operatorname{fQuantile}(d_1, d_2, p) = \frac{d_2\, y}{\operatorname{guard}\big(d_1 (1 - y)\big)} .\]

The F critical value. The F CDF is an incomplete beta, so the beta quantile y of \((d_1/2, d_2/2)\) is transformed back as \(d_2\, y / \big(d_1 (1 - y)\big)\).

Intuition The F critical value, with no iteration of its own: an F ratio \(x\) corresponds to the beta variable \(d_1 x/(d_1 x + d_2)\), so inverting the beta and mapping the answer back inverts the F.

Use it when you need the rejection threshold of an ANOVA or regression F test at level α, fQuantile(d1, d2, 1 − α), or confidence bounds for the ratio of two normal variances. It inherits the cost of betaQuantile, so compute it once per test.

hypothesis-testingcomparing-variancesconfidence-intervals

def fQuantile(d1, d2, p) ≔ d2 · betaQuantile(d1 / 2, d2 / 2, p) / guard(d1 · (1 − betaQuantile(d1 / 2, d2 / 2, p)))

Building blocks

• tTurn(v, p, x) — Let \(v > 0\), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{tTurn}(v, p, x) = x - \big(\operatorname{studentCdf}(v, x) - p\big) / \operatorname{guard}\big(\operatorname{studentPdf}(v, x)\big)\). One Newton step on the Student CDF.
• tLoop(k, v, p, x) — Let \(k \in \mathbb{N}\), \(v > 0\), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{tLoop}(0, v, p, x) = x\) and \(\operatorname{tLoop}(k, v, p, x) = \operatorname{tLoop}(k - 1, v, p, \operatorname{tTurn}(v, p, x))\). Repeats tTurn k times; used by studentQuantile.
• wilson(v, p) — Let \(v > 0\), \(p \in (0, 1)\) and \(z = \operatorname{standardQuantile}(p)\). Then \(\operatorname{wilson}(v, p) = v \Big(1 - \frac{2}{9v} + z \sqrt{\frac{2}{9v}}\Big)^{3}\). The Wilson–Hilferty approximation to the χ² quantile; the starting point for χ² and gamma.
• chiTurn(v, p, x) — Let \(v > 0\), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{chiTurn}(v, p, x) = \max\big(x - (\operatorname{chiCdf}(v, x) - p) / \operatorname{guard}(\operatorname{chiPdf}(v, x)),\; 10^{-12}\big)\). One Newton step on the χ² CDF, kept positive.
• chiLoop(k, v, p, x) — Let \(k \in \mathbb{N}\), \(v > 0\), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{chiLoop}(0, v, p, x) = x\) and \(\operatorname{chiLoop}(k, v, p, x) = \operatorname{chiLoop}(k - 1, v, p, \operatorname{chiTurn}(v, p, x))\). Repeats chiTurn k times; used by chiQuantile.
• gammaTurn(s, r, p, x) — Let \(s > 0\), \(\theta > 0\) (the argument r), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{gammaTurn}(s, \theta, p, x) = \max\big(x - (\operatorname{gammaCdf}(s, \theta, x) - p) / \operatorname{guard}(\operatorname{gammaPdf}(s, \theta, x)),\; 10^{-12}\big)\). One Newton step on the gamma CDF, kept positive.
• gammaLoop(k, s, r, p, x) — Let \(k \in \mathbb{N}\), \(s > 0\), \(\theta > 0\) (the argument r), \(p \in (0, 1)\) and \(x \in \mathbb{R}\). Then \(\operatorname{gammaLoop}(0, s, \theta, p, x) = x\) and \(\operatorname{gammaLoop}(k, s, \theta, p, x) = \operatorname{gammaLoop}(k - 1, s, \theta, p, \operatorname{gammaTurn}(s, \theta, p, x))\). Repeats gammaTurn k times; used by gammaQuantile.
• betaMid(lo, hi) — Let \(\ell, u \in \mathbb{R}\) (the arguments lo, hi). Then \(\operatorname{betaMid}(\ell, u) = (\ell + u)/2\). The midpoint of the bisection interval.
• betaUnder(a, b, p, lo, hi) — Let \(a, b > 0\), \(p \in [0, 1]\), \(\ell, u \in [0, 1]\) (the arguments lo, hi) and \(\mu = (\ell + u)/2\). Then \(\operatorname{betaUnder}(a, b, p, \ell, u) = [\operatorname{betaCdf}(a, b, \mu) < p]\). 1 when the beta CDF at the midpoint is below p, else 0.
• betaLow(a, b, p, lo, hi) — Let \(a, b, p, \ell, u, \mu\) be as in betaUnder and \(\omega = \operatorname{betaUnder}(a, b, p, \ell, u)\). Then \(\operatorname{betaLow}(a, b, p, \ell, u) = \omega\,\mu + (1 - \omega)\,\ell\). The new lower end of the interval.
• betaHigh(a, b, p, lo, hi) — Let \(a, b, p, \ell, u, \mu\) and \(\omega\) be as in betaLow. Then \(\operatorname{betaHigh}(a, b, p, \ell, u) = \omega\,u + (1 - \omega)\,\mu\). The new upper end of the interval.
• betaLoop(k, a, b, p, lo, hi) — Let \(k \in \mathbb{N}\) and \(a, b, p, \ell, u\) be as in betaUnder. Then \(\operatorname{betaLoop}(0, a, b, p, \ell, u) = (\ell + u)/2\) and \(\operatorname{betaLoop}(k, a, b, p, \ell, u) = \operatorname{betaLoop}\big(k - 1, a, b, p, \operatorname{betaLow}(a, b, p, \ell, u), \operatorname{betaHigh}(a, b, p, \ell, u)\big)\). Runs k bisection steps and returns the final midpoint; used by betaQuantile.

Quantiles of counts

A discrete CDF is a step function, so Newton’s method cannot invert it. The quantile of a count is the smallest k whose CDF reaches p, found by counting upwards from 0 and stopping on the data. The cost therefore grows with the answer.

poissonQuantile

poissonQuantile(l, p) · l mean λ, positive · p probability in (0, 1) → whole number scalar

Let \(\lambda > 0\) (the argument l) and \(0 < p < 1\). Then

\[\operatorname{poissonQuantile}(\lambda, p) = \operatorname{poissonFrom}(0, \lambda, p) = \min\{\, k \in \{0, 1, 2, \dots\} : \operatorname{poissonCdf}(\lambda, k) \ge p \,\}.\]

The smallest k with \(P(X \le k) \ge p\) for a Poisson variable with mean λ.

Intuition The smallest count \(k\) that covers probability \(p\): stock \(k\) units and Poisson demand exceeds the stock with probability at most \(1 - p\). It is found by trying \(k = 0, 1, 2, \dots\) until the CDF reaches \(p\).

Use it when you size capacity or inventory for a service level, set an alert threshold on event counts, or draw Poisson counts from uniforms. The cost grows with the answer (about λ steps for central \(p\), each an incomplete-gamma evaluation), so large λ is slow; there normalQuantile(l, sqrt(l), p) is a close approximation. A \(p\) very close to 1 makes the search long.

capacity-planningcount-datarandom-variate-generation

def poissonQuantile(l, p)     ≔ poissonFrom(0, l, p)

binomialQuantile

binomialQuantile(n, q, p) · n number of trials · q probability of success · p probability in (0, 1] → whole number scalar

Let \(n \in \mathbb{N}\), \(q \in [0, 1]\) and \(0 < p \le 1\). Then

\[\operatorname{binomialQuantile}(n, q, p) = \operatorname{binomialFrom}(0, n, q, p) = \min\{\, k \in \{0, 1, 2, \dots\} : \operatorname{binomialCdf}(n, q, k) \ge p \,\}.\]

The smallest k with \(P(X \le k) \ge p\) for a binomial variable; never more than n, because \(\operatorname{binomialCdf}(n, q, n) = 1\).

Intuition The smallest number of successes \(k\) out of \(n\) whose cumulative probability reaches \(p\), found by counting up from 0; it never passes \(n\).

Use it when you set an acceptance threshold (reject a batch whose defects exceed binomialQuantile(n, q, 0.95)), plan overbooking (how many of \(n\) booked customers show up, at 99% confidence), or draw binomial counts from uniforms. The cost grows with the answer, up to \(n + 1\) CDF evaluations.

quality-controlcapacity-planningrandom-variate-generation

def binomialQuantile(n, q, p) ≔ binomialFrom(0, n, q, p)

Building blocks

• poissonFrom(k, l, p) — Let \(k \in \mathbb{N}\), \(\lambda > 0\) (the argument l) and \(p \in (0, 1)\). Then \(\operatorname{poissonFrom}(k, \lambda, p) = k\) if \(\operatorname{poissonCdf}(\lambda, k) \ge p\), and \(\operatorname{poissonFrom}(k, \lambda, p) = \operatorname{poissonFrom}(k + 1, \lambda, p)\) otherwise. Returns k if the Poisson CDF at k has reached p, otherwise tries k + 1.
• binomialFrom(k, n, q, p) — Let \(k, n \in \mathbb{N}\), \(q \in [0, 1]\) and \(p \in (0, 1]\). Then \(\operatorname{binomialFrom}(k, n, q, p) = k\) if \(\operatorname{binomialCdf}(n, q, k) \ge p\), and \(\operatorname{binomialFrom}(k, n, q, p) = \operatorname{binomialFrom}(k + 1, n, q, p)\) otherwise. Returns k if the binomial CDF at k has reached p, otherwise tries k + 1.

Examples

All examples below are compiled from one file, Probability.hb, so the C# class is Probability. Each snippet shows only the doorways that example uses.

Probabilities and cut-offs of the normal

What fraction of a normal population with mean 100 and standard deviation 15 lies below 130, and which score marks the top 10%?

open Prelude.Distributions

fn belowNormal(m, s, x) ↦ normalCdf(m, s, x)
fn normalCut(m, s, p) ↦ normalQuantile(m, s, p)

using Examples;

double below = Probability.BelowNormal(100.0, 15.0, 130.0);
double top10 = Probability.NormalCut(100.0, 15.0, 0.9);
double z975  = Probability.NormalCut(0.0, 1.0, 0.975);

Output (verified, rounded to ten significant digits):

normalCdf(100, 15, 130)      = 0.9772498681
normalQuantile(100, 15, 0.9) = 119.2232734832
normalQuantile(0, 1, 0.975)  = 1.9599639845

130 is two standard deviations above the mean, so the answer is \(\Phi(2) \approx 0.97725\); the familiar 1.96 is the 97.5% point of the standard normal.

Turn uniform draws into samples

Sampling is a quantile applied to uniform draws. Here the draws are chosen by hand so the results can be checked: 0.5 is the median, 0.8413447 is \(\Phi(1)\) and 0.0227501 is \(\Phi(-2)\).

open Prelude.Distributions

fn normalDraws(m, s, u) ↦ normalSample(m, s, u)
fn exponentialDraws(r, u) ↦ exponentialQuantile(r, u)

var normals = Probability.NormalDraws(10.0, 2.0, new[] { 0.5, 0.8413447, 0.0227501, 0.975 });
var waits   = Probability.ExponentialDraws(2.0, new[] { 0.0, 0.5, 0.75 });

Output (verified, rounded to six decimal places):

normalSample(10, 2, u)       = [10, 12, 5.999999, 13.919928]
exponentialQuantile(2, u)    = [0, 0.346574, 0.693147]

The normal draws are the mean, one and minus two standard deviations away (the 5.999999 reflects the rounded input), and the 97.5% point. The exponential draws with rate 2 are \(-\ln(1-u)/2\): 0, \(\ln 2 / 2\) and \(\ln 4 / 2\).

A t-test p-value and critical value

A t statistic of 2.228 on 10 degrees of freedom, and the critical values for a two-sided 5% test at 10 and a one-sided 5% test at 3 degrees of freedom.

open Prelude.Distributions

fn tTwoSided(v, t) ↦ 2 · (1 − studentCdf(v, abs(t)))
fn tCritical(v, p) ↦ studentQuantile(v, p)

double pValue = Probability.TTwoSided(10.0, 2.228);
double crit10 = Probability.TCritical(10.0, 0.975);
double crit3  = Probability.TCritical(3.0, 0.95);

Output (verified, rounded to six significant digits):

two-sided p (v = 10, t = 2.228) = 0.0500118
studentQuantile(10, 0.975)      = 2.22814
studentQuantile(3, 0.95)        = 2.35336

These match the printed t tables (2.228 and 2.353).

A χ² p-value and critical value

open Prelude.Distributions

fn chiPValue(v, x) ↦ 1 − chiCdf(v, x)
fn chiCritical(v, p) ↦ chiQuantile(v, p)

double p    = Probability.ChiPValue(2.0, 5.991);
double crit = Probability.ChiCritical(4.0, 0.95);

Output (verified, rounded to six significant digits):

1 − chiCdf(2, 5.991)   = 0.0500116
chiQuantile(4, 0.95)   = 9.48773

With 2 degrees of freedom the χ² upper tail is \(e^{-x/2}\), and \(e^{-5.991/2} \approx 0.05\); 9.488 is the tabulated 5% point for 4 degrees of freedom.

Counts: Poisson and binomial

The chance of at most 2 events when 3 are expected, the 90% point of that Poisson, and binomial probabilities for fair coin tosses.

open Prelude.Distributions

fn poissonAtMost(l, k) ↦ poissonCdf(l, k)
fn poissonCut(l, p) ↦ poissonQuantile(l, p)
fn binomialAtMost(n, q, k) ↦ binomialCdf(n, q, k)
fn binomialMass(n, q, k) ↦ binomialPdf(n, q, k)

double atMost2 = Probability.PoissonAtMost(3.0, 2.0);
double cut90   = Probability.PoissonCut(3.0, 0.9);
double heads3  = Probability.BinomialAtMost(10.0, 0.5, 3.0);
var    masses  = Probability.BinomialMass(4.0, 0.5, new[] { 0d, 1d, 2d, 3d, 4d });

Output (verified, rounded to six significant digits):

poissonCdf(3, 2)             = 0.423190
poissonQuantile(3, 0.9)      = 5
binomialCdf(10, 0.5, 3)      = 0.171875
binomialPdf(4, 0.5, 0..4)    = [0.0625, 0.25, 0.375, 0.25, 0.0625]

By hand: \(e^{-3}(1 + 3 + 4.5) \approx 0.42319\); at most 3 heads in 10 tosses is \((1 + 10 + 45 + 120)/1024 = 0.171875\); and four tosses give \(\tfrac{1}{16}(1, 4, 6, 4, 1)\).

Special functions and a U-shaped beta

open Prelude.Distributions

fn errorFunction(x) ↦ erf(x)
fn logGamma(a) ↦ lgamma(a)
fn betaCut(a, b, p) ↦ betaQuantile(a, b, p)

double erf1   = Probability.ErrorFunction(1.0);
double lg5    = Probability.LogGamma(5.0);
double lgHalf = Probability.LogGamma(0.5);
double q99    = Probability.BetaCut(0.5, 0.5, 0.99);

Output (verified, rounded to ten significant digits):

erf(1)                     = 0.8427007929
lgamma(5)                  = 3.1780538303
lgamma(0.5)                = 0.5723649429
betaQuantile(0.5, 0.5, 0.99) = 0.9997532802

\(\ln\Gamma(5) = \ln 24\) and \(\ln\Gamma(\tfrac12) = \ln\sqrt\pi\). Beta(½, ½) is the arcsine distribution, whose quantile is \(\sin^2(\pi p / 2)\); at \(p = 0.99\) that is 0.99975328.