Core Distributions¶

heavy_tails ¶

Cauchy `dataclass` ¶

Cauchy(x0=0.0, gamma=1.0)

Bases: InverseTransformSampling

Cauchy(location x0, scale gamma>0). PDF: f(x) = [1/πgamma] * [1 / (1 + ((x-x0)/gamma)^2)] CDF: F(x) = 0.5 + (1/π) * arctan((x - x0)/gamma) PPF: x = x0 + gamma * tan(π(u - 0.5))

Examples:

Symmetric, and heavy enough to have no mean at all -- the sample mean of Cauchy draws is itself Cauchy, so more data does not help.

>>> cauchy = Cauchy(x0=0.0, gamma=1.0)
>>> cauchy.cdf(0.0)
0.5
>>> round(cauchy.ppf(0.75), 10)
1.0
>>> round(cauchy.pdf(0.0), 6)
0.31831

The tail decays like 1/x, so the survival function times x tends to 1/pi:

>>> round(cauchy.sf(1e6) * 1e6, 5)
0.31831

cdf ¶

cdf(x)

Distribution function.

For very negative z the arctangent approaches -pi/2 and adding one half cancels it away, leaving about seven digits of a probability that is the whole point of the left tail. arctan(-1/z) is the same number with its argument near zero, where it is exact.

Source code in heavytails/heavy_tails.py

def cdf(self, x: ArrayLike) -> Any:
    """Distribution function.

    For very negative ``z`` the arctangent approaches ``-pi/2`` and adding
    one half cancels it away, leaving about seven digits of a probability
    that is the whole point of the left tail. ``arctan(-1/z)`` is the same
    number with its argument near zero, where it is exact.
    """
    values, scalar = as_array(x)
    z = (values - self.x0) / self.gamma
    with np.errstate(divide="ignore", invalid="ignore"):
        lower = np.arctan(-1.0 / z) / math.pi
        upper = 1.0 - np.arctan(1.0 / z) / math.pi
    middle = 0.5 + np.arctan(z) / math.pi
    return restore(select(z < -1.0, lower, select(z > 1.0, upper, middle)), scalar)

pdf ¶

pdf(x)

Density.

Source code in heavytails/heavy_tails.py

def pdf(self, x: ArrayLike) -> Any:
    """Density."""
    values, scalar = as_array(x)
    z = (values - self.x0) / self.gamma
    return restore(1.0 / (math.pi * self.gamma * (1.0 + z * z)), scalar)

ppf ¶

ppf(u)

Quantile function.

Returns inf when the quantile exceeds the float range. At extreme parameters the true value is genuinely not representable, and reporting it is more useful than raising, which aborts a parameter sweep at the first point that overflows.

