Extra Distributions¶

extra_distributions ¶

BetaPrime `dataclass` ¶

BetaPrime(a, b, s=1.0)

Bases: Samplable

Beta-Prime (a.k.a. Inverse-Beta, Pearson Type VI) with shapes a>0, b>0 and scale s>0.

PDF: f(x) = 1 / (s * B(a,b)) * (x/s)^(a-1) * (1 + x/s)^(-a-b), x>0 CDF: F(x) = I_{ y }(a,b) with y = x / (x + s) (regularized incomplete beta) PPF: No closed form in general -> monotone numeric inversion. Sampling: If U~Gamma(a,1), V~Gamma(b,1), then X = s * U/V ~ BetaPrime(a,b,s).

Examples:

A ratio of two gamma variates, with tail index b:

>>> beta_prime = BetaPrime(a=2.0, b=3.0, s=1.0)
>>> round(beta_prime.cdf(1.0), 10)
0.6875
>>> round(beta_prime.ppf(0.5), 6)
0.627942

The survival function stays accurate where one minus the distribution function would have run out of digits:

>>> f"{beta_prime.sf(1e6):.6e}"
'3.999985e-18'

pdf ¶

pdf(x)

Density.

Elementary, so it is one NumPy expression. The probabilities below need the incomplete beta, which NumPy does not have, and go one element at a time through :func:~heavytails._array.elementwise.

Source code in heavytails/extra_distributions.py

def pdf(self, x: ArrayLike) -> Any:
    """Density.

    Elementary, so it is one NumPy expression. The probabilities below
    need the incomplete beta, which NumPy does not have, and go one
    element at a time through :func:`~heavytails._array.elementwise`.
    """
    values, scalar = as_array(x)
    a, b, s = self.a, self.b, self.s
    z = select(values > 0.0, values, 1.0) / s
    with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
        density = (
            math.exp(-(math.log(s) + _log_beta(a, b)))
            * z ** (a - 1.0)
            * (1.0 + z) ** (-(a + b))
        )
    return restore(select(values > 0.0, density, 0.0), scalar)

BurrXII `dataclass` ¶

BurrXII(c, k, s=1.0)

Bases: InverseTransformSampling

Burr Type XII with shapes c>0, k>0 and scale s>0.

CDF: F(x) = 1 - (1 + (x/s)^c)(-k), x > 0 PDF: f(x) = (ck/s) * (x/s)^(c-1) * (1 + (x/s)^c)(-k-1) PPF: x = s * ( (1 - u)^(-1/k) - 1 )^(1/c)

Examples:

The tail index is the product c * k, so shape and scale can be traded against each other:

>>> burr = BurrXII(c=2.0, k=1.0, s=1.0)
>>> burr.cdf(1.0)
0.5
>>> round(burr.sf(3.0), 10)
0.1
>>> round(BurrXII(c=1.0, k=2.0, s=1.0).sf(3.0), 6)
0.0625

GeneralizedPareto `dataclass` ¶

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

Bases: InverseTransformSampling

Generalized Pareto Distribution (GPD) with shape xi, scale sigma>0, location mu.

Support

x >= mu if xi >= 0 (heavy-tailed when xi>0) mu <= x <= mu - sigma/xi if xi < 0 (bounded tail; not heavy)

CDF

F(x) = 1 - (1 + xi (x-mu)/sigma)^(-1/xi), valid where bracket > 0

PDF: f(x) = (1/sigma) * (1 + xi z)^(-1/xi - 1), z=(x-mu)/sigma PPF: x = mu + (sigma/xi) * ( (1-u)^(-xi) - 1 ) if xi != 0 x = mu - sigma * ln(1-u) if xi = 0 (exponential limit)

Examples:

The limit law for exceedances over a high threshold, which is what makes peaks-over-threshold work at all:

>>> gpd = GeneralizedPareto(xi=0.5, sigma=1.0, mu=0.0)
>>> round(gpd.sf(1.0), 6)
0.444444
>>> round(gpd.ppf(0.5), 6)
0.828427

Positive xi gives a Pareto tail with index 1 / xi; negative xi bounds it above at mu - sigma / xi, past which nothing falls:

