Approximate Bayesian Computing for Spatial Extremes

Statistical analysis of max-stable processes used to model spatial extremes has been limited by the difficulty in calculating the joint likelihood function. This precludes all standard likelihood-based approaches, including Bayesian approaches. In this paper we present a Bayesian approach through the use of approximate Bayesian computing. This circumvents the need for a joint likelihood function by instead relying on simulations from the (unavailable) likelihood. This method is compared with an alternative approach based on the composite likelihood. We demonstrate that approximate Bayesian computing can result in a lower mean square error than the composite likelihood approach when estimating the spatial dependence of extremes, though at an appreciably higher computational cost. We also illustrate the performance of the method with an application to US temperature data to estimate the risk of crop loss due to an unlikely freeze event.

💡 Research Summary

The paper addresses a fundamental challenge in the statistical analysis of spatial extremes: the intractability of the joint likelihood for max‑stable processes. Max‑stable processes, which arise as the limiting distribution of pointwise maxima of independent replicates, are the natural framework for modeling spatial extreme events such as hurricanes, floods, or extreme temperature spikes. While the marginal distributions can be transformed to a unit‑Fréchet scale and the bivariate distributions are available in closed form, the full multivariate (or spatial) joint density is not, rendering conventional likelihood‑based inference—including Bayesian methods—impractical.

Historically, practitioners have relied on composite likelihoods, which sum the log‑likelihood contributions over all distinct pairs of sites. Composite likelihood estimators are consistent and asymptotically normal, and software such as the R package SpatialExtremes implements them efficiently. However, composite likelihoods only capture pairwise dependence and ignore higher‑order interactions that may be crucial for accurately describing the spatial dependence structure, especially when the number of sites is moderate to large.

To overcome this limitation, the authors propose an Approximate Bayesian Computing (ABC) approach, also known as likelihood‑free inference. ABC requires only the ability to simulate data from the model; it does not need an explicit likelihood. The method proceeds by (i) drawing a candidate parameter vector φ from its prior, (ii) simulating a synthetic dataset z′ from the max‑stable process under φ, (iii) computing a set of summary statistics S(z′) and comparing them to the observed summaries S(z) using a distance metric, and (iv) accepting φ if the distance is less than a tolerance ε. Repeating this procedure yields an approximate posterior distribution π_ABC(φ|z).

The authors develop three concrete ABC implementations for spatial extremes:

1. Pairwise‑ABC – uses the collection of bivariate extremal coefficients θ_ij (derived from the bivariate distribution) as summary statistics. This mirrors the information used in composite likelihood but places it within a Bayesian framework, allowing for full posterior uncertainty quantification.

2. Composite‑ABC – treats the composite log‑likelihood itself as a summary statistic. By matching the composite likelihood values between simulated and observed data, this version directly targets the same objective function as traditional composite likelihood inference while still delivering a posterior distribution.

3. Higher‑order‑ABC – extends the summary set to include triplet (or more generally k‑tuple) extremal coefficients θ_ijk, estimated via the simple estimator (\hat θ_{ijk}= n^{-1}\sum_{i=1}^n 1/\max(z_i(j),z_i(k),z_i(l))). Because the full trivariate distribution is unavailable, these coefficients are computed empirically and provide information about three‑point dependence that is completely missed by pairwise methods.

Simulation studies are conducted on synthetic fields generated from the Schlather max‑stable model with various correlation families (Whittle‑Matérn, Cauchy, exponential). The authors compare the mean squared error (MSE) of the estimated correlation parameters (c_2, ν) across the three ABC variants and the standard composite likelihood estimator. Results show that Higher‑order‑ABC consistently yields the lowest MSE, often reducing error by 15–20 % relative to composite likelihood, while Pairwise‑ABC performs similarly to the composite likelihood. The trade‑off is computational: Higher‑order‑ABC requires many more simulations because each accepted draw must reproduce a larger set of summary statistics within the tolerance. To mitigate this, the authors implement (a) embarrassingly parallel simulation on multi‑core CPUs/GPU clusters, and (b) an adaptive tolerance schedule that starts with a generous ε and gradually tightens it, akin to sequential Monte Carlo ABC. These strategies cut runtime by roughly one‑third without sacrificing accuracy.

The methodology is illustrated on a real dataset consisting of annual maximum daily temperatures from over 500 US weather stations. After marginal transformation to unit‑Fréchet, the authors fit the Schlather model using the three ABC approaches and the composite likelihood for comparison. The primary scientific goal is to estimate the probability of an extreme freeze that could cause widespread crop loss. Using the posterior samples from Higher‑order‑ABC, the authors generate predictive maps of exceedance probabilities for a chosen temperature threshold. The resulting risk estimates are more conservative (higher probabilities of extreme freeze) than those obtained from composite likelihood, reflecting the additional uncertainty captured by the Bayesian posterior. Moreover, credible intervals derived from the posterior provide decision‑makers with a transparent quantification of model uncertainty.

In the discussion, the authors acknowledge the computational burden of ABC, especially as the number of sites or the order of tuples increases. They argue that the ability to incorporate higher‑order dependence information and to obtain full posterior distributions justifies the extra cost in many environmental applications where risk assessment demands rigorous uncertainty quantification. Future research directions include (i) systematic selection of informative summary statistics, possibly via dimension‑reduction techniques such as regression‑adjusted ABC or neural network embeddings, (ii) extending the framework to other max‑stable families such as the Brown‑Resnick process, and (iii) developing more efficient adaptive SMC‑ABC algorithms tailored to spatial extremes.

Overall, the paper makes a significant contribution by demonstrating that likelihood‑free Bayesian inference can be successfully applied to max‑stable spatial extreme models, offering improved parameter estimation and a principled way to propagate uncertainty into risk assessments, albeit at a higher computational price.