Notebooks · lib/stats/special

The special functions under lib/stats's distributions

Each function takes and returns one number. The methods, and where each is held to mpmath at 50 digits, are in lib/stats/accuracy.py (the grid) and lib/stats/accuracy.json (the worst errors found, beside SciPy's and R's).

log1p(x), expm1(x): log(1 + x) and exp(x) - 1 by Kahan's corrections, log1pmx(t): log(1 + t) - t, by the series in u = t / (2 + t) for |t| < 1/2.
lgamma(x): log |Gamma(x)|. Taylor series about 1 and 2 with the coefficients zeta(k) - 1 (lib/stats/tables.viv), the recurrence down to [3/2, 5/2) below 10, Stirling's series from 10, reflection below 0.
digamma(x): psi(x). The Taylor series about its positive zero within 1/4 of it, the series about 1 and 2, the recurrence, the asymptotic series from 10, reflection below 0.
stirlerr(x): log Gamma(x + 1) - (x + 1/2) log x + x - log sqrt(2 pi), Stirling's series from 10, and below it the recurrence stirlerr(x) = stirlerr(x + 1) + sum_k u^2k / (2k + 1), u = 1 / (2x + 1).
bd0(x, m): x log(x / m) + m - x, Loader's deviance term.
dpois_raw(k, m): m^k e^-m / Gamma(k + 1) for real k, Loader's saddle-point form exp(-stirlerr(k) - bd0(k, m)) / sqrt(2 pi k).
dbinom_raw(k, n, p, q): the binomial term for real k and n, Loader's form (R's dbinom_raw).
erfc(x), erfcx(x) = exp(x^2) erfc(x): Cody's rational function on [1/2, 4] (W. J. Cody, Math. Comp. 23 (1969), as in netlib's specfun CALERF), the continued fraction above 4, and exp(-x^2) split so that the square is exact. erf below 1/2 is its Taylor series (stats_erf).
ndtr(x), ndtri(p): the standard normal distribution function and its inverse (Acklam's approximation refined by Halley steps on ndtr).
erfinv(y), erfcinv(q): inverses by Halley steps on erf or erfc.
gammainc(a, x), gammaincc(a, x): the regularized incomplete gamma functions P and Q. The series for P below x = a + 1, the continued fraction for Q above it, and for a < 1, x <= 1.1 both from one alternating series (P = c (1 - a s), Q = -expm1(log c) + c a s, c = x^a / Gamma(1 + a)). The factor x^a e^-x / Gamma(a) is a dpois_raw.
gammaincinv(a, p), gammainccinv(a, q): their inverses in x, from Wilson-Hilferty's or the small-x start, by safeguarded Halley steps on the smaller tail.
betainc(a, b, x), betaincc(a, b, x): the regularized incomplete beta function and its complement by the continued fraction of DLMF 8.17.22 (modified Lentz, Thompson and Barnett 1986) on the side of the mean where it converges, with the factor x^a (1 - x)^b / B(a, b) in Stirling form.
betaincinv(a, b, p), betainccinv(a, b, q): inverses in x, from the starting value of algorithm AS 109 (Cran, Martin and Thomas 1977), by safeguarded Newton steps.
lbeta(a, b): log B(a, b).

Used with from stats.special import ....

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 54–1058

from stats.tables import ZETA_M1, PSI_ROOT_HI, PSI_ROOT_LO, PSI_ROOT_TAYLOR, EULER_GAMMA, LN_SQRT_2PI, SQRT_PI, SQRT_2PI, SQRT_2

EPS = 2.220446049250313e-16
TINY = 1.0e-300
# iterations a series or continued fraction may take before it stops
MAX_TERMS = 1000000
# a root iteration's steps
MAX_STEPS = 100
# where the series of log1pmx is used
LOG1PMX_SERIES = 0.5
# Stirling's series from here; below, the recurrences
STIRLING_MIN = 10.0
# B_2k / (2k (2k - 1)), k = 1 to 10: Stirling's series of log Gamma
STIRLING_COEF = [0.08333333333333333, -0.002777777777777778, 0.0007936507936507937, -0.0005952380952380953, 0.0008417508417508417, -0.0019175269175269176, 0.00641025641025641, -0.029550653594771242, 0.17964437236883057, -1.3924322169059011]
# B_2k / (2k), k = 1 to 8: the asymptotic series of psi
PSI_ASYM_COEF = [0.08333333333333333, -0.008333333333333333, 0.003968253968253968, -0.004166666666666667, 0.007575757575757576, -0.021092796092796094, 0.08333333333333333, -0.4432598039215686]
# the Taylor series of psi about its zero is used within this distance
PSI_ROOT_RADIUS = 0.25
# exp(-x^2) splits x at a multiple of 2^-20, whose square is exact below 2^6
SPLIT_SCALE = 1048576.0
SPLIT_MAX = 64.0
# erfc: the builtin erf below, Cody's rational function to, the continued fraction above
ERF_DIRECT = 0.5
CODY_MAX = 4.0
# terms of the continued fraction above CODY_MAX (11 reach full precision at x = 4)
ERFC_CF_TERMS = 40.0
CODY_P = [0.564188496988670089, 8.88314979438837594, 66.1191906371416295, 298.635138197400131, 881.952221241769090, 1712.04761263407058, 2051.07837782607147, 1230.33935479799725, 2.15311535474403846e-8]
CODY_Q = [15.7449261107098347, 117.693950891312499, 537.181101862009858, 1621.38957456669019, 3290.79923573345963, 4362.61909014324716, 3439.36767414372164, 1230.33935480374942]
# Acklam's rational approximation of the normal quantile (relative error 1.15e-9)
ACKLAM_A = [-39.69683028665376, 220.9460984245205, -275.9285104469687, 138.3577518672690, -30.66479806614716, 2.506628277459239]
ACKLAM_B = [-54.47609879822406, 161.5858368580409, -155.6989798598866, 66.80131188771972, -13.28068155288572]
ACKLAM_C = [-0.007784894002430293, -0.3223964580411365, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783]
ACKLAM_D = [0.007784695709041462, 0.3224671290700398, 2.445134137142996, 3.754408661907416]
ACKLAM_LOW = 0.02425
NDTRI_STEPS = 3
# the incomplete gamma function's small-a form is used for a < 1, x <= this
GAMMA_SMALL_X = 1.1
# Wilson-Hilferty's start is used from this a
WH_MIN_A = 1.0
NAN = inf - inf

# ------------------------------------------------------------ elementary

#| log(1 + x), accurate for small x (Kahan: log(u) x / (u - 1), u = 1 + x).
def log1p(x) {
    if x == inf {
        return inf
    }
    u = 1.0 + x
    if u == 1.0 {
        return x
    }
    return log(u) * x / (u - 1.0)
}

#| exp(x) - 1, accurate for small x (Kahan: (u - 1) x / log(u), u = e^x).
def expm1(x) {
    u = exp(x)
    if u == 1.0 {
        return x

from stats.tables import ZETA_M1, PSI_ROOT_HI, PSI_ROOT_LO, PSI_ROOT_TAYLOR, EULER_GAMMA, LN_SQRT_2PI, SQRT_PI, SQRT_2PI, SQRT_2

EPS = 2.220446049250313e-16
TINY = 1.0e-300
# iterations a series or continued fraction may take before it stops
MAX_TERMS = 1000000
# a root iteration's steps
MAX_STEPS = 100
# where the series of log1pmx is used
LOG1PMX_SERIES = 0.5
# Stirling's series from here; below, the recurrences
STIRLING_MIN = 10.0
# B_2k / (2k (2k - 1)), k = 1 to 10: Stirling's series of log Gamma
STIRLING_COEF = [0.08333333333333333, -0.002777777777777778, 0.0007936507936507937, -0.0005952380952380953, 0.0008417508417508417, -0.0019175269175269176, 0.00641025641025641, -0.029550653594771242, 0.17964437236883057, -1.3924322169059011]
# B_2k / (2k), k = 1 to 8: the asymptotic series of psi
PSI_ASYM_COEF = [0.08333333333333333, -0.008333333333333333, 0.003968253968253968, -0.004166666666666667, 0.007575757575757576, -0.021092796092796094, 0.08333333333333333, -0.4432598039215686]
# the Taylor series of psi about its zero is used within this distance
PSI_ROOT_RADIUS = 0.25
# exp(-x^2) splits x at a multiple of 2^-20, whose square is exact below 2^6
SPLIT_SCALE = 1048576.0
SPLIT_MAX = 64.0
# erfc: the builtin erf below, Cody's rational function to, the continued fraction above
ERF_DIRECT = 0.5
CODY_MAX = 4.0
# terms of the continued fraction above CODY_MAX (11 reach full precision at x = 4)
ERFC_CF_TERMS = 40.0
CODY_P = [0.564188496988670089, 8.88314979438837594, 66.1191906371416295, 298.635138197400131, 881.952221241769090, 1712.04761263407058, 2051.07837782607147, 1230.33935479799725, 2.15311535474403846e-8]
CODY_Q = [15.7449261107098347, 117.693950891312499, 537.181101862009858, 1621.38957456669019, 3290.79923573345963, 4362.61909014324716, 3439.36767414372164, 1230.33935480374942]
# Acklam's rational approximation of the normal quantile (relative error 1.15e-9)
ACKLAM_A = [-39.69683028665376, 220.9460984245205, -275.9285104469687, 138.3577518672690, -30.66479806614716, 2.506628277459239]
ACKLAM_B = [-54.47609879822406, 161.5858368580409, -155.6989798598866, 66.80131188771972, -13.28068155288572]
ACKLAM_C = [-0.007784894002430293, -0.3223964580411365, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783]
ACKLAM_D = [0.007784695709041462, 0.3224671290700398, 2.445134137142996, 3.754408661907416]
ACKLAM_LOW = 0.02425
NDTRI_STEPS = 3
# the incomplete gamma function's small-a form is used for a < 1, x <= this
GAMMA_SMALL_X = 1.1
# Wilson-Hilferty's start is used from this a
WH_MIN_A = 1.0
NAN = inf - inf

# ------------------------------------------------------------ elementary

#| log(1 + x), accurate for small x (Kahan: log(u) x / (u - 1), u = 1 + x).
def log1p(x) {
    if x == inf {
        return inf
    }
    u = 1.0 + x
    if u == 1.0 {
        return x
    }
    return log(u) * x / (u - 1.0)
}

#| exp(x) - 1, accurate for small x (Kahan: (u - 1) x / log(u), u = e^x).
def expm1(x) {
    u = exp(x)
    if u == 1.0 {
        return x