Uses the cotangent form in the tails. The textbook tan(pi (u - 1/2)) puts its argument next to +-pi/2, where the tangent is arbitrarily steep, so a rounding error of one ulp in the argument becomes an enormous error in the result: at u = 1e-9 it returns -318309868.8 where the answer is -318309886.2, wrong in the eighth digit. cot(pi u) has its argument near zero instead, where it is computed exactly, and the tail is the only part of a Cauchy anyone cares about.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` when the quantile exceeds the float range. At extreme
    parameters the true value is genuinely not representable, and reporting
    it is more useful than raising, which aborts a parameter sweep at the
    first point that overflows.

    Uses the cotangent form in the tails. The textbook
    ``tan(pi (u - 1/2))`` puts its argument next to ``+-pi/2``, where the
    tangent is arbitrarily steep, so a rounding error of one ulp in the
    argument becomes an enormous error in the result: at ``u = 1e-9`` it
    returns -318309868.8 where the answer is -318309886.2, wrong in the
    eighth digit. ``cot(pi u)`` has its argument near zero instead, where
    it is computed exactly, and the tail is the only part of a Cauchy
    anyone cares about.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    with np.errstate(divide="ignore", invalid="ignore"):
        low = self.x0 - self.gamma / np.tan(math.pi * values)
        high = self.x0 + self.gamma / np.tan(math.pi * (1.0 - values))
    middle = self.x0 + self.gamma * np.tan(math.pi * (values - 0.5))
    return restore(
        select(values < 0.25, low, select(values > 0.75, high, middle)),
        scalar,
    )

sf ¶

sf(x)

Survival function 1 - CDF.

Computed as atan(1/z)/pi in the upper tail rather than as 1 - cdf(x). For large z the latter subtracts two numbers that agree to every displayed digit and collapses to exactly zero, whereas atan(1/z) -> 1/z keeps full relative precision.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function 1 - CDF.

    Computed as ``atan(1/z)/pi`` in the upper tail rather than as
    ``1 - cdf(x)``. For large z the latter subtracts two numbers that agree
    to every displayed digit and collapses to exactly zero, whereas
    ``atan(1/z) -> 1/z`` keeps full relative precision.
    """
    values, scalar = as_array(x)
    z = (values - self.x0) / self.gamma
    with np.errstate(divide="ignore", invalid="ignore"):
        positive = np.arctan(1.0 / z) / math.pi
    return restore(select(z > 0.0, positive, 0.5 - np.arctan(z) / math.pi), scalar)

Frechet `dataclass` ¶

Frechet(alpha, s=1.0, m=0.0)

Bases: InverseTransformSampling

Fréchet(alpha, s, m): heavy-tailed extreme-value distribution. Support x > m. alpha>0 (shape), s>0 (scale), m (location). CDF: F(x) = exp( - ((x - m)/s)^(-alpha) ), x>m PDF: f(x) = (alpha/s) * ((x - m)/s)^(-alpha-1) * exp( - ((x - m)/s)^(-alpha) ), x>m PPF: x = m + s * [ -ln(u) ]^{-1/alpha}

Examples:

The limit law for maxima of Pareto-tailed data, so cdf(m + s) is exp(-1) whatever the shape:

>>> round(Frechet(alpha=2.0, s=1.0, m=0.0).cdf(1.0), 6)
0.367879
>>> round(Frechet(alpha=5.0, s=1.0, m=0.0).cdf(1.0), 6)
0.367879
>>> round(Frechet(alpha=2.0, s=1.0, m=0.0).ppf(0.5), 6)
1.201122

ppf ¶

ppf(u)

Quantile function.

Returns inf when the quantile exceeds the float range. At extreme parameters the true value is genuinely not representable, and reporting it is more useful than raising, which aborts a parameter sweep at the first point that overflows.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` when the quantile exceeds the float range. At extreme
    parameters the true value is genuinely not representable, and reporting
    it is more useful than raising, which aborts a parameter sweep at the
    first point that overflows.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    with np.errstate(over="ignore"):
        quantile = self.m + self.s * (-np.log(values)) ** (-1.0 / self.alpha)
    return restore(quantile, scalar)

sf ¶

sf(x)

Survival function 1 - CDF, via -expm1 for tail accuracy.

In the upper tail the exponent tends to zero and exp of it rounds to exactly 1.0, so 1 - cdf(x) would yield 0. -expm1(t) is accurate for small t and preserves the true decay.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function 1 - CDF, via ``-expm1`` for tail accuracy.

    In the upper tail the exponent tends to zero and ``exp`` of it rounds to
    exactly 1.0, so ``1 - cdf(x)`` would yield 0. ``-expm1(t)`` is accurate
    for small t and preserves the true decay.
    """
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        above = -np.expm1(-(((values - self.m) / self.s) ** (-self.alpha)))
    return restore(select(values <= self.m, 1.0, above), scalar)

GEV_Frechet `dataclass` ¶

GEV_Frechet(xi, mu=0.0, sigma=1.0)

Bases: InverseTransformSampling

Generalized Extreme Value (Fréchet-type) with xi>0, mu (loc), sigma>0 (scale). Heavy-tailed when xi>0.

CDF: F(x) = exp( -[1 + xi ( (x-mu)/sigma )]^(-1/xi) ), for 1 + xi (x-mu)/sigma > 0 PDF: f(x) = (1/sigma) * [1 + xi z]^(-1/xi - 1) * exp( -[1 + xi z]^(-1/xi) ), z=(x-mu)/sigma PPF: x = mu + (sigma/xi) * ( (-ln u)^(-xi) - 1 )

Examples:

The generalized extreme value distribution in its heavy-tailed branch, parameterised by xi = 1 / alpha:

>>> gev = GEV_Frechet(xi=0.5, mu=0.0, sigma=1.0)
>>> round(gev.cdf(1.0), 6)
0.64118
>>> round(gev.ppf(0.5), 6)
0.402245

Its support starts at mu - sigma / xi, and nothing falls below it:

>>> gev.cdf(-2.0)
0.0

ppf ¶

ppf(u)

Quantile function.

Returns inf when the quantile exceeds the float range. At extreme parameters the true value is genuinely not representable, and reporting it is more useful than raising, which aborts a parameter sweep at the first point that overflows.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` when the quantile exceeds the float range. At extreme
    parameters the true value is genuinely not representable, and reporting
    it is more useful than raising, which aborts a parameter sweep at the
    first point that overflows.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    with np.errstate(over="ignore"):
        quantile = self.mu + (self.sigma / self.xi) * (
            (-np.log(values)) ** (-self.xi) - 1.0
        )
    return restore(quantile, scalar)

sf ¶

sf(x)

Survival function 1 - CDF, via -expm1 for tail accuracy.

See :meth:Frechet.sf: forming 1 - exp(-t) for small t loses every significant digit, while -expm1(-t) does not.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function 1 - CDF, via ``-expm1`` for tail accuracy.

    See :meth:`Frechet.sf`: forming ``1 - exp(-t)`` for small t loses every
    significant digit, while ``-expm1(-t)`` does not.
    """
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        t = 1.0 + self.xi * (values - self.mu) / self.sigma
        above = -np.expm1(-(t ** (-1.0 / self.xi)))
    return restore(select(t > 0.0, above, 1.0), scalar)

