Distributions

NumPyro-native extreme value distributions. Every class subclasses numpyro.distributions.Distribution, so it plugs directly into MCMC / SVI / Predictive and is differentiable end to end. Beyond sample / log_prob, each carries the EVT toolkit: cdf, icdf, survival_function, mean, variance, entropy, and return_level.

All of them delegate their math to the primitives — reach for those directly if you want a density without a distribution object.

Generalized Extreme Value

The distribution for block maxima (annual temperature maxima, seasonal peak discharge). The concentration \(\xi\) selects the tail type continuously, so inference can learn which regime the data is in instead of you committing up front — this is the default choice unless you have reason to fix \(\xi\).

NumPyro-compatible extreme value distributions.

GeneralizedExtremeValueDistribution(loc: ArrayLike = 0.0, scale: ArrayLike = 1.0, concentration: ArrayLike | None = None, shape: ArrayLike | None = None, validate_args: bool | None = None)

Bases: Distribution

Generalized Extreme Value Distribution (GEVD) for NumPyro.

The GEVD is a family of continuous probability distributions developed within extreme value theory. It encompasses three types of extreme value distributions:

• Type I (Gumbel): ξ = 0, exponential tails
• Type II (Fréchet): ξ > 0, heavy polynomial tails
• Type III (Weibull): ξ < 0, bounded upper tail

Probability Density Function:

For ξ ≠ 0: f(x) = (1/σ) * t(x)^(-(1/ξ + 1)) * exp(-t(x)^(-1/ξ)) where t(x) = 1 + ξ(x - μ)/σ

For ξ = 0 (Gumbel limit): f(x) = (1/σ) * exp(-(x - μ)/σ) * exp(-exp(-(x - μ)/σ))

Cumulative Distribution Function:

For ξ ≠ 0: F(x) = exp(-t(x)^(-1/ξ)) where t(x) = 1 + ξ(x - μ)/σ

For ξ = 0: F(x) = exp(-exp(-(x - μ)/σ))

Support:

• ξ > 0: x ≥ μ - σ/ξ (lower bounded)
• ξ = 0: x ∈ ℝ (unbounded)
• ξ < 0: x ≤ μ - σ/ξ (upper bounded)

Parameters:

Name	Type	Description	Default
loc	float	Location parameter μ ∈ ℝ	0.0
scale	float	Scale parameter σ > 0	1.0
shape	float	Shape parameter ξ ∈ ℝ	None

Examples:

>>> import jax
>>> import jax.numpy as jnp
>>>
>>> # Create GEVD instances for each type
>>> gumbel = GeneralizedExtremeValueDistribution(loc=0, scale=1, shape=0)
>>> frechet = GeneralizedExtremeValueDistribution(loc=0, scale=1, shape=0.2)
>>> weibull = GeneralizedExtremeValueDistribution(loc=0, scale=1, shape=-0.2)
>>>
>>> # Sample and evaluate
>>> key = jax.random.PRNGKey(42)
>>> samples = frechet.sample(key, sample_shape=(1000,))
>>> log_probs = frechet.log_prob(samples)

Initialize the Generalized Extreme Value Distribution.

Parameters:

Name	Type	Description	Default
loc	ArrayLike	Location parameter μ (real number)	0.0
scale	ArrayLike	Scale parameter σ (positive real number)	1.0
concentration	ArrayLike | None	Shape parameter ξ (real number)	None
shape	ArrayLike | None	Deprecated backward-compatible alias for concentration	None
validate_args	bool | None	Whether to validate input arguments	None

Raises:

Type	Description
ValueError	If scale <= 0

Source code in src/xtremax/distributions/gevd.py

def __init__(
    self,
    loc: ArrayLike = 0.0,
    scale: ArrayLike = 1.0,
    concentration: ArrayLike | None = None,
    shape: ArrayLike | None = None,
    validate_args: bool | None = None,
):
    """
    Initialize the Generalized Extreme Value Distribution.

    Args:
        loc: Location parameter μ (real number)
        scale: Scale parameter σ (positive real number)
        concentration: Shape parameter ξ (real number)
        shape: Deprecated backward-compatible alias for ``concentration``
        validate_args: Whether to validate input arguments

    Raises:
        ValueError: If scale <= 0
    """
    if shape is not None:
        if concentration is not None:
            raise ValueError(
                "Pass only one of 'concentration' or the deprecated 'shape' alias."
            )
        warnings.warn(
            "'shape' is deprecated; use 'concentration' instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        concentration = shape

    if concentration is None:
        concentration = 0.0

    self.loc, self.scale, self.concentration = promote_shapes(
        loc, scale, concentration
    )

    # Determine batch shape from broadcasted parameters
    batch_shape = lax.broadcast_shapes(
        jnp.shape(self.loc), jnp.shape(self.scale), jnp.shape(self.concentration)
    )

    super().__init__(batch_shape=batch_shape, validate_args=validate_args)

support: constraints.Constraint property

Return the support of the distribution.

Returns:

Type	Description
Constraint	Constraint object reflecting the shape-dependent support:
Constraint	ξ > 0: [μ - σ/ξ, +∞) → greater_than_eq(lower_bound) (closed, matching the old interval membership; the density vanishes at the endpoint so log_prob is -inf)
Constraint	ξ = 0: (-∞, +∞) → real
Constraint	ξ < 0: (-∞, μ - σ/ξ] → less_than_eq(upper_bound) (closed: the endpoint density is finite at ξ = -1 and diverges for ξ < -1, so the endpoint is a genuine support point)

The constraint kind is dispatched on the statically known sign of the concentration so that biject_to(support) returns a bijector that genuinely handles the half-unbounded geometry — the previous interval(lo, hi) construction mapped to Sigmoid ∘ Affine which yields inf/nan when a bound is infinite, breaking every NUTS/SVI model that used this distribution as a latent site.

When the concentration is traced (inside jit/inference) or a mixed-sign batch, the sign is not statically known and the support falls back to constraints.real: membership then accepts points outside the true parameter-dependent support, but log_prob returns -inf there, so inference remains correct (init simply rejects such draws).

mean: jnp.ndarray property

Mean. Thin wrapper for :func:~xtremax.primitives.gev.gev_mean.

mode: jnp.ndarray property

Compute the mode of the distribution.

For ξ ≠ 0: mode = μ + (σ/ξ) * ((1+ξ)^(-ξ) - 1)

This interior stationary point only exists for ξ > -1. For ξ ≤ -1 (Weibull branch), the density increases all the way to the finite upper endpoint, so the mode is the upper bound itself.

For ξ = 0: mode = μ

Returns:

Type	Description
ndarray	Mode of the distribution

variance: jnp.ndarray property

Compute the variance of the distribution.

The variance exists when ξ < 1/2:

For ξ ≠ 0, ξ < 1/2: Var[X] = (σ²/ξ²) * (Γ(1-2ξ) - Γ²(1-ξ))

For ξ = 0: Var[X] = σ² * π²/6

Thin wrapper for :func:~xtremax.primitives._moments.gev_variance, which evaluates the Γ-difference through a cancellation-free reduced form — the direct subtraction above loses every significant digit for small ξ.

Returns:

Type	Description
ndarray	Variance, or +∞ when it doesn't exist (ξ ≥ 1/2)

sample(key: jnp.ndarray, sample_shape: tuple = (), *, shape: tuple | None = None) -> jnp.ndarray

Generate samples from the GEVD using inverse transform sampling.

The sampling uses the quantile function (inverse CDF):

For ξ ≠ 0: Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

For ξ = 0: Q(p) = μ - σ * ln(-ln(p))

Parameters:

Name	Type	Description	Default
key	ndarray	JAX random key for sampling	required
sample_shape	tuple	Shape of samples to generate	()
shape	tuple | None	Keyword-only alias for sample_shape matching the pipekit_cycle.ObservationNoise protocol's parameter name (sample(key, shape)), so callers typed against that protocol can pass it by keyword. Passing both a non-empty sample_shape and shape raises ValueError.	None

Returns:

Type	Description
ndarray	Array of samples from the GEVD

Source code in src/xtremax/distributions/gevd.py

def sample(
    self,
    key: jnp.ndarray,
    sample_shape: tuple = (),
    *,
    shape: tuple | None = None,
) -> jnp.ndarray:
    """
    Generate samples from the GEVD using inverse transform sampling.

    The sampling uses the quantile function (inverse CDF):

    For ξ ≠ 0:
        Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

    For ξ = 0:
        Q(p) = μ - σ * ln(-ln(p))

    Args:
        key: JAX random key for sampling
        sample_shape: Shape of samples to generate
        shape: Keyword-only alias for ``sample_shape`` matching the
            ``pipekit_cycle.ObservationNoise`` protocol's parameter
            name (``sample(key, shape)``), so callers typed against
            that protocol can pass it by keyword. Passing both a
            non-empty ``sample_shape`` and ``shape`` raises
            ``ValueError``.

    Returns:
        Array of samples from the GEVD
    """
    if shape is not None:
        if sample_shape != ():
            raise ValueError(
                "Pass only one of 'sample_shape' (NumPyro spelling) or "
                "'shape' (ObservationNoise protocol spelling), not both."
            )
        sample_shape = shape
    check_prng_key(key)
    extended_shape = sample_shape + self.batch_shape

    # JAX's Uniform(0, 1) sampler can emit exact 0 or 1 at the
    # endpoints; passing those to icdf yields -inf/+inf and poisons
    # downstream computations. Clamp away from the endpoints.
    uniform_samples = dist.Uniform(0.0, 1.0).sample(key, extended_shape)
    eps = jnp.finfo(uniform_samples.dtype).eps
    uniform_samples = jnp.clip(uniform_samples, eps, 1.0 - eps)

    # Apply inverse CDF transformation
    return self.icdf(uniform_samples)

log_prob(value: jnp.ndarray) -> jnp.ndarray

Log PDF. Thin wrapper for :func:~xtremax.primitives.gev.gev_log_prob.

Source code in src/xtremax/distributions/gevd.py

@validate_sample
def log_prob(self, value: jnp.ndarray) -> jnp.ndarray:
    """Log PDF. Thin wrapper for :func:`~xtremax.primitives.gev.gev_log_prob`."""
    return gev_log_prob(value, self.loc, self.scale, self.concentration)

cdf(value: jnp.ndarray) -> jnp.ndarray

CDF. Thin wrapper for :func:~xtremax.primitives.gev.gev_cdf.

Source code in src/xtremax/distributions/gevd.py

def cdf(self, value: jnp.ndarray) -> jnp.ndarray:
    """CDF. Thin wrapper for :func:`~xtremax.primitives.gev.gev_cdf`."""
    return gev_cdf(value, self.loc, self.scale, self.concentration)

icdf(q: jnp.ndarray) -> jnp.ndarray

Quantile function. Thin wrapper for gev_icdf.

Source code in src/xtremax/distributions/gevd.py

def icdf(self, q: jnp.ndarray) -> jnp.ndarray:
    """Quantile function. Thin wrapper for ``gev_icdf``."""
    return gev_icdf(q, self.loc, self.scale, self.concentration)

upper_bound() -> jnp.ndarray

Compute the upper bound of the support.

Returns:

Type	Description
ndarray	Upper bound: μ - σ/ξ for ξ < 0, +∞ otherwise

Source code in src/xtremax/distributions/gevd.py

def upper_bound(self) -> jnp.ndarray:
    """
    Compute the upper bound of the support.

    Returns:
        Upper bound: μ - σ/ξ for ξ < 0, +∞ otherwise
    """
    loc = jnp.asarray(self.loc)
    scale = jnp.asarray(self.scale)
    shape = jnp.asarray(self.concentration)
    # Replace ξ=0 with a safe placeholder (1.0) inside the division
    # so plain Python inputs don't raise `ZeroDivisionError` before
    # `jnp.where` selects the +∞ branch. The non-finite branch value
    # is masked out.
    safe_shape = jnp.where(shape == 0.0, 1.0, shape)
    bounded = loc - scale / safe_shape
    return jnp.where(shape < 0, bounded, jnp.inf)

lower_bound() -> jnp.ndarray

Compute the lower bound of the support.

Returns:

Type	Description
ndarray	Lower bound: μ - σ/ξ for ξ > 0, -∞ otherwise

Source code in src/xtremax/distributions/gevd.py

def lower_bound(self) -> jnp.ndarray:
    """
    Compute the lower bound of the support.

    Returns:
        Lower bound: μ - σ/ξ for ξ > 0, -∞ otherwise
    """
    loc = jnp.asarray(self.loc)
    scale = jnp.asarray(self.scale)
    shape = jnp.asarray(self.concentration)
    safe_shape = jnp.where(shape == 0.0, 1.0, shape)
    bounded = loc - scale / safe_shape
    return jnp.where(shape > 0, bounded, -jnp.inf)

covariance() -> jnp.ndarray

Marginal variance broadcast to the batch shape.

Structural seam for pipekit_cycle.ObservationNoise (which requires covariance() and sample(key, shape)): together with the existing :meth:sample, this lets the distribution act as a non-Gaussian observation-error model in data-assimilation experiments — without xtremax importing pipekit (see docs/interop.md). The observation errors are treated as independent per batch element, so the "covariance" is the marginal variance vector rather than a dense matrix; consuming analysis steps interpret it as a diagonal.

Returns:

Type	Description
ndarray	Variance broadcast to batch_shape (+inf where the
ndarray	variance does not exist, i.e. ξ >= 1/2).

Source code in src/xtremax/distributions/gevd.py

def covariance(self) -> jnp.ndarray:
    """Marginal variance broadcast to the batch shape.

    Structural seam for ``pipekit_cycle.ObservationNoise`` (which
    requires ``covariance()`` and ``sample(key, shape)``): together
    with the existing :meth:`sample`, this lets the distribution act
    as a non-Gaussian observation-error model in data-assimilation
    experiments — without xtremax importing pipekit (see
    ``docs/interop.md``). The observation errors are treated as
    independent per batch element, so the "covariance" is the
    marginal variance vector rather than a dense matrix; consuming
    analysis steps interpret it as a diagonal.

    Returns:
        Variance broadcast to ``batch_shape`` (``+inf`` where the
        variance does not exist, i.e. ξ >= 1/2).
    """
    return jnp.broadcast_to(jnp.asarray(self.variance), self.batch_shape)

kurtosis() -> jnp.ndarray

Compute the excess kurtosis of the distribution.

Excess kurtosis exists when ξ < 1/4. The formula involves fourth-order moments and gamma functions:

κ = μ₄/σ⁴ - 3

where μ₄ is the fourth central moment.

Thin wrapper for :func:~xtremax.primitives._moments.gev_excess_kurtosis, which evaluates the fourth-moment gamma combination through a cancellation-free reduced form (it loses two leading orders in the naive one) and reduces to the Gumbel value 12/5 at ξ = 0.

Returns:

Type	Description
ndarray	Excess kurtosis, or NaN when it doesn't exist (ξ ≥ 1/4)

Source code in src/xtremax/distributions/gevd.py

def kurtosis(self) -> jnp.ndarray:
    """
    Compute the excess kurtosis of the distribution.

    Excess kurtosis exists when ξ < 1/4. The formula involves fourth-order
    moments and gamma functions:

    κ = μ₄/σ⁴ - 3

    where μ₄ is the fourth central moment.

    Thin wrapper for
    :func:`~xtremax.primitives._moments.gev_excess_kurtosis`, which
    evaluates the fourth-moment gamma combination through a
    cancellation-free reduced form (it loses *two* leading orders in
    the naive one) and reduces to the Gumbel value 12/5 at ξ = 0.

    Returns:
        Excess kurtosis, or NaN when it doesn't exist (ξ ≥ 1/4)
    """
    return gev_excess_kurtosis(self.concentration)

skew() -> jnp.ndarray

Compute the skewness of the distribution.

Skewness exists when ξ < 1/3:

For Gumbel (ξ = 0): skew ≈ 1.1396 For general GEVD: involves third-order moments and gamma functions

Thin wrapper for :func:~xtremax.primitives._moments.gev_skewness, which evaluates the third-moment gamma combination through a cancellation-free reduced form. The reduced ratio carries the sign of ξ on its own, so no explicit sign(ξ) factor is needed.

Returns:

Type	Description
ndarray	Skewness, or NaN when it doesn't exist (ξ ≥ 1/3)

Source code in src/xtremax/distributions/gevd.py

def skew(self) -> jnp.ndarray:
    """
    Compute the skewness of the distribution.

    Skewness exists when ξ < 1/3:

    For Gumbel (ξ = 0): skew ≈ 1.1396
    For general GEVD: involves third-order moments and gamma functions

    Thin wrapper for :func:`~xtremax.primitives._moments.gev_skewness`,
    which evaluates the third-moment gamma combination through a
    cancellation-free reduced form. The reduced ratio carries the sign
    of ξ on its own, so no explicit ``sign(ξ)`` factor is needed.

    Returns:
        Skewness, or NaN when it doesn't exist (ξ ≥ 1/3)
    """
    return gev_skewness(self.concentration)

entropy() -> jnp.ndarray

Differential entropy of the GEV distribution (in nats).

For any :math:\xi (including the Gumbel limit),

.. math:: H = \log\sigma + 1 + \gamma(1 + \xi)

where :math:\gamma is the Euler–Mascheroni constant. At :math:\xi = 0 this correctly reduces to :math:\log\sigma + 1 + \gamma (the Gumbel entropy), so no separate branch is needed.

Source code in src/xtremax/distributions/gevd.py

def entropy(self) -> jnp.ndarray:
    r"""Differential entropy of the GEV distribution (in nats).

    For any :math:`\xi` (including the Gumbel limit),

    .. math:: H = \log\sigma + 1 + \gamma(1 + \xi)

    where :math:`\gamma` is the Euler–Mascheroni constant. At
    :math:`\xi = 0` this correctly reduces to :math:`\log\sigma + 1 + \gamma`
    (the Gumbel entropy), so no separate branch is needed.
    """
    scale, shape = self.scale, self.concentration
    euler_gamma = 0.5772156649015329
    return jnp.log(scale) + 1.0 + euler_gamma * (1.0 + shape)

survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the survival function S(x) = 1 - F(x).

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the survival function	required

Returns:

Type	Description
ndarray	Survival probabilities

Source code in src/xtremax/distributions/gevd.py

def survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the survival function S(x) = 1 - F(x).

    Args:
        value: Points at which to evaluate the survival function

    Returns:
        Survival probabilities
    """
    # Delegate to the stable primitive so the class and functional APIs
    # share the smooth ξ→0 handling (no threshold to drift against).
    return gev_survival(value, self.loc, self.scale, self.concentration)

log_survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the log survival function log S(x) = log( 1 - F(x) ).

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the survival function	required

Returns:

Type	Description
ndarray	Survival probabilities

Source code in src/xtremax/distributions/gevd.py

def log_survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the log survival function log S(x) = log( 1 - F(x) ).

    Args:
        value: Points at which to evaluate the survival function

    Returns:
        Survival probabilities
    """
    # Delegate to the stable primitive (smooth ξ→0, matches the class cdf).
    return gev_log_survival(value, self.loc, self.scale, self.concentration)

hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the hazard rate h(x) = f(x) / S(x).

The hazard rate represents the instantaneous rate of failure given survival up to time x.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the hazard rate	required

Returns:

Type	Description
ndarray	Hazard rate values

Source code in src/xtremax/distributions/gevd.py

def hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the hazard rate h(x) = f(x) / S(x).

    The hazard rate represents the instantaneous rate of failure
    given survival up to time x.

    Args:
        value: Points at which to evaluate the hazard rate

    Returns:
        Hazard rate values
    """
    return jnp.exp(self.log_prob(value) - self.log_survival_function(value))

cumulative_hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the cumulative hazard rate Λ(x) = -log(S(x)).

This represents the accumulated hazard up to time x.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the cumulative hazard rate	required

Returns:

Type	Description
ndarray	Cumulative hazard rate values

Source code in src/xtremax/distributions/gevd.py

def cumulative_hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the cumulative hazard rate Λ(x) = -log(S(x)).

    This represents the accumulated hazard up to time x.

    Args:
        value: Points at which to evaluate the cumulative hazard rate

    Returns:
        Cumulative hazard rate values
    """
    return -self.log_survival_function(value)

return_level(return_period: float | jnp.ndarray) -> jnp.ndarray

Return level. Thin wrapper for gev_return_level.

Source code in src/xtremax/distributions/gevd.py

def return_level(self, return_period: float | jnp.ndarray) -> jnp.ndarray:
    """Return level. Thin wrapper for ``gev_return_level``."""
    return gev_return_level(return_period, self.loc, self.scale, self.concentration)

tail_index() -> jnp.ndarray

Compute the tail index α = 1/ξ for ξ > 0.

For Fréchet distributions (ξ > 0), the tail follows a power law: P(X > x) ~ x^(-1/ξ) as x → ∞

Returns:

Type	Description
ndarray	Tail index (1/ξ for ξ > 0, ∞ for ξ ≤ 0)

Source code in src/xtremax/distributions/gevd.py

def tail_index(self) -> jnp.ndarray:
    """
    Compute the tail index α = 1/ξ for ξ > 0.

    For Fréchet distributions (ξ > 0), the tail follows a power law:
    P(X > x) ~ x^(-1/ξ) as x → ∞

    Returns:
        Tail index (1/ξ for ξ > 0, ∞ for ξ ≤ 0)
    """
    shape = jnp.asarray(self.concentration)
    # Safe divisor: with Python-float parameters the division happens
    # eagerly in both jnp.where branches, so ξ = 0 raised
    # ZeroDivisionError before the mask applied.
    safe_shape = jnp.where(shape > 0, shape, 1.0)
    return jnp.where(shape > 0, 1.0 / safe_shape, jnp.inf)

exceedance_probability(threshold: jnp.ndarray) -> jnp.ndarray

Compute probability of exceeding a threshold: P(X > threshold).

This is equivalent to the survival function at the threshold.

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value	required

Returns:

Type	Description
ndarray	Exceedance probabilities

Source code in src/xtremax/distributions/gevd.py

def exceedance_probability(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute probability of exceeding a threshold: P(X > threshold).

    This is equivalent to the survival function at the threshold.

    Args:
        threshold: Threshold value

    Returns:
        Exceedance probabilities
    """
    return self.survival_function(threshold)

conditional_excess_mean(threshold: jnp.ndarray) -> jnp.ndarray

Mean excess :math:E[X - u \mid X > u] for the GEV distribution.

Computed via quantile-space quadrature:

.. math:: E[X \mid X > u] = \frac{\int_{F(u)}^{1} F^{-1}(p)\, dp}{1 - F(u)}

using a per-threshold trapezoidal rule on a grid from :math:F(u) to :math:1 - \epsilon. This is exact up to quadrature error across all three GEV branches (Gumbel, Fréchet, Weibull).

The previous implementation used the GPD linear form :math:(\sigma + \xi(u - \mu)) / (1 - \xi), which is the Pickands–Balkema–de Haan POT limit — correct only as :math:u \to upper endpoint, and notably returning a constant :math:\sigma for all :math:u when :math:\xi = 0 even though the true Gumbel mean excess depends on :math:u.

Returns NaN where the mean does not exist (:math:\xi \ge 1).

Source code in src/xtremax/distributions/gevd.py

def conditional_excess_mean(self, threshold: jnp.ndarray) -> jnp.ndarray:
    r"""Mean excess :math:`E[X - u \mid X > u]` for the GEV distribution.

    Computed via quantile-space quadrature:

    .. math:: E[X \mid X > u]
        = \frac{\int_{F(u)}^{1} F^{-1}(p)\, dp}{1 - F(u)}

    using a per-threshold trapezoidal rule on a grid from
    :math:`F(u)` to :math:`1 - \epsilon`. This is exact up to
    quadrature error across all three GEV branches (Gumbel,
    Fréchet, Weibull).

    The previous implementation used the GPD linear form
    :math:`(\sigma + \xi(u - \mu)) / (1 - \xi)`, which is the
    *Pickands–Balkema–de Haan POT limit* — correct only as
    :math:`u \to` upper endpoint, and notably returning a constant
    :math:`\sigma` for all :math:`u` when :math:`\xi = 0` even
    though the true Gumbel mean excess depends on :math:`u`.

    Returns NaN where the mean does not exist (:math:`\xi \ge 1`).
    """
    threshold_arr = jnp.asarray(threshold)
    shape = self.concentration

    # Log-space tail-probability quadrature:
    #   ∫_{F(u)}^1 F⁻¹(p) dp = ∫_0^{S(u)} F⁻¹(1-y) dy
    #                        = ∫_{-∞}^{log S(u)} F⁻¹(1-eᵛ) · eᵛ dv
    # with v = log y. We truncate the lower end at log(1e-12),
    # which is negligible for any ξ < 1. This replaces the earlier
    # linear-p grid that had two problems: (a) the normaliser used
    # the truncated mass `(1-ε) - F(u)` instead of `S(u)`, and (b)
    # the fixed linear grid systematically underresolved heavy
    # tails where most of the conditional mass lives near p = 1.
    n_grid = 1024
    p0 = self.cdf(threshold_arr)
    # Joint batch of threshold and distribution parameters — the grid
    # axis is appended after this, so every per-batch quantity must be
    # broadcast to it before gaining the trailing axis.
    batch = jnp.shape(p0)
    s_u = 1.0 - p0
    # Floor S(u) so log is finite; the NaN mask below drops results
    # whose threshold is effectively at or beyond the upper support.
    # Also cap just below 1.0 so the integration grid stays strictly
    # inside the open interval (0, 1): ``p_grid`` would otherwise
    # reach 0 when ``S(u) = 1`` (e.g. u far below the location),
    # and ``icdf(0)`` diverges to -∞ for unbounded-below branches.
    s_u_safe = jnp.clip(s_u, 1e-12, 1.0 - 1e-6)
    log_s_u = jnp.log(s_u_safe)
    # Lower endpoint of the log-y integration range: step ~20 e-folds
    # below log_s_u so the grid always captures the bulk of the
    # tail-integrand mass AND stays strictly ascending (``v_grid``
    # dx > 0). Without this, a fixed floor like log(1e-6) would run
    # backward when S(u) < 1e-6 and trapezoid would sign-flip the
    # result. The mass below the truncation point — which behaves
    # like y^{1-ξ} and reached ~14% of the total at ξ = 0.9 — is
    # added back analytically below rather than by widening the
    # span: a 1/(1-ξ)-scaled span would push ``exp(v_grid)`` under
    # the float32 exp underflow (~e⁻¹⁰³) and produce 0·inf = NaN in
    # the integrand. The ``-jnp.log1p(-y_grid)`` computation below
    # is accurate down to float32 subnormals.
    span = jnp.asarray(20.0, dtype=log_s_u.dtype)
    log_y_min = log_s_u - span

    unit = jnp.linspace(0.0, 1.0, n_grid)
    log_s_u_exp = jnp.expand_dims(log_s_u, axis=-1)
    log_y_min_exp = jnp.expand_dims(jnp.broadcast_to(log_y_min, batch), axis=-1)
    # v grid ascending from log_y_min up to log S(u). If log_y_min
    # >= log_s_u (untrustworthy regime; mask below returns NaN),
    # the grid is degenerate but numerically bounded.
    v_grid = log_y_min_exp * (1.0 - unit) + log_s_u_exp * unit
    y_grid = jnp.exp(v_grid)
    # Compute ``x = F⁻¹(1 - y)`` directly from ``y`` using
    # ``-log1p(-y)`` instead of going through ``self.icdf(1 - y)``.
    # For ``y`` near ``eps(float32)`` the round-tripped ``p = 1 - y``
    # in float32 only retains ~1 ulp of precision, and the
    # downstream ``log(p)`` in ``gev_icdf`` loses that — icdf comes
    # out finite but wildly wrong and the quadrature sign-flips.
    # ``log1p(-y)`` is designed for exactly this small-y regime.
    # Parameters are broadcast to the joint batch and given the same
    # trailing grid axis as ``neg_log_q`` — evaluating batched
    # (..., n_grid) grids against unexpanded (batch,) parameters
    # raised a broadcasting error.
    neg_log_q = -jnp.log1p(-y_grid)
    is_gumbel = jnp.abs(shape) < GUMBEL_THRESHOLD
    xi = jnp.where(is_gumbel, 1.0, shape)
    loc_e = jnp.expand_dims(jnp.broadcast_to(self.loc, batch), axis=-1)
    scale_e = jnp.expand_dims(jnp.broadcast_to(self.scale, batch), axis=-1)
    xi_e = jnp.expand_dims(jnp.broadcast_to(xi, batch), axis=-1)
    gum_e = jnp.expand_dims(jnp.broadcast_to(is_gumbel, batch), axis=-1)
    gumbel_x = loc_e - scale_e * jnp.log(neg_log_q)
    gev_x = loc_e + (scale_e / xi_e) * (jnp.power(neg_log_q, -xi_e) - 1.0)
    x_grid = jnp.where(gum_e, gumbel_x, gev_x)

    integrand = x_grid * y_grid  # F⁻¹(1-y) · y, multiplied by dv=dy/y
    numerator = jnp.trapezoid(integrand, x=v_grid, axis=-1)

    # Closed-form remainder for the truncated (0, y_c] tail with
    # y_c = exp(log_y_min) = S(u)·e⁻²⁰. There -log(1-y) ≈ y to ~y_c
    # relative accuracy, so
    #   ξ ≠ 0: x(y) ≈ (μ - σ/ξ) + (σ/ξ)·y^{-ξ}
    #          → ∫₀^{y_c} x dy = (μ - σ/ξ)·y_c + (σ/ξ)·y_c^{1-ξ}/(1-ξ)
    #   ξ = 0: x(y) ≈ μ - σ·log y
    #          → ∫₀^{y_c} x dy = y_c·(μ - σ·(log y_c - 1))
    # For light tails the remainder is negligible; at ξ = 0.9 it is
    # ~14% of the numerator — dropping it was the truncation bias.
    loc_b = jnp.broadcast_to(self.loc, batch)
    scale_b = jnp.broadcast_to(self.scale, batch)
    xi_b = jnp.broadcast_to(xi, batch)
    gum_b = jnp.broadcast_to(is_gumbel, batch)
    shape_b = jnp.broadcast_to(shape, batch)
    y_c = jnp.exp(log_y_min)
    # Sanitize the discarded branch: 1-ξ ≤ 0 (ξ ≥ 1, masked to NaN
    # below) and the Gumbel branch would otherwise divide by zero.
    one_minus_xi = jnp.where(gum_b | (shape_b >= 1.0), 1.0, 1.0 - xi_b)
    gev_rem = (loc_b - scale_b / xi_b) * y_c + (scale_b / xi_b) * jnp.power(
        y_c, one_minus_xi
    ) / one_minus_xi
    gumbel_rem = y_c * (loc_b - scale_b * (log_y_min - 1.0))
    numerator = numerator + jnp.where(gum_b, gumbel_rem, gev_rem)

    mean_conditional = numerator / s_u_safe
    mean_excess = mean_conditional - threshold_arr

    mean_exists = shape < 1.0
    valid = mean_exists & (s_u > 1e-12)
    return jnp.where(valid, mean_excess, jnp.nan)

reliability_function(time: jnp.ndarray) -> jnp.ndarray

Compute the reliability function R(t) = P(X > t).

This is equivalent to the survival function but emphasizes the reliability engineering interpretation for Weibull Type III.

Parameters:

Name	Type	Description	Default
time	ndarray	Time points for reliability evaluation	required

Returns:

Type	Description
ndarray	Reliability values (probability of surviving past time t)

Source code in src/xtremax/distributions/gevd.py

def reliability_function(self, time: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the reliability function R(t) = P(X > t).

    This is equivalent to the survival function but emphasizes
    the reliability engineering interpretation for Weibull Type III.

    Args:
        time: Time points for reliability evaluation

    Returns:
        Reliability values (probability of surviving past time t)
    """
    return self.survival_function(time)

mean_residual_life(time: jnp.ndarray) -> jnp.ndarray

Compute the mean residual life function: E[X - t | X > t].

For reliability applications, this represents the expected additional lifetime given survival to time t.

Parameters:

Name	Type	Description	Default
time	ndarray	Current time/age	required

Returns:

Type	Description
ndarray	Mean residual life values

Source code in src/xtremax/distributions/gevd.py

def mean_residual_life(self, time: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the mean residual life function: E[X - t | X > t].

    For reliability applications, this represents the expected
    additional lifetime given survival to time t.

    Args:
        time: Current time/age

    Returns:
        Mean residual life values
    """
    return self.conditional_excess_mean(time)

percentile_residual_life(time: jnp.ndarray, percentile: float) -> jnp.ndarray

Compute percentile residual life: (Q(p|X>t) - t).

This gives the additional time until the p-th percentile failure, conditional on survival to time t.

Parameters:

Name	Type	Description	Default
time	ndarray	Current time/age	required
percentile	float	Percentile level (0 < percentile < 1)	required

Returns:

Type	Description
ndarray	Percentile residual life values

Source code in src/xtremax/distributions/gevd.py

def percentile_residual_life(
    self, time: jnp.ndarray, percentile: float
) -> jnp.ndarray:
    """
    Compute percentile residual life: (Q(p|X>t) - t).

    This gives the additional time until the p-th percentile failure,
    conditional on survival to time t.

    Args:
        time: Current time/age
        percentile: Percentile level (0 < percentile < 1)

    Returns:
        Percentile residual life values
    """
    # Conditional CDF: F(x | X > t) = p ⟺ F(x) = 1 - (1 - p) S(t).
    # (Previously used `1 - p S(t)`, which mapped p=0 to F^{-1}(1),
    # i.e. the upper endpoint, instead of returning zero residual.)
    survival_prob = self.survival_function(time)
    total_prob = 1.0 - (1.0 - percentile) * survival_prob
    return self.icdf(total_prob) - time

expand(batch_shape: tuple[int, ...]) -> dist.Distribution

Expand to batch_shape by reconstructing via __init__.

Going through the constructor (rather than an ExpandedDistribution wrapper) keeps every EVT method available on the returned instance.

Source code in src/xtremax/distributions/gevd.py

def expand(self, batch_shape: tuple[int, ...]) -> dist.Distribution:
    """Expand to ``batch_shape`` by reconstructing via ``__init__``.

    Going through the constructor (rather than an ``ExpandedDistribution``
    wrapper) keeps every EVT method available on the returned instance.
    """
    batch_shape = tuple(batch_shape)
    if batch_shape == self.batch_shape:
        return self
    return type(self)(
        loc=jnp.broadcast_to(self.loc, batch_shape),
        scale=jnp.broadcast_to(self.scale, batch_shape),
        concentration=jnp.broadcast_to(self.concentration, batch_shape),
        validate_args=self._validate_args,
    )

Ragged records: masking the likelihood

Station networks rarely share a calendar. Annual maxima over \(S\) stations and \(T\) years arrive as an \((S, T)\) grid with gaps where a station had no record that year, and the likelihood should only see the observed station-years.

Use NumPyro's masking distribution — fill the gaps with any in-support value (the per-station mean works) and let the mask discard them:

import jax.numpy as jnp
import numpyro
from xtremax import GeneralizedExtremeValueDistribution as GEV

filled = jnp.where(mask, maxima, station_mean[:, None])  # (S, T)
gev = GEV(loc=mu[:, None], scale=sigma[:, None], concentration=xi[:, None])
numpyro.sample("obs", gev.mask(mask), obs=filled)

numpyro.handlers.mask(mask=mask) around a plain sample is equivalent, as is summing log_prob yourself into a numpyro.factor:

lp = gev.log_prob(filled)
numpyro.factor("obs", jnp.where(mask, lp, 0.0).sum())

All three score identically. Prefer the first: factor does register an observed site, but an auxiliary one carrying Unit(log_factor) rather than the GEV, so everything downstream sees a summed scalar where the other two keep the distribution. log_likelihood returns one total per posterior draw instead of the per-observation matrix WAIC and PSIS-LOO need, and Predictive has no distribution left to draw the site from.

Why the fill value does not bias the fit

Masked entries contribute neither density nor gradient, so what you write into the gaps cannot move the posterior. This holds even when the filler falls outside the support: the GEV endpoint is parameter-dependent (\(1 + \xi(y-\mu)/\sigma > 0\)), so a per-station mean can sit past a \(\xi < 0\) upper bound at the parameters the sampler happens to be holding.

Genuinely observed maxima are a different matter — one outside the support at a trial parameter draw scores \(-\infty\), and find_valid_initial_params rejects that draw (it wants a finite potential and finite gradients) and tries another, up to 100 times. nan fails the same isfinite test, so the distinction is not which one gets retried but whether retrying can help: an out-of-support observation is \(-\infty\) only at some parameters, so a later draw lands somewhere valid. A nan leaking out of a masked-out gap is nan at every parameter, and no number of retries recovers from that.

Generalized Pareto

The distribution for threshold exceedances (peaks-over-threshold). Pair it with the extraction threshold and declustering utilities: POT uses more of the data than block maxima, at the cost of a threshold choice.

NumPyro-compatible extreme value distributions.

GeneralizedParetoDistribution(scale: ArrayLike = 1.0, concentration: ArrayLike | None = None, shape: ArrayLike | None = None, validate_args: bool | None = None)

Bases: Distribution

Generalized Pareto Distribution (GPD) for NumPyro.

The Generalized Pareto Distribution is fundamental in extreme value theory for modeling threshold exceedances in the Peaks-Over-Threshold (POT) framework. It emerges as the limiting distribution of scaled excesses over high thresholds for a wide class of underlying distributions.

Key Characteristics: - Models excesses above a threshold: Y = X - u | X > u - Three families based on shape parameter ξ: * ξ > 0: Pareto-type (heavy tails, power-law decay) * ξ = 0: Exponential-type (exponential tails) * ξ < 0: Beta-type (bounded support, light tails) - Direct connection to GEVD via POT theory - Foundation for threshold-based extreme value modeling

Probability Density Function:

For ξ ≠ 0: f(x) = (1/σ) * (1 + ξx/σ)^(-(1/ξ + 1))

For ξ = 0 (exponential limit): f(x) = (1/σ) * exp(-x/σ)

Cumulative Distribution Function:

For ξ ≠ 0: F(x) = 1 - (1 + ξx/σ)^(-1/ξ)

For ξ = 0: F(x) = 1 - exp(-x/σ)

Support:

• ξ ≥ 0: x ≥ 0 (non-negative)
• ξ < 0: 0 ≤ x ≤ -σ/ξ (bounded above)

Parameters:

• scale (σ): Scale parameter > 0
• shape (ξ): Shape parameter ∈ ℝ

Connection to GEVD:

If block maxima follow GEVD(μ, σ*, ξ), then threshold excesses follow GPD(σ + ξ(u - μ), ξ) where u is the threshold.

Applications:

• Financial risk: Value-at-Risk, Expected Shortfall modeling
• Hydrology: Flood frequency analysis above design levels
• Insurance: Large claim modeling, catastrophe reinsurance
• Engineering: Structural reliability, extreme load analysis
• Environmental: Pollution exceedances, extreme weather events
• Telecommunications: Network traffic bursts, service failures

Examples:

>>> import jax
>>> import jax.numpy as jnp
>>>
>>> # Heavy-tailed excesses (financial losses)
>>> heavy_tail_gpd = GeneralizedParetoDistribution(scale=1000, shape=0.25)
>>>
>>> # Exponential excesses (reliability applications)
>>> exponential_gpd = GeneralizedParetoDistribution(scale=500, shape=0.0)
>>>
>>> # Bounded excesses (physical constraints)
>>> bounded_gpd = GeneralizedParetoDistribution(scale=100, shape=-0.2)
>>>
>>> # Key properties
>>> print(f"Mean (heavy tail): {heavy_tail_gpd.mean}")
>>> print(f"Upper bound (bounded): {bounded_gpd.upper_bound()}")
>>> print(f"Tail index (heavy): {heavy_tail_gpd.tail_index()}")
>>>
>>> # Sample and analyze
>>> key = jax.random.PRNGKey(42)
>>> samples = heavy_tail_gpd.sample(key, sample_shape=(1000,))
>>> log_probs = heavy_tail_gpd.log_prob(samples)

Initialize the Generalized Pareto Distribution.

Parameters:

Name	Type	Description	Default
scale	ArrayLike	Scale parameter σ > 0	1.0
concentration	ArrayLike | None	Shape parameter ξ ∈ ℝ * ξ > 0: Heavy tails (Pareto-type) * ξ = 0: Exponential tails * ξ < 0: Light tails, bounded support (Beta-type)	None
shape	ArrayLike | None	Deprecated backward-compatible alias for concentration. Stored as self.concentration; the NumPyro-inherited Distribution.shape() method is kept callable.	None
validate_args	bool | None	Whether to validate input arguments	None

Raises:

Type	Description
ValueError	If scale <= 0

Source code in src/xtremax/distributions/gpd.py

def __init__(
    self,
    scale: ArrayLike = 1.0,
    concentration: ArrayLike | None = None,
    shape: ArrayLike | None = None,
    validate_args: bool | None = None,
):
    """
    Initialize the Generalized Pareto Distribution.

    Args:
        scale: Scale parameter σ > 0
        concentration: Shape parameter ξ ∈ ℝ
              * ξ > 0: Heavy tails (Pareto-type)
              * ξ = 0: Exponential tails
              * ξ < 0: Light tails, bounded support (Beta-type)
        shape: Deprecated backward-compatible alias for ``concentration``.
            Stored as ``self.concentration``; the NumPyro-inherited
            ``Distribution.shape()`` method is kept callable.
        validate_args: Whether to validate input arguments

    Raises:
        ValueError: If scale <= 0
    """
    if shape is not None:
        if concentration is not None:
            raise ValueError(
                "Pass only one of 'concentration' or the deprecated 'shape' alias."
            )
        warnings.warn(
            "'shape' is deprecated; use 'concentration' instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        concentration = shape

    if concentration is None:
        concentration = 0.0

    self.scale, self.concentration = promote_shapes(scale, concentration)

    # Determine batch shape from broadcasted parameters
    batch_shape = lax.broadcast_shapes(
        jnp.shape(self.scale), jnp.shape(self.concentration)
    )

    super().__init__(batch_shape=batch_shape, validate_args=validate_args)

support: constraints.Constraint property

Return the support constraint for GPD.

Returns:

Type	Description
Constraint	Constraint object reflecting the shape-dependent support:
Constraint	ξ ≥ 0: [0, +∞) → nonnegative
Constraint	ξ < 0: [0, -σ/ξ] → interval(0, upper_bound)

Only the ξ < 0 case uses interval — its bounds are finite, so biject_to maps it to a well-behaved Sigmoid ∘ Affine. For ξ ≥ 0 the old interval(0, inf) construction made biject_to return inf, breaking NUTS/SVI models with a GPD latent site; constraints.nonnegative has a registered exp-bijector and accepts the closed endpoint x = 0.

When the concentration is traced (inside jit/inference) or a mixed-sign batch, the support falls back to nonnegative — the lower bound is parameter-independent, and points above a finite ξ < 0 upper endpoint are handled by log_prob returning -inf.

mean: jnp.ndarray property

Mean. Thin wrapper for :func:~xtremax.primitives.gpd.gpd_mean.

mode: jnp.ndarray property

Compute the mode of the GPD.

For ξ > -1 the density is monotone decreasing from x = 0, so the mode is 0 (the threshold). For ξ < -1 the density increases toward the finite upper endpoint -σ/ξ, so the mode is the upper endpoint itself (at ξ = -1, the uniform case, every point is modal and 0 is returned by convention).

Returns:

Name	Type	Description
Mode	ndarray	0 for ξ ≥ -1, the upper endpoint -σ/ξ for ξ < -1

variance: jnp.ndarray property

Compute the variance of the GPD.

Variance exists when ξ < 1/2:

For ξ < 1/2: Var[X] = σ² / ((1-ξ)²(1-2ξ))

For ξ ≥ 1/2: Var[X] = +∞ (infinite variance)

Returns:

Type	Description
ndarray	Variance or +∞ when it doesn't exist

sample(key: jnp.ndarray, sample_shape: tuple = (), *, shape: tuple | None = None) -> jnp.ndarray

Generate samples from the GPD using inverse transform sampling.

The sampling uses the quantile function:

For ξ ≠ 0: Q(p) = (σ/ξ) * ((1-p)^(-ξ) - 1)

For ξ = 0: Q(p) = -σ * ln(1-p)

Parameters:

Name	Type	Description	Default
key	ndarray	JAX random key for sampling	required
sample_shape	tuple	Shape of samples to generate	()
shape	tuple | None	Keyword-only alias for sample_shape matching the pipekit_cycle.ObservationNoise protocol's parameter name (sample(key, shape)), so callers typed against that protocol can pass it by keyword. Passing both a non-empty sample_shape and shape raises ValueError.	None

Returns:

Type	Description
ndarray	Array of samples from the GPD (all within support)

Source code in src/xtremax/distributions/gpd.py

def sample(
    self,
    key: jnp.ndarray,
    sample_shape: tuple = (),
    *,
    shape: tuple | None = None,
) -> jnp.ndarray:
    """
    Generate samples from the GPD using inverse transform sampling.

    The sampling uses the quantile function:

    For ξ ≠ 0:
        Q(p) = (σ/ξ) * ((1-p)^(-ξ) - 1)

    For ξ = 0:
        Q(p) = -σ * ln(1-p)

    Args:
        key: JAX random key for sampling
        sample_shape: Shape of samples to generate
        shape: Keyword-only alias for ``sample_shape`` matching the
            ``pipekit_cycle.ObservationNoise`` protocol's parameter
            name (``sample(key, shape)``), so callers typed against
            that protocol can pass it by keyword. Passing both a
            non-empty ``sample_shape`` and ``shape`` raises
            ``ValueError``.

    Returns:
        Array of samples from the GPD (all within support)
    """
    if shape is not None:
        if sample_shape != ():
            raise ValueError(
                "Pass only one of 'sample_shape' (NumPyro spelling) or "
                "'shape' (ObservationNoise protocol spelling), not both."
            )
        sample_shape = shape
    check_prng_key(key)
    extended_shape = sample_shape + self.batch_shape

    # JAX's Uniform(0, 1) sampler can emit exact 0 or 1 at the
    # endpoints; passing those to icdf yields -inf/+inf and poisons
    # downstream computations. Clamp away from the endpoints.
    uniform_samples = dist.Uniform(0.0, 1.0).sample(key, extended_shape)
    eps = jnp.finfo(uniform_samples.dtype).eps
    uniform_samples = jnp.clip(uniform_samples, eps, 1.0 - eps)

    # Apply inverse CDF transformation
    return self.icdf(uniform_samples)

log_prob(value: jnp.ndarray) -> jnp.ndarray

Log PDF. Thin wrapper for :func:~xtremax.primitives.gpd.gpd_log_prob.

Source code in src/xtremax/distributions/gpd.py

@validate_sample
def log_prob(self, value: jnp.ndarray) -> jnp.ndarray:
    """Log PDF. Thin wrapper for :func:`~xtremax.primitives.gpd.gpd_log_prob`."""
    return gpd_log_prob(value, self.scale, self.concentration)

cdf(value: jnp.ndarray) -> jnp.ndarray

CDF. Thin wrapper for :func:~xtremax.primitives.gpd.gpd_cdf.

Source code in src/xtremax/distributions/gpd.py

def cdf(self, value: jnp.ndarray) -> jnp.ndarray:
    """CDF. Thin wrapper for :func:`~xtremax.primitives.gpd.gpd_cdf`."""
    return gpd_cdf(value, self.scale, self.concentration)

icdf(q: jnp.ndarray) -> jnp.ndarray

Quantile function. Thin wrapper for gpd_icdf.

Source code in src/xtremax/distributions/gpd.py

def icdf(self, q: jnp.ndarray) -> jnp.ndarray:
    """Quantile function. Thin wrapper for ``gpd_icdf``."""
    return gpd_icdf(q, self.scale, self.concentration)

upper_bound() -> jnp.ndarray

Compute the upper bound of the support.

Returns:

Type	Description
ndarray	Upper bound:
ndarray	ξ ≥ 0: +∞
ndarray	ξ < 0: -σ/ξ

Source code in src/xtremax/distributions/gpd.py

def upper_bound(self) -> jnp.ndarray:
    """
    Compute the upper bound of the support.

    Returns:
        Upper bound:
        - ξ ≥ 0: +∞
        - ξ < 0: -σ/ξ
    """
    scale = jnp.asarray(self.scale)
    shape = jnp.asarray(self.concentration)
    # Replace ξ=0 with a safe placeholder inside the division so Python
    # scalars don't raise ZeroDivisionError before jnp.where selects
    # the +∞ branch.
    safe_shape = jnp.where(shape == 0.0, 1.0, shape)
    return jnp.where(shape < 0, -scale / safe_shape, jnp.inf)

lower_bound() -> jnp.ndarray

Compute the lower bound of the support.

Returns:

Type	Description
ndarray	Lower bound: 0 (always non-negative for GPD)

Source code in src/xtremax/distributions/gpd.py

def lower_bound(self) -> jnp.ndarray:
    """
    Compute the lower bound of the support.

    Returns:
        Lower bound: 0 (always non-negative for GPD)
    """
    return jnp.zeros_like(self.scale)

covariance() -> jnp.ndarray

Marginal variance broadcast to the batch shape.

Structural seam for pipekit_cycle.ObservationNoise (which requires covariance() and sample(key, shape)): together with the existing :meth:sample, this lets the distribution act as a non-Gaussian observation-error model in data-assimilation experiments — without xtremax importing pipekit (see docs/interop.md). The observation errors are treated as independent per batch element, so the "covariance" is the marginal variance vector rather than a dense matrix; consuming analysis steps interpret it as a diagonal.

Returns:

Type	Description
ndarray	Variance broadcast to batch_shape (+inf where the
ndarray	variance does not exist, i.e. ξ >= 1/2).

Source code in src/xtremax/distributions/gpd.py

def covariance(self) -> jnp.ndarray:
    """Marginal variance broadcast to the batch shape.

    Structural seam for ``pipekit_cycle.ObservationNoise`` (which
    requires ``covariance()`` and ``sample(key, shape)``): together
    with the existing :meth:`sample`, this lets the distribution act
    as a non-Gaussian observation-error model in data-assimilation
    experiments — without xtremax importing pipekit (see
    ``docs/interop.md``). The observation errors are treated as
    independent per batch element, so the "covariance" is the
    marginal variance vector rather than a dense matrix; consuming
    analysis steps interpret it as a diagonal.

    Returns:
        Variance broadcast to ``batch_shape`` (``+inf`` where the
        variance does not exist, i.e. ξ >= 1/2).
    """
    return jnp.broadcast_to(jnp.asarray(self.variance), self.batch_shape)

kurtosis() -> jnp.ndarray

Compute the excess kurtosis of the GPD.

Excess kurtosis exists when ξ < 1/4:

κ = 3(1-2ξ)(2ξ²+ξ+3) / ((1-3ξ)(1-4ξ)) - 3

Returns:

Type	Description
ndarray	Excess kurtosis or +∞ when it doesn't exist

Source code in src/xtremax/distributions/gpd.py

def kurtosis(self) -> jnp.ndarray:
    """
    Compute the excess kurtosis of the GPD.

    Excess kurtosis exists when ξ < 1/4:

    κ = 3(1-2ξ)(2ξ²+ξ+3) / ((1-3ξ)(1-4ξ)) - 3

    Returns:
        Excess kurtosis or +∞ when it doesn't exist
    """
    shape = self.concentration

    # Kurtosis exists for ξ < 1/4
    kurt_exists = shape < 0.25

    # The rational formula is exact and smooth at ξ = 0 (evaluating to
    # the exponential value 6), so no separate branch or threshold is
    # needed — the old class-local 1e-8 threshold is gone.
    numerator = 3.0 * (1.0 - 2.0 * shape) * (2.0 * shape**2 + shape + 3.0)
    denominator = (1.0 - 3.0 * shape) * (1.0 - 4.0 * shape)
    kurtosis_val = numerator / denominator - 3.0

    # NaN, not inf: beyond ξ = 1/4 the standardized fourth moment is
    # undefined, unlike the variance which genuinely diverges to +inf.
    return jnp.where(kurt_exists, kurtosis_val, jnp.nan)

skew() -> jnp.ndarray

Compute the skewness of the GPD.

Skewness exists when ξ < 1/3:

γ₃ = 2(1+ξ)√(1-2ξ) / (1-3ξ)

Returns:

Type	Description
ndarray	Skewness or +∞ when it doesn't exist

Source code in src/xtremax/distributions/gpd.py

def skew(self) -> jnp.ndarray:
    """
    Compute the skewness of the GPD.

    Skewness exists when ξ < 1/3:

    γ₃ = 2(1+ξ)√(1-2ξ) / (1-3ξ)

    Returns:
        Skewness or +∞ when it doesn't exist
    """
    shape = self.concentration

    # Skewness exists for ξ < 1/3
    skew_exists = shape < 1.0 / 3.0

    # The closed-form expression is exact and smooth at ξ = 0
    # (evaluating to the exponential value 2), so no separate branch or
    # threshold is needed — the old class-local 1e-8 threshold is gone.
    numerator = 2.0 * (1.0 + shape) * jnp.sqrt(1.0 - 2.0 * shape)
    denominator = 1.0 - 3.0 * shape
    skewness_val = numerator / denominator

    # NaN, not inf: beyond ξ = 1/3 the standardized third moment is
    # undefined, unlike the variance which genuinely diverges to +inf.
    return jnp.where(skew_exists, skewness_val, jnp.nan)

entropy() -> jnp.ndarray

Compute the differential entropy of the GPD.

For the GPD: H = log(σ) + ξ + 1

Returns:

Type	Description
ndarray	Differential entropy in nats

Source code in src/xtremax/distributions/gpd.py

def entropy(self) -> jnp.ndarray:
    """
    Compute the differential entropy of the GPD.

    For the GPD:
    H = log(σ) + ξ + 1

    Returns:
        Differential entropy in nats
    """
    return jnp.log(self.scale) + self.concentration + 1.0

survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the survival function S(x) = 1 - F(x).

For ξ ≠ 0: S(x) = (1 + ξx/σ)^(-1/ξ)

For ξ = 0: S(x) = exp(-x/σ)

This has direct interpretation in POT models as the probability of observing an exceedance larger than x.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate survival function	required

Returns:

Type	Description
ndarray	Survival probabilities

Source code in src/xtremax/distributions/gpd.py

def survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the survival function S(x) = 1 - F(x).

    For ξ ≠ 0:
        S(x) = (1 + ξx/σ)^(-1/ξ)

    For ξ = 0:
        S(x) = exp(-x/σ)

    This has direct interpretation in POT models as the probability
    of observing an exceedance larger than x.

    Args:
        value: Points at which to evaluate survival function

    Returns:
        Survival probabilities
    """
    # Delegate to the stable primitive so the class and functional APIs
    # share the same smooth ξ→0 handling (no threshold to drift against).
    return gpd_survival(value, self.scale, self.concentration)

hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the hazard rate h(x) = f(x) / S(x).

For GPD: h(x) = 1 / (σ + ξx)

This represents the instantaneous failure rate in reliability applications or the rate of threshold exceedance.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate hazard rate	required

Returns:

Type	Description
ndarray	Hazard rate values

Source code in src/xtremax/distributions/gpd.py

def hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the hazard rate h(x) = f(x) / S(x).

    For GPD:
    h(x) = 1 / (σ + ξx)

    This represents the instantaneous failure rate in reliability
    applications or the rate of threshold exceedance.

    Args:
        value: Points at which to evaluate hazard rate

    Returns:
        Hazard rate values
    """
    scale, shape = self.scale, self.concentration

    # Hazard rate: h(x) = 1 / (σ + ξx). Enforce BOTH the lower support
    # bound (GPD is defined for x ≥ 0; f(x)=0, S(x)=1 ⇒ h(x)=0 below)
    # and the upper bound in the bounded ξ < 0 case (σ + ξx > 0).
    denominator = scale + shape * value
    valid = (denominator > 0.0) & (value >= 0.0)
    hazard_val = 1.0 / jnp.where(valid, denominator, 1.0)

    return jnp.where(valid, hazard_val, 0.0)

cumulative_hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the cumulative hazard rate Λ(x) = -log(S(x)).

For ξ ≠ 0: Λ(x) = (1/ξ) * log(1 + ξx/σ)

For ξ = 0: Λ(x) = x/σ

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate cumulative hazard rate	required

Returns:

Type	Description
ndarray	Cumulative hazard rate values

Source code in src/xtremax/distributions/gpd.py

def cumulative_hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the cumulative hazard rate Λ(x) = -log(S(x)).

    For ξ ≠ 0:
        Λ(x) = (1/ξ) * log(1 + ξx/σ)

    For ξ = 0:
        Λ(x) = x/σ

    Args:
        value: Points at which to evaluate cumulative hazard rate

    Returns:
        Cumulative hazard rate values
    """
    # Λ(x) = -log S(x); delegate to the stable log-survival primitive.
    return -gpd_log_survival(value, self.scale, self.concentration)

return_level(return_period: float | jnp.ndarray) -> jnp.ndarray

Return level. Thin wrapper for gpd_return_level.

Source code in src/xtremax/distributions/gpd.py

def return_level(self, return_period: float | jnp.ndarray) -> jnp.ndarray:
    """Return level. Thin wrapper for ``gpd_return_level``."""
    return gpd_return_level(return_period, self.scale, self.concentration)

tail_index() -> jnp.ndarray

Compute the tail index for GPD.

For GPD with ξ > 0: tail index α = 1/ξ For ξ ≤ 0: tail index is not defined (no power-law tail)

Returns:

Type	Description
ndarray	Tail index (1/ξ for ξ > 0, ∞ otherwise)

Source code in src/xtremax/distributions/gpd.py

def tail_index(self) -> jnp.ndarray:
    """
    Compute the tail index for GPD.

    For GPD with ξ > 0: tail index α = 1/ξ
    For ξ ≤ 0: tail index is not defined (no power-law tail)

    Returns:
        Tail index (1/ξ for ξ > 0, ∞ otherwise)
    """
    return jnp.where(self.concentration > 0, 1.0 / self.concentration, jnp.inf)

exceedance_probability(threshold: jnp.ndarray) -> jnp.ndarray

Compute probability of exceeding a threshold: P(X > threshold).

This is the survival function and fundamental for POT analysis.

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value	required

Returns:

Type	Description
ndarray	Exceedance probabilities

Source code in src/xtremax/distributions/gpd.py

def exceedance_probability(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute probability of exceeding a threshold: P(X > threshold).

    This is the survival function and fundamental for POT analysis.

    Args:
        threshold: Threshold value

    Returns:
        Exceedance probabilities
    """
    return self.survival_function(threshold)

conditional_excess_mean(threshold: jnp.ndarray) -> jnp.ndarray

Compute the mean excess function: E[X - u | X > u].

For GPD, this has the simple form: E[X - u | X > u] = (σ + ξu) / (1 - ξ)

This linear relationship is a key property of GPD and forms the basis for threshold selection in POT modeling.

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value u	required

Returns:

Type	Description
ndarray	Conditional excess mean values

Source code in src/xtremax/distributions/gpd.py

def conditional_excess_mean(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the mean excess function: E[X - u | X > u].

    For GPD, this has the simple form:
    E[X - u | X > u] = (σ + ξu) / (1 - ξ)

    This linear relationship is a key property of GPD and forms
    the basis for threshold selection in POT modeling.

    Args:
        threshold: Threshold value u

    Returns:
        Conditional excess mean values
    """
    scale, shape = self.scale, self.concentration

    # Only valid for ξ < 1 and threshold within support
    mean_exists = shape < 1.0
    within_support = (threshold >= 0.0) & (threshold < self.upper_bound())
    valid = mean_exists & within_support

    # Mean excess: (σ + ξu) / (1 - ξ)
    excess_mean = (scale + shape * threshold) / (1.0 - shape)

    return jnp.where(valid, excess_mean, jnp.inf)

threshold_stability_plot_data(thresholds: jnp.ndarray) -> dict

Generate data for threshold stability plots in POT analysis.

Returns modified scale parameters σ* = σ + ξ(u - u₀) for different thresholds, which should be approximately constant if GPD fits well.

Parameters:

Name	Type	Description	Default
thresholds	ndarray	Array of threshold values	required

Returns:

Type	Description
dict	Dictionary with threshold stability metrics

Source code in src/xtremax/distributions/gpd.py

def threshold_stability_plot_data(self, thresholds: jnp.ndarray) -> dict:
    """
    Generate data for threshold stability plots in POT analysis.

    Returns modified scale parameters σ* = σ + ξ(u - u₀) for different
    thresholds, which should be approximately constant if GPD fits well.

    Args:
        thresholds: Array of threshold values

    Returns:
        Dictionary with threshold stability metrics
    """
    scale, shape = self.scale, self.concentration

    # Reference threshold (typically the lowest)
    u0 = thresholds[0]

    # Modified scale parameters for each threshold
    modified_scales = scale + shape * (thresholds - u0)

    # Shape parameters (should remain constant)
    shapes = jnp.full_like(thresholds, shape)

    return {
        "thresholds": thresholds,
        "modified_scales": modified_scales,
        "shapes": shapes,
        "reference_threshold": u0,
    }

hill_plot_data(order_statistics: jnp.ndarray, k_values: jnp.ndarray) -> dict

Generate data for Hill plots (tail index estimation).

Computes Hill estimator: α̂k = (1/k) * Σᵢ₌₁ᵏ log(X)} / X_{n-k,n

Vectorized over k_values via a cumulative sum of log order statistics, so it is safe under jax.jit (the previous Python loop branched on traced values and raised a TracerBoolConversion error). Invalid k (≤ 0 or ≥ n) yield NaN.

Parameters:

Name	Type	Description	Default
order_statistics	ndarray	Sorted sample in descending order	required
k_values	ndarray	Numbers of upper order statistics to use	required

Returns:

Type	Description
dict	Dictionary with Hill plot data

Source code in src/xtremax/distributions/gpd.py

def hill_plot_data(
    self, order_statistics: jnp.ndarray, k_values: jnp.ndarray
) -> dict:
    """
    Generate data for Hill plots (tail index estimation).

    Computes Hill estimator: α̂_k = (1/k) * Σᵢ₌₁ᵏ log(X_{n-i+1,n} / X_{n-k,n})

    Vectorized over ``k_values`` via a cumulative sum of log order
    statistics, so it is safe under ``jax.jit`` (the previous Python
    loop branched on traced values and raised a ``TracerBoolConversion``
    error). Invalid ``k`` (≤ 0 or ≥ n) yield NaN.

    Args:
        order_statistics: Sorted sample in descending order
        k_values: Numbers of upper order statistics to use

    Returns:
        Dictionary with Hill plot data
    """
    order_statistics = jnp.asarray(order_statistics)
    k_arr = jnp.asarray(k_values)
    n = order_statistics.shape[0]

    # (1/k)·Σ_{i<k} log X_i − log X_k, for every k at once.
    logs = jnp.log(order_statistics)
    cumsum = jnp.cumsum(logs)
    valid = (k_arr > 0) & (k_arr < n)
    k_safe = jnp.clip(k_arr, 1, max(n - 1, 1))
    mean_log_top = cumsum[k_safe - 1] / k_safe
    hill_est = mean_log_top - logs[k_safe]
    estimates = jnp.where(hill_est > 0, 1.0 / hill_est, jnp.inf)
    estimates = jnp.where(valid, estimates, jnp.nan)

    return {
        "k_values": k_values,
        "hill_estimates": estimates,
        "theoretical_tail_index": self.tail_index(),
    }

expand(batch_shape: tuple[int, ...]) -> dist.Distribution

Expand to batch_shape by reconstructing via __init__.

We deliberately go through the constructor (rather than an ExpandedDistribution wrapper) so every EVT method stays available on the returned distribution.

Source code in src/xtremax/distributions/gpd.py

def expand(self, batch_shape: tuple[int, ...]) -> dist.Distribution:
    """Expand to ``batch_shape`` by reconstructing via ``__init__``.

    We deliberately go through the constructor (rather than an
    ``ExpandedDistribution`` wrapper) so every EVT method stays
    available on the returned distribution.
    """
    batch_shape = tuple(batch_shape)
    if batch_shape == self.batch_shape:
        return self
    return type(self)(
        scale=jnp.broadcast_to(self.scale, batch_shape),
        concentration=jnp.broadcast_to(self.concentration, batch_shape),
        validate_args=self._validate_args,
    )

Fixed-type GEV families

The three GEV types as standalone distributions, for when theory or prior evidence pins the tail behaviour and you would rather not estimate \(\xi\). Fixing the type removes a parameter and usually tightens return-level intervals — at the cost of an assumption the data can no longer overturn.

Class	Type	Tail
GumbelType1GEVD	I (\(\xi = 0\))	Exponential
FrechetType2GEVD	II (\(\xi > 0\))	Heavy polynomial, lower-bounded
WeibullType3GEVD	III (\(\xi < 0\))	Bounded above

NumPyro-compatible extreme value distributions.

GumbelType1GEVD(loc: float = 0.0, scale: float = 1.0, validate_args: bool | None = None)

Bases: Gumbel

Gumbel Type I Generalized Extreme Value Distribution for NumPyro.

The distribution core — sample/log_prob/cdf/icdf/ mean/variance/support — is inherited directly from :class:numpyro.distributions.Gumbel (parameter-for-parameter identical), so sampling uses numpyro's native tail-accurate jax.random.gumbel path and the class inherits numpyro's batching and KL registrations for free. This subclass adds the EVT sugar: return levels, moments beyond variance, hazard/survival methods, the characteristic function, and probability-paper utilities.

The Gumbel Type I is the GEVD with shape parameter ξ = 0, characterized by exponential tails and unbounded support. This distribution is the limiting case of the GEVD and represents the "double exponential" distribution for extreme values. It's the most commonly used extreme value distribution in meteorology, hydrology, and engineering applications.

Key Characteristics: - Unbounded support: x ∈ (-∞, +∞) - Exponential tails: P(X > x) ~ exp(-exp((x-μ)/σ)) for large x - Moderate extreme behavior (neither heavy nor light tails) - Asymmetric with positive skewness ≈ 1.14 - Often emerges from block maxima of exponential-type distributions

Probability Density Function:

f(x) = (1/σ) * exp(-z - exp(-z)) where z = (x - μ)/σ

Cumulative Distribution Function:

F(x) = exp(-exp(-z)) where z = (x - μ)/σ

Quantile Function:

Q(p) = μ - σ * ln(-ln(p))

Support:

x ∈ (-∞, +∞) (unbounded)

Tail Behavior:

Upper tail: P(X > x) ~ exp(-exp((x-μ)/σ)) (exponential decay) Lower tail: P(X < x) ~ exp(exp((μ-x)/σ)) (double exponential growth)

Parameters:

• loc (μ): Location parameter ∈ ℝ (mode of the distribution)
• scale (σ): Scale parameter > 0 (controls spread)

Statistical Properties:

• Mean: E[X] = μ + σγ (where γ ≈ 0.5772 is Euler-Mascheroni constant)
• Mode: μ (location parameter)
• Variance: Var[X] = σ²π²/6
• Skewness: γ₃ ≈ 1.1396 (constant, independent of parameters)
• Kurtosis: γ₄ = 12/5 = 2.4 (excess kurtosis, constant)

Applications in Geosciences:

• Daily maximum/minimum temperatures (temperate climates)
• Wind speed extremes (non-hurricane conditions)
• River discharge peaks (regulated systems)
• Wave height extremes (moderate sea states)
• Precipitation extremes (stable atmospheric conditions)
• Seismic background noise levels
• Atmospheric pressure extremes

Examples:

>>> import jax
>>> import jax.numpy as jnp
>>>
>>> # Create Gumbel distribution for temperature extremes
>>> temp_extremes = GumbelType1GEVD(loc=35.0, scale=5.0)  # °C
>>>
>>> # Key properties
>>> print(f"Mean temperature: {temp_extremes.mean:.2f}°C")
>>> print(f"Mode (most likely): {temp_extremes.mode:.2f}°C")
>>> print(f"Standard deviation: {jnp.sqrt(temp_extremes.variance):.2f}°C")
>>> print(f"Skewness: {temp_extremes.skew():.3f}")
>>>
>>> # Sample and evaluate
>>> key = jax.random.PRNGKey(42)
>>> samples = temp_extremes.sample(key, sample_shape=(1000,))
>>> log_probs = temp_extremes.log_prob(samples)
>>>
>>> # Return level analysis
>>> return_50yr = temp_extremes.return_level(50)
>>> print(f"50-year return level: {return_50yr:.2f}°C")

Initialize the Gumbel Type I GEVD.

The core distribution machinery is delegated to :class:numpyro.distributions.Gumbel; only the mathematical constants used by the EVT sugar methods are set up here.

Parameters:

Name	Type	Description	Default
loc	float	Location parameter μ (mode of distribution)	0.0
scale	float	Scale parameter σ (positive real number)	1.0
validate_args	bool | None	Whether to validate input arguments	None

Raises:

Type	Description
ValueError	If scale <= 0

Source code in src/xtremax/distributions/gumbel.py

def __init__(
    self, loc: float = 0.0, scale: float = 1.0, validate_args: bool | None = None
):
    """
    Initialize the Gumbel Type I GEVD.

    The core distribution machinery is delegated to
    :class:`numpyro.distributions.Gumbel`; only the mathematical
    constants used by the EVT sugar methods are set up here.

    Args:
        loc: Location parameter μ (mode of distribution)
        scale: Scale parameter σ (positive real number)
        validate_args: Whether to validate input arguments

    Raises:
        ValueError: If scale <= 0
    """
    # Mathematical constants for Gumbel distribution
    self._euler_gamma = 0.5772156649015329  # Euler-Mascheroni constant
    self._pi_squared_over_six = (jnp.pi**2) / 6.0  # π²/6 for variance
    self._gumbel_skewness = 1.1395470994046486  # 12*ζ(3)/π² where ζ(3) ≈ 1.202
    self._gumbel_kurtosis = 12.0 / 5.0  # Excess kurtosis = 12/5

    super().__init__(loc, scale, validate_args=validate_args)

mode: jnp.ndarray property

Compute the mode of the Gumbel Type I GEVD.

The mode is simply the location parameter: mode = μ

This is the most likely extreme value and represents the peak of the probability density function.

Returns:

Type	Description
ndarray	Mode of the distribution (equals location parameter)

sample(key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray

Sample via numpyro's native jax.random.gumbel path.

Only adds xtremax's up-front PRNG-key validation (a clear TypeError instead of numpyro's bare AssertionError) before delegating to :meth:numpyro.distributions.Gumbel.sample.

Source code in src/xtremax/distributions/gumbel.py

def sample(self, key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray:
    """Sample via numpyro's native ``jax.random.gumbel`` path.

    Only adds xtremax's up-front PRNG-key validation (a clear
    ``TypeError`` instead of numpyro's bare ``AssertionError``) before
    delegating to :meth:`numpyro.distributions.Gumbel.sample`.
    """
    check_prng_key(key)
    return super().sample(key, sample_shape)

upper_bound() -> jnp.ndarray

Compute the upper bound of the support.

For Gumbel Type I: upper bound = +∞

Returns:

Type	Description
ndarray	Upper bound (+∞ for Gumbel)

Source code in src/xtremax/distributions/gumbel.py

def upper_bound(self) -> jnp.ndarray:
    """
    Compute the upper bound of the support.

    For Gumbel Type I: upper bound = +∞

    Returns:
        Upper bound (+∞ for Gumbel)
    """
    return jnp.full_like(self.loc, jnp.inf)

lower_bound() -> jnp.ndarray

Compute the lower bound of the support.

For Gumbel Type I: lower bound = -∞

Returns:

Type	Description
ndarray	Lower bound (-∞ for Gumbel)

Source code in src/xtremax/distributions/gumbel.py

def lower_bound(self) -> jnp.ndarray:
    """
    Compute the lower bound of the support.

    For Gumbel Type I: lower bound = -∞

    Returns:
        Lower bound (-∞ for Gumbel)
    """
    return jnp.full_like(self.loc, -jnp.inf)

kurtosis() -> jnp.ndarray

Compute the excess kurtosis of the Gumbel Type I GEVD.

The excess kurtosis is constant: κ = 12/5 = 2.4

This value is independent of the distribution parameters and indicates heavier tails than the normal distribution (which has κ = 0).

Returns:

Type	Description
ndarray	Excess kurtosis (constant value 12/5)

Source code in src/xtremax/distributions/gumbel.py

def kurtosis(self) -> jnp.ndarray:
    """
    Compute the excess kurtosis of the Gumbel Type I GEVD.

    The excess kurtosis is constant:
    κ = 12/5 = 2.4

    This value is independent of the distribution parameters and indicates
    heavier tails than the normal distribution (which has κ = 0).

    Returns:
        Excess kurtosis (constant value 12/5)
    """
    return jnp.full_like(self.loc, self._gumbel_kurtosis)

skew() -> jnp.ndarray

Compute the skewness of the Gumbel Type I GEVD.

The skewness is constant: γ₃ ≈ 1.1396470994

This positive value indicates right-skewness (longer right tail), which is characteristic of extreme value distributions for maxima. The value is independent of distribution parameters.

Returns:

Type	Description
ndarray	Skewness (constant positive value ≈ 1.14)

Source code in src/xtremax/distributions/gumbel.py

def skew(self) -> jnp.ndarray:
    """
    Compute the skewness of the Gumbel Type I GEVD.

    The skewness is constant:
    γ₃ ≈ 1.1396470994

    This positive value indicates right-skewness (longer right tail),
    which is characteristic of extreme value distributions for maxima.
    The value is independent of distribution parameters.

    Returns:
        Skewness (constant positive value ≈ 1.14)
    """
    return jnp.full_like(self.loc, self._gumbel_skewness)

entropy() -> jnp.ndarray

Compute the differential entropy of the Gumbel Type I GEVD.

The entropy is: H = log(σ) + γ + 1 where γ is the Euler-Mascheroni constant.

This measures the uncertainty in the distribution and increases with the scale parameter.

Returns:

Type	Description
ndarray	Differential entropy in nats

Source code in src/xtremax/distributions/gumbel.py

def entropy(self) -> jnp.ndarray:
    """
    Compute the differential entropy of the Gumbel Type I GEVD.

    The entropy is:
    H = log(σ) + γ + 1
    where γ is the Euler-Mascheroni constant.

    This measures the uncertainty in the distribution and increases
    with the scale parameter.

    Returns:
        Differential entropy in nats
    """
    return jnp.log(self.scale) + self._euler_gamma + 1.0

survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the survival function S(x) = 1 - F(x).

For Gumbel Type I: S(x) = 1 - exp(-exp(-z)) where z = (x - μ)/σ

This represents the probability of exceeding a given value, which is fundamental for risk analysis and return level calculations.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the survival function	required

Returns:

Type	Description
ndarray	Survival probabilities

Source code in src/xtremax/distributions/gumbel.py

def survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the survival function S(x) = 1 - F(x).

    For Gumbel Type I:
    S(x) = 1 - exp(-exp(-z)) where z = (x - μ)/σ

    This represents the probability of exceeding a given value,
    which is fundamental for risk analysis and return level calculations.

    Args:
        value: Points at which to evaluate the survival function

    Returns:
        Survival probabilities
    """
    # Computed as `-expm1(-exp(-z))` rather than `1 - cdf` so the upper
    # tail stays numerically resolvable (1 - exp(-small) catastrophically
    # cancels to 0 once S(x) is below ~1e-7).
    z = (value - self.loc) / self.scale
    return -jnp.expm1(-jnp.exp(-z))

log_survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the log survival function log S(x) = log(1 - F(x)).

For extreme quantiles, this provides better numerical stability than computing log(1 - F(x)) directly.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate log survival function	required

Returns:

Type	Description
ndarray	Log survival probabilities

Source code in src/xtremax/distributions/gumbel.py

def log_survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the log survival function log S(x) = log(1 - F(x)).

    For extreme quantiles, this provides better numerical stability
    than computing log(1 - F(x)) directly.

    Args:
        value: Points at which to evaluate log survival function

    Returns:
        Log survival probabilities
    """
    loc, scale = self.loc, self.scale
    z = (value - loc) / scale
    # S(x) = 1 - F(x) = -expm1(-exp(-z)).
    # log S(x) = log(-expm1(-exp(-z))) — numerically stable across both
    # tails. For large z: exp(-z) → 0, -expm1(-u) → u, so log S → -z
    # (the correct Gumbel upper-tail asymptote). The previous branch
    # returned ≈ -exp(-z) at large z, which tends to 0 instead of the
    # large negative log-survival.
    return jnp.log(-jnp.expm1(-jnp.exp(-z)))

hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Hazard rate :math:h(x) = f(x) / S(x) for the Gumbel.

With :math:F(x) = \exp(-e^{-z}) and :math:z = (x - \mu)/\sigma, :math:S(x) = 1 - F(x) = -\mathrm{expm1}(-e^{-z}), so the hazard does not simplify to :math:e^{z}/\sigma — the upper-tail limit is :math:1/\sigma, not diverging.

Source code in src/xtremax/distributions/gumbel.py

def hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    r"""Hazard rate :math:`h(x) = f(x) / S(x)` for the Gumbel.

    With :math:`F(x) = \exp(-e^{-z})` and :math:`z = (x - \mu)/\sigma`,
    :math:`S(x) = 1 - F(x) = -\mathrm{expm1}(-e^{-z})`, so the hazard
    does *not* simplify to :math:`e^{z}/\sigma` — the upper-tail limit
    is :math:`1/\sigma`, not diverging.
    """
    z = (value - self.loc) / self.scale
    log_pdf = gumbel_log_prob(value, self.loc, self.scale)
    # log S(x) = log(-expm1(-exp(-z))) — numerically stable across both tails.
    log_surv = jnp.log(-jnp.expm1(-jnp.exp(-z)))
    return jnp.exp(log_pdf - log_surv)

cumulative_hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Cumulative hazard :math:\Lambda(x) = -\log S(x) for the Gumbel.

With :math:S(x) = -\mathrm{expm1}(-e^{-z}), :math:\Lambda(x) = -\log(-\mathrm{expm1}(-e^{-z})).

Source code in src/xtremax/distributions/gumbel.py

def cumulative_hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    r"""Cumulative hazard :math:`\Lambda(x) = -\log S(x)` for the Gumbel.

    With :math:`S(x) = -\mathrm{expm1}(-e^{-z})`,
    :math:`\Lambda(x) = -\log(-\mathrm{expm1}(-e^{-z}))`.
    """
    z = (value - self.loc) / self.scale
    return -jnp.log(-jnp.expm1(-jnp.exp(-z)))

return_level(return_period: float | jnp.ndarray) -> jnp.ndarray

Return level. Thin wrapper for gumbel_return_level.

Source code in src/xtremax/distributions/gumbel.py

def return_level(self, return_period: float | jnp.ndarray) -> jnp.ndarray:
    """Return level. Thin wrapper for ``gumbel_return_level``."""
    return gumbel_return_level(return_period, self.loc, self.scale)

tail_index() -> jnp.ndarray

Compute the tail index for Gumbel Type I GEVD.

For Gumbel distributions, there is no power-law tail behavior. The tail index concept (from extreme value theory) is not applicable since Gumbel has exponential (not polynomial) tails.

Returns:

Type	Description
ndarray	NaN (tail index not defined for exponential tails)

Source code in src/xtremax/distributions/gumbel.py

def tail_index(self) -> jnp.ndarray:
    """
    Compute the tail index for Gumbel Type I GEVD.

    For Gumbel distributions, there is no power-law tail behavior.
    The tail index concept (from extreme value theory) is not applicable
    since Gumbel has exponential (not polynomial) tails.

    Returns:
        NaN (tail index not defined for exponential tails)
    """
    return jnp.full_like(self.loc, jnp.nan)

exceedance_probability(threshold: jnp.ndarray) -> jnp.ndarray

Compute probability of exceeding a threshold: P(X > threshold).

This is equivalent to the survival function and follows the double exponential decay characteristic of Gumbel distributions.

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value	required

Returns:

Type	Description
ndarray	Exceedance probabilities

Source code in src/xtremax/distributions/gumbel.py

def exceedance_probability(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute probability of exceeding a threshold: P(X > threshold).

    This is equivalent to the survival function and follows the
    double exponential decay characteristic of Gumbel distributions.

    Args:
        threshold: Threshold value

    Returns:
        Exceedance probabilities
    """
    return self.survival_function(threshold)

conditional_excess_mean(threshold: jnp.ndarray) -> jnp.ndarray

Mean excess :math:E[X - u \mid X > u] for the Gumbel.

Computed from the identity

.. math:: E[X - u \mid X > u] \, S(u) = \int_u^\infty S(x)\,dx

via a trapezoidal quadrature over a grid from u to u + 50σ (well past the Gumbel tail). Works for scalar, batched, and broadcasted threshold / scale — the grid axis is always placed last so trapezoidal integration stays on the quadrature axis rather than the batch axis.

Returns NaN where the survival probability is effectively zero.

Source code in src/xtremax/distributions/gumbel.py

def conditional_excess_mean(self, threshold: jnp.ndarray) -> jnp.ndarray:
    r"""Mean excess :math:`E[X - u \mid X > u]` for the Gumbel.

    Computed from the identity

    .. math:: E[X - u \mid X > u] \, S(u) = \int_u^\infty S(x)\,dx

    via a trapezoidal quadrature over a grid from ``u`` to ``u + 50σ``
    (well past the Gumbel tail). Works for scalar, batched, and
    broadcasted ``threshold`` / ``scale`` — the grid axis is always
    placed last so trapezoidal integration stays on the quadrature
    axis rather than the batch axis.

    Returns ``NaN`` where the survival probability is effectively zero.
    """
    threshold_arr = jnp.asarray(threshold)
    scale_arr = jnp.asarray(self.scale)

    # Trapezoidal quadrature of the IBP form
    #   E[X - u | X > u] = (1/S(u)) ∫_u^∞ S(x) dx
    # with an *adaptive* upper cap chosen as the further of
    # ``u + 50σ`` (enough for exponential-tail decay from u) and
    # ``icdf(1 - 1e-6)`` (absolute far-tail quantile). The fixed
    # ``u + 50σ`` cap used previously truncated substantial mass
    # when u sat far below the location (e.g. u = μ - 1000σ),
    # biasing the mean excess low by O(|u|). The adaptive cap
    # handles both far-left and far-right thresholds correctly.
    n_grid = 1024
    q_top = jnp.asarray(1.0 - 1e-6, dtype=scale_arr.dtype)
    x_top = self.icdf(q_top)
    upper = jnp.maximum(threshold_arr + 50.0 * scale_arr, x_top)

    # Joint batch of threshold and parameters; per-batch quantities are
    # broadcast to it before gaining the trailing grid axis — batched
    # loc/scale evaluated against a (..., n_grid) grid otherwise raise
    # a broadcasting error inside survival_function.
    batch = jnp.shape(upper)
    unit = jnp.linspace(0.0, 1.0, n_grid)
    t_exp = jnp.expand_dims(jnp.broadcast_to(threshold_arr, batch), axis=-1)
    u_exp = jnp.expand_dims(upper, axis=-1)
    x_grid = t_exp + (u_exp - t_exp) * unit  # (..., n_grid)

    loc_e = jnp.expand_dims(jnp.broadcast_to(self.loc, batch), axis=-1)
    scale_e = jnp.expand_dims(jnp.broadcast_to(scale_arr, batch), axis=-1)
    z_grid = (x_grid - loc_e) / scale_e
    integrand = -jnp.expm1(-jnp.exp(-z_grid))  # S(x) on the grid
    integral = jnp.trapezoid(integrand, x=x_grid, axis=-1)
    s_u = self.survival_function(threshold_arr)
    return jnp.where(s_u > 1e-15, integral / s_u, jnp.nan)

median() -> jnp.ndarray

Compute the median (50th percentile) of the Gumbel Type I GEVD.

The median is: Q(0.5) = μ - σ * ln(-ln(0.5)) = μ - σ * ln(ln(2))

Returns:

Type	Description
ndarray	Median value of the distribution

Source code in src/xtremax/distributions/gumbel.py

def median(self) -> jnp.ndarray:
    """
    Compute the median (50th percentile) of the Gumbel Type I GEVD.

    The median is:
    Q(0.5) = μ - σ * ln(-ln(0.5)) = μ - σ * ln(ln(2))

    Returns:
        Median value of the distribution
    """
    return self.icdf(0.5)

interquartile_range() -> jnp.ndarray

Compute the interquartile range (IQR) of the Gumbel Type I GEVD.

IQR = Q(0.75) - Q(0.25)

Returns:

Type	Description
ndarray	Interquartile range

Source code in src/xtremax/distributions/gumbel.py

def interquartile_range(self) -> jnp.ndarray:
    """
    Compute the interquartile range (IQR) of the Gumbel Type I GEVD.

    IQR = Q(0.75) - Q(0.25)

    Returns:
        Interquartile range
    """
    q25 = self.icdf(0.25)
    q75 = self.icdf(0.75)
    return q75 - q25

coefficient_of_variation() -> jnp.ndarray

Compute the coefficient of variation (CV = σ/μ) for risk assessment.

For Gumbel distributions: CV = (σπ/√6) / (μ + σγ)

Returns:

Type	Description
ndarray	Coefficient of variation

Source code in src/xtremax/distributions/gumbel.py

def coefficient_of_variation(self) -> jnp.ndarray:
    """
    Compute the coefficient of variation (CV = σ/μ) for risk assessment.

    For Gumbel distributions:
    CV = (σπ/√6) / (μ + σγ)

    Returns:
        Coefficient of variation
    """
    std_dev = jnp.sqrt(self.variance)
    return std_dev / jnp.abs(self.mean)

moment_generating_function(t: jnp.ndarray) -> jnp.ndarray

Compute the moment generating function M(t) = E[exp(tX)].

For Gumbel Type I with t < 1/σ: M(t) = exp(μt) * Γ(1 - σt) where Γ is the gamma function.

Parameters:

Name	Type	Description	Default
t	ndarray	Values at which to evaluate MGF (must satisfy t < 1/σ)	required

Returns:

Type	Description
ndarray	Moment generating function values

Source code in src/xtremax/distributions/gumbel.py

def moment_generating_function(self, t: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the moment generating function M(t) = E[exp(tX)].

    For Gumbel Type I with t < 1/σ:
    M(t) = exp(μt) * Γ(1 - σt)
    where Γ is the gamma function.

    Args:
        t: Values at which to evaluate MGF (must satisfy t < 1/σ)

    Returns:
        Moment generating function values
    """
    loc, scale = self.loc, self.scale

    # MGF exists for t < 1/σ
    valid = t < (1.0 / scale)

    # M(t) = exp(μt) * Γ(1 - σt)
    mgf_val = jnp.exp(loc * t) * jnp.exp(gammaln(1.0 - scale * t))

    return jnp.where(valid, mgf_val, jnp.inf)

characteristic_function(t: jnp.ndarray) -> jnp.ndarray

Compute the characteristic function φ(t) = E[exp(itX)].

For Gumbel Type I: φ(t) = exp(iμt) * Γ(1 - iσt)

jax.scipy.special.gammaln is real-only, so the complex argument 1 - iσt is evaluated via jax.pure_callback to scipy.special.loggamma. This is host-hop but numerically correct — a naive JAX call would have failed or returned garbage.

Parameters:

Name	Type	Description	Default
t	ndarray	Real values at which to evaluate the characteristic function	required

Returns:

Type	Description
ndarray	Complex-valued characteristic function

Source code in src/xtremax/distributions/gumbel.py

def characteristic_function(self, t: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the characteristic function φ(t) = E[exp(itX)].

    For Gumbel Type I:
    φ(t) = exp(iμt) * Γ(1 - iσt)

    ``jax.scipy.special.gammaln`` is real-only, so the complex argument
    ``1 - iσt`` is evaluated via ``jax.pure_callback`` to
    ``scipy.special.loggamma``. This is host-hop but numerically
    correct — a naive JAX call would have failed or returned garbage.

    Args:
        t: Real values at which to evaluate the characteristic function

    Returns:
        Complex-valued characteristic function
    """
    loc, scale = self.loc, self.scale
    t_arr = jnp.asarray(t)
    scale_arr = jnp.asarray(scale)
    loc_arr = jnp.asarray(loc)

    complex_arg = 1.0 - 1j * scale_arr * t_arr
    log_gamma = jax.pure_callback(
        _host_complex_loggamma,
        jax.ShapeDtypeStruct(complex_arg.shape, complex_arg.dtype),
        complex_arg,
    )
    return jnp.exp(1j * loc_arr * t_arr) * jnp.exp(log_gamma)

gumbel_probability_paper_coordinates(value: jnp.ndarray) -> tuple[jnp.ndarray, jnp.ndarray]

Convert values to Gumbel probability paper coordinates for graphical analysis.

Transforms (x, F(x)) to coordinates where the Gumbel CDF appears linear: y-coordinate: -ln(-ln(F(x))) x-coordinate: x (unchanged)

Parameters:

Name	Type	Description	Default
value	ndarray	Data values	required

Returns:

Type	Description
tuple[ndarray, ndarray]	Tuple of (x-coordinates, y-coordinates) for probability paper

Source code in src/xtremax/distributions/gumbel.py

def gumbel_probability_paper_coordinates(
    self, value: jnp.ndarray
) -> tuple[jnp.ndarray, jnp.ndarray]:
    """
    Convert values to Gumbel probability paper coordinates for graphical analysis.

    Transforms (x, F(x)) to coordinates where the Gumbel CDF appears linear:
    y-coordinate: -ln(-ln(F(x)))
    x-coordinate: x (unchanged)

    Args:
        value: Data values

    Returns:
        Tuple of (x-coordinates, y-coordinates) for probability paper
    """
    cdf_vals = self.cdf(value)

    # Avoid numerical issues at boundaries
    cdf_vals = jnp.clip(cdf_vals, 1e-15, 1 - 1e-15)

    # Gumbel probability paper transformation: y = -ln(-ln(F))
    y_coords = -jnp.log(-jnp.log(cdf_vals))

    return value, y_coords

expand(batch_shape: tuple[int, ...]) -> dist.Distribution

Expand to batch_shape by reconstructing via __init__.

Going through the constructor re-populates all cached constants (_pi_squared_over_six, _gumbel_skewness, _gumbel_kurtosis, _euler_gamma) on the returned instance. Bypassing __init__ broke variance/skew/kurtosis/ entropy on expanded distributions.

Source code in src/xtremax/distributions/gumbel.py

def expand(self, batch_shape: tuple[int, ...]) -> dist.Distribution:
    """Expand to ``batch_shape`` by reconstructing via ``__init__``.

    Going through the constructor re-populates all cached constants
    (``_pi_squared_over_six``, ``_gumbel_skewness``,
    ``_gumbel_kurtosis``, ``_euler_gamma``) on the returned instance.
    Bypassing ``__init__`` broke ``variance``/``skew``/``kurtosis``/
    ``entropy`` on expanded distributions.
    """
    batch_shape = tuple(batch_shape)
    if batch_shape == self.batch_shape:
        return self
    return type(self)(
        loc=jnp.broadcast_to(self.loc, batch_shape),
        scale=jnp.broadcast_to(self.scale, batch_shape),
        validate_args=self._validate_args,
    )

FrechetType2GEVD(loc: ArrayLike = 0.0, scale: ArrayLike = 1.0, concentration: ArrayLike | None = None, shape: ArrayLike | None = None, validate_args: bool | None = None)

Bases: Distribution

Fréchet Type II Generalized Extreme Value Distribution for NumPyro.

The Fréchet Type II is the GEVD with shape parameter ξ > 0, characterized by heavy polynomial tails and lower-bounded support. This distribution is fundamental in extreme value theory for modeling phenomena with power-law tail behavior, such as financial extreme losses, natural disaster magnitudes, and network traffic bursts.

Key Characteristics: - Lower bounded support: x ≥ μ - σ/ξ (since ξ > 0) - Heavy polynomial tails: P(X > x) ~ x^(-1/ξ) as x → ∞ - Power-law tail index α = 1/ξ - Infinite moments when ξ ≥ 1/k for k-th moment - Models phenomena with "fat tails" and extreme outliers

Probability Density Function:

f(x) = (1/σ) * t(x)^(-(1/ξ + 1)) * exp(-t(x)^(-1/ξ)) where t(x) = 1 + ξ(x - μ)/σ and ξ > 0

Cumulative Distribution Function:

F(x) = exp(-t(x)^(-1/ξ)) for t(x) > 0 F(x) = 0 for t(x) ≤ 0 (below lower bound) where t(x) = 1 + ξ(x - μ)/σ

Support:

x ∈ [μ - σ/ξ, +∞) where μ - σ/ξ is the lower bound

Tail Behavior:

P(X > x) ~ C * x^(-1/ξ) as x → ∞, where C is a constant This gives tail index α = 1/ξ, with smaller ξ meaning heavier tails

Parameters:

• loc (μ): Location parameter ∈ ℝ
• scale (σ): Scale parameter > 0
• shape (ξ): Shape parameter > 0 (enforced for Type II)

Statistical Properties:

• Mean exists for ξ < 1: E[X] = μ + (σ/ξ)(Γ(1-ξ) - 1)
• Variance exists for ξ < 1/2: Var[X] = (σ²/ξ²)(Γ(1-2ξ) - Γ²(1-ξ))
• k-th moment exists for ξ < 1/k
• Heavy right tail makes it suitable for extreme event modeling

Examples:

>>> import jax
>>> import jax.numpy as jnp
>>>
>>> # Create Fréchet Type II GEVD with different tail heaviness
>>> light_tail = FrechetType2GEVD(loc=0.0, scale=1.0, shape=0.1)  # α = 10
>>> heavy_tail = FrechetType2GEVD(loc=0.0, scale=1.0, shape=0.8)  # α = 1.25
>>>
>>> # Key properties
>>> print(f"Lower bound: {heavy_tail.lower_bound()}")
>>> print(f"Tail index: {heavy_tail.tail_index()}")
>>> print(f"Mean exists: {heavy_tail.concentration < 1}")
>>>
>>> # Sample and evaluate
>>> key = jax.random.PRNGKey(42)
>>> samples = heavy_tail.sample(key, sample_shape=(1000,))
>>> log_probs = heavy_tail.log_prob(samples)

Initialize the Fréchet Type II GEVD.

Parameters:

Name	Type	Description	Default
loc	ArrayLike	Location parameter μ (real number)	0.0
scale	ArrayLike	Scale parameter σ (positive real number)	1.0
concentration	ArrayLike | None	Shape parameter ξ (positive real number, ξ > 0)	None
shape	ArrayLike | None	Deprecated backward-compatible alias for concentration. Stored as self.concentration; the NumPyro-inherited Distribution.shape() method is kept callable.	None
validate_args	bool | None	Whether to validate input arguments	None

Raises:

Type	Description
ValueError	If validate_args=True and concentration <= 0 (enforced by the arg_constraints entry). The previous Python branch if jnp.any(shape <= 0) was tracer-unsafe: building the distribution inside a jit-compiled NumPyro model would raise a tracer concretization error before sampling/log-prob could run. arg_constraints already expresses the domain.

Source code in src/xtremax/distributions/frechet.py

def __init__(
    self,
    loc: ArrayLike = 0.0,
    scale: ArrayLike = 1.0,
    concentration: ArrayLike | None = None,
    shape: ArrayLike | None = None,
    validate_args: bool | None = None,
):
    """
    Initialize the Fréchet Type II GEVD.

    Args:
        loc: Location parameter μ (real number)
        scale: Scale parameter σ (positive real number)
        concentration: Shape parameter ξ (positive real number, ξ > 0)
        shape: Deprecated backward-compatible alias for ``concentration``.
            Stored as ``self.concentration``; the NumPyro-inherited
            ``Distribution.shape()`` method is kept callable.
        validate_args: Whether to validate input arguments

    Raises:
        ValueError: If ``validate_args=True`` and ``concentration <= 0``
            (enforced by the ``arg_constraints`` entry). The
            previous Python branch ``if jnp.any(shape <= 0)`` was
            tracer-unsafe: building the distribution inside a
            ``jit``-compiled NumPyro model would raise a tracer
            concretization error before sampling/log-prob could
            run. ``arg_constraints`` already expresses the domain.
    """
    if shape is not None:
        if concentration is not None:
            raise ValueError(
                "Pass only one of 'concentration' or the deprecated 'shape' alias."
            )
        warnings.warn(
            "'shape' is deprecated; use 'concentration' instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        concentration = shape

    if concentration is None:
        concentration = 0.1

    self.loc, self.scale, self.concentration = promote_shapes(
        loc, scale, concentration
    )

    # Determine batch shape from broadcasted parameters
    batch_shape = lax.broadcast_shapes(
        jnp.shape(self.loc), jnp.shape(self.scale), jnp.shape(self.concentration)
    )

    super().__init__(batch_shape=batch_shape, validate_args=validate_args)

support: constraints.Constraint property

Return the support constraint for Fréchet Type II GEVD.

Returns:

Type	Description
Constraint	greater_than_eq(μ - σ/ξ) — lower-bounded, closed at the
Constraint	endpoint (matching the closed interval membership this
Constraint	class always had; the density vanishes there so log_prob
Constraint	is -inf). The old interval(lo, +inf) construction made
Constraint	biject_to(support) return inf, breaking NUTS/SVI models
Constraint	with a Fréchet latent site; greater_than_eq has a
Constraint	registered bijector that handles the unbounded tail.

mean: jnp.ndarray property

Mean. Thin wrapper for :func:~xtremax.primitives.frechet.frechet_mean.

mode: jnp.ndarray property

Compute the mode of the Fréchet Type II GEVD.

The mode is: mode = μ + (σ/ξ) * ((1+ξ)^(-ξ) - 1)

This represents the most likely value in the heavy-tailed distribution.

Returns:

Type	Description
ndarray	Mode of the distribution

variance: jnp.ndarray property

Compute the variance of the Fréchet Type II GEVD.

The variance exists when ξ < 1/2: Var[X] = (σ²/ξ²) * (Γ(1-2ξ) - Γ²(1-ξ))

For very heavy tails (ξ ≥ 1/2), the variance is infinite.

Returns:

Type	Description
ndarray	Variance or +∞ when it doesn't exist (ξ ≥ 1/2)

sample(key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray

Generate samples from the Fréchet Type II GEVD using inverse transform sampling.

The sampling uses the quantile function: Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

Since ξ > 0, this generates samples with heavy right tails.

Parameters:

Name	Type	Description	Default
key	ndarray	JAX random key for sampling	required
sample_shape	tuple	Shape of samples to generate	()

Returns:

Type	Description
ndarray	Array of samples from the Fréchet Type II GEVD, all above lower bound

Source code in src/xtremax/distributions/frechet.py

def sample(self, key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray:
    """
    Generate samples from the Fréchet Type II GEVD using inverse transform sampling.

    The sampling uses the quantile function:
    Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

    Since ξ > 0, this generates samples with heavy right tails.

    Args:
        key: JAX random key for sampling
        sample_shape: Shape of samples to generate

    Returns:
        Array of samples from the Fréchet Type II GEVD, all above lower bound
    """
    check_prng_key(key)
    shape = sample_shape + self.batch_shape

    # JAX's Uniform(0, 1) sampler can emit exact 0 or 1 at the
    # endpoints; passing those to icdf yields -inf/+inf and poisons
    # downstream computations. Clamp away from the endpoints.
    uniform_samples = dist.Uniform(0.0, 1.0).sample(key, shape)
    eps = jnp.finfo(uniform_samples.dtype).eps
    uniform_samples = jnp.clip(uniform_samples, eps, 1.0 - eps)

    # Apply inverse CDF transformation
    return self.icdf(uniform_samples)

log_prob(value: jnp.ndarray) -> jnp.ndarray

Log PDF. Thin wrapper for frechet_log_prob.

Source code in src/xtremax/distributions/frechet.py

@validate_sample
def log_prob(self, value: jnp.ndarray) -> jnp.ndarray:
    """Log PDF. Thin wrapper for ``frechet_log_prob``."""
    return frechet_log_prob(value, self.loc, self.scale, self.concentration)

cdf(value: jnp.ndarray) -> jnp.ndarray

CDF. Thin wrapper for :func:~xtremax.primitives.frechet.frechet_cdf.

Source code in src/xtremax/distributions/frechet.py

def cdf(self, value: jnp.ndarray) -> jnp.ndarray:
    """CDF. Thin wrapper for :func:`~xtremax.primitives.frechet.frechet_cdf`."""
    return frechet_cdf(value, self.loc, self.scale, self.concentration)

icdf(q: jnp.ndarray) -> jnp.ndarray

Quantile function. Thin wrapper for frechet_icdf.

Source code in src/xtremax/distributions/frechet.py

def icdf(self, q: jnp.ndarray) -> jnp.ndarray:
    """Quantile function. Thin wrapper for ``frechet_icdf``."""
    return frechet_icdf(q, self.loc, self.scale, self.concentration)

upper_bound() -> jnp.ndarray

Compute the upper bound of the support.

For Fréchet Type II: upper bound = +∞

Returns:

Type	Description
ndarray	Upper bound (+∞ for Type II)

Source code in src/xtremax/distributions/frechet.py

def upper_bound(self) -> jnp.ndarray:
    """
    Compute the upper bound of the support.

    For Fréchet Type II: upper bound = +∞

    Returns:
        Upper bound (+∞ for Type II)
    """
    return jnp.full_like(self.loc, jnp.inf)

lower_bound() -> jnp.ndarray

Compute the lower bound of the support.

For Fréchet Type II (ξ > 0): lower bound = μ - σ/ξ

Returns:

Type	Description
ndarray	Lower bound of the distribution support

Source code in src/xtremax/distributions/frechet.py

def lower_bound(self) -> jnp.ndarray:
    """
    Compute the lower bound of the support.

    For Fréchet Type II (ξ > 0): lower bound = μ - σ/ξ

    Returns:
        Lower bound of the distribution support
    """
    return self.loc - self.scale / self.concentration

kurtosis() -> jnp.ndarray

Compute the excess kurtosis of the Fréchet Type II GEVD.

Excess kurtosis exists when ξ < 1/4 and involves fourth-order moments. For heavy-tailed distributions, kurtosis captures the tail heaviness.

Returns:

Type	Description
ndarray	Excess kurtosis, or NaN when it doesn't exist (ξ ≥ 1/4)

Source code in src/xtremax/distributions/frechet.py

def kurtosis(self) -> jnp.ndarray:
    """
    Compute the excess kurtosis of the Fréchet Type II GEVD.

    Excess kurtosis exists when ξ < 1/4 and involves fourth-order moments.
    For heavy-tailed distributions, kurtosis captures the tail heaviness.

    Returns:
        Excess kurtosis, or NaN when it doesn't exist (ξ ≥ 1/4)
    """
    return gev_excess_kurtosis(self.concentration)

skew() -> jnp.ndarray

Compute the skewness of the Fréchet Type II GEVD.

Skewness exists when ξ < 1/3 and involves third-order moments.

For Fréchet Type II, skewness is typically positive due to the heavy right tail creating right-skewed distributions.

Returns:

Type	Description
ndarray	Skewness, or NaN when it doesn't exist (ξ ≥ 1/3)

Source code in src/xtremax/distributions/frechet.py

def skew(self) -> jnp.ndarray:
    """
    Compute the skewness of the Fréchet Type II GEVD.

    Skewness exists when ξ < 1/3 and involves third-order moments.

    For Fréchet Type II, skewness is typically positive due to the
    heavy right tail creating right-skewed distributions.

    Returns:
        Skewness, or NaN when it doesn't exist (ξ ≥ 1/3)
    """
    return gev_skewness(self.concentration)

entropy() -> jnp.ndarray

Compute the differential entropy of the Fréchet Type II GEVD.

The entropy is: H = log(σ) + 1 + γ * (1 + ξ)

where γ is the Euler-Mascheroni constant ≈ 0.5772. This is the standard GEV entropy formula evaluated at the Fréchet branch (ξ > 0); at ξ = 0 it reduces to the Gumbel entropy log σ + 1 + γ.

Returns:

Type	Description
ndarray	Differential entropy in nats

Source code in src/xtremax/distributions/frechet.py

def entropy(self) -> jnp.ndarray:
    """
    Compute the differential entropy of the Fréchet Type II GEVD.

    The entropy is:
    H = log(σ) + 1 + γ * (1 + ξ)

    where γ is the Euler-Mascheroni constant ≈ 0.5772. This is the
    standard GEV entropy formula evaluated at the Fréchet branch
    (ξ > 0); at ξ = 0 it reduces to the Gumbel entropy `log σ + 1 + γ`.

    Returns:
        Differential entropy in nats
    """
    scale, shape = self.scale, self.concentration
    euler_gamma = 0.5772156649015329

    return jnp.log(scale) + 1.0 + euler_gamma * (1.0 + shape)

survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the survival function S(x) = 1 - F(x).

For Fréchet Type II, this exhibits power-law decay: S(x) ~ C * x^(-1/ξ) for large x, where C is a constant.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the survival function	required

Returns:

Type	Description
ndarray	Survival probabilities (power-law tail behavior)

Source code in src/xtremax/distributions/frechet.py

def survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the survival function S(x) = 1 - F(x).

    For Fréchet Type II, this exhibits power-law decay:
    S(x) ~ C * x^(-1/ξ) for large x, where C is a constant.

    Args:
        value: Points at which to evaluate the survival function

    Returns:
        Survival probabilities (power-law tail behavior)
    """
    # Computed as ``-expm1(log F)`` rather than ``1 - F`` so the
    # far right tail stays resolvable: once ``F ≈ 1 - 1e-8``, float32
    # rounds ``1 - F`` to zero while ``-expm1(-t^(-1/ξ))`` preserves
    # the tail mass down to subnormals, keeping
    # ``exceedance_probability`` / ``hazard_rate`` finite.
    loc, scale, shape = self.loc, self.scale, self.concentration
    z = (value - loc) / scale
    t = 1.0 + shape * z
    valid = t > 0.0
    log_cdf = -jnp.power(jnp.where(valid, t, 1.0), -1.0 / shape)
    surv_inside = -jnp.expm1(log_cdf)
    # Below the lower support bound (t ≤ 0 with ξ > 0), F(x) = 0 so
    # S(x) = 1.
    return jnp.where(valid, surv_inside, 1.0)

log_survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the log survival function log S(x) = log(1 - F(x)).

For heavy-tailed distributions, this provides better numerical stability for extreme quantiles.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate log survival function	required

Returns:

Type	Description
ndarray	Log survival probabilities

Source code in src/xtremax/distributions/frechet.py

def log_survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the log survival function log S(x) = log(1 - F(x)).

    For heavy-tailed distributions, this provides better numerical
    stability for extreme quantiles.

    Args:
        value: Points at which to evaluate log survival function

    Returns:
        Log survival probabilities
    """
    loc, scale, shape = self.loc, self.scale, self.concentration
    z = (value - loc) / scale
    t = 1.0 + shape * z

    valid = t > 0.0
    # With F(x) = exp(-t^(-1/ξ)), log F(x) = -t^(-1/ξ). The survival
    # is 1 - F(x) = -expm1(-t^(-1/ξ)), and log S = log(-expm1(...)).
    log_f = -jnp.power(jnp.where(valid, t, 1.0), -1.0 / shape)
    log_surv_inside = jnp.log(-jnp.expm1(log_f))
    # Below the lower bound (t <= 0), F(x) = 0 so log S = 0.
    return jnp.where(valid, log_surv_inside, 0.0)

hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the hazard rate h(x) = f(x) / S(x).

For Fréchet Type II, the hazard rate approaches zero as x → ∞, indicating decreasing failure rates for extreme values.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the hazard rate	required

Returns:

Type	Description
ndarray	Hazard rate values

Source code in src/xtremax/distributions/frechet.py

def hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the hazard rate h(x) = f(x) / S(x).

    For Fréchet Type II, the hazard rate approaches zero as x → ∞,
    indicating decreasing failure rates for extreme values.

    Args:
        value: Points at which to evaluate the hazard rate

    Returns:
        Hazard rate values
    """
    log_hazard = self.log_prob(value) - self.log_survival_function(value)
    return jnp.exp(log_hazard)

cumulative_hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the cumulative hazard rate Λ(x) = -log(S(x)).

This represents the accumulated hazard up to time x. For Fréchet Type II, this grows slower than exponentially.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the cumulative hazard rate	required

Returns:

Type	Description
ndarray	Cumulative hazard rate values

Source code in src/xtremax/distributions/frechet.py

def cumulative_hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the cumulative hazard rate Λ(x) = -log(S(x)).

    This represents the accumulated hazard up to time x.
    For Fréchet Type II, this grows slower than exponentially.

    Args:
        value: Points at which to evaluate the cumulative hazard rate

    Returns:
        Cumulative hazard rate values
    """
    return -self.log_survival_function(value)

return_level(return_period: float | jnp.ndarray) -> jnp.ndarray

Return level. Thin wrapper for frechet_return_level.

Source code in src/xtremax/distributions/frechet.py

def return_level(self, return_period: float | jnp.ndarray) -> jnp.ndarray:
    """Return level. Thin wrapper for ``frechet_return_level``."""
    return frechet_return_level(
        return_period, self.loc, self.scale, self.concentration
    )

tail_index() -> jnp.ndarray

Compute the tail index α = 1/ξ for Fréchet Type II GEVD.

The tail index characterizes the power-law tail behavior: P(X > x) ~ x^(-α) as x → ∞

Lower values of α (higher ξ) indicate heavier tails.

Returns:

Type	Description
ndarray	Tail index α = 1/ξ

Source code in src/xtremax/distributions/frechet.py

def tail_index(self) -> jnp.ndarray:
    """
    Compute the tail index α = 1/ξ for Fréchet Type II GEVD.

    The tail index characterizes the power-law tail behavior:
    P(X > x) ~ x^(-α) as x → ∞

    Lower values of α (higher ξ) indicate heavier tails.

    Returns:
        Tail index α = 1/ξ
    """
    return 1.0 / self.concentration

exceedance_probability(threshold: jnp.ndarray) -> jnp.ndarray

Compute probability of exceeding a threshold: P(X > threshold).

For large thresholds, this exhibits power-law decay: P(X > x) ~ C * x^(-1/ξ)

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value	required

Returns:

Type	Description
ndarray	Exceedance probabilities (power-law tail behavior)

Source code in src/xtremax/distributions/frechet.py

def exceedance_probability(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute probability of exceeding a threshold: P(X > threshold).

    For large thresholds, this exhibits power-law decay:
    P(X > x) ~ C * x^(-1/ξ)

    Args:
        threshold: Threshold value

    Returns:
        Exceedance probabilities (power-law tail behavior)
    """
    return self.survival_function(threshold)

conditional_excess_mean(threshold: jnp.ndarray) -> jnp.ndarray

Mean excess :math:E[X - u \mid X > u] for Fréchet Type II GEVD.

Computed via quantile-space quadrature, same machinery as the general GEV class. The previous GPD linear approximation :math:(\sigma + \xi(u - \mu)) / (1 - \xi) is only the POT asymptote, not the finite-threshold mean excess.

Returns NaN for :math:\xi \ge 1 (heavy tails where the mean does not exist).

Source code in src/xtremax/distributions/frechet.py

def conditional_excess_mean(self, threshold: jnp.ndarray) -> jnp.ndarray:
    r"""Mean excess :math:`E[X - u \mid X > u]` for Fréchet Type II GEVD.

    Computed via quantile-space quadrature, same machinery as the
    general GEV class. The previous GPD linear approximation
    :math:`(\sigma + \xi(u - \mu)) / (1 - \xi)` is only the POT
    asymptote, not the finite-threshold mean excess.

    Returns NaN for :math:`\xi \ge 1` (heavy tails where the mean
    does not exist).
    """
    threshold_arr = jnp.asarray(threshold)
    shape = self.concentration

    # Log-space tail-probability quadrature (see GEVD implementation
    # for the derivation). This replaces the earlier linear-p grid
    # whose (a) normaliser used truncated mass instead of S(u), and
    # (b) linear spacing underresolved heavy Fréchet tails.
    n_grid = 1024
    p0 = self.cdf(threshold_arr)
    # Joint batch of threshold and parameters; per-batch quantities are
    # broadcast to it before gaining the trailing grid axis.
    batch = jnp.shape(p0)
    s_u = 1.0 - p0
    # Cap away from 1.0 so p_grid stays strictly in (0, 1). See
    # GEVD.conditional_excess_mean.
    s_u_safe = jnp.clip(s_u, 1e-12, 1.0 - 1e-6)
    log_s_u = jnp.log(s_u_safe)
    # Lower log-y endpoint: step ~20 e-folds below log_s_u so the
    # grid always captures the bulk of the tail-integrand mass AND
    # stays strictly ascending. The truncated mass below y_c —
    # which behaves like y^{1-ξ} and reached ~14% of the total at
    # ξ = 0.9 — is added back analytically below rather than by
    # widening the span (which would underflow float32 exp). The
    # ``-jnp.log1p(-y)`` step below is accurate to subnormals in
    # float32. See GEVD.conditional_excess_mean for details.
    span = jnp.asarray(20.0, dtype=log_s_u.dtype)
    log_y_min = log_s_u - span

    unit = jnp.linspace(0.0, 1.0, n_grid)
    log_s_u_exp = jnp.expand_dims(log_s_u, axis=-1)
    log_y_min_exp = jnp.expand_dims(jnp.broadcast_to(log_y_min, batch), axis=-1)
    v_grid = log_y_min_exp * (1.0 - unit) + log_s_u_exp * unit
    y_grid = jnp.exp(v_grid)
    # Compute ``x = F⁻¹(1 - y)`` directly from ``y`` using
    # ``-log1p(-y)``. For ``y`` near ``eps(float32)`` the
    # round-tripped ``p = 1 - y`` keeps only ~1 ulp of precision and
    # ``log(p)`` inside ``frechet_icdf`` loses that entirely, so the
    # quadrature sign-flips at very high thresholds. See
    # GEVD.conditional_excess_mean for details. ξ > 0 is enforced
    # by ``arg_constraints``, so the Gumbel branch is unreachable.
    # Parameters are broadcast to the joint batch and given the same
    # trailing grid axis — batched (batch,) parameters against a
    # (..., n_grid) grid otherwise raise a broadcasting error.
    neg_log_q = -jnp.log1p(-y_grid)
    loc_e = jnp.expand_dims(jnp.broadcast_to(self.loc, batch), axis=-1)
    scale_e = jnp.expand_dims(jnp.broadcast_to(self.scale, batch), axis=-1)
    shape_e = jnp.expand_dims(jnp.broadcast_to(shape, batch), axis=-1)
    x_grid = loc_e + (scale_e / shape_e) * (jnp.power(neg_log_q, -shape_e) - 1.0)

    integrand = x_grid * y_grid
    numerator = jnp.trapezoid(integrand, x=v_grid, axis=-1)

    # Closed-form remainder for the truncated (0, y_c] tail with
    # y_c = S(u)·e⁻²⁰, where x(y) ≈ (μ - σ/ξ) + (σ/ξ)·y^{-ξ}:
    #   ∫₀^{y_c} x dy = (μ - σ/ξ)·y_c + (σ/ξ)·y_c^{1-ξ}/(1-ξ).
    # Dropping this was the ~14%-at-ξ=0.9 truncation bias. ξ ≥ 1 is
    # masked to NaN below; sanitize its 1-ξ ≤ 0 divisor.
    loc_b = jnp.broadcast_to(self.loc, batch)
    scale_b = jnp.broadcast_to(self.scale, batch)
    shape_b = jnp.broadcast_to(shape, batch)
    y_c = jnp.exp(log_y_min)
    one_minus_xi = jnp.where(shape_b >= 1.0, 1.0, 1.0 - shape_b)
    remainder = (loc_b - scale_b / shape_b) * y_c + (scale_b / shape_b) * jnp.power(
        y_c, one_minus_xi
    ) / one_minus_xi
    numerator = numerator + remainder

    mean_conditional = numerator / s_u_safe
    mean_excess = mean_conditional - threshold_arr

    mean_exists = shape < 1.0
    valid = mean_exists & (s_u > 1e-12)
    return jnp.where(valid, mean_excess, jnp.nan)

hill_estimator_threshold(quantile: float = 0.95) -> jnp.ndarray

Suggest threshold for Hill tail index estimation.

For Fréchet distributions, the threshold should be high enough to be in the tail region but low enough to have sufficient data.

Parameters:

Name	Type	Description	Default
quantile	float	Quantile level for threshold selection	0.95

Returns:

Type	Description
ndarray	Suggested threshold for tail analysis

Source code in src/xtremax/distributions/frechet.py

def hill_estimator_threshold(self, quantile: float = 0.95) -> jnp.ndarray:
    """
    Suggest threshold for Hill tail index estimation.

    For Fréchet distributions, the threshold should be high enough
    to be in the tail region but low enough to have sufficient data.

    Args:
        quantile: Quantile level for threshold selection

    Returns:
        Suggested threshold for tail analysis
    """
    return self.icdf(quantile)

pareto_tail_approx(threshold: jnp.ndarray) -> tuple[jnp.ndarray, jnp.ndarray]

Generalized Pareto approximation for excesses over threshold.

For Fréchet GEVD, excesses over high thresholds approximately follow GPD(σ, ξ) where: σ = σ + ξ(u - μ) (modified scale) ξ remains the same (shape parameter)

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value for POT modeling	required

Returns:

Type	Description
tuple[ndarray, ndarray]	Tuple of (GPD scale parameter, GPD shape parameter)

Source code in src/xtremax/distributions/frechet.py

def pareto_tail_approx(
    self, threshold: jnp.ndarray
) -> tuple[jnp.ndarray, jnp.ndarray]:
    """
    Generalized Pareto approximation for excesses over threshold.

    For Fréchet GEVD, excesses over high thresholds approximately
    follow GPD(σ*, ξ) where:
    σ* = σ + ξ(u - μ) (modified scale)
    ξ remains the same (shape parameter)

    Args:
        threshold: Threshold value for POT modeling

    Returns:
        Tuple of (GPD scale parameter, GPD shape parameter)
    """
    loc, scale, shape = self.loc, self.scale, self.concentration

    # GPD parameters for excesses over threshold
    gpd_scale = scale + shape * (threshold - loc)
    gpd_shape = shape

    return gpd_scale, gpd_shape

value_at_risk(alpha: float) -> jnp.ndarray

Compute Value-at-Risk (VaR) at confidence level α.

VaR_α = Q(α) where Q is the quantile function. For heavy-tailed distributions, VaR can be very large.

Parameters:

Name	Type	Description	Default
alpha	float	Confidence level (e.g., 0.95, 0.99, 0.999)	required

Returns:

Type	Description
ndarray	Value-at-Risk at confidence level α

Source code in src/xtremax/distributions/frechet.py

def value_at_risk(self, alpha: float) -> jnp.ndarray:
    """
    Compute Value-at-Risk (VaR) at confidence level α.

    VaR_α = Q(α) where Q is the quantile function.
    For heavy-tailed distributions, VaR can be very large.

    Args:
        alpha: Confidence level (e.g., 0.95, 0.99, 0.999)

    Returns:
        Value-at-Risk at confidence level α
    """
    return self.icdf(alpha)

expected_shortfall(alpha: float) -> jnp.ndarray

Compute Expected Shortfall (Conditional VaR) at confidence level α.

ES_α = E[X | X > VaR_α] for Fréchet distributions. This captures the expected loss beyond the VaR threshold.

Parameters:

Name	Type	Description	Default
alpha	float	Confidence level (e.g., 0.95, 0.99, 0.999)	required

Returns:

Type	Description
ndarray	Expected Shortfall at confidence level α

Source code in src/xtremax/distributions/frechet.py

def expected_shortfall(self, alpha: float) -> jnp.ndarray:
    """
    Compute Expected Shortfall (Conditional VaR) at confidence level α.

    ES_α = E[X | X > VaR_α] for Fréchet distributions.
    This captures the expected loss beyond the VaR threshold.

    Args:
        alpha: Confidence level (e.g., 0.95, 0.99, 0.999)

    Returns:
        Expected Shortfall at confidence level α
    """
    var_alpha = self.value_at_risk(alpha)
    excess_mean = self.conditional_excess_mean(var_alpha)

    return var_alpha + excess_mean

power_law_scaling(scale_factor: float) -> jnp.ndarray

Compute scaling relationship for power-law tails.

For Fréchet distributions: if X ~ Fréchet(μ,σ,ξ), then aX has tail behavior ~ (ax)^(-1/ξ) = a^(-1/ξ) * x^(-1/ξ)

Parameters:

Name	Type	Description	Default
scale_factor	float	Scaling factor a > 0	required

Returns:

Type	Description
ndarray	Relative change in tail probabilities under scaling

Source code in src/xtremax/distributions/frechet.py

def power_law_scaling(self, scale_factor: float) -> jnp.ndarray:
    """
    Compute scaling relationship for power-law tails.

    For Fréchet distributions: if X ~ Fréchet(μ,σ,ξ),
    then aX has tail behavior ~ (ax)^(-1/ξ) = a^(-1/ξ) * x^(-1/ξ)

    Args:
        scale_factor: Scaling factor a > 0

    Returns:
        Relative change in tail probabilities under scaling
    """
    return jnp.power(scale_factor, -1.0 / self.concentration)

expand(batch_shape: tuple[int, ...]) -> dist.Distribution

Expand to batch_shape by reconstructing via __init__.

Source code in src/xtremax/distributions/frechet.py

def expand(self, batch_shape: tuple[int, ...]) -> dist.Distribution:
    """Expand to ``batch_shape`` by reconstructing via ``__init__``."""
    batch_shape = tuple(batch_shape)
    if batch_shape == self.batch_shape:
        return self
    return type(self)(
        loc=jnp.broadcast_to(self.loc, batch_shape),
        scale=jnp.broadcast_to(self.scale, batch_shape),
        concentration=jnp.broadcast_to(self.concentration, batch_shape),
        validate_args=self._validate_args,
    )

WeibullType3GEVD(loc: ArrayLike = 0.0, scale: ArrayLike = 1.0, concentration: ArrayLike | None = None, shape: ArrayLike | None = None, validate_args: bool | None = None)

Bases: Distribution

Weibull Type III Generalized Extreme Value Distribution for NumPyro.

The Weibull Type III (reversed Weibull) is the GEVD with shape parameter ξ < 0, characterized by a bounded upper tail. This distribution is fundamental in extreme value theory for modeling bounded maxima and reliability applications with upper failure limits.

Key Characteristics: - Bounded upper tail: x ≤ μ - σ/ξ (since ξ < 0) - Light tail behavior with exponential decay - Often used for modeling material strength limits, maximum system capacity - Opposite tail behavior to Fréchet (Type II)

Probability Density Function:

f(x) = (1/σ) * t(x)^(-(1/ξ + 1)) * exp(-t(x)^(-1/ξ)) where t(x) = 1 + ξ(x - μ)/σ and ξ < 0

Cumulative Distribution Function:

F(x) = exp(-t(x)^(-1/ξ)) for t(x) > 0 F(x) = 1 for t(x) ≤ 0 (above upper bound) where t(x) = 1 + ξ(x - μ)/σ

Support:

x ∈ (-∞, μ - σ/ξ] where μ - σ/ξ is the upper bound

Parameters:

• loc (μ): Location parameter ∈ ℝ
• scale (σ): Scale parameter > 0
• shape (ξ): Shape parameter < 0 (enforced for Type III)

Statistical Properties:

• Mean exists for ξ > -1: E[X] = μ + (σ/ξ)(Γ(1-ξ) - 1)
• Variance exists for ξ > -1/2: Var[X] = (σ²/ξ²)(Γ(1-2ξ) - Γ²(1-ξ))
• Upper bounded support makes it suitable for reliability analysis

Examples:

>>> import jax
>>> import jax.numpy as jnp
>>>
>>> # Create Weibull Type III GEVD
>>> weibull_gevd = WeibullType3GEVD(loc=10.0, scale=2.0, shape=-0.3)
>>>
>>> # Key properties
>>> print(f"Upper bound: {weibull_gevd.upper_bound()}")
>>> print(f"Mean: {weibull_gevd.mean}")
>>> print(f"Mode: {weibull_gevd.mode}")
>>>
>>> # Sample and evaluate
>>> key = jax.random.PRNGKey(42)
>>> samples = weibull_gevd.sample(key, sample_shape=(1000,))
>>> log_probs = weibull_gevd.log_prob(samples)

Initialize the Weibull Type III GEVD.

Parameters:

Name	Type	Description	Default
loc	ArrayLike	Location parameter μ (real number)	0.0
scale	ArrayLike	Scale parameter σ (positive real number)	1.0
concentration	ArrayLike | None	Shape parameter ξ (negative real number, ξ < 0)	None
shape	ArrayLike | None	Deprecated backward-compatible alias for concentration. Stored as self.concentration; the NumPyro-inherited Distribution.shape() method is kept callable.	None
validate_args	bool | None	Whether to validate input arguments	None

Raises:

Type	Description
ValueError	If validate_args=True and concentration >= 0 (enforced by the arg_constraints entry). The previous Python branch if jnp.any(shape >= 0) was tracer-unsafe: building the distribution inside a jit-compiled NumPyro model would raise a tracer concretization error before sampling/log-prob could run. arg_constraints already expresses the domain.

Source code in src/xtremax/distributions/weibull.py

def __init__(
    self,
    loc: ArrayLike = 0.0,
    scale: ArrayLike = 1.0,
    concentration: ArrayLike | None = None,
    shape: ArrayLike | None = None,
    validate_args: bool | None = None,
):
    """
    Initialize the Weibull Type III GEVD.

    Args:
        loc: Location parameter μ (real number)
        scale: Scale parameter σ (positive real number)
        concentration: Shape parameter ξ (negative real number, ξ < 0)
        shape: Deprecated backward-compatible alias for ``concentration``.
            Stored as ``self.concentration``; the NumPyro-inherited
            ``Distribution.shape()`` method is kept callable.
        validate_args: Whether to validate input arguments

    Raises:
        ValueError: If ``validate_args=True`` and ``concentration >= 0``
            (enforced by the ``arg_constraints`` entry). The
            previous Python branch ``if jnp.any(shape >= 0)`` was
            tracer-unsafe: building the distribution inside a
            ``jit``-compiled NumPyro model would raise a tracer
            concretization error before sampling/log-prob could
            run. ``arg_constraints`` already expresses the domain.
    """
    if shape is not None:
        if concentration is not None:
            raise ValueError(
                "Pass only one of 'concentration' or the deprecated 'shape' alias."
            )
        warnings.warn(
            "'shape' is deprecated; use 'concentration' instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        concentration = shape

    if concentration is None:
        concentration = -0.1

    self.loc, self.scale, self.concentration = promote_shapes(
        loc, scale, concentration
    )

    # Determine batch shape from broadcasted parameters
    batch_shape = lax.broadcast_shapes(
        jnp.shape(self.loc), jnp.shape(self.scale), jnp.shape(self.concentration)
    )

    super().__init__(batch_shape=batch_shape, validate_args=validate_args)

support: constraints.Constraint property

Return the support constraint for Weibull Type III GEVD.

Returns:

Type	Description
Constraint	less_than_eq(μ - σ/ξ) — upper-bounded, closed at the
Constraint	endpoint (the endpoint density is finite at ξ = -1 and diverges
Constraint	for ξ < -1, so the endpoint is a genuine support point). The old
Constraint	interval(-inf, hi) construction made biject_to(support)
Constraint	return nan/inf, breaking NUTS/SVI models with a Weibull latent
Constraint	site; less_than_eq has a registered bijector that handles
Constraint	the unbounded lower tail.

mean: jnp.ndarray property

Mean. Thin wrapper for :func:~xtremax.primitives.weibull.weibull_mean.

mode: jnp.ndarray property

Compute the mode of the Weibull Type III GEVD.

The mode is: mode = μ + (σ/ξ) * ((1+ξ)^(-ξ) - 1)

This interior stationary point only exists for -1 < ξ < 0. For ξ ≤ -1, the density increases all the way to the finite upper endpoint, so the mode is the upper bound itself.

Returns:

Type	Description
ndarray	Mode of the distribution

variance: jnp.ndarray property

Compute the variance of the Weibull Type III GEVD.

Var[X] = (σ²/ξ²) * (Γ(1-2ξ) - Γ²(1-ξ))

For Weibull Type III, ξ < 0 means the distribution is bounded above, so all moments exist for every valid ξ in the Type III range (no guard is needed). The earlier guard ξ > -1/2 was mistakenly inherited from the Fréchet moment-existence condition reflected to negative ξ; it silently returned NaN for valid shapes like ξ = -1.

Returns:

Type	Description
ndarray	Variance (always finite on the Type III domain).

sample(key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray

Generate samples from the Weibull Type III GEVD via inverse transform sampling.

The sampling uses the quantile function: Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

Since ξ < 0, this becomes: Q(p) = μ + (σ/|ξ|) * (1 - (-ln(p))^(|ξ|))

Parameters:

Name	Type	Description	Default
key	ndarray	JAX random key for sampling	required
sample_shape	tuple	Shape of samples to generate	()

Returns:

Type	Description
ndarray	Array of samples from the Weibull Type III GEVD, all bounded above

Source code in src/xtremax/distributions/weibull.py

def sample(self, key: jnp.ndarray, sample_shape: tuple = ()) -> jnp.ndarray:
    """
    Generate samples from the Weibull Type III GEVD via inverse transform sampling.

    The sampling uses the quantile function:
    Q(p) = μ + (σ/ξ) * ((-ln(p))^(-ξ) - 1)

    Since ξ < 0, this becomes:
    Q(p) = μ + (σ/|ξ|) * (1 - (-ln(p))^(|ξ|))

    Args:
        key: JAX random key for sampling
        sample_shape: Shape of samples to generate

    Returns:
        Array of samples from the Weibull Type III GEVD, all bounded above
    """
    check_prng_key(key)
    shape = sample_shape + self.batch_shape

    # JAX's Uniform(0, 1) sampler can emit exact 0 or 1 at the
    # endpoints; passing those to icdf yields -inf/+inf and poisons
    # downstream computations. Clamp away from the endpoints.
    uniform_samples = dist.Uniform(0.0, 1.0).sample(key, shape)
    eps = jnp.finfo(uniform_samples.dtype).eps
    uniform_samples = jnp.clip(uniform_samples, eps, 1.0 - eps)

    # Apply inverse CDF transformation
    return self.icdf(uniform_samples)

log_prob(value: jnp.ndarray) -> jnp.ndarray

Log PDF. Thin wrapper for weibull_log_prob.

Source code in src/xtremax/distributions/weibull.py

@validate_sample
def log_prob(self, value: jnp.ndarray) -> jnp.ndarray:
    """Log PDF. Thin wrapper for ``weibull_log_prob``."""
    return weibull_log_prob(value, self.loc, self.scale, self.concentration)

cdf(value: jnp.ndarray) -> jnp.ndarray

CDF. Thin wrapper for :func:~xtremax.primitives.weibull.weibull_cdf.

Source code in src/xtremax/distributions/weibull.py

def cdf(self, value: jnp.ndarray) -> jnp.ndarray:
    """CDF. Thin wrapper for :func:`~xtremax.primitives.weibull.weibull_cdf`."""
    return weibull_cdf(value, self.loc, self.scale, self.concentration)

icdf(q: jnp.ndarray) -> jnp.ndarray

Quantile function. Thin wrapper for weibull_icdf.

Source code in src/xtremax/distributions/weibull.py

def icdf(self, q: jnp.ndarray) -> jnp.ndarray:
    """Quantile function. Thin wrapper for ``weibull_icdf``."""
    return weibull_icdf(q, self.loc, self.scale, self.concentration)

upper_bound() -> jnp.ndarray

Compute the upper bound of the support.

For Weibull Type III (ξ < 0): upper bound = μ - σ/ξ

Returns:

Type	Description
ndarray	Upper bound of the distribution support

Source code in src/xtremax/distributions/weibull.py

def upper_bound(self) -> jnp.ndarray:
    """
    Compute the upper bound of the support.

    For Weibull Type III (ξ < 0): upper bound = μ - σ/ξ

    Returns:
        Upper bound of the distribution support
    """
    # Since ξ < 0, -σ/ξ = σ/|ξ| > 0
    return self.loc - self.scale / self.concentration

lower_bound() -> jnp.ndarray

Compute the lower bound of the support.

For Weibull Type III: lower bound = -∞

Returns:

Type	Description
ndarray	Lower bound (-∞ for Type III)

Source code in src/xtremax/distributions/weibull.py

def lower_bound(self) -> jnp.ndarray:
    """
    Compute the lower bound of the support.

    For Weibull Type III: lower bound = -∞

    Returns:
        Lower bound (-∞ for Type III)
    """
    return jnp.full_like(self.loc, -jnp.inf)

kurtosis() -> jnp.ndarray

Compute the excess kurtosis of the Weibull Type III GEVD.

κ = μ₄/σ⁴ - 3

As with variance, every moment exists for every ξ < 0 (bounded support), so no ξ > -1/4 guard is applied.

Returns:

Type	Description
ndarray	Excess kurtosis (always finite on the Type III domain).

Source code in src/xtremax/distributions/weibull.py

def kurtosis(self) -> jnp.ndarray:
    """
    Compute the excess kurtosis of the Weibull Type III GEVD.

    κ = μ₄/σ⁴ - 3

    As with variance, every moment exists for every ξ < 0 (bounded
    support), so no ``ξ > -1/4`` guard is applied.

    Returns:
        Excess kurtosis (always finite on the Type III domain).
    """
    return gev_excess_kurtosis(self.concentration)

skew() -> jnp.ndarray

Compute the skewness of the Weibull Type III GEVD.

For Weibull Type III, skewness is typically negative due to the upper bound creating left-skewed distributions. All moments exist across the full valid ξ < 0 range.

Returns:

Type	Description
ndarray	Skewness (always finite on the Type III domain).

Source code in src/xtremax/distributions/weibull.py

def skew(self) -> jnp.ndarray:
    """
    Compute the skewness of the Weibull Type III GEVD.

    For Weibull Type III, skewness is typically negative due to the
    upper bound creating left-skewed distributions. All moments
    exist across the full valid ξ < 0 range.

    Returns:
        Skewness (always finite on the Type III domain).
    """
    # X = μ + (σ/ξ)(W − 1) with ξ < 0: the (σ/ξ)³ factor in the third
    # central moment flips the sign of the standardized skew. The
    # reduced form carries that sign on its own, so the Type III branch
    # needs no static ``-`` here.
    return gev_skewness(self.concentration)

entropy() -> jnp.ndarray

Compute the differential entropy of the Weibull Type III GEVD.

The entropy is: H = log(σ) + 1 + γ * (1 + ξ)

where γ is the Euler-Mascheroni constant ≈ 0.5772. This is the standard GEV entropy formula evaluated at the Weibull branch (ξ < 0); at ξ = 0 it reduces to the Gumbel entropy log σ + 1 + γ.

Returns:

Type	Description
ndarray	Differential entropy in nats

Source code in src/xtremax/distributions/weibull.py

def entropy(self) -> jnp.ndarray:
    """
    Compute the differential entropy of the Weibull Type III GEVD.

    The entropy is:
    H = log(σ) + 1 + γ * (1 + ξ)

    where γ is the Euler-Mascheroni constant ≈ 0.5772. This is the
    standard GEV entropy formula evaluated at the Weibull branch
    (ξ < 0); at ξ = 0 it reduces to the Gumbel entropy
    ``log σ + 1 + γ``.

    Returns:
        Differential entropy in nats
    """
    scale, shape = self.scale, self.concentration
    euler_gamma = 0.5772156649015329

    return jnp.log(scale) + 1.0 + euler_gamma * (1.0 + shape)

survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the survival function S(x) = 1 - F(x).

For Weibull Type III, this represents the probability of not exceeding the value x, which is relevant for reliability analysis.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the survival function	required

Returns:

Type	Description
ndarray	Survival probabilities

Source code in src/xtremax/distributions/weibull.py

def survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the survival function S(x) = 1 - F(x).

    For Weibull Type III, this represents the probability of not
    exceeding the value x, which is relevant for reliability analysis.

    Args:
        value: Points at which to evaluate the survival function

    Returns:
        Survival probabilities
    """
    # Computed as ``-expm1(log F)`` rather than ``1 - F`` so the
    # near-upper-endpoint regime stays resolvable: once
    # ``F ≈ 1 - 1e-8``, float32 rounds ``1 - F`` to zero and
    # downstream ``log(S)`` in hazard_rate/cumulative_hazard_rate
    # returns ``-inf``/``nan`` even for valid in-support inputs.
    loc, scale, shape = self.loc, self.scale, self.concentration
    z = (value - loc) / scale
    t = 1.0 + shape * z
    valid = t > 0.0
    log_cdf = -jnp.power(jnp.where(valid, t, 1.0), -1.0 / shape)
    surv_inside = -jnp.expm1(log_cdf)
    # Above the upper support bound (t ≤ 0 with ξ < 0), F(x) = 1 so
    # S(x) = 0.
    return jnp.where(valid, surv_inside, 0.0)

log_survival_function(value: jnp.ndarray) -> jnp.ndarray

Compute the log survival function log S(x) = log(1 - F(x)).

Directly evaluates log(-expm1(log F)) to keep the near-upper- endpoint regime numerically stable for hazard/cumulative-hazard consumers.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate log survival function	required

Returns:

Type	Description
ndarray	Log survival probabilities

Source code in src/xtremax/distributions/weibull.py

def log_survival_function(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the log survival function log S(x) = log(1 - F(x)).

    Directly evaluates ``log(-expm1(log F))`` to keep the near-upper-
    endpoint regime numerically stable for hazard/cumulative-hazard
    consumers.

    Args:
        value: Points at which to evaluate log survival function

    Returns:
        Log survival probabilities
    """
    loc, scale, shape = self.loc, self.scale, self.concentration
    z = (value - loc) / scale
    t = 1.0 + shape * z
    valid = t > 0.0
    log_cdf = -jnp.power(jnp.where(valid, t, 1.0), -1.0 / shape)
    log_surv_inside = jnp.log(-jnp.expm1(log_cdf))
    # Above the upper support bound, S(x) = 0 ⇒ log S = -∞.
    return jnp.where(valid, log_surv_inside, -jnp.inf)

hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the hazard rate h(x) = f(x) / S(x).

The hazard rate represents the instantaneous failure rate given survival up to time x. For Weibull Type III, the hazard rate approaches infinity as x approaches the upper bound.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the hazard rate	required

Returns:

Type	Description
ndarray	Hazard rate values

Source code in src/xtremax/distributions/weibull.py

def hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the hazard rate h(x) = f(x) / S(x).

    The hazard rate represents the instantaneous failure rate
    given survival up to time x. For Weibull Type III, the
    hazard rate approaches infinity as x approaches the upper bound.

    Args:
        value: Points at which to evaluate the hazard rate

    Returns:
        Hazard rate values
    """
    log_hazard = self.log_prob(value) - self.log_survival_function(value)
    return jnp.exp(log_hazard)

cumulative_hazard_rate(value: jnp.ndarray) -> jnp.ndarray

Compute the cumulative hazard rate Λ(x) = -log(S(x)).

This represents the accumulated hazard up to time x. For Weibull Type III, this is finite at the upper bound.

Parameters:

Name	Type	Description	Default
value	ndarray	Points at which to evaluate the cumulative hazard rate	required

Returns:

Type	Description
ndarray	Cumulative hazard rate values

Source code in src/xtremax/distributions/weibull.py

def cumulative_hazard_rate(self, value: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the cumulative hazard rate Λ(x) = -log(S(x)).

    This represents the accumulated hazard up to time x.
    For Weibull Type III, this is finite at the upper bound.

    Args:
        value: Points at which to evaluate the cumulative hazard rate

    Returns:
        Cumulative hazard rate values
    """
    return -self.log_survival_function(value)

return_level(return_period: float | jnp.ndarray) -> jnp.ndarray

Return level. Thin wrapper for weibull_return_level.

Source code in src/xtremax/distributions/weibull.py

def return_level(self, return_period: float | jnp.ndarray) -> jnp.ndarray:
    """Return level. Thin wrapper for ``weibull_return_level``."""
    return weibull_return_level(
        return_period, self.loc, self.scale, self.concentration
    )

tail_index() -> jnp.ndarray

Compute the tail index for Weibull Type III GEVD.

For Weibull Type III (ξ < 0), there is no power law tail behavior. Instead, the tail is exponentially bounded. The concept of tail index (1/ξ) is not meaningful for bounded distributions.

Returns:

Type	Description
ndarray	NaN (tail index not defined for bounded distributions)

Source code in src/xtremax/distributions/weibull.py

def tail_index(self) -> jnp.ndarray:
    """
    Compute the tail index for Weibull Type III GEVD.

    For Weibull Type III (ξ < 0), there is no power law tail behavior.
    Instead, the tail is exponentially bounded. The concept of tail
    index (1/ξ) is not meaningful for bounded distributions.

    Returns:
        NaN (tail index not defined for bounded distributions)
    """
    return jnp.full_like(self.concentration, jnp.nan)

exceedance_probability(threshold: jnp.ndarray) -> jnp.ndarray

Compute probability of exceeding a threshold: P(X > threshold).

For Weibull Type III, this is 0 when threshold ≥ upper bound.

Parameters:

Name	Type	Description	Default
threshold	ndarray	Threshold value	required

Returns:

Type	Description
ndarray	Exceedance probabilities (0 above upper bound)

Source code in src/xtremax/distributions/weibull.py

def exceedance_probability(self, threshold: jnp.ndarray) -> jnp.ndarray:
    """
    Compute probability of exceeding a threshold: P(X > threshold).

    For Weibull Type III, this is 0 when threshold ≥ upper bound.

    Args:
        threshold: Threshold value

    Returns:
        Exceedance probabilities (0 above upper bound)
    """
    return self.survival_function(threshold)

conditional_excess_mean(threshold: jnp.ndarray) -> jnp.ndarray

Mean excess :math:E[X - u \mid X > u] for Weibull Type III GEVD.

Computed via quantile-space quadrature, same as the general GEV class. The previous GPD linear approximation :math:(\sigma + \xi(u - \mu)) / (1 - \xi) is the POT asymptote and is wrong for finite thresholds. As :math:u approaches the upper endpoint :math:\mu - \sigma/\xi, the quadrature correctly decays to zero.

Source code in src/xtremax/distributions/weibull.py

def conditional_excess_mean(self, threshold: jnp.ndarray) -> jnp.ndarray:
    r"""Mean excess :math:`E[X - u \mid X > u]` for Weibull Type III GEVD.

    Computed via quantile-space quadrature, same as the general GEV
    class. The previous GPD linear approximation
    :math:`(\sigma + \xi(u - \mu)) / (1 - \xi)` is the POT asymptote
    and is wrong for finite thresholds. As :math:`u` approaches the
    upper endpoint :math:`\mu - \sigma/\xi`, the quadrature
    correctly decays to zero.
    """
    threshold_arr = jnp.asarray(threshold)

    # Log-space tail-probability quadrature (see GEVD implementation
    # for the derivation). Correctly normalises by S(u) rather than
    # the truncated mass `(1-ε) - F(u)`.
    n_grid = 1024
    p0 = self.cdf(threshold_arr)
    # Joint batch of threshold and parameters; per-batch quantities are
    # broadcast to it before gaining the trailing grid axis.
    batch = jnp.shape(p0)
    s_u = 1.0 - p0
    # Cap away from 1.0 so p_grid stays strictly in (0, 1). See
    # GEVD.conditional_excess_mean.
    s_u_safe = jnp.clip(s_u, 1e-12, 1.0 - 1e-6)
    log_s_u = jnp.log(s_u_safe)
    # Lower log-y endpoint: step ~20 e-folds below log_s_u so the
    # grid always captures the bulk of the tail-integrand mass AND
    # stays strictly ascending. A fixed floor like log(1e-6) would
    # run backward when S(u) < 1e-6 (near the upper endpoint) and
    # trapezoid would sign-flip. The tiny truncated remainder is
    # added back analytically below (negligible for ξ < 0 but kept
    # for structural parity with the heavy-tailed classes). The
    # ``-jnp.log1p(-y)`` step below is accurate to subnormals in
    # float32. See GEVD.conditional_excess_mean.
    shape = self.concentration
    span = jnp.asarray(20.0, dtype=log_s_u.dtype)
    log_y_min = log_s_u - span

    unit = jnp.linspace(0.0, 1.0, n_grid)
    log_s_u_exp = jnp.expand_dims(log_s_u, axis=-1)
    log_y_min_exp = jnp.expand_dims(jnp.broadcast_to(log_y_min, batch), axis=-1)
    v_grid = log_y_min_exp * (1.0 - unit) + log_s_u_exp * unit
    y_grid = jnp.exp(v_grid)
    # Compute ``x = F⁻¹(1 - y)`` directly from ``y`` using
    # ``-log1p(-y)``. For ``y`` near ``eps(float32)`` the
    # round-tripped ``p = 1 - y`` keeps only ~1 ulp of precision and
    # ``log(p)`` inside ``weibull_icdf`` loses that entirely. See
    # GEVD.conditional_excess_mean for details. ξ < 0 is enforced
    # by ``arg_constraints``, so the Gumbel branch is unreachable.
    # Parameters are broadcast to the joint batch and given the same
    # trailing grid axis — batched (batch,) parameters against a
    # (..., n_grid) grid otherwise raise a broadcasting error.
    neg_log_q = -jnp.log1p(-y_grid)
    loc_e = jnp.expand_dims(jnp.broadcast_to(self.loc, batch), axis=-1)
    scale_e = jnp.expand_dims(jnp.broadcast_to(self.scale, batch), axis=-1)
    shape_e = jnp.expand_dims(jnp.broadcast_to(shape, batch), axis=-1)
    x_grid = loc_e + (scale_e / shape_e) * (jnp.power(neg_log_q, -shape_e) - 1.0)

    integrand = x_grid * y_grid
    numerator = jnp.trapezoid(integrand, x=v_grid, axis=-1)

    # Closed-form remainder for the truncated (0, y_c] tail with
    # y_c = S(u)·e⁻²⁰, where x(y) ≈ (μ - σ/ξ) + (σ/ξ)·y^{-ξ}:
    #   ∫₀^{y_c} x dy = (μ - σ/ξ)·y_c + (σ/ξ)·y_c^{1-ξ}/(1-ξ).
    # Negligible for ξ < 0 (1-ξ > 1 makes it decay faster than y_c)
    # but kept for parity with the heavy-tailed classes.
    loc_b = jnp.broadcast_to(self.loc, batch)
    scale_b = jnp.broadcast_to(self.scale, batch)
    shape_b = jnp.broadcast_to(shape, batch)
    y_c = jnp.exp(log_y_min)
    remainder = (loc_b - scale_b / shape_b) * y_c + (scale_b / shape_b) * jnp.power(
        y_c, 1.0 - shape_b
    ) / (1.0 - shape_b)
    numerator = numerator + remainder

    mean_conditional = numerator / s_u_safe
    mean_excess = mean_conditional - threshold_arr

    return jnp.where(s_u > 1e-12, mean_excess, jnp.nan)

reliability_function(time: jnp.ndarray) -> jnp.ndarray

Compute the reliability function R(t) = P(X > t).

This is equivalent to the survival function but emphasizes the reliability engineering interpretation for Weibull Type III.

Parameters:

Name	Type	Description	Default
time	ndarray	Time points for reliability evaluation	required

Returns:

Type	Description
ndarray	Reliability values (probability of surviving past time t)

Source code in src/xtremax/distributions/weibull.py

def reliability_function(self, time: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the reliability function R(t) = P(X > t).

    This is equivalent to the survival function but emphasizes
    the reliability engineering interpretation for Weibull Type III.

    Args:
        time: Time points for reliability evaluation

    Returns:
        Reliability values (probability of surviving past time t)
    """
    return self.survival_function(time)

mean_residual_life(time: jnp.ndarray) -> jnp.ndarray

Compute the mean residual life function: E[X - t | X > t].

For reliability applications, this represents the expected additional lifetime given survival to time t.

Parameters:

Name	Type	Description	Default
time	ndarray	Current time/age	required

Returns:

Type	Description
ndarray	Mean residual life values

Source code in src/xtremax/distributions/weibull.py

def mean_residual_life(self, time: jnp.ndarray) -> jnp.ndarray:
    """
    Compute the mean residual life function: E[X - t | X > t].

    For reliability applications, this represents the expected
    additional lifetime given survival to time t.

    Args:
        time: Current time/age

    Returns:
        Mean residual life values
    """
    return self.conditional_excess_mean(time)

percentile_residual_life(time: jnp.ndarray, percentile: float) -> jnp.ndarray

Compute percentile residual life: (Q(p|X>t) - t).

This gives the additional time until the p-th percentile failure, conditional on survival to time t.

Parameters:

Name	Type	Description	Default
time	ndarray	Current time/age	required
percentile	float	Percentile level (0 < percentile < 1)	required

Returns:

Type	Description
ndarray	Percentile residual life values

Source code in src/xtremax/distributions/weibull.py

def percentile_residual_life(
    self, time: jnp.ndarray, percentile: float
) -> jnp.ndarray:
    """
    Compute percentile residual life: (Q(p|X>t) - t).

    This gives the additional time until the p-th percentile failure,
    conditional on survival to time t.

    Args:
        time: Current time/age
        percentile: Percentile level (0 < percentile < 1)

    Returns:
        Percentile residual life values
    """
    # Conditional CDF: F(x | X > t) = p ⟺ F(x) = 1 - (1 - p) S(t).
    # (Previously used `1 - p S(t)`, which mapped p=0 to F^{-1}(1),
    # i.e. the upper endpoint, instead of returning zero residual.)
    survival_prob = self.survival_function(time)
    total_prob = 1.0 - (1.0 - percentile) * survival_prob
    return self.icdf(total_prob) - time

expand(batch_shape: tuple[int, ...]) -> dist.Distribution

Expand to batch_shape by reconstructing via __init__.

Source code in src/xtremax/distributions/weibull.py

def expand(self, batch_shape: tuple[int, ...]) -> dist.Distribution:
    """Expand to ``batch_shape`` by reconstructing via ``__init__``."""
    batch_shape = tuple(batch_shape)
    if batch_shape == self.batch_shape:
        return self
    return type(self)(
        loc=jnp.broadcast_to(self.loc, batch_shape),
        scale=jnp.broadcast_to(self.scale, batch_shape),
        concentration=jnp.broadcast_to(self.concentration, batch_shape),
        validate_args=self._validate_args,
    )