📊 statistics

Nonparametric inference for nonstationary spatial point processes

This paper proposes a discretization-free MCMC algorithm for a spatial Cox process model that utilizes a random Voronoi partition to flexibly capture nonstationary features like abrupt intensity changes and hotspots in point pattern data while ensuring exact inference and computational scalability.

Original authors: Izabel Nolau, Flávio B. Gonçalves, Dani Gamerman

Published 2026-06-25

📖 4 min read☕ Coffee break read

CC0 1.0

Original authors: Izabel Nolau, Flávio B. Gonçalves, Dani Gamerman

Original paper dedicated to the public domain under CC0 1.0 (http://creativecommons.org/publicdomain/zero/1.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are looking at a map of a forest, and you want to understand where fires are most likely to start. You see some spots with huge clusters of fire (hotspots) and other areas that are completely empty. You also notice that the rules seem to change abruptly: one side of the river might be a fire-prone dry grassland, while the other side is a wet, fire-resistant swamp.

Traditional statistical tools for mapping these patterns are like a blurry camera. They try to smooth out the data to make a pretty picture. If there's a sharp line between a fire zone and a safe zone, a blurry camera just paints a muddy gray transition. It misses the sudden changes and the distinct "hot" and "cold" spots.

This paper introduces a new, sharper tool called the NSPP model (Nonstationary Spatial Point Process). Here is how it works, using simple analogies:

1. The "Random Jigsaw Puzzle" Approach

Instead of trying to force the whole map into one smooth, continuous picture, this new method cuts the map into a random jigsaw puzzle.

The Partition: Imagine dropping a handful of pins on a map. The area around each pin becomes its own unique "neighborhood" (a Voronoi cell).
The Rules: Inside each neighborhood, the fire rules are smooth and consistent (like a calm lake). But the moment you cross the border into the next neighborhood, the rules can change completely. One neighborhood might be a fire zone; the next might be a safe zone.
The Magic: The computer doesn't just guess where the puzzle pieces go; it figures out the best puzzle shape while it figures out the fire rules. This allows it to capture sharp, jagged edges where the fire risk suddenly spikes or drops, which older models miss.

2. The "Ghost Points" Trick (Data Augmentation)

Calculating the exact probability of these fire patterns is usually a mathematical nightmare because it involves infinite possibilities. It's like trying to count every single grain of sand on a beach to predict a storm.

The Solution: The authors invented a clever trick. They imagine "ghost points" (invisible fire starts) that didn't actually happen but could have.
Why it helps: By adding these ghost points to the math, the impossible calculation suddenly becomes solvable. It's like turning a complex riddle into a simple puzzle. This allows them to get the exact answer without having to chop the map into a tiny, rigid grid (which usually introduces errors).

3. Why It's Better Than the Old Ways

No Blurring: Old models smooth out the data, making a sharp cliff look like a gentle hill. This new model keeps the cliff sharp.
No Grid Errors: Many methods force the map into a grid (like a chessboard). If a fire hotspot sits right between two squares, the model gets confused. This new method flows naturally over the land without being trapped in a grid.
Robustness: Even if you tell the computer to cut the map into 20 pieces when there are really only 2 distinct zones, the model is smart enough to realize the extra pieces are unnecessary. It doesn't break; it just ignores the extra noise.

Real-World Tests

The authors tested this on two things:

Fake Data: They created computer simulations with known sharp boundaries and hotspots. The new model found them perfectly, while the old "blurry" models failed to see the sharp lines.
Real Data:
- Trees in Panama: They mapped where a specific tree species grows. The model correctly identified that the trees love steep slopes and avoid flat areas, capturing the sharp drop-off in density where the terrain changes.
- Fires in Brazil: They mapped fires in Mato Grosso. The model successfully identified specific high-risk zones in the north and low-risk zones in the south, drawing a clear line between them that matched real-world land use changes (like deforestation).

The Bottom Line

This paper presents a new way to map random events (like fires, trees, or disease cases) that respects the reality that nature often changes abruptly. It uses a "random puzzle" strategy to find sharp edges and a "ghost point" trick to do the math exactly. The result is a map that is sharper, more accurate, and better at showing where the real danger zones are, without the blurring or grid errors of older methods.

Technical Summary: Nonparametric Inference for Nonstationary Spatial Point Processes

Problem Statement
Spatial point pattern data frequently exhibit complex features such as abrupt changes, hotspots, and spatially varying dependence in local intensity. Under a Poisson process framework, these phenomena correspond to discontinuities and nonstationarity in the underlying intensity function (IF). Standard modeling approaches, including stationary models and many existing nonstationary formulations, struggle to capture these features. Stationary models fail to represent heterogeneous regimes, while standard nonstationary models (e.g., those driven by a single latent Gaussian process) tend to oversmooth, blurring the sharp transitions that characterize the data.

The core difficulty lies in estimating flexible nonstationary covariance structures when the driving Gaussian process (GP) is latent and enters the model nonlinearly through the IF of a Poisson process. Furthermore, the likelihood of a GP-driven Cox process is intractable, requiring the integration of an infinite-dimensional GP against a Poisson likelihood with no closed form. Existing literature typically circumvents this via approximations, such as discretizing the domain into grids or using Laplace approximations. These methods introduce approximation errors that are difficult to quantify, often scale cubically with the number of grid cells, and can distort structural features like sharp intensity boundaries.

Methodology
The paper proposes a Nonstationary Spatial Point Process (NSPP) model, a Cox process where nonstationarity is induced through a random partition of the spatial domain. The key components of the methodology are:

Random Partition Framework: The spatial domain $S$ is partitioned into $L$ regions ( $S_1, \dots, S_L$ ) via a Voronoi tessellation generated by a set of unknown generator points $U$ . Within each region $S_\ell$ , the intensity is modeled as $\lambda^*_\ell F(\beta_\ell(s))$ , where:
- $\lambda^*_\ell$ is a region-specific upper bound for the intensity.
- $F$ is a known monotone link function (e.g., standard Gaussian CDF).
- $\beta_\ell$ is a conditionally independent, stationary Gaussian Process specific to region $S_\ell$ .
- This construction allows for heterogeneous local behavior and sharp discontinuities at region boundaries while retaining the tractable machinery of stationary GPs within regions.
Exact Inference via Data Augmentation: To address the intractable likelihood, the authors extend the Poisson thinning scheme to the random partition setting. They introduce latent processes $\tilde{Y}_\ell$ (thinned events) and $Z_\ell$ (events outside the region) alongside the observed process $Y_\ell$ . This yields an augmented likelihood that is analytically tractable, enabling exact Bayesian inference without discretizing the spatial domain or approximating the integral of the intensity function.
MCMC Algorithm: A discretization-free Markov Chain Monte Carlo (MCMC) algorithm is developed to target the exact infinite-dimensional posterior distribution. The algorithm employs a Gibbs sampling scheme with five blocks:
- Updates for latent processes $(\tilde{Y}, Z)$ .
- Joint updates for the partition $U$ and process labels $(Y, \tilde{Y}, Z)$ using a Metropolis-Hastings step that respects the repulsive prior on generator points.
- Updates for the latent GPs $\beta$ (using skew-normal sampling or Metropolis-Hastings).
- Updates for intensity bounds $\lambda^*$ .
- Updates for GP hyperparameters $\theta$ and the infinite-dimensional remainder of $\beta$ via retrospective sampling.
The use of retrospective sampling ensures that the algorithm targets the exact posterior without finite-dimensional approximation. Additionally, the random partitioning localizes the cubic computational cost of GP operations to subdomains, and the framework can be scaled to large datasets using sparse GP constructions (e.g., Vecchia-type approximations) while preserving exactness.

Key Contributions

Exact Inference for Nonstationary Cox Processes: The paper presents the first GP-driven Cox process that combines flexible nonstationary intensity estimation with exact, discretization-error-free Bayesian inference.
Handling Sharp Discontinuities: By utilizing a random Voronoi partition with conditionally independent GPs, the model captures abrupt changes and sharp transitions in intensity that are typically smoothed out by standard nonstationary models.
Computational Efficiency: The partitioning strategy reduces the computational burden of GP inference by confining calculations to subdomains. The method scales to large datasets (e.g., thousands of points) while maintaining exact inference properties.
Robustness to Model Specification: The methodology demonstrates robustness to the choice of the number of partitions $L$ , performing well even when $L$ is larger than the true number of underlying regimes.

Results
The methodology was evaluated through synthetic examples and real-world applications:

Synthetic Data: In scenarios with discontinuous intensity functions and localized hotspots, the NSPP model outperformed stationary models, kernel-based estimators, and other nonstationary alternatives. It accurately recovered intensity surfaces and sharp boundaries, whereas competing models either oversmoothed the data or failed to capture the discontinuities. The model remained robust even when the number of partitions $L$ was misspecified.
Real Data (Beilschmiedia pendula): Applied to tree locations in Panama, the NSPP captured high- and low-density regions more accurately than a Log-Gaussian Cox Process (LGCP) and an XGBoost-based nonstationary model. The LGCP produced an overly smooth surface, while XGBoostPP generated unreliable local estimates. The NSPP naturally identified low-intensity zones without prior specification of such classes.
Real Data (Mato Grosso Fires): Applied to fire occurrences in Brazil, the model successfully identified high-risk regions (hotspots) in the central-northern and northwestern parts of the state, aligning with known land-use dynamics and deforestation patterns. It provided spatially resolved intensity estimates that standard stationary models could not achieve.

Significance and Claims
The paper claims that the proposed NSPP model offers a substantial methodological advance for point process analysis under nonstationarity. It addresses the long-standing trade-off between flexibility, exactness, and robust uncertainty quantification. Unlike existing approaches that sacrifice one of these properties, the NSPP model:

Provides coherent and robust uncertainty quantification for all structural features (intensities, correlations, and partitions) within a fully Bayesian framework.
Avoids discretization errors by operating directly on the continuous domain.
Effectively captures sharp spatial changes and heterogeneous regimes that are characteristic of real-world phenomena like fire outbreaks and ecological distributions.

The authors conclude that the method is competitive with, and often superior to, stationary and nonstationary alternatives in various scenarios, offering a principled framework for analyzing complex spatial point patterns. Future work is suggested to include non-spatial covariates and further computational improvements.