mixle.stats.univariate.continuous.generalized_pareto module

Generalized Pareto distribution (GPD): the peaks-over-threshold law of extreme exceedances.

By the Pickands-Balkema-de Haan theorem the distribution of exceedances over a high threshold of almost any distribution converges to a GPD, which makes it the workhorse for modelling tail risk (hydrology, finance, reliability). With threshold loc = mu, scale sigma > 0 and shape xi, for y = x - mu >= 0:

f(x) = (1/sigma) (1 + xi * y / sigma) ** (-1/xi - 1) (xi != 0), f(x) = (1/sigma) exp(-y / sigma) (xi == 0, the exponential tail).

xi > 0 is a heavy (Pareto) tail, xi = 0 exponential, xi < 0 a tail with a finite upper endpoint at mu - sigma/xi. The threshold mu is treated as a fixed, known level (chosen, not fit – the standard peaks-over-threshold setup); sigma and xi are fit by method of moments, which is closed-form: xi = (1 - m^2/v)/2 and sigma = m (1 - xi) from the exceedance mean m and variance v (valid for xi < 1/2, where the variance is finite).

Reference: Pickands, ‘Statistical inference using extreme order statistics’, Ann. Statist. (1975).

class GeneralizedParetoSuffStat(sum_, sum2_, count_, anchored=None, max_x=None, max_unverified=False)

Bases: tuple

A (sum, sum2, count) sufficient statistic that also carries two side payloads.

Behaves exactly like the plain 3-tuple everywhere it is indexed, unpacked, or iterated (it is one); anchored is the shift-anchored moment payload (anchor, sum_i w_i*(x_i - anchor), sum_i w_i*(x_i - anchor)^2) that GeneralizedParetoAccumulator maintains alongside the raw moments so the moment inversion survives a threshold far from zero. max_x is the largest observation actually folded into these moments with positive weight (None when no such max is known to be trustworthy – see GeneralizedParetoAccumulator.combine()); GeneralizedParetoEstimator.estimate() uses it to keep a shape<0 fit’s implied support consistent with the data it was fit from (T1-01). max_unverified disambiguates WHY max_x is None for a prior-carrying accumulator: it is True when a foreign/pre-payload combine() left real observations with no tracked max, so the T1-01 protection could not be applied though it is supposed to be unconditional there, versus the ordinary False case where the max was never relevant to begin with. estimate reads it to disclose that ambiguity via numerical_repairs() rather than silently falling back to the pre-T1-01 unclamped behaviour (T1-02). Code that doesn’t know about any of these payloads (generic scale_suff_stat, engine kernels, …) sees an ordinary tuple and the estimate falls back to the historical raw path.

Parameters:

• sum_ (float)

• sum2_ (float)

• count_ (float)

• anchored (tuple[float, float, float] | None)

• max_x (float | None)

• max_unverified (bool)

class GeneralizedParetoPriorMoments(mean, second_moment, variance=None, mean_offset=None)

Bases: tuple

The (mean, second_moment) prior pair, carrying the variance and mean offset it was built from.

This family is the only one of its siblings whose pseudo_count prior is encoded as RAW MOMENTS: Gumbel, Student-t and Logistic all take (loc, scale) parameters, from which the prior variance is computed exactly at any magnitude. Recovering it here as second_moment - mean**2 is the same cancellation this module repairs everywhere else, and it is worse than the data path because it is a two-number encoding with no second chance: a prior at a threshold of 1.7e9 recovers a variance of 0.0 from a true 10.4167, and at 1e8 recovers 10.0 – 4% low.

mean has the same disease in the OTHER direction: it is loc plus an O(1) exceedance offset, and forming that sum rounds the offset away to loc’s own ULP once loc is large enough – a threshold of 1e17 already loses the entire offset, well before the variance’s cancellation (which destroys a SQUARED quantity, not an additive one) would even register. mean_offset hangs the exact, loc-invariant mean - loc beside the pair for exactly this reason: captured once while it is still precise (at construction, or carried forward UNCHANGED across a loc retarget rather than re-derived by subtracting a large loc back out of an already-rounded mean), a consumer that needs the prior’s exceedance mean reads it directly instead of differencing.

Keeping the documented pair AS the tuple and hanging the exact variance and mean offset beside it fixes the library’s own GeneralizedParetoDistribution.estimator() path without changing the shape any caller passes: a plain 2-tuple still means exactly what it always meant.

Parameters:

• mean (float)

• second_moment (float)

• variance (float | None)

• mean_offset (float | None)

class GeneralizedParetoDistribution(scale, shape, loc=0.0, name=None, keys=None)

Bases: SequenceEncodableProbabilityDistribution

Generalized Pareto distribution with threshold loc, scale > 0 and shape xi.

Parameters:

• scale (float)

• shape (float)

• loc (float)

• name (str | None)

• keys (str | None)

density(x)

Return the probability density at a single observation.

Parameters:

x (float)

Return type:

float

log_density(x)

Return the log-density at a single observation (-inf outside the support).

Parameters:

x (float)

Return type:

float

seq_log_density(x)

Return vectorized log-density values for sequence-encoded observations.

