Notebooks · lib/stats/dist

Probability distributions in standard form

Fourteen continuous and five discrete families, each in SciPy's standard form (loc 0, scale 1; vp.stats applies loc and scale in Python) with SciPy's shape parameters. dist_value(fam, what, x, s1, s2, s3) evaluates one of

PDF  LOGPDF  CDF  SF  PPF  ISF

of family fam at one point; the family and the function are the named codes below. Each tail is computed directly (SF is never 1 - CDF where the CDF is near 1), through the functions of lib/stats/special.viv:

normal and lognormal: ndtr, ndtri;
Student t, F, beta: betainc, betaincc and their inverses, with the argument (x^2 / (df + x^2) or df / (df + x^2) for t) chosen so that the smaller one is passed;
chi-squared and gamma: gammainc, gammaincc and their inverses; the gamma density is Loader's dpois_raw (as R's dgamma);
exponential, Weibull, Cauchy, logistic, uniform, Laplace, Pareto: closed forms with log1p and expm1 where they lose digits otherwise;
binomial, Poisson, negative binomial: Loader's saddle-point terms for the mass (R's dbinom, dpois, dnbinom) and the incomplete beta or gamma function for the tails; geometric in closed form; hypergeometric as R's dhyper (three binomial terms) with the tails summed from the far end. Discrete quantiles are the smallest k with CDF(k) >= p, searched from the normal approximation.

dist_sample(fam, s1, s2, s3, n, seed): n draws by inversion, PPF of the uniforms np.random.default_rng(seed).random(n) gives (lib/nprandom's PCG64, bit for bit), so a seed gives the same draws in every engine and in Python.

Used with from stats.dist 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 36–722

from stats.tables import LN_SQRT_2PI, SQRT_2PI
from stats.special import NAN, log1p, expm1, lgamma, lbeta, ndtr, ndtri, exp_mhalfsq, gammainc, gammaincc, gammaincinv, gammainccinv, betainc, betaincc, betaincinv, betainccinv, beta_prefix, beta_delta, dpois_raw, dbinom_raw, stirlerr, bd0
from nprandom.pcg64 import npr_pcg_from_entropy, npr_pcg_double
from nprandom.seedseq import npr_ss_words

# the families
NORM = 0
LOGNORM = 1
T = 2
CHI2 = 3
F = 4
GAMMA = 5
BETA = 6
EXPON = 7
WEIBULL = 8
CAUCHY = 9
LOGISTIC = 10
UNIFORM = 11
LAPLACE = 12
PARETO = 13
BINOM = 20
POISSON = 21
NBINOM = 22
GEOM = 23
HYPERGEOM = 24
# the first discrete family's code
DISCRETE = 20
# the functions
PDF = 0
LOGPDF = 1
CDF = 2
SF = 3
PPF = 4
ISF = 5
# Student t: the argument x^2 / (df + x^2) is passed below this, df / (df + x^2) above
T_SWITCH = 0.5
# a discrete quantile search's steps before it gives up
MAX_SEARCH = 100000000

# ------------------------------------------------------------ continuous

def norm_value(what, x) {
    if what == PDF {
        return exp_mhalfsq(x) / SQRT_2PI
    }
    if what == LOGPDF {
        return -0.5 * x * x - LN_SQRT_2PI
    }
    if what == CDF {
        return ndtr(x)
    }
    if what == SF {
        return ndtr(-x)
    }
    if what == PPF {
        return ndtri(x)
    }
    return -ndtri(x)
}

from stats.tables import LN_SQRT_2PI, SQRT_2PI
from stats.special import NAN, log1p, expm1, lgamma, lbeta, ndtr, ndtri, exp_mhalfsq, gammainc, gammaincc, gammaincinv, gammainccinv, betainc, betaincc, betaincinv, betainccinv, beta_prefix, beta_delta, dpois_raw, dbinom_raw, stirlerr, bd0
from nprandom.pcg64 import npr_pcg_from_entropy, npr_pcg_double
from nprandom.seedseq import npr_ss_words

# the families
NORM = 0
LOGNORM = 1
T = 2
CHI2 = 3
F = 4
GAMMA = 5
BETA = 6
EXPON = 7
WEIBULL = 8
CAUCHY = 9
LOGISTIC = 10
UNIFORM = 11
LAPLACE = 12
PARETO = 13
BINOM = 20
POISSON = 21
NBINOM = 22
GEOM = 23
HYPERGEOM = 24
# the first discrete family's code
DISCRETE = 20
# the functions
PDF = 0
LOGPDF = 1
CDF = 2
SF = 3
PPF = 4
ISF = 5
# Student t: the argument x^2 / (df + x^2) is passed below this, df / (df + x^2) above
T_SWITCH = 0.5
# a discrete quantile search's steps before it gives up
MAX_SEARCH = 100000000

# ------------------------------------------------------------ continuous

def norm_value(what, x) {
    if what == PDF {
        return exp_mhalfsq(x) / SQRT_2PI
    }
    if what == LOGPDF {
        return -0.5 * x * x - LN_SQRT_2PI
    }
    if what == CDF {
        return ndtr(x)
    }
    if what == SF {
        return ndtr(-x)
    }
    if what == PPF {
        return ndtri(x)
    }
    return -ndtri(x)
}