Distributions and special functions#

Probability distributions: cumulative distribution functions, survival functions, their logarithms and quantiles, compiled with numba.

Names and parametrisations follow scipy.stats: norm, t, f, chi2, gamma (shape and scale), beta and hypergeom (population size M, number of successes n in the population, number of draws N). The survival function sf is P(X > x); logsf and logcdf stay finite far beyond the point where the plain probabilities underflow, which p-values of strong effects routinely reach.

Every public function is vectorised: its arguments broadcast with numpy’s rules and the result is a float64 array of the broadcast shape, or a Python float when every argument is a scalar. NaN arguments give NaN, and so do parameters outside their domain; the functions never raise for numeric input. Arguments that are not real numbers raise TypeError, and so do numpy masked arrays, whose mask numpy would drop: NaN marks a missing value instead.

The compiled scalar routines behind the public functions are named <function>_scalar (for example t_sf_scalar) and can be called from other compiled code, such as per-gene loops.

Methods

  • Normal: Q(x) = erfc(x / sqrt 2) / 2 from the C library’s erfc, with the rounding error of x / sqrt 2 corrected to first order (Dekker’s exact product, Dekker 1971), so the tails keep full relative accuracy; beyond x = 37, log Q(x) from Laplace’s continued fraction for Mills’ ratio (DLMF 7.9.2). Quantiles by Wichura’s (1988) algorithm AS 241.

  • Student t (Abramowitz and Stegun 1964, 26.7): P(|T| > |x|) = I_w(df/2, 1/2) with w = df / (df + x^2); df may be non-integer and df = inf gives the normal distribution. Quantiles: closed forms for df = 1, 2 and inf; otherwise Newton’s method in log t on log P(T > t), or on log P(|T| <= t) for p > 1/4 where the quantile is small, started from the better of the Cornish-Fisher expansion (Abramowitz and Stegun 1964, 26.7.5) and the power-law tail, and kept inside a bisection bracket.

  • F (Abramowitz and Stegun 1964, 26.6; Johnson, Kotz and Balakrishnan 1995, chapter 27): P(F > x) = I_w(dfd/2, dfn/2) with w = dfd / (dfd + dfn x); dfd = inf gives dfn F ~ chi-squared(dfn), and dfn = inf gives dfd / F ~ chi-squared(dfd). Quantiles: Newton’s method in log x on the logarithm of the tail that holds min(p, 1 - p), started from the better of Paulson’s (1942) normal approximation to F^(1/3) and the power law of that tail, and kept inside a bisection bracket; dfn = inf uses 1 / F ~ F(dfd, dfn).

  • Chi-squared and gamma (Abramowitz and Stegun 1964, 26.4): the regularized incomplete gamma function; chi-squared quantiles are those of X / df ~ F(df, inf).

  • Beta: the regularized incomplete beta function.

  • Hypergeometric: each probability from Loader’s (2000) saddle-point form of the binomial probabilities, P(X = k) = b(k; n, p) b(N - k; M - n, p) / b(N; M, p) with p = N / M; the tail on the far side of the mode is summed from its largest term by the ratio of successive probabilities, with compensated addition, until the rest of the tail (bounded through the log-concavity of the probabilities) is negligible, and the other tail is its complement in log space.

The incomplete gamma and beta functions are called with their arguments formed from the distribution’s own inputs (for example lambda = a w’ - b w for the beta function is formed as df (x^2 - 1) / (2 (df + x^2)) for the t distribution), so no accuracy is lost to rounding a transformed argument; see cytepool.stats.special.

References