Parameters:

x (ndarray)

Return type:

ndarray

classmethod compute_capabilities()

Describe backend support for generated generalized-Pareto kernels.

classmethod compute_declaration()

Return the structured compute declaration for generalized Pareto distributions.

static backend_legacy_sufficient_statistics(x, params, engine)

Per-row GPD moment sums in accumulator order (sum, sum2, count).

Parameters:

• x (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, …]

static backend_log_density_from_params(x, scale, shape, loc, engine)

Engine-neutral GPD log-density; the abs(xi) < tol exponential limit is selected per element.

Parameters:

• x (Any)

• scale (Any)

• shape (Any)

• loc (Any)

• engine (Any)

Return type:

Any

backend_seq_log_density(x, engine)

Engine-neutral vectorized log-density for encoded data.

Parameters:

• x (Any)

• engine (Any)

Return type:

Any

classmethod backend_stacked_params(dists, engine)

Stacked GPD parameters for a homogeneous mixture kernel.

Parameters:

• dists (Sequence[GeneralizedParetoDistribution])

• engine (Any)

Return type:

dict[str, Any]

classmethod backend_stacked_log_density(x, params, engine)

Return an (n, k) matrix of GPD log densities.

Parameters:

• x (ndarray)

• params (dict[str, Any])

• engine (Any)

Return type:

Any

classmethod backend_stacked_sufficient_statistics(x, weights, params, engine)

Stacked GPD moment sums (sum, sum2, count) using engine-resident arrays.

Parameters:

• x (ndarray)

• weights (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, Any, Any]

cdf(x)

Cumulative distribution function P(X <= x) (exact).

Parameters:

x (float)

Return type:

float

quantile(q)

Inverse CDF F^{-1}(q).

Parameters:

q (float)

Return type:

float

mean()

Mean loc + scale/(1-xi) for xi < 1, else inf.

Return type:

float

variance()

Variance scale^2 / ((1-xi)^2 (1-2xi)) for xi < 1/2, else inf.

Return type:

float

entropy()

Differential entropy log(scale) + xi + 1.

Return type:

float

sampler(seed=None)

Return a sampler for drawing observations from this distribution.

Parameters:

seed (int | None)

Return type:

GeneralizedParetoSampler

estimator(pseudo_count=None)

Return a method-of-moments estimator for scale and shape at the fixed threshold loc.

Parameters:

pseudo_count (float | None)

Return type:

GeneralizedParetoEstimator

dist_to_encoder()

Return the data encoder used by this distribution for vectorized methods.

Return type:

GeneralizedParetoDataEncoder

class GeneralizedParetoSampler(dist, seed=None)

Bases: DistributionSampler

Draw iid GPD observations by inverse-CDF transform.

Parameters:

• dist (GeneralizedParetoDistribution)

• seed (int | None)

sample(size=None, *, batched=True)

Draw one sample or an array of iid samples by inverse CDF.

Parameters:

• size (int | None)

• batched (bool)

Return type:

float | ndarray

class GeneralizedParetoAccumulator(loc=0.0, name=None, keys=None, has_prior=False)

Bases: AnchoredMomentTrack, SequenceEncodableStatisticAccumulator

Accumulate weighted first and second moments for GPD estimation.

Alongside the raw (sum, sum2, count) the accumulator keeps the conditioning-gated shift-anchored moment track of AnchoredMomentTrack. Peaks-over-threshold data is the case that needs it most: BOTH moments the inversion uses are differences taken at data magnitude – the exceedance mean sum_x/count - loc and the variance E[x^2]-E[x]^2 – so a threshold at epoch seconds turned a shape of +0.16 into a flat 0 and the exponential limit, silently. Well-conditioned data never activates the track and accumulates bit-identically to the historical single-pass path.

Parameters:

• loc (float)

• name (str | None)

• keys (str | None)

• has_prior (bool)

update(x, weight, estimate)

Accumulate weighted first and second moments for one observation.

Parameters:

• x (float)

• weight (float)

• estimate (GeneralizedParetoDistribution | None)

Return type:

None

initialize(x, weight, rng)

Initialize statistics from one observation.

Parameters:

• x (float)

• weight (float)

• rng (RandomState | None)

Return type:

None

seq_update(x, weights, estimate)

Accumulate weighted first and second moments from encoded data.

Parameters:

• x (ndarray)

• weights (ndarray)

• estimate (GeneralizedParetoDistribution | None)

Return type:

None

seq_initialize(x, weights, rng)

Initialize statistics from encoded observations.

Parameters:

• x (ndarray)

• weights (ndarray)

• rng (RandomState | None)

Return type:

None

combine(suff_stat)

Merge another generalized-Pareto sufficient-statistic tuple.

Parameters:

suff_stat (tuple[float, float, float])

Return type:

GeneralizedParetoAccumulator

value()

Return accumulated sum, second moment sum, and count.

A plain 3-tuple for every consumer that treats it as one; once the shift-anchored moment track is live, or a trustworthy observation max is known AND relevant, it is a GeneralizedParetoSuffStat additionally carrying the anchored moments in its .anchored attribute and/or the max in .max_x, so combine() can fold them in and GeneralizedParetoEstimator.estimate() can invert threshold-invariant moments and keep a shape<0 fit’s implied support consistent with the data (T1-01). “Relevant” is always true for a prior-carrying estimator (see _max_x_relevant()); for one without a prior it is the same cheap raw-moment pre-check as before, so ordinary well-conditioned data still returns the historical plain tuple, byte for byte.

For a prior-carrying estimator that protection is supposed to be unconditional, but a foreign (pre-payload) combine() can still leave real observations with no trustworthy max (self._max_tainted); rather than silently degrade to the plain tuple that hides this from estimate, that case is still reported as a GeneralizedParetoSuffStat with max_x=None and max_unverified=True, so estimate can disclose the ambiguity instead of reproducing the T1-01 crash class with zero signal (T1-02).

Return type:

tuple[float, float, float]

from_value(x)

Replace accumulator contents from a sufficient-statistic tuple.

Parameters:

x (tuple[float, float, float])

Return type:

GeneralizedParetoAccumulator

scale(c)

Scale the accumulated statistics in place by c, anchored track included.

The structural default round-trips through value()/from_value() and scale_suff_stat rebuilds a PLAIN tuple, which would drop the payload and undo the repair; loc is a fixed threshold, not a statistic, and is left alone. The tracked max is likewise left alone: uniform weight scaling reweights the aggregate contribution of the same observations, it does not change which values were actually observed.

Parameters:

c (float)

Return type:

GeneralizedParetoAccumulator

acc_to_encoder()

Return the encoder used by this accumulator.

Return type:

GeneralizedParetoDataEncoder

class GeneralizedParetoAccumulatorFactory(loc=0.0, name=None, keys=None, has_prior=False)

Bases: StatisticAccumulatorFactory

Factory for GeneralizedParetoAccumulator.

Parameters:

• loc (float)

• name (str | None)

• keys (str | None)

• has_prior (bool)

make()

Create a fresh generalized-Pareto accumulator.

Return type:

GeneralizedParetoAccumulator

class GeneralizedParetoEstimator(loc=0.0, pseudo_count=None, suff_stat=None, min_scale=1.0e-12, xi_max=0.5 - 1.0e-6, xi_min=-10.0, name=None, keys=None)

Bases: ParameterEstimator

Method-of-moments estimator for GPD scale and shape at a fixed threshold loc.

Both inverted moments – the exceedance mean sum_x/count - loc and the variance E[x^2]-E[x]^2 – are differences taken at DATA magnitude, so a threshold far from zero destroys them: peaks over a threshold at epoch seconds fitted shape=0 (the exponential limit) instead of +0.16, silently. Whenever the accumulated statistics carry the shift-anchored payload (see GeneralizedParetoAccumulator) both moments are formed about a data anchor instead, which makes the fit THRESHOLD-EQUIVARIANT: estimating on x + c with the threshold at loc + c returns the same scale and shape, for any constant c the data can carry. With a plain raw tuple – statistics the conditioning gate never needed to anchor, or ones that arrived already reduced from an engine kernel or an older serialization – the historical raw path is used bit-identically, and estimate warns rather than returning moments it cannot stand behind.

The pseudo_count prior is carried as the raw pair suff_stat = (mean0, second0); a plain 2-tuple has its variance and exceedance-mean offset recovered by differencing (second0 - mean0**2 and mean0 - loc respectively), which is the prior’s own encoding and is unaffected by the anchored track – a prior whose mean sits far from zero loses precision in both reconstructions, so express such a prior in threshold-relative coordinates, or build it with GeneralizedParetoDistribution.estimator(), which carries both exactly (see GeneralizedParetoPriorMoments).

loc is a property, not a bare attribute: reassigning it re-targets this SAME estimator at a new threshold (the “threshold stability” / mean-residual-life re-thresholding workflow) and, when a pseudo_count prior is configured, re-anchors its frozen (mean0, second0) moments to the new origin so the prior keeps meaning what it meant when built instead of being silently reinterpreted relative to the wrong loc (T1-01).

Parameters:

• loc (float)

• pseudo_count (float | None)

• suff_stat (tuple[float, float] | None)

• min_scale (float)

• xi_max (float)

• xi_min (float)

• name (str | None)

• keys (str | None)

property loc: float

The fixed threshold mu observations are exceedances over.

accumulator_factory()

Return an accumulator factory for generalized-Pareto moments.

Return type:

GeneralizedParetoAccumulatorFactory

estimate(nobs, suff_stat)

Estimate scale and shape from exceedance moments at the fixed location.

Parameters:

• nobs (float | None)

• suff_stat (tuple[float, float, float])

Return type:

GeneralizedParetoDistribution

class GeneralizedParetoDataEncoder(loc=0.0, scale=None, shape=None)

Bases: DataSequenceEncoder

Encode GPD observations as a float array.

Parameters:

• loc (float)

• scale (float | None)

• shape (float | None)

seq_encode(x)

Encode observations as a floating-point array.

Parameters:

x (Sequence[float])

Return type:

ndarray