>>> GeneralizedPareto(xi=-0.5, sigma=1.0, mu=0.0).cdf(2.0)
1.0
>>> gpd.cdf(-1.0)
0.0

sf ¶

sf(x)

Survival function, from the power itself rather than 1 - cdf.

The distribution function approaches 1 in the tail, so subtracting it from 1 keeps only the digits the tail has already lost. This is the tail the distribution exists to describe, so it is computed directly.

Source code in heavytails/extra_distributions.py

def sf(self, x: ArrayLike) -> Any:
    """Survival function, from the power itself rather than ``1 - cdf``.

    The distribution function approaches 1 in the tail, so subtracting it
    from 1 keeps only the digits the tail has already lost. This is the
    tail the distribution exists to describe, so it is computed directly.
    """
    values, scalar = as_array(x)
    z = (values - self.mu) / self.sigma
    with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
        if self.xi == 0.0:
            survival = np.exp(-z)
        else:
            survival = np.exp(-np.log1p(self._inner(z)) / self.xi)
    outside = select(values < self.mu, 1.0, 0.0)
    return restore(select(self._valid(values), survival, outside), scalar)

InverseGamma `dataclass` ¶

InverseGamma(alpha, beta)

Bases: Samplable

Inverse-Gamma with shape alpha>0 and scale β>0 (support x>0). PDF: f(x) = β^alpha / Gamma(alpha) * x^{-alpha-1} * exp(-β/x) CDF: F(x) = Q(alpha, β/x) = Gamma(alpha, β/x) / Gamma(alpha) (regularized upper gamma) where Q = 1 - P and P is the regularized lower gamma. Sampling: If G ~ Gamma(alpha, scale=1), then X = β / G has InvGamma(alpha, β).

Examples:

Its tail index is alpha, and its lower tail is the interesting numerical case: the probability there is far too small to reach by subtracting from one.

>>> inverse_gamma = InverseGamma(alpha=2.0, beta=1.0)
>>> round(inverse_gamma.ppf(0.5), 6)
0.595824
>>> f"{inverse_gamma.cdf(0.02):.6e}"
'9.836624e-21'
>>> round(inverse_gamma.sf(1.0), 6)
0.264241

pdf ¶

pdf(x)

Density.

Elementary, so it is a single NumPy expression. The probabilities below are not -- they need the incomplete gamma, which NumPy does not have -- and they go one element at a time through :func:~heavytails._array.elementwise.

Source code in heavytails/extra_distributions.py

def pdf(self, x: ArrayLike) -> Any:
    """Density.

    Elementary, so it is a single NumPy expression. The probabilities
    below are not -- they need the incomplete gamma, which NumPy does not
    have -- and they go one element at a time through
    :func:`~heavytails._array.elementwise`.
    """
    values, scalar = as_array(x)
    a, b = self.alpha, self.beta
    positive = select(values > 0.0, values, 1.0)
    with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
        density = (
            (b**a / math.exp(math.lgamma(a)))
            * positive ** (-a - 1.0)
            * np.exp(-b / positive)
        )
    return restore(select(values > 0.0, density, 0.0), scalar)

LogLogistic `dataclass` ¶

LogLogistic(kappa, lam=1.0)

Bases: InverseTransformSampling

Log-Logistic (Fisk) with shape kappa>0 and scale lambda_>0 (support x>0). CDF: F(x) = 1 / (1 + (lambda_/x)^kappa) = (x^kappa) / (x^kappa + lambda_^kappa) PDF: f(x) = (kappa/lambda_) (x/lambda_)^(kappa-1) / (1 + (x/lambda_)^kappa)2 PPF: x = lambda_ * (u/(1-u))^(1/kappa)

Examples:

The median is lam exactly, whatever the shape:

>>> loglogistic = LogLogistic(kappa=2.0, lam=3.0)
>>> loglogistic.cdf(3.0)
0.5
>>> round(loglogistic.ppf(0.5), 10)
3.0

The tail index is kappa:

>>> round(loglogistic.sf(30.0), 6)
0.009901

Extra Distributions¶