Abramowitz M, Stegun IA, eds (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Applied Mathematics Series 55, National Bureau of Standards, Washington.

Dekker TJ (1971). A floating-point technique for extending the available precision. Numerische Mathematik 18, 224-242. doi:10.1007/BF01397083

Johnson NL, Kotz S, Balakrishnan N (1995). Continuous Univariate Distributions, Volume 2, 2nd edition. Wiley, New York.

Loader C (2000). Fast and accurate computation of binomial probabilities. Technical note, Lucent Technologies, Murray Hill; the version of 5 February 2002 is available at https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf

Paulson E (1942). An approximate normalization of the analysis of variance distribution. Annals of Mathematical Statistics 13, 233-235. doi:10.1214/aoms/1177731608

Wichura MJ (1988). Algorithm AS 241: the percentage points of the normal distribution. Journal of the Royal Statistical Society C 37, 477-484. doi:10.2307/2347330

Olver FWJ, Olde Daalhuis AB, Lozier DW, Schneider BI, Boisvert RF, Clark CW, Miller BR, Saunders BV, Cohl HS, McClain MA, eds (2026). NIST Digital Library of Mathematical Functions, release 1.2.8 of 2026-09-15. https://dlmf.nist.gov/

cytepool.stats.distributions.beta_cdf(x, a, b)#

Cumulative distribution function of the beta distribution, P(X <= x).

Parameters:
  • x – Quantile; x <= 0 gives 0 and x >= 1 gives 1.

  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

Returns:

P(X <= x), as a float or a float64 array of the broadcast shape. Shapes outside their domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.beta_sf(x, a, b)#

Survival function of the beta distribution, P(X > x).

Parameters:
  • x – Quantile; x <= 0 gives 1 and x >= 1 gives 0.

  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

Returns:

P(X > x), as a float or a float64 array of the broadcast shape. Shapes outside their domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.chi2_cdf(x, df)#

Cumulative distribution function of the chi-squared distribution, P(X <= x).

Parameters:
  • x – Quantile; x <= 0 gives 0.

  • df – Degrees of freedom, 0 < df < inf, not necessarily an integer.

Returns:

P(X <= x), as a float or a float64 array of the broadcast shape. df outside (0, inf) gives NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.chi2_isf(p, df)#

Inverse survival function of the chi-squared distribution: P(X > x) = p.

Accurate for p down to the smallest doubles.

Parameters:
  • p – Probability in [0, 1]; 0 gives +inf, 1 gives 0, values outside give NaN.

  • df – Degrees of freedom, 0 < df < inf, not necessarily an integer.

Returns:

The quantile, as a float or a float64 array of the broadcast shape. df outside (0, inf) gives NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.chi2_logsf(x, df)#

Logarithm of the survival function of the chi-squared distribution, log P(X > x).

Parameters:
  • x – Quantile; x <= 0 gives 0.

  • df – Degrees of freedom, 0 < df < inf, not necessarily an integer.

Returns:

log P(X > x), as a float or a float64 array of the broadcast shape. df outside (0, inf) gives NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.chi2_ppf(p, df)#

Quantile function of the chi-squared distribution: the x with P(X <= x) = p.

Parameters:
  • p – Probability in [0, 1]; 0 gives 0, 1 gives +inf, values outside give NaN.

  • df – Degrees of freedom, 0 < df < inf, not necessarily an integer.

Returns:

The quantile, as a float or a float64 array of the broadcast shape. df outside (0, inf) gives NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.chi2_sf(x, df)#

Survival function of the chi-squared distribution, P(X > x).

Parameters:
  • x – Quantile; x <= 0 gives 1.

  • df – Degrees of freedom, 0 < df < inf, not necessarily an integer.

Returns:

P(X > x), as a float or a float64 array of the broadcast shape. df outside (0, inf) gives NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.f_cdf(x, dfn, dfd)#

Cumulative distribution function of the F distribution, P(F <= x).

Parameters:
  • x – Quantile; x <= 0 gives 0.

  • dfn – Numerator degrees of freedom, dfn > 0; inf gives dfd / F ~ chi-squared(dfd).

  • dfd – Denominator degrees of freedom, dfd > 0; inf gives dfn F ~ chi-squared(dfn).

Returns:

P(F <= x), as a float or a float64 array of the broadcast shape. Degrees of freedom <= 0 give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.f_isf(p, dfn, dfd)#

Inverse survival function of the F distribution: the x with P(F > x) = p.

Accurate for p down to the smallest doubles; use this rather than f_ppf(1 - p) for upper-tail quantiles.

Parameters:
  • p – Probability in [0, 1]; 0 gives +inf, 1 gives 0, values outside give NaN.

  • dfn – Numerator degrees of freedom, dfn > 0; inf gives dfd / F ~ chi-squared(dfd).

  • dfd – Denominator degrees of freedom, dfd > 0; inf gives dfn F ~ chi-squared(dfn).

Returns:

The quantile, as a float or a float64 array of the broadcast shape. Degrees of freedom <= 0 give NaN; a quantile beyond the largest double gives +inf.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.f_logsf(x, dfn, dfd)#

Logarithm of the survival function of the F distribution, log P(F > x).

Parameters:
  • x – Quantile; x <= 0 gives 0.

  • dfn – Numerator degrees of freedom, dfn > 0; inf gives dfd / F ~ chi-squared(dfd).

  • dfd – Denominator degrees of freedom, dfd > 0; inf gives dfn F ~ chi-squared(dfn).

Returns:

log P(F > x), as a float or a float64 array of the broadcast shape. Degrees of freedom <= 0 give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.f_ppf(p, dfn, dfd)#

Quantile function of the F distribution: the x with P(F <= x) = p.

Accurate in both tails: for p above 1/2 the upper tail 1 - p, which is exact, is solved instead.

Parameters:
  • p – Probability in [0, 1]; 0 gives 0, 1 gives +inf, values outside give NaN.

  • dfn – Numerator degrees of freedom, dfn > 0; inf gives dfd / F ~ chi-squared(dfd).

  • dfd – Denominator degrees of freedom, dfd > 0; inf gives dfn F ~ chi-squared(dfn).

Returns:

The quantile, as a float or a float64 array of the broadcast shape. Degrees of freedom <= 0 give NaN; a quantile beyond the largest double gives +inf.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.f_sf(x, dfn, dfd)#

Survival function of the F distribution, P(F > x).

Parameters:
  • x – Quantile; x <= 0 gives 1.

  • dfn – Numerator degrees of freedom, dfn > 0; inf gives dfd / F ~ chi-squared(dfd).

  • dfd – Denominator degrees of freedom, dfd > 0; inf gives dfn F ~ chi-squared(dfn).

Returns:

P(F > x), as a float or a float64 array of the broadcast shape. Degrees of freedom <= 0 give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.gamma_cdf(x, shape, scale=1.0)#

Cumulative distribution function of the gamma distribution, P(X <= x).

The density is x^(shape - 1) exp(-x / scale) / (Gamma(shape) scale^shape).

Parameters:
  • x – Quantile; x <= 0 gives 0.

  • shape – Shape, 0 < shape < inf.

  • scale – Scale, 0 < scale < inf.

Returns:

P(X <= x), as a float or a float64 array of the broadcast shape. Parameters outside their domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.gamma_sf(x, shape, scale=1.0)#

Survival function of the gamma distribution, P(X > x).

Parameters:
  • x – Quantile; x <= 0 gives 1.

  • shape – Shape, 0 < shape < inf.

  • scale – Scale, 0 < scale < inf.

Returns:

P(X > x), as a float or a float64 array of the broadcast shape. Parameters outside their domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.hypergeom_logsf(k, M, n, N)#

Logarithm of the survival function of the hypergeometric distribution, log P(X > k).

Parameters:
  • k – Quantile; a non-integer k gives log P(X > floor(k)).

  • M – Population size, a non-negative integer below 2^53.

  • n – Number of successes in the population, an integer with 0 <= n <= M.

  • N – Number of draws, an integer with 0 <= N <= M.

Returns:

log P(X > k), as a float or a float64 array of the broadcast shape. Parameters that are not integers in their ranges give NaN, and so does a tail that needs more than 1e7 terms (near the mode of populations above about 2e13).

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.hypergeom_sf(k, M, n, N)#

Survival function of the hypergeometric distribution, P(X > k).

X counts the successes among N draws without replacement from a population of M items of which n are successes (scipy’s convention).

Parameters:
  • k – Quantile; a non-integer k gives P(X > floor(k)).

  • M – Population size, a non-negative integer below 2^53.

  • n – Number of successes in the population, an integer with 0 <= n <= M.

  • N – Number of draws, an integer with 0 <= N <= M.

Returns:

P(X > k), as a float or a float64 array of the broadcast shape. Parameters that are not integers in their ranges give NaN, and so does a tail that needs more than 1e7 terms (near the mode of populations above about 2e13).

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.norm_cdf(x)#

Cumulative distribution function of the standard normal distribution, P(Z <= x).

Parameters:

x – Quantile.

Returns:

P(Z <= x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.distributions.norm_isf(p)#

Inverse survival function of the standard normal distribution, P(Z > z) = p.

Parameters:

p – Probability in [0, 1]; 0 gives +inf, 1 gives -inf, values outside give NaN.

Returns:

The quantile, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If p is not real-valued.

cytepool.stats.distributions.norm_logcdf(x)#

Logarithm of the standard normal distribution function, log P(Z <= x).

Parameters:

x – Quantile; finite down to x of about -1.9e154.

Returns:

log P(Z <= x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.distributions.norm_logsf(x)#

Logarithm of the standard normal survival function, log P(Z > x).

Parameters:

x – Quantile; finite up to x of about 1.9e154.

Returns:

log P(Z > x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.distributions.norm_ppf(p)#

Quantile function of the standard normal distribution: the z with P(Z <= z) = p.

Parameters:

p – Probability in [0, 1]; 0 gives -inf, 1 gives +inf, values outside give NaN.

Returns:

The quantile, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If p is not real-valued.

cytepool.stats.distributions.norm_sf(x)#

Survival function of the standard normal distribution, P(Z > x).

Parameters:

x – Quantile.

Returns:

P(Z > x), with full relative accuracy in the upper tail, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.distributions.t_cdf(x, df)#

Cumulative distribution function of Student’s t distribution, P(T <= x).

Parameters:
  • x – Quantile.

  • df – Degrees of freedom, df > 0, not necessarily an integer; df = inf gives the standard normal distribution. df <= 0 gives NaN.

Returns:

P(T <= x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.t_isf(p, df)#

Inverse survival function of Student’s t distribution: the t with P(T > t) = p.

Accurate for p down to the smallest doubles; use this rather than t_ppf(1 - p) for upper-tail quantiles. With few degrees of freedom the quantile of a small p can exceed the largest double; the result is then +inf.

Parameters:
  • p – Probability in [0, 1]; 0 gives +inf, 1 gives -inf, values outside give NaN.

  • df – Degrees of freedom, df > 0, not necessarily an integer; df = inf gives the standard normal distribution. df <= 0 gives NaN.

Returns:

The quantile, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.t_logsf(x, df)#

Logarithm of the survival function of Student’s t distribution, log P(T > x).

Parameters:
  • x – Quantile.

  • df – Degrees of freedom, df > 0, not necessarily an integer; df = inf gives the standard normal distribution. df <= 0 gives NaN.

Returns:

log P(T > x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.t_ppf(p, df)#

Quantile function of Student’s t distribution: the t with P(T <= t) = p.

With few degrees of freedom the quantile of a small p can be below minus the largest double; the result is then -inf.

Parameters:
  • p – Probability in [0, 1]; 0 gives -inf, 1 gives +inf, values outside give NaN.

  • df – Degrees of freedom, df > 0, not necessarily an integer; df = inf gives the standard normal distribution. df <= 0 gives NaN.

Returns:

The quantile, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.distributions.t_sf(x, df)#

Survival function of Student’s t distribution, P(T > x).

Parameters:
  • x – Quantile.

  • df – Degrees of freedom, df > 0, not necessarily an integer; df = inf gives the standard normal distribution. df <= 0 gives NaN.

Returns:

P(T > x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

Explained in: Empirical Bayes moderation and TREAT

Special functions for statistics, compiled with numba.

The log-gamma, digamma and trigamma functions, the inverse of the trigamma function, the log-beta function, and the regularized incomplete gamma and beta functions with their complements and logarithms.

Every public function is vectorised: its arguments broadcast with numpy’s rules and the result is a float64 array of the broadcast shape, or a Python float when every argument is a scalar. NaN arguments give NaN, and so do arguments outside a function’s domain; the functions never raise for numeric input. Arguments that are not real numbers (strings, complex numbers, objects) raise TypeError, and so do numpy masked arrays, whose mask numpy would drop: NaN marks a missing value instead.

The compiled scalar routines behind the public functions are named <function>_scalar and can be called from other compiled code.

Methods

  • lgamma is the C library’s lgamma.

  • digamma: a Taylor series about its positive root on [1, 2], which keeps the relative accuracy near the root; the recurrence psi(x + 1) = psi(x) + 1/x maps (0, 10) onto [1, 2]; the asymptotic expansion in Bernoulli numbers for x >= 10; reflection for x < 0 (DLMF 5.5, 5.11, 5.15).

  • trigamma: recurrence up to x >= 12, then the asymptotic expansion; reflection for x < 0 (DLMF 5.15).

  • trigamma_inverse: Newton’s method on 1/trigamma(x), which is convex and close to linear (Smyth 2004). Both starting values lie to the right of the root, so the iteration decreases monotonically to it.

  • log_beta and the prefactors of the incomplete functions use Stirling’s series with Loader’s (2000) deviance form for large arguments, so that no large logarithms cancel.

  • Incomplete gamma: the power series for P (DLMF 8.7.1) and Legendre’s continued fraction for Q (DLMF 8.9.2), each used where it converges fast and computing the smaller tail; the larger tail is its complement in log space. For shapes below 1 and x <= 1.5, Q is computed directly from the alternating series of the lower function, since 1 - P would lose accuracy there. For shapes from 1e3 near the center (0.58 a < x < 1.59 a), Temme’s uniform asymptotic expansion (DLMF 8.12) replaces the series and the fraction, whose step counts grow with the shape.

  • Incomplete beta: the continued fraction DLMF 8.17.22 in its even contraction, evaluated with the modified Lentz method (Lentz 1976; Thompson and Barnett 1986) for I_x(a, b) when x < (a + 1)/(a + b + 2), and for 1 - I_x(a, b) = I_{1-x}(b, a) otherwise. The contraction is written in lambda = a(1 - x) - b x, formed without cancellation, so that it stays accurate when x is close to 1 and a is large (the t and F distributions with many degrees of freedom). For a shape below 1 the tail given by the fraction can be close to 1; the other tail then comes from a term-by-term series of its own rather than from 1 minus it (the analogue of the incomplete gamma case).

Limits:

Accuracy is checked against 50-digit mpmath values (validation/distributions.py) for shapes from 1e-3 to 1e6, incomplete gamma shapes up to 1e18 above the center, and the t and F distributions with degrees of freedom from 0.5 to 1e12 and infinity: relative errors stay below 1e-12 down to probabilities of 1e-300. In deep tails they are of order eps |log p| (about 2e-13 at p = 1e-300), because a probability is evaluated through its logarithm. Larger shapes keep working: the incomplete beta fraction needs a few thousand steps at shapes of 1e12, and the uniform expansion serves the incomplete gamma function near the center for every shape from 1e3. A series or continued fraction that has not converged after 1e7 steps gives NaN rather than an inaccurate value.

References

Olver FWJ, Olde Daalhuis AB, Lozier DW, Schneider BI, Boisvert RF, Clark CW, Miller BR, Saunders BV, Cohl HS, McClain MA, eds (2026). NIST Digital Library of Mathematical Functions, release 1.2.8 of 2026-09-15. https://dlmf.nist.gov/

Dekker TJ (1971). A floating-point technique for extending the available precision. Numerische Mathematik 18, 224-242. doi:10.1007/BF01397083

Lentz WJ (1976). Generating Bessel functions in Mie scattering calculations using continued fractions. Applied Optics 15, 668-671. doi:10.1364/AO.15.000668

Thompson IJ, Barnett AR (1986). Coulomb and Bessel functions of complex arguments and order. Journal of Computational Physics 64, 490-509. doi:10.1016/0021-9991(86)90046-X

Loader C (2000). Fast and accurate computation of binomial probabilities. Technical note, Lucent Technologies, Murray Hill; the version of 5 February 2002 is available at https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf

Maechler M (2012). Accurately computing log(1 - exp(-|a|)) assessed by the Rmpfr package. Vignette of the R package Rmpfr, https://cran.r-project.org/package=Rmpfr

Temme NM (1979). The asymptotic expansion of the incomplete gamma functions. SIAM Journal on Mathematical Analysis 10, 757-766. doi:10.1137/0510071

Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027

cytepool.stats.special.betainc(a, b, x)#

Regularized incomplete beta function I_x(a, b).

Parameters:
  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

  • x – Argument, 0 <= x <= 1.

Returns:

I_x(a, b), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.betaincc(a, b, x)#

Complement of the regularized incomplete beta function, 1 - I_x(a, b).

Equal to I_{1-x}(b, a) and computed directly, so it keeps full relative accuracy when it is small.

Parameters:
  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

  • x – Argument, 0 <= x <= 1.

Returns:

1 - I_x(a, b), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.digamma(x)#

Digamma function psi(x) = d/dx log Gamma(x).

Parameters:

x – Real argument. The poles x = 0, -1, -2, …, where psi has no limit, and x = -inf give NaN; x = +inf gives +inf.

Returns:

psi(x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.special.gammainc(a, x)#

Regularized lower incomplete gamma function P(a, x) = gamma(a, x) / Gamma(a).

Parameters:
  • a – Shape, 0 < a < inf.

  • x – Argument, x >= 0 (x = +inf gives 1).

Returns:

P(a, x), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.gammaincc(a, x)#

Regularized upper incomplete gamma function Q(a, x) = 1 - P(a, x).

Computed directly, so it keeps full relative accuracy when it is small.

Parameters:
  • a – Shape, 0 < a < inf.

  • x – Argument, x >= 0 (x = +inf gives 0).

Returns:

Q(a, x), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.lgamma(x)#

Logarithm of the absolute value of the gamma function, log|Gamma(x)|.

Parameters:

x – Real argument. At the poles x = 0, -1, -2, … the result is +inf; at x = -inf it is NaN.

Returns:

log|Gamma(x)|, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.special.log_beta(a, b)#

Logarithm of the beta function, log B(a, b).

log B(a, b) = log Gamma(a) + log Gamma(b) - log Gamma(a + b), computed without the cancellation of the three log-gamma values for large arguments.

Parameters:
  • a – First parameter, a >= 0.

  • b – Second parameter, b >= 0. If a or b is 0 the result is +inf (NaN if the other is +inf); if one is +inf and the other positive it is -inf. Negative values give NaN.

Returns:

log B(a, b), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.log_betainc(a, b, x)#

Logarithm of the regularized incomplete beta function, log I_x(a, b).

Parameters:
  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

  • x – Argument, 0 <= x <= 1 (x = 0 gives -inf).

Returns:

log I_x(a, b), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.log_betaincc(a, b, x)#

Logarithm of the complement of the incomplete beta function, log(1 - I_x(a, b)).

Parameters:
  • a – First shape, 0 < a < inf.

  • b – Second shape, 0 < b < inf.

  • x – Argument, 0 <= x <= 1 (x = 1 gives -inf).

Returns:

log(1 - I_x(a, b)), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.log_gammainc(a, x)#

Logarithm of the regularized lower incomplete gamma function, log P(a, x).

Finite far beyond the point where P(a, x) underflows.

Parameters:
  • a – Shape, 0 < a < inf.

  • x – Argument, x >= 0 (x = 0 gives -inf).

Returns:

log P(a, x), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.log_gammaincc(a, x)#

Logarithm of the regularized upper incomplete gamma function, log Q(a, x).

Finite far beyond the point where Q(a, x) underflows.

Parameters:
  • a – Shape, 0 < a < inf.

  • x – Argument, x >= 0 (x = +inf gives -inf).

Returns:

log Q(a, x), as a float or a float64 array of the broadcast shape. Arguments outside the domain give NaN.

Raises:

TypeError – If an argument is not real-valued.

cytepool.stats.special.trigamma(x)#

Trigamma function psi’(x), the second derivative of log Gamma(x).

Parameters:

x – Real argument. The poles x = 0, -1, -2, … give +inf (the limit from both sides), x = +inf gives 0 and x = -inf gives NaN.

Returns:

psi’(x), as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If x is not real-valued.

cytepool.stats.special.trigamma_inverse(y)#

Inverse of the trigamma function on x > 0: the x with psi’(x) = y.

Solved by Newton’s method to full double precision (Smyth 2004).

Parameters:

y – Value of the trigamma function, y > 0. y = 0 gives +inf, y = +inf gives 0 and y < 0 gives NaN.

Returns:

x > 0 with trigamma(x) = y, as a float or a float64 array of the broadcast shape.

Raises:

TypeError – If y is not real-valued.