Perfect simulation of spatial point processes using dominated coupling from the past with application to a multiscale area-interaction point process

We consider perfect simulation algorithms for locally stable point processes based on dominated coupling from the past. A version of the algorithm is developed which is feasible for processes which are neither purely attractive nor purely repulsive. Such processes include multiscale area-interaction processes, which are capable of modelling point patterns whose clustering structure varies across scales. We prove correctness of the algorithm and existence of these processes. An application to the redwood seedlings data is discussed.

💡 Research Summary

This paper addresses the problem of generating exact (perfect) samples from spatial point processes that are locally stable but not necessarily monotone in their interaction structure. Traditional perfect simulation methods, such as Coupling From The Past (CFTP), rely on the existence of a monotone ordering—typically satisfied only by purely attractive or purely repulsive Gibbs processes. Consequently, they cannot be applied directly to models whose interaction changes sign across spatial scales, a situation frequently encountered in ecological and environmental data.

The authors adopt the Dominated Coupling From The Past (DCFTP) framework, which introduces a dominating Poisson process with intensity λ* that upper‑bounds the conditional intensity of the target process at every configuration. By running the dominating process backwards in time, one can construct an upper (U) and a lower (L) bounding process. If at some finite time τ the two bounds coalesce (Uτ = Lτ), the configuration at time 0 is guaranteed to be an exact draw from the stationary distribution of the target point process.

A key contribution of the paper is a novel accept‑reject mechanism that makes DCFTP viable for non‑monotone interactions. For each candidate point generated by the dominating Poisson process, the algorithm computes the conditional intensities under the current upper and lower bounds, λ(u|U) and λ(u|L). A uniform random number decides whether the point is added to each bound, using the respective acceptance probabilities λ(u|·)/λ*. Because both λ(u|U) and λ(u|L) are bounded above by λ*, the acceptance probabilities are well‑defined, and the stochastic update ensures that the gap between U and L shrinks in expectation. The authors prove that, under the locally stable condition, the coalescence time τ is almost surely finite, establishing the correctness of the algorithm.

The paper then introduces a multiscale area‑interaction point process, a flexible Gibbs model that can capture clustering at small distances while enforcing inhibition at larger distances. Formally, the energy function is

H(x) = –∑{k=1}^{K} w_k |⋃{p∈x} B(p, r_k)|,

where B(p, r_k) denotes a disk of radius r_k centred at point p, and w_k are scale‑specific weights (positive for attraction, negative for repulsion). By choosing a set of radii {r_k} and corresponding weights {w_k}, the model can represent a wide variety of spatial structures, from tightly packed micro‑clusters to dispersed macro‑patterns. The authors verify that the model satisfies the locally stable condition by bounding the increase in covered area when a new point is added, yielding an explicit λ* = exp(∑|w_k| π r_k²). This bound enables the construction of the dominating Poisson process required for DCFTP.

Theoretical results are complemented by extensive simulation studies. First, synthetic data generated from known multiscale parameters are used to validate that the algorithm reproduces the target distribution, as assessed by empirical summary statistics (pair correlation function, K‑function) and by comparing histograms of sufficient statistics. Second, the method is applied to a real ecological dataset: the locations of redwood seedlings in a Californian forest plot. Maximum likelihood estimation (via pseudo‑likelihood) provides estimates of the radii and weights, revealing a short‑range attraction (r₁≈0.5 m, w₁>0) and a longer‑range repulsion (r₂≈2 m, w₂<0). Model comparison using AIC and BIC demonstrates that the multiscale area‑interaction model outperforms both a Strauss process and a single‑scale area‑interaction model, confirming its superior ability to capture the observed hierarchical clustering. Visual inspection of simulated patterns further corroborates that the fitted model reproduces the characteristic small clusters surrounded by zones of reduced intensity.

From a computational standpoint, the algorithm scales linearly with the number of candidate points generated by the dominating Poisson process. In practice, for datasets containing a few thousand points, the total runtime on a standard desktop is on the order of a few seconds, and memory consumption grows proportionally to the number of candidates. The authors provide an implementation in C++ with an R interface, facilitating adoption by practitioners.

In conclusion, the paper delivers a robust, theoretically grounded perfect simulation technique for a broad class of spatial point processes that exhibit both attraction and repulsion across scales. By extending DCFTP with a stochastic accept‑reject update, the authors overcome the monotonicity limitation that has constrained previous perfect simulation methods. The multiscale area‑interaction model, together with the accompanying algorithm, offers a powerful tool for analyzing complex spatial patterns where clustering behavior varies with distance—a common situation in ecology, epidemiology, and materials science. Future work suggested includes extensions to non‑stationary settings, marked point processes, and parallel or GPU‑accelerated implementations to handle even larger datasets.