mixle.analysis.max_stable module

Max-stable processes for spatial extremes: the Smith (Gaussian-storm) model.

Block maxima of a spatial field (annual flood peaks, peak seismic amplitude, extreme porosity) are spatially dependent, and that dependence has its own limit law – a max-stable process – which the ordinary GEV/GPD (treated independently per site) misses. The Smith model is the canonical one: Z(s) = max_i xi_i * phi_Sigma(s - U_i) over a Poisson storm process, giving unit-Frechet margins and a closed-form pairwise dependence. The extremal coefficient theta(h) in [1, 2] summarizes it: 1 = full dependence (extremes always co-occur), 2 = independence.

extremal_coefficient and bivariate_cdf are closed-form and exact. The storm process itself is, in principle, an infinite Poisson process over an unbounded domain; SmithMaxStableSampler.sample() necessarily truncates both the storm count and the domain, so simulated margins are only approximately unit Frechet – see that method’s docstring for exactly what is and is not rigorously controlled, and SmithMaxStableSampler.approximation_diagnostic() to check it empirically.

class SmithMaxStable(sigma)[source]

Bases: object

The Smith max-stable process with Gaussian storm-profile covariance sigma (d x d, SPD).

A spatial process, not an i.i.d. leaf distribution: its full likelihood is intractable, so it is not a SequenceEncodableProbabilityDistribution – it exposes the things that do have closed forms (extremal_coefficient, bivariate_cdf) plus a sampler, and is fitted by the module-level fit_smith_maxstable() (composite/madogram estimation), mirroring the functional fit style of the other non-leaf spatial models. Margins are unit Frechet; spatial dependence grows with sigma.

Parameters:

sigma (np.ndarray)

property sigma: ndarray

Owned read-only covariance used consistently by dependence formulas and simulation.

extremal_coefficient(h)[source]

theta(h) = 2 * Phi(a/2) with a the Mahalanobis lag length – 1 at h=0 (full dependence) rising to 2 as the lag grows (independence).

Parameters:

h (ndarray)

Return type:

float

bivariate_cdf(z1, z2, h)[source]

P(Z(s) <= z1, Z(s+h) <= z2) = exp(-V(z1, z2)) – the Smith bivariate distribution.

z1/z2 are thresholds on unit-Frechet margins and so must be strictly positive (the unit-Frechet support is (0, inf)); the closed form below divides by z1/z2 and takes log(z2/z1), which silently produces nan/inf rather than a real probability at or below 0.

Parameters:
Return type:

float

sampler(locations, seed=None)[source]

Return a sampler over the requested spatial locations.

Parameters:
Return type:

SmithMaxStableSampler

class SmithMaxStableFit(model, n_locations, n_replicates, n_pairs, residual, status='ok')[source]

Bases: object

A fitted isotropic Smith max-stable process, plus a diagnostic of whether the pairwise dependence-scale match actually identified an interior optimum.

status is one of:

  • "ok": the fitted scale is an interior point of the search bracket – a genuine stationary point of the F-madogram least-squares objective.

  • "boundary": the fit landed at (or numerically indistinguishable from) the edge of the search bracket – the signature of an objective with no interior minimum. This happens, for example, when replicates carry almost no independent information (every location’s rank collapses together), which pushes the empirical extremal coefficient to the same value for every pair and drives the least-squares match monotonically toward one edge instead of a genuine optimum. model is still populated with the boundary value, but it is a search artifact, not a dependence-scale estimate, and should not be trusted.

Parameters:

  • model (SmithMaxStable)

  • n_locations (int)

  • n_replicates (int)

  • n_pairs (int)

  • residual (float)

  • status (str)

model

the fitted SmithMaxStable.

Type:

mixle.analysis.max_stable.SmithMaxStable

n_locations

distinct spatial locations used.

Type:

int

n_replicates

replicated fields used.

Type:

int

n_pairs

