Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Layer 0 — Primitives

Pure JAX functions. Stateless, differentiable, composable. No classes, no NumPyro dependency, no xarray — just (Array, ...) → Array.

All functions are compatible with jax.jit, jax.grad, and jax.vmap. They implement the mathematical core of extreme value theory; the NumPyro Distribution wrappers at Layer 1 compose these into the standard probabilistic interface.

GEV Functions (xtremax._src.gev)

The Generalized Extreme Value distribution unifies three classical types under one parameterization.

CDF:

F(x;μ,σ,ξ)=exp⁡ ⁣{−[1+ξx−μσ]−1/ξ}F(x; \mu, \sigma, \xi) = \exp\!\left\{-\left[1 + \xi\frac{x - \mu}{\sigma}\right]^{-1/\xi}\right\}
(1)

with the convention ξ=0⇒F(x)=exp⁡{−exp⁡(−(x−μ)/σ)}\xi = 0 \Rightarrow F(x) = \exp\{-\exp(-(x-\mu)/\sigma)\} (Gumbel limit).

Support: {x:1+ξ(x−μ)/σ>0}\{x : 1 + \xi(x-\mu)/\sigma > 0\}.

gev_log_prob(x, loc, scale, shape)

Mathematical definition:

log⁡f(x;μ,σ,ξ)=−log⁡σ−(1+1/ξ)log⁡t−t−1/ξ\log f(x; \mu, \sigma, \xi) = -\log\sigma - (1 + 1/\xi)\log t - t^{-1/\xi}
(2)

where t=1+ξ(x−μ)/σt = 1 + \xi(x - \mu)/\sigma. For ξ=0\xi = 0: log⁡f=−log⁡σ−s−e−s\log f = -\log\sigma - s - e^{-s} where s=(x−μ)/σs = (x-\mu)/\sigma.

Numerical method: Uses jnp.where to branch between ξ≠0\xi \neq 0 and ξ≈0\xi \approx 0 (Gumbel) to avoid division by zero. Gradient-safe via JAX autodiff.

Complexity: O(N)O(N) element-wise.

def gev_log_prob(
    x: Float[Array, "..."],
    loc: Float[Array, "..."],
    scale: Float[Array, "..."],
    shape: Float[Array, "..."],
) -> Float[Array, "..."]:

gev_cdf(x, loc, scale, shape)

F(x)=exp⁡ ⁣{−t−1/ξ},t=1+ξ(x−μ)/σF(x) = \exp\!\left\{-t^{-1/\xi}\right\}, \quad t = 1 + \xi(x - \mu)/\sigma
(3)
def gev_cdf(x, loc, scale, shape) -> Array:

gev_icdf(q, loc, scale, shape)

Quantile function (inverse CDF):

F−1(q)=μ+σξ[(−log⁡q)−ξ−1]F^{-1}(q) = \mu + \frac{\sigma}{\xi}\left[(-\log q)^{-\xi} - 1\right]
(4)

For ξ=0\xi = 0: F−1(q)=μ−σlog⁡(−log⁡q)F^{-1}(q) = \mu - \sigma \log(-\log q).

def gev_icdf(q, loc, scale, shape) -> Array:

gev_return_level(period, loc, scale, shape)

TT-year return level:

zT=F−1(1−1/T)=μ+σξ[(−log⁡(1−1/T))−ξ−1]z_T = F^{-1}(1 - 1/T) = \mu + \frac{\sigma}{\xi}\left[(-\log(1 - 1/T))^{-\xi} - 1\right]
(5)

This is the level expected to be exceeded on average once every TT years.

def gev_return_level(period, loc, scale, shape) -> Array:

gev_mean(loc, scale, shape)

E[X]=μ+σΓ(1−ξ)−1ξ,ξ<1\mathbb{E}[X] = \mu + \sigma \frac{\Gamma(1 - \xi) - 1}{\xi}, \quad \xi < 1
(6)