InverseTransformSampling ¶

Bases: Samplable

Sampling by feeding uniforms through the quantile function.

Draws the uniforms in a loop -- the generator is Python's and each draw depends on the one before, so there is nothing to vectorise there -- and then inverts all of them in a single call.

That second half is where the time was going. Inverting 50,000 uniforms one at a time costs about 0.30s, because it is 50,000 separate dispatches through :func:~heavytails._array.as_array and back; inverting them in one call costs about 0.8ms.

The uniforms are the same either way, so a seeded sample is reproducible and its distribution is unchanged. The variates are almost the same: NumPy's vectorised loops for log and pow round differently from its scalar ones, so a few draws land one ULP away from what the scalar path gave. For Weibull(k=0.5) that is 14 values in 20,000, worst case a relative 1.99e-16, with the sample mean unchanged to twelve digits. Tests that pin seeded variates to the last bit will see it; nothing measuring a distribution will.

Examples:

>>> Pareto(alpha=2.0, xm=1.0).rvs(3, seed=7) == Pareto(
...     alpha=2.0, xm=1.0
... ).rvs(3, seed=7)
True

ppf ¶

ppf(u)

The quantile function, which every family mixing this in defines.

Declared here because the batched sampler below is written in terms of it: the mixin is not merely a helper, it is a statement that the family samples by inversion and therefore has a quantile function to invert.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """The quantile function, which every family mixing this in defines.

    Declared here because the batched sampler below is written in terms of
    it: the mixin is not merely a helper, it is a statement that the family
    samples by inversion and therefore has a quantile function to invert.
    """
    raise NotImplementedError

LogNormal `dataclass` ¶

LogNormal(mu=0.0, sigma=1.0)

Bases: Samplable

LogNormal with underlying Normal(mu, sigma^2), sigma>0. PDF: f(x) = [1/(x sigma sqrt(2π))] * exp( -(ln x - mu)^2 / (2sigma^2) ), x>0 CDF: F(x) = 0.5 * [1 + erf( (ln x - mu) / (sigma sqrt(2)) )], x>0

Examples:

The median is exp(mu) exactly:

>>> lognormal = LogNormal(mu=0.0, sigma=1.0)
>>> lognormal.ppf(0.5)
1.0
>>> round(lognormal.cdf(1.0), 10)
0.5

Heavy-tailed but not regularly varying: the tail decays faster than any power, so a tail index estimator applied to it returns something that does not mean what it usually means.

>>> round(lognormal.sf(10.0), 6)
0.010651

cdf ¶

cdf(x)

Distribution function.

Evaluated one point at a time: it needs the complementary error function and NumPy has none, so the loop is the implementation rather than a stand-in. :meth:pdf is elementary and does vectorise.

Source code in heavytails/heavy_tails.py

def cdf(self, x: ArrayLike) -> Any:
    """Distribution function.

    Evaluated one point at a time: it needs the complementary error
    function and NumPy has none, so the loop is the implementation rather
    than a stand-in. :meth:`pdf` is elementary and does vectorise.
    """
    values, scalar = as_array(x)
    return restore(elementwise(self._cdf_one, values), scalar)

pdf ¶

pdf(x)

Density. Elementary, so this one vectorises.

Source code in heavytails/heavy_tails.py

def pdf(self, x: ArrayLike) -> Any:
    """Density. Elementary, so this one vectorises."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        z = (np.log(values) - self.mu) / self.sigma
        density = np.exp(-0.5 * z * z) / (
            values * self.sigma * math.sqrt(2.0 * math.pi)
        )
    return restore(select(values <= 0.0, 0.0, density), scalar)

ppf ¶

ppf(u)

Quantile function.

Returns inf when the quantile exceeds the float range rather than raising. For large mu the answer is genuinely not representable -- the median of LogNormal(mu=1000) is exp(1000) -- and inf is the correct value to report for it.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` when the quantile exceeds the float range rather than
    raising. For large ``mu`` the answer is genuinely not representable --
    the median of ``LogNormal(mu=1000)`` is ``exp(1000)`` -- and ``inf`` is
    the correct value to report for it.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    return restore(elementwise(self._ppf_one, values), scalar)

sf ¶

sf(x)

Survival function 1 - CDF, computed with erfc for tail accuracy.

1 - cdf(x) collapses to exactly zero once cdf(x) rounds to 1.0, which happens well inside the range of interest. erfc is accurate for large arguments and keeps the true decay.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function 1 - CDF, computed with ``erfc`` for tail accuracy.

    ``1 - cdf(x)`` collapses to exactly zero once ``cdf(x)`` rounds to 1.0,
    which happens well inside the range of interest. ``erfc`` is accurate
    for large arguments and keeps the true decay.
    """
    values, scalar = as_array(x)
    return restore(elementwise(self._sf_one, values), scalar)

ParameterError ¶

Bases: ValueError

Raised when distribution parameters are invalid.

Pareto `dataclass` ¶

Pareto(alpha, xm=1.0)

Bases: InverseTransformSampling

Pareto Type I with scale xm>0 and shape alpha>0. PDF: f(x) = alpha x_m^alpha / x^{alpha+1}, x >= x_m CDF: F(x) = 1 - (x_m / x)^alpha PPF: F^{-1}(u) = x_m * (1 - u)^{-1/alpha}

Examples:

The defining property is a straight line on a log-log tail plot: the survival function falls by a factor of 10 ** alpha per decade.

>>> pareto = Pareto(alpha=2.0, xm=1.0)
>>> round(pareto.sf(10.0), 10)
0.01
>>> round(pareto.sf(100.0), 10)
0.0001
>>> pareto.cdf(2.0)
0.75
>>> round(pareto.ppf(0.5), 6)
1.414214

cdf ¶

cdf(x)

Distribution function.

Source code in heavytails/heavy_tails.py

def cdf(self, x: ArrayLike) -> Any:
    """Distribution function."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        below = 1.0 - (self.xm / values) ** self.alpha
    return restore(select(values < self.xm, 0.0, below), scalar)

pdf ¶

pdf(x)

Density. A number in gives a number out, a sequence gives an array.

Source code in heavytails/heavy_tails.py

def pdf(self, x: ArrayLike) -> Any:
    """Density. A number in gives a number out, a sequence gives an array."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        density = self.alpha * self.xm**self.alpha / values ** (self.alpha + 1.0)
    return restore(select(values < self.xm, 0.0, density), scalar)

ppf ¶

ppf(u)

Quantile function.

Returns inf where the quantile exceeds the float range. At extreme parameters the true value is genuinely not representable, and reporting it is more useful than raising, which aborts a parameter sweep at the first point that overflows.

Raises:

Type	Description
`ValueError`	If any probability is outside (0, 1). The whole input is checked before any work starts, so the message can name the offending value instead of failing part way through.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` where the quantile exceeds the float range. At extreme
    parameters the true value is genuinely not representable, and reporting
    it is more useful than raising, which aborts a parameter sweep at the
    first point that overflows.

    Raises:
        ValueError: If any probability is outside (0, 1). The whole input
            is checked before any work starts, so the message can name the
            offending value instead of failing part way through.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    with np.errstate(over="ignore"):
        quantile = self.xm * (1.0 - values) ** (-1.0 / self.alpha)
    return restore(quantile, scalar)

sf ¶

sf(x)

Survival function.

Computed directly rather than as 1 - cdf, which returns exactly zero once the survival probability drops below the spacing of one.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function.

    Computed directly rather than as ``1 - cdf``, which returns exactly
    zero once the survival probability drops below the spacing of one.
    """
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        above = (self.xm / values) ** self.alpha
    return restore(select(values < self.xm, 1.0, above), scalar)

RNG ¶

RNG(seed=None)

Thin wrapper around random.Random for reproducibility and isolation.

Attributes¶

rng : random.Random Underlying random number generator.

Examples:

Seeded, which is what makes a simulation study checkable rather than merely plausible:

>>> RNG(42).uniform_0_1() == RNG(42).uniform_0_1()
True

Source code in heavytails/heavy_tails.py

def __init__(self, seed: int | None = None) -> None:
    self.rng = random.Random(seed)

chisquare ¶

chisquare(df)

χ²(df) via Gamma(k=df/2, θ=2).

Source code in heavytails/heavy_tails.py

def chisquare(self, df: float) -> float:
    """χ²(df) via Gamma(k=df/2, θ=2)."""
    if df <= 0:
        raise ParameterError("Chi-square requires df > 0.")
    return self.gamma(shape_k=df / 2.0, scale_theta=2.0)

gamma ¶

gamma(shape_k, scale_theta=1.0)

X ~ Gamma(k, θ) with k>0, θ>0 using Marsaglia-Tsang (2000). Works for all k>0 (uses boost for k<1).

References¶

G. Marsaglia and W. W. Tsang (2000). A Simple Method for Generating Gamma Variables. ACM Transactions on Mathematical Software 26(3):363-372.

Source code in heavytails/heavy_tails.py

def gamma(self, shape_k: float, scale_theta: float = 1.0) -> float:
    """
    X ~ Gamma(k, θ) with k>0, θ>0 using Marsaglia-Tsang (2000).
    Works for all k>0 (uses boost for k<1).

    References
    ----------
    G. Marsaglia and W. W. Tsang (2000). A Simple Method for Generating Gamma Variables.
    ACM Transactions on Mathematical Software 26(3):363-372.
    """
    if not (shape_k > 0 and scale_theta > 0):
        raise ParameterError("Gamma requires shape k>0 and scale θ>0.")

    k = shape_k
    if k < 1.0:
        # Boost: sample from Gamma(k+1, 1) then * U^(1/k)
        x = self._gamma_mt(k + 1.0)
        u = self.uniform_0_1()
        return float(scale_theta * (x * (u ** (1.0 / k))))
    else:
        return float(scale_theta * self._gamma_mt(k))

standard_normal ¶

standard_normal()

Z ~ N(0,1). Uses Python stdlib Box-Muller via random.gauss (Ziggurat internally).

Source code in heavytails/heavy_tails.py

def standard_normal(self) -> float:
    """Z ~ N(0,1). Uses Python stdlib Box-Muller via random.gauss (Ziggurat internally)."""
    return self.rng.gauss(0.0, 1.0)

uniform_0_1 ¶

uniform_0_1()

U ~ Uniform(0,1) in (0,1), clipped away from exact 0 and 1 for log/ppf stability.

Source code in heavytails/heavy_tails.py

def uniform_0_1(self) -> float:
    """U ~ Uniform(0,1) in (0,1), clipped away from exact 0 and 1 for log/ppf stability."""
    # Avoid 0 and 1 to prevent log(0) or tan(pi*(U-0.5)) exploding from exactly 0.5
    u = self.rng.random()
    eps = 1e-16
    return min(max(u, eps), 1.0 - eps)

Samplable ¶

Mixin to provide vectorized sampling with a given RNG.

Examples:

Every distribution here draws through it, so the same seed gives the same sample:

>>> Pareto(alpha=2.0, xm=1.0).rvs(3, seed=1) == Pareto(
...     alpha=2.0, xm=1.0
... ).rvs(3, seed=1)
True

rvs ¶

rvs(n, seed=None)

Draw n IID variates. Subclasses implement ._rvs_one(rng), or ._rvs_many(rng, n) when the whole sample can be drawn at once.

Source code in heavytails/heavy_tails.py

def rvs(self, n: int, seed: int | None = None) -> list[float]:
    """
    Draw n IID variates. Subclasses implement ._rvs_one(rng), or
    ._rvs_many(rng, n) when the whole sample can be drawn at once.
    """
    if not isinstance(n, int) or n <= 0:
        raise ValueError("n must be a positive integer.")
    rng = RNG(seed)
    return self._rvs_many(rng, n)

StudentT `dataclass` ¶

StudentT(nu)

Bases: Samplable

Student's t with degrees of freedom nu>0. PDF: f(x) = Gamma((nu+1)/2) / [ sqrt(nuπ) Gamma(nu/2) ] * (1 + x^2/nu)(-(nu+1)/2) Sampling: X = Z / sqrt(Y/nu) with Z~N(0,1), Y~χ²(nu)

The CDF, survival function and quantile function are expressed through the regularized incomplete beta function in heavytails._special, so no third-party dependency is required.

Examples:

The tail index is nu, so moments below it exist and the rest do not. One degree of freedom is the Cauchy exactly:

>>> round(StudentT(nu=1.0).cdf(1.0), 10)
0.75
>>> round(Cauchy(x0=0.0, gamma=1.0).cdf(1.0), 10)
0.75

More degrees of freedom means a lighter tail:

>>> round(StudentT(nu=1.0).sf(10.0), 6)
0.031726
>>> round(StudentT(nu=4.0).sf(10.0), 6)
0.000281

cdf ¶

cdf(x)

Distribution function, via the regularized incomplete beta.

Evaluated one point at a time. NumPy has no incomplete beta, so there is no expression to vectorise and the loop is the implementation rather than a stand-in for one. :meth:pdf is elementary and does vectorise.

Source code in heavytails/heavy_tails.py

def cdf(self, x: ArrayLike) -> Any:
    """Distribution function, via the regularized incomplete beta.

    Evaluated one point at a time. NumPy has no incomplete beta, so there
    is no expression to vectorise and the loop is the implementation
    rather than a stand-in for one. :meth:`pdf` is elementary and does
    vectorise.
    """
    values, scalar = as_array(x)
    return restore(elementwise(self._cdf_one, values), scalar)

pdf ¶

pdf(x)

Density. Elementary, so this one vectorises.

Source code in heavytails/heavy_tails.py

def pdf(self, x: ArrayLike) -> Any:
    """Density. Elementary, so this one vectorises."""
    values, scalar = as_array(x)
    nu = self.nu
    constant = math.gamma((nu + 1.0) / 2.0) / (
        math.sqrt(nu * math.pi) * math.gamma(nu / 2.0)
    )
    density = constant * (1.0 + values * values / nu) ** (-(nu + 1.0) / 2.0)
    return restore(density, scalar)

ppf ¶

ppf(u)

Quantile function, obtained by inverting the incomplete beta.

Uses the symmetry of the Student-t about zero so that only the upper half is solved, then inverts I_y(nu/2, 1/2) = 2(1-u) directly. Inverting in y rather than in x is what keeps the far tail accurate: a quantile solver working to an absolute tolerance on the CDF loses most of its digits where the density is vanishingly small.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function, obtained by inverting the incomplete beta.

    Uses the symmetry of the Student-t about zero so that only the upper
    half is solved, then inverts ``I_y(nu/2, 1/2) = 2(1-u)`` directly.
    Inverting in ``y`` rather than in ``x`` is what keeps the far tail
    accurate: a quantile solver working to an absolute tolerance on the CDF
    loses most of its digits where the density is vanishingly small.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    return restore(elementwise(self._ppf_one, values), scalar)

sf ¶

sf(x)

Survival function, computed directly for tail accuracy.

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function, computed directly for tail accuracy."""
    values, scalar = as_array(x)
    return restore(elementwise(self._sf_one, values), scalar)

Weibull `dataclass` ¶

Weibull(k, lam=1.0)

Bases: InverseTransformSampling

Weibull(k, lambda_) with shape k>0 and scale lambda_>0. PDF: f(x) = (k/lambda_) (x/lambda_)^{k-1} exp(-(x/lambda_)^k), x>=0 CDF: F(x) = 1 - exp(-(x/lambda_)^k), x>=0 PPF: x = lambda_ * (-ln(1-u))^{1/k} Heavy-tailed for k in (0,1) (subexponential, slower than exponential decay).

Examples:

Shape one is the exponential distribution, whose survival function at the scale is 1/e:

>>> round(Weibull(k=1.0, lam=1.0).sf(1.0), 6)
0.367879

Below one it is heavy-tailed, and the density is unbounded at the origin:

>>> Weibull(k=0.5, lam=1.0).pdf(0.0)
inf
>>> round(Weibull(k=0.5, lam=1.0).sf(4.0), 6)
0.135335

cdf ¶

cdf(x)

Distribution function, via -expm1 so the lower tail survives.

Source code in heavytails/heavy_tails.py

def cdf(self, x: ArrayLike) -> Any:
    """Distribution function, via ``-expm1`` so the lower tail survives."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        below = -np.expm1(-((values / self.lam) ** self.k))
    return restore(select(values < 0.0, 0.0, below), scalar)

pdf ¶

pdf(x)

Density, which diverges at the origin for shape below one.

Source code in heavytails/heavy_tails.py

def pdf(self, x: ArrayLike) -> Any:
    """Density, which diverges at the origin for shape below one."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        z = (values / self.lam) ** self.k
        density = (
            (self.k / self.lam) * (values / self.lam) ** (self.k - 1.0) * np.exp(-z)
        )
    if self.k < 1.0:
        # The limit at the origin. The expression above is 0 ** negative
        # there, which is an error rather than a statement about the
        # density.
        density = select(values == 0.0, np.inf, density)
    return restore(select(values < 0.0, 0.0, density), scalar)

ppf ¶

ppf(u)

Quantile function.

Returns inf when the quantile exceeds the float range. At extreme parameters the true value is genuinely not representable, and reporting it is more useful than raising, which aborts a parameter sweep at the first point that overflows.

Source code in heavytails/heavy_tails.py

def ppf(self, u: ArrayLike) -> Any:
    """Quantile function.

    Returns ``inf`` when the quantile exceeds the float range. At extreme
    parameters the true value is genuinely not representable, and reporting
    it is more useful than raising, which aborts a parameter sweep at the
    first point that overflows.
    """
    values, scalar = as_array(u)
    check_probabilities(values)
    with np.errstate(over="ignore"):
        quantile = self.lam * (-np.log1p(-values)) ** (1.0 / self.k)
    return restore(quantile, scalar)

sf ¶

sf(x)

Survival function: 1 - CDF(x).

Source code in heavytails/heavy_tails.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function: 1 - CDF(x)."""
    values, scalar = as_array(x)
    with np.errstate(divide="ignore", invalid="ignore"):
        above = np.exp(-((values / self.lam) ** self.k))
    return restore(select(values < 0.0, 1.0, above), scalar)

Core Distributions¶

heavy_tails ¶

Cauchy dataclass ¶

cdf ¶

pdf ¶

ppf ¶

sf ¶

Frechet dataclass ¶

ppf ¶

sf ¶

GEV_Frechet dataclass ¶

ppf ¶

sf ¶

InverseTransformSampling ¶

ppf ¶

LogNormal dataclass ¶

cdf ¶

pdf ¶

ppf ¶

sf ¶

ParameterError ¶

Pareto dataclass ¶

cdf ¶

pdf ¶

ppf ¶

sf ¶

RNG ¶

Attributes¶

chisquare ¶

gamma ¶

References¶

standard_normal ¶

uniform_0_1 ¶

Samplable ¶

rvs ¶

StudentT dataclass ¶

cdf ¶

pdf ¶

ppf ¶

sf ¶

Weibull dataclass ¶

cdf ¶

pdf ¶

ppf ¶

sf ¶

Cauchy `dataclass` ¶

Frechet `dataclass` ¶

GEV_Frechet `dataclass` ¶

LogNormal `dataclass` ¶

Pareto `dataclass` ¶

StudentT `dataclass` ¶

Weibull `dataclass` ¶