location pairs the least-squares match was computed over (n_locations choose 2).

Type:

int

residual

the F-madogram objective’s value at the fitted scale (mean squared error between the empirical and model extremal coefficients across pairs; exactly 0 for a 2-location fit, which has one pair and one parameter).

Type:

float

status

see above.

Type:

str

Construction validates the record against its own documented domain (MXR-080-1900). It previously validated nothing at all, so SmithMaxStableFit(model=None, n_locations=-1, n_replicates=0, n_pairs=-5, residual=nan) constructed cleanly and reported converged is True – a record claiming a converged fit that had no model, no data and no objective value. The invariants are exactly what fit_smith_maxstable() produces: a real SmithMaxStable; at least two locations (a dependence scale needs a pair) and one replicate; n_pairs the n_locations choose 2 the docstring above states; a finite, non-negative residual (it is a mean of squares); and a status this class actually defines.

Deliberately NOT changed: the record stays a MUTABLE dataclass. Its frozen siblings in this package were converted under a separate finding, and freezing this one is a broader API change than the shape invariants this validates – the fields can still be rebound after construction and escape the checks below.

property converged: bool

True only when the fit is a genuine interior optimum, not a search-bound artifact.

class SmithMaxStableSampler(dist, locations, seed=None)[source]

Bases: object

Sampler for a fitted Smith max-stable process at fixed locations.

sample() truncates the (in principle unbounded, infinite-storm) Smith spectral representation to a finite spatial box and a tolerance-controlled number of storms; see its docstring for exactly what that truncation does and does not rigorously control, and approximation_diagnostic() for an empirical check of the result.

Parameters:

  • dist (SmithMaxStable)

  • locations (np.ndarray)

  • seed (int | None)

sample(size=None, *, n_storms=20_000, tol=1e-3, box_sigma=5.0)[source]

Draw max-stable field(s) at the locations via a boxed, tolerance-truncated storm sum.

The true Smith representation is a max over an infinite Poisson process of storms (xi_i, U_i) spread over the whole of R^d; only then are margins exactly unit Frechet. This sampler approximates that in two distinct ways, one rigorously controlled and one not:

  1. Storm-count truncation (rigorously controlled by ``tol``). Storms are drawn in decreasing order of their Poisson intensity mark xi (xi_i = 1/Gamma_i, the natural arrival order of this construction). Because the Gaussian storm kernel is bounded by its peak density phi_Sigma(0), a storm with mark xi can contribute at most xi * vol * phi_Sigma(0) to any site, no matter where its center U_i lands (vol the sampling box’s volume). Since xi is strictly decreasing, once that bound drops below tol for one storm it has also dropped below tol for every later one – so it is exact, not heuristic, to stop there: every dropped storm’s true contribution is provably below tol at every site, which bounds the total storm-count truncation error (relative to the boxed field’s own true value) by tol. n_storms is a safety cap on iterations, not the primary control; a UserWarning is raised if that cap binds before tol is reached for some replicate, since then the tol guarantee was not actually achieved. Every storm drawn is used before its own tol check is applied, so at least one storm always contributes to the field regardless of how loose tol is – an all-zero field (like the previous n_storms=0 behavior) is just as invalid a “draw” as a negative one, since 0 is outside unit Frechet’s (0, inf) support.

  2. Spatial truncation (NOT rigorously controlled). Storm centers are drawn only from a finite box around the query locations, expanded by box_sigma storm-kernel standard deviations in each dimension, rather than the true unbounded domain. This introduces a bias that empirically shrinks as box_sigma grows, but is not exactly bounded here – use approximation_diagnostic() to check it for your particular locations/sigma.

Net effect: simulated margins are only approximately, not exactly, unit Frechet. The storm-count half of that approximation is controlled to within tol (when n_storms doesn’t bind first); the box half is not formally bounded and should be checked empirically.