Undefined for ξ≥1\xi \geq 1 (infinite mean — heavy tail).

gev_variance(loc, scale, shape)

Var⁡(X)=σ2ξ2[Γ(1−2ξ)−Γ2(1−ξ)],ξ<1/2\operatorname{Var}(X) = \frac{\sigma^2}{\xi^2}\left[\Gamma(1 - 2\xi) - \Gamma^2(1 - \xi)\right], \quad \xi < 1/2
(7)

GPD Functions (xtremax._src.gpd)

The Generalized Pareto Distribution models threshold exceedances Y=X−u∣X>uY = X - u \mid X > u.

CDF:

G(y;σ,ξ)=1−(1+ξyσ)−1/ξ,y>0G(y; \sigma, \xi) = 1 - \left(1 + \xi \frac{y}{\sigma}\right)^{-1/\xi}, \quad y > 0
(8)

For ξ=0\xi = 0: G(y)=1−exp⁡(−y/σ)G(y) = 1 - \exp(-y/\sigma) (exponential).

gpd_log_prob(y, scale, shape)

log⁡g(y;σ,ξ)=−log⁡σ−(1+1/ξ)log⁡(1+ξy/σ)\log g(y; \sigma, \xi) = -\log\sigma - (1 + 1/\xi)\log\left(1 + \xi y/\sigma\right)
(9)
def gpd_log_prob(y, scale, shape) -> Array:

gpd_cdf(y, scale, shape)

gpd_survival(y, scale, shape)

Survival function: Gˉ(y)=(1+ξy/σ)−1/ξ\bar{G}(y) = (1 + \xi y/\sigma)^{-1/\xi}

gpd_return_level(period, scale, shape, rate)

TT-year return level for POT:

zT=u+σξ[(λT)ξ−1]z_T = u + \frac{\sigma}{\xi}\left[(\lambda T)^{\xi} - 1\right]
(10)

where uu is the threshold and λ\lambda is the exceedance rate (expected exceedances per year).

def gpd_return_level(period, scale, shape, rate) -> Array:

gpd_mean_excess(scale, shape)

Mean excess function: E[Y−u∣Y>u]=σ+ξu1−ξ\mathbb{E}[Y - u \mid Y > u] = \frac{\sigma + \xi u}{1 - \xi} for ξ<1\xi < 1.

Point Process Functions (xtremax._src.point_process)

poisson_log_likelihood(event_times, rate, T)

Homogeneous Poisson process log-likelihood:

log⁡L=∑i=1nlog⁡λ−λT\log L = \sum_{i=1}^n \log \lambda - \lambda T
(11)

where λ\lambda is the rate and TT is the observation window.

inhomogeneous_poisson_log_likelihood(event_times, log_intensity_fn, T)

Inhomogeneous Poisson process:

log⁡L=∑i=1nlog⁡λ(ti)−∫0Tλ(t) dt\log L = \sum_{i=1}^n \log \lambda(t_i) - \int_0^T \lambda(t)\,dt
(12)

The integral is approximated numerically (trapezoidal or quadrature).

hawkes_log_likelihood(event_times, background_rate, alpha, beta)

Hawkes self-exciting process:

λ(t)=μ+α∑ti<tβe−β(t−ti)\lambda(t) = \mu + \alpha \sum_{t_i < t} \beta e^{-\beta(t - t_i)}
(13)
log⁡L=∑i=1nlog⁡λ(ti)−μT−αβ∑i=1n(1−e−β(T−ti))\log L = \sum_{i=1}^n \log \lambda(t_i) - \mu T - \frac{\alpha}{\beta}\sum_{i=1}^n \left(1 - e^{-\beta(T - t_i)}\right)
(14)

extremal_index(exceedances, threshold, run_length)

Extremal index θ∈(0,1]\theta \in (0, 1] — measures clustering of threshold exceedances. θ=1\theta = 1 means no clustering (independent exceedances). Estimated via the runs method:

θ^=number of clustersnumber of exceedances\hat{\theta} = \frac{\text{number of clusters}}{\text{number of exceedances}}
(15)

Max-Stable Functions (xtremax._src.max_stable)

brown_resnick_extremal_coeff(h, variogram_fn)

Brown-Resnick extremal coefficient:

θ(h)=2Φ ⁣(γ(h)/2)\theta(h) = 2\Phi\!\left(\sqrt{\gamma(h)/2}\right)
(16)

where γ(h)\gamma(h) is the semivariogram and Φ\Phi is the standard normal CDF.

smith_extremal_coeff(h, cov_matrix)

θ(h)=2Φ ⁣(hTΣ−1h/2)\theta(h) = 2\Phi\!\left(\sqrt{h^T \Sigma^{-1} h / 2}\right)
(17)

schlather_extremal_coeff(h, correlation_fn)

θ(h)=1+(1−ρ(h))/2\theta(h) = 1 + \sqrt{(1 - \rho(h))/2}
(18)

extremal_t_extremal_coeff(h, correlation_fn, df)

θ(h)=2 Tν+1 ⁣((ν+1)(1−ρ(h))/(1+ρ(h)))\theta(h) = 2\,T_{\nu+1}\!\left(\sqrt{(\nu+1)(1 - \rho(h))/(1 + \rho(h))}\right)
(19)

pairwise_log_likelihood(z, extremal_coeff_fn, coordinates)

Composite pairwise log-likelihood for max-stable processes (Padoan et al. 2010):

ℓc=∑i<jlog⁡fij(zi,zj)\ell_c = \sum_{i < j} \log f_{ij}(z_i, z_j)
(20)

where fijf_{ij} is the bivariate density derived from the bivariate CDF:

Fij(zi,zj)=exp⁡ ⁣(−1ziA ⁣(log⁡zj/zilog⁡zj/zi+log⁡zi/zj)−1zjA ⁣(log⁡zi/zjlog⁡zj/zi+log⁡zi/zj))F_{ij}(z_i, z_j) = \exp\!\left(-\frac{1}{z_i}A\!\left(\frac{\log z_j / z_i}{\log z_j / z_i + \log z_i / z_j}\right) - \frac{1}{z_j}A\!\left(\frac{\log z_i / z_j}{\log z_j / z_i + \log z_i / z_j}\right)\right)
(21)

and AA is the Pickands dependence function determined by θ\theta.

madogram(data, coordinates, n_bins)

F-madogram for nonparametric extremal coefficient estimation:

ν^(h)=12mean⁡∣Fn(Z(s1))−Fn(Z(s2))∣,∥s1−s2∥≈h\hat{\nu}(h) = \frac{1}{2}\operatorname{mean}\left|F_n(Z(s_1)) - F_n(Z(s_2))\right|, \quad \|s_1 - s_2\| \approx h
(22)
θ^(h)=1+2ν^(h)1−2ν^(h)\hat{\theta}(h) = \frac{1 + 2\hat{\nu}(h)}{1 - 2\hat{\nu}(h)}
(23)

Variogram Functions (xtremax._src.variogram)

Used by max-stable models and spatial GP layers.

power_variogram(h, range, smoothness)

γ(h)=(∥h∥/λ)α,α∈(0,2]\gamma(h) = \left(\|h\| / \lambda\right)^\alpha, \quad \alpha \in (0, 2]
(24)

matern_variogram(h, range, smoothness, sill)

γ(h)=σ2(1−21−νΓ(ν)(∥h∥λ)νKν ⁣(∥h∥λ))\gamma(h) = \sigma^2 \left(1 - \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\|h\|}{\lambda}\right)^\nu K_\nu\!\left(\frac{\|h\|}{\lambda}\right)\right)
(25)

where KνK_\nu is the modified Bessel function of the second kind.

Design
API Overview
Design
Components