Parameters:

  • size (int | None) – number of independent replicate fields (None for a single field).

  • n_storms (int) – safety cap on storms drawn per replicate; must be a positive integer. Storms stop earlier, as soon as tol is satisfied. Must be a positive integer – with n_storms=0 no storms are drawn at all, which cannot produce a valid unit-Frechet draw (0 is outside its (0, inf) support), so that case is rejected outright rather than silently returned as if it were one.

  • tol (float) – stop once every remaining storm’s max-possible contribution to every site is below this (absolute, field-scale) tolerance. Must be finite and strictly positive.

  • box_sigma (float) – half-width, in storm-kernel standard deviations, of the box storm centers are drawn from beyond the query locations’ own extent. Must be finite and strictly positive.

Returns:

Field(s) at the locations. Shape (len(locations),) if size is None, else (size, len(locations)).

Return type:

ndarray

size (when supplied) and n_storms are exact non-Boolean positive integers; zero does not represent a valid request for a field or a storm process.

approximation_diagnostic(n=2000, probability_levels=(0.25, 0.5, 0.75, 0.9), **sample_kwargs)[source]

Empirically estimate how close this sampler’s truncated simulation is to true unit-Frechet margins.

Draws n replicate fields (with sample_kwargs forwarded to sample(), so a specific n_storms/tol/box_sigma setting can be checked) and compares their empirical CDF, at each site and each requested probability level, against the true unit-Frechet CDF – see FrechetApproximationDiagnostic. n is an exact non-Boolean positive integer and probability_levels is a nonempty finite one-dimensional sequence strictly inside (0, 1).

Parameters:
Return type:

FrechetApproximationDiagnostic

class FrechetApproximationDiagnostic(n_replicates, probability_levels, max_abs_error, mean_abs_error, per_site_max_abs_error)[source]

Bases: object

Empirical check of how close a truncated SmithMaxStableSampler draw is to true unit Frechet margins.

At each site and each requested probability level p, compares the empirical CDF of simulated draws (evaluated at the true unit-Frechet quantile for p) against p itself – a QQ/PP-plot-style discrepancy. Both summary statistics should shrink toward 0 as the sampler’s truncation is loosened (smaller tol, larger box_sigma/n_storms) and should NOT shrink just from raising n alone (n only reduces Monte Carlo noise in the estimate of a possibly nonzero bias).

Parameters:
n_replicates

number of simulated draws the diagnostic was computed from.

Type:

int

probability_levels

the probability levels checked.

Type:

tuple[float, …]

max_abs_error

worst-case |empirical_cdf(q_p) - p| over every (site, level) pair.

Type:

float

mean_abs_error

mean |empirical_cdf(q_p) - p| over every (site, level) pair.

Type:

float

per_site_max_abs_error

worst-case error at each site, shape (n_locations,).

Type:

numpy.ndarray

fit_smith_maxstable(locations, fields)[source]

Fit an isotropic Smith max-stable process (sigma = s^2 I) to replicated spatial extremes.

locations is (n_locations, d) and fields is (n_replicates, n_locations) of block maxima. Estimation matches the binned empirical extremal coefficient (from the F-madogram) to the model 2 Phi(|h| / (2 s)).

This is a single-parameter (s) fit estimated from pairwise statistics, so it is only identifiable given enough structure to compute those statistics from:

  • at least 2 locations, so there is at least one pair (with at least one nonzero lag between some pair) to compute an empirical extremal coefficient from at all;

  • at least 2 replicates, so the per-location rank transform is not degenerate – with exactly 1 replicate every location’s rank trivially collapses to a single point, making the empirical extremal coefficient identically 1 (full dependence) for every pair regardless of distance, which is not a real dependence signal.

Returns a SmithMaxStableFit (the fitted model plus fit-quality diagnostics); raises ValueError if the inputs are malformed or too small to attempt a fit at all.

Parameters:
Return type:

SmithMaxStableFit