Approximation intensity for pairwise interaction Gibbs point processes using determinantal point processes

The intensity of a Gibbs point process is usually an intractable function of the model parameters. For repulsive pairwise interaction point processes, this intensity can be expressed as the Laplace transform of some particular function. Baddeley and Nair (2002) developped the Poisson-saddlepoint approximation which consists, for basic models, in calculating this Laplace transform with respect to a homogeneous Poisson point process. In this paper, we develop an approximation which consists in calculating the same Laplace transform with respect to a specific determinantal point process. This new approximation is efficiently implemented and turns out to be more accurate than the Poisson-saddlepoint approximation, as demonstrated by some numerical examples.


Introduction
Due to their simple interpretation, Gibbs point processes and in particular pairwise interaction point processes play a central role in the analysis of spatial point patterns (see van Lieshout (2000); Møller and Waagepetersen (2004); Baddeley et al. (2015)).In a nutshell, such models (in the homogeneous case) are defined in a bounded domain by a density with respect to the unit rate Poisson point process which takes the form where x is a finite configuration of points, where β > 0 represents the activity parameter, |x| is the number of elements of x and where g : R d → R + is the pairwise interaction function.
However, many important theoretical properties of these models are in general intractable, like for instance the simplest one, the intensity λ ∈ R + , representing the mean number of points per unit volume.It is known (see e.g.Section 2.2) that Such an expectation is in general intractable.As clearly outlined by Baddeley and Nair (2012), this intractability constitutes a severe drawback.For example, simulating a Gibbs point process with a prescribed value of λ cannot be done beforehand even for simple models such as Strauss models.Baddeley and Nair (2012) suggest to evaluate the expectation with respect to a homogeneous Poisson point process with intensity λ.This results in the Poisson-saddlepoint approximation, denoted by λ ps , obtained as the solution of ) du (provided this integral is finite).The general idea of the present paper is to evaluate the same expectaction with respect to a determinantal point process (with intensity λ).Determinantal point processes (DPP), see e.g.Lavancier et al. (2015), are a class of repulsive models which is more tractable than Gibbs models.For example all moments are explicit.If g ≤ 1 and has a finite range R > 0, our approximation denoted by λ dpp is the solution of where (1 − g) 2 |B(0, R)| , imsart-ejs ver.2014/10/16 file: appIntDPP.texdate: November 13, 2018 Coeurjolly and Lavancier/Intensity approximation using DPP

3
|A| denotes the volume of some bounded domain A ⊂ R d , B(0, ρ) is the Euclidean ball centered at 0 with radius ρ and δ ≥ 0 is some possible hard-core distance.
Both approximations λ dpp and λ ps can be obtained very quickly with a unitroot search algorithm.Figure 1 reports λ dpp and λ ps as well as the true intensity λ (obtained by Monte-Carlo simulations) for Strauss models in terms of the interaction parameter γ 1 ∈ [0, 1].This setting is considered by Baddeley and Nair (2012).The DPP approximation outperforms the Poisson-saddlepoint approximation especially when γ 1 is close to zero, i.e. for very repulsive point processes.More numerical illustrations are displayed in Section 4. The rest of the paper is organized as follows.Section 2 provides necessary notation and background material on point processes, Gibbs point processes and determinantal point processes.Intensity approximations are discussed in detail in Section 3. Finally, Section 4 presents numerical experiments for several classes of pairwise interaction point processes.(2003) and Møller and Waagepetersen (2004).

Background and Poisson point processes
A spatial point process is said to have an nth order intensity function ρ (n) if for any nonnegative measurable function h : (R d ) n → R + , the following formula referred to as Campbell-Mecke formula holds where the sign = over the sum means that u 1 , . . ., u n are pairwise distinct.Then, ρ (n) (u 1 , . . ., u n ) du 1 • • • du n can be interpreted as the approximate probability for X having a point in each of infinitesimally small regions around u 1 , . . ., u n of volumes du 1 , . . .du n , respectively.We also write ρ(u) for the intensity function ρ (1) (u).A spatial point process X in R d is said to be stationary (respectively isotropic) if its distribution is invariant under translations (respectively under rotations).When X is stationary, the intensity function reduces to a constant denoted by λ in the rest of this paper.As a matter of fact, λ measures the mean number of points per unit volume.
The Poisson point process, often defined as follows (see e.g.Møller and Waagepetersen (2004)), serves as the reference model.Definition 2.1.Let ρ be a locally integrable function on S, for S ⊆ R d .A point process X satisfying the following statements is called the Poisson point process on S with intensity function ρ: • for any m ≥ 1, and for any disjoint and bounded B 1 , . . ., B m ⊂ S, the random variables X B1 , . . ., X Bm are independent; • N (B) follows a Poisson distribution with parameter B ρ(u) du for any bounded B ⊂ S.
Among the many properties of Poisson point processes, it is to be noticed that the nth order intensity function writes ρ Let Z be a unit rate Poisson point process on S, which means that its intensity is constant and equal to one.Assume, first, that S is bounded (|S| < ∞).We say that a spatial point process X has a density f if the distribution of X is absolutely continuous with respect to the one of Z and with density f .Thus, for any nonnegative measurable function h defined on N , E h(X) = E(f (Z)h(Z)).Now, suppose that f is hereditary, i.e., for any pairwise distinct u 0 , u 1 , . . ., u n ∈ S, f ({u 1 , . . ., u n }) > 0 whenever f ({u 0 , u 1 , . . ., u n }) > 0. We can then define the so-called Papangelou conditional intensity by for any u ∈ S and x ∈ N , setting 0/0 = 0.By the interpretation of f , λ(u, x) du can be considered as the conditional probability of observing one event in a small ball, say B, centered at u with volume du, given that X outside B agrees with x.When f is hereditary, there is a one-to-one correspondence between f and λ.
Because the notion of density for Z when S = R d makes no sense, the Papangelou conditional intensity cannot be defined through a ratio of densities in R d .But it still makes sense as the Papangelou conditional intensity can actually be defined at the Radon-Nykodym derivative of P !u the reduced Palm distribution of X with respect to P, the distribution of X (see Daley and Vere-Jones (2003)).We do not want to enter in too much detail here and prefer to refer the interested reader to Coeurjolly et al. (2017).
Finally, we mention the celebrated Georgii-Nguyen-Zessin formula (see Georgii, 1976;Nguyen and Zessin, 1979), which states that for any h : R d ×N → R (such that the following expectations are finite) (2.3) By identification of (2.1) and ( 2.3), we see a link between the intensity function of a point process and the Papangelou conditional intensity: for any which in the stationary case reduces to λ = E (λ(0, X)) . (2.4)

Gibbs point processes
For a recent and detailed presentation, we refer to Dereudre (2017).Gibbs processes are characterized by an energy function H (or Hamiltonian) that maps any finite point configuration to R ∪ {∞}.Specifically, if |S| < ∞, a Gibbs point process on S associated to H and with activity β > 0 admits the following density with respect to the unit rate Poisson process: where ∝ means "proportional to".This definition makes sense under some regularity conditions on H, typically non degeneracy (H(∅) < ∞) and stability (there exists A ∈ R such that H(x) ≥ A|x| for any x ∈ N ).Consequently, configurations x having a small energy H(x) are more likely to be generated by a Gibbs point process than by a Poisson point process, and conversely for imsart-ejs ver.2014/10/16 file: appIntDPP.texdate: November 13, 2018 configurations having a high energy.In the extreme case where H(x) = ∞, then x cannot, almost surely, be the realization of a Gibbs point process associated to H.
In this paper, we focus on pairwise interaction point processes.To be close to the original paper by Baddeley and Nair (2012) the present contribution is based on, we use their notation: a Gibbs point process in S is said to be a pairwise interaction point process with pairwise interaction function g : If |S| = ∞, this definition and more generally Definition (2.5) do not make sense since H(x) can be infinite or even undefined if |x| = ∞.In this case, Gibbs point processes have to be defined via their conditional specifications and for pairwise interactions Gibbs point processes, restrictions on g have to be imposed for existence (see again Dereudre (2017) and the references therein for details).Nonetheless, as mentioned in the previous section, the concept of Papangelou conditional intensity applies whenever |S| < ∞ or |S| = ∞, and in either case it has the explicit form for any u ∈ S. Note that when S = R d , a pairwise interaction Gibbs point process is stationary if g is symmetric and it is further isotropic if g(u − v) depends simply on u − v .From (2.4), we deduce that the intensity parameter of a stationary pairwise interaction process writes (2.7) Let us give a few examples (which are in particular well-defined in R d ).Many other examples can be found e.g. in the recent monograph by Baddeley et al. (2015).
Let us note that a Strauss model with γ = 0 and radius R is actually a hard-core model with radius R. The Diggle-Graton potential can be found in Baddeley et al. (2015) in a slightly different parameterization.The one chosen here makes comparisons with the Strauss model easier.For instance, when γ = 0 the model reduces to a Strauss model with γ = 0 and radius R. When γ = 1, the function g grows linearly from 0 to 1. Figure 2 depicts the form of some of the pairwise interaction functions presented above.A Gibbs point process has a finite range R if for any u ∈ R d and x ∈ N , λ(u, x) = λ(u, x ∩ B(u, R)).For pairwise interaction point processes, this property translates to g(u) = 1 for any u ∈ R d such that u > R. All previous models have a finite range R < ∞.An example of infinite range pairwise interaction point process which will not be considered in this paper is the Lennard-Jones model (see e.g.Ruelle (1969); Baddeley et al. (2015)).

Determinantal point processes
Determinantal point processes (DPPs) are models for inhibitive point patterns.We refer to Lavancier et al. (2015) for their main statistical properties.They are defined through a kernel function K which is a function from S × S to C. A point process is a DPP on S with kernel K, denoted by DPP(K), if for any n, its nth order intensity function takes the form (2.9) for every (u 1 , . . ., u n ) ∈ S n , where [K](u 1 , . . ., u n ) denotes the matrix with entries K(u i , u j ), 1 ≤ i, j ≤ n.In particular, the intensity function of DPP(K) is K(u, u).
Conditions on the kernel K are required to ensure the existence of DPP(K).For our purpose, we will only consider DPPs on a compact set.So let us assume that S is compact and suppose that K is a continuous real-valued covariance function on S × S. In this setting, by the Mercer theorem (see Riesz and Nagy (1990)), K admits the spectral expansion where {φ i } i≥1 is an orthonormal basis of L 2 (S) and where λ i , i ≥ 1, are referred to as the eigenvalues of K.Under the above assumptions, DPP(K) exists if and only if λ i ≤ 1 for all i.Due to their tractability, DPPs have many interesting properties.Many of them have been obtained by Shirai and Takahashi (2003), from which we derive the following key-equation used by our intensity approximation.
Proposition 2.2.Let X be a DPP on a compact set S with kernel K. Assume that K is a continuous real-valued covariance function on S × S whose all eigenvalues are not greater than 1.For any function g where λi , for i ≥ 1, are the eigenvalues of the kernel K : S × S → R given by Proof.Note that where L X denotes the Laplace transform of X. From Theorem 1.2 in Shirai and Takahashi (2003), for any nonnegative measurable function f on S L X (f ) = Det(I − K) where Det denotes the Fredholm determinant of an operator and K is the integral operator associated to the kernel On the other hand, see for instance (2.10) in Shirai and Takahashi (2003), where Tr denotes the trace operator.The result follows from the fact that for any n ≥ 1

Poisson-saddlepoint approximation
We remind that the intensity parameter of a Gibbs point process, and in particular a pairwise interaction point process satisfies (2.7).The expectaction in (2.7) is to be regarded with respect to P the distribution of the Gibbs point process X. Baddeley and Nair (2012) suggest to replace P by a simpler distribution, say Q, for which the right-hand-side of (2.7) becomes tractable.The Poissonsaddlepoint approximation consists in choosing Π(λ), the Poisson distribution with parameter λ, as distribution Q.As a result, the Poisson-saddlepoint approximation consists in resolving the equation with the convention that log 0 = −∞ and where, to avoid any ambiguity, we denote by Y a Poisson point process with intensity λ defined on R d and stress also this by indexing the E with the distribution Π(λ).It turns out that if g(u) ∈ [0, 1] for any u ∈ R d , the right-hand side of (3.1) is the Laplace transform of some Poisson functional and equals β exp(−λG) where G = R d (1−g(u)) du, see e.g.Møller and Waagepetersen (2004, Proposition 3.3).As noticed in Baddeley and Nair (2012), this formula extends to more general functions g, provided G > −∞.Hence, the Poisson-saddlepoint approximation, denoted by λ PS in this paper, is defined as the solution of imsart-ejs ver.2014/10/16 file: appIntDPP.texdate: November 13, 2018 Coeurjolly and Lavancier/Intensity approximation using DPP where W is the inverse function of x → x exp(x).
For stationary pairwise Gibbs models with finite range R, and such that λ(u, x) ≤ β (or equivalently such that g ≤ 1), then 0 ≤ G ≤ |B(0, R)|.In this case, Baddeley and Nair (2012) prove, among other properties, that λ ps exists uniquely and is an increasing function of β.From a numerical point of view, λ PS can be very efficiently and quickly estimated using root-finding algorithms.

DPP approximation
Following the same idea as the Poisson-saddlepoint approximation, for a repulsive stationary pairwise interaction point process with pairwise interaction function g ≤ 1 having a finite range R, we suggest to substitute the measure P involved in the expectation (2.7) by the measure Q corresponding to a DPP defined on B(0, R) with some kernel K (to be chosen) and intensity λ, i.e.K(u, u) = λ.Similarly to the previous section, by letting DPP(K; λ) denote the distribution of such a DPP and Y ∼ DPP(K; λ), the DPP approximation of the intensity λ is the solution of (3.3) From Proposition 2.2 and in particular from (2.11), this yields the estimating equation where the eigenvalues λi of K are related to λ by the relation To complete this approximation, the eigenvalues λi need to be specified.
In the following we choose the eigenvalues λi to be zero except a finite number N of them that are all equal.Given that this means that for some N ≥ λG, λi = λG N , for i = 1, . . ., N (3.4) and λi = 0 for i ≥ N + 1.With this choice, the integer N remains the single parameter to choose in our approximation.Note that N ≥ λG is a necessary condition to ensure λi ≤ 1 and so the existence of a DPP with kernel K, but it is in general not sufficient to ensure the existence of the relation between K and K where K defines a DPP.This will be clearly illustrated below when g is the Strauss interaction function.For the choice (3.4), the DPP approximation of the intensity, denoted by λ dpp , becomes the solution of To motivate (3.4) and how we should set N , assume for a moment that g is the interaction function of a Strauss model with range R and interaction parameter γ ∈ [0, 1], see (2.8).In this case K(u, v) = (1 − γ)K(u, v) for any u, v ∈ B(0, R) and the eigenvalues λ i of K satisfy λi = (1−γ)λ i .In the approximation (3.3), we start by choosing a kernel K with a finite number of non-vanishing eigenvalues λ i that are all equal.In view of λ i = K(u, u)du = λb, where b denotes the volume of B(0, R), this leads to λ i = λb/N for i = 1, . . ., N and N ≥ λb.Note that the latter inequality is necessary to ensure the existence of DPP(K).Going back to λi , this means that (3.4) follows with the necessary and sufficient condition N ≥ λb = λG/(1 − γ) which is greater than λG.
In order to set N precisely for the Strauss model, remember that a homogeneous DPP is more repulsive when its eigenvalues are close to 1, see Lavancier et al. (2015); Biscio and Lavancier (2016), and at the opposite a DPP is close to a Poisson point process when its eigenvalues are all close to 0. This suggests that in order to make the approximation (3.3) efficient, we should choose λ i close to 1 when the Gibbs process we want to approximate is very repulsive, that is when γ is close to 0. Moreover the eigenvalues should decrease to 0 when γ increases to 1.If λ i = λb/N , this is equivalent to choosing N an integer that increases from λb to infinity when γ increases from 0 to 1.A natural option is thus to choose N as the smallest integer larger than λb/(1 − γ).Our final choice for the Strauss model is therefore N = λb/(1 − γ) , where .denotes the ceiling function, which we may write, for later purposes, N = λG/(1 − γ) 2 .
However, with the latter choice, the function in the right-hand side of equation (3.5) is not continuous in λ, which may lead to none or several solutions to this equation.As a last step in our approximation, we therefore consider the upper convex envelope of this function, ensuring a unique solution to (3.5).This finally leads for the Strauss interaction process to the approximation λ dpp defined as the solution of Let us now discuss the case of a general pairwise interaction function g.In this setting, it is in general not possible to relate the eigenvalues λ i of K with the eigenvalues λi of K. Motivated by the Strauss case, we choose λi as in (3.4) where N = λG/κ and κ ∈ [0, 1] is a parameter that takes into account the repulsiveness encoded in g.In general κ must be close to 0 when g is close to 1 (the Poisson case), and close to 1 when g is close to a pure hard-core interaction.We decide to quantify the repulsiveness of the model by b −1 (1 − g) 2 , in agreement with our choice for the Strauss model for which κ = (1 − γ) 2 .Note that for a pairwise interaction g with range R and involving a possible hard-core distance δ, we have |B(0, δ)| ≤ (1 − g) 2 ≤ |B(0, R)| where the left or right equality occurs for a pure hard-core interaction (if δ > 0), a situation where κ must be 1.This leads us to the choice Plugging N = λG/κ into (3.5) and considering the upper convex envelope to ensure the existence of a unique solution, we finally end up with our general DPP approximation being the solution of where κ is given by (3.6) and where W κ is the inverse function of . In view of (3.7)-(3.8)and similarly to λ ps , the approximation λ dpp can be very efficiently implemented using root-finding algorithms.We further have the following properties.
Theorem 3.1.Consider a stationary pairwise interaction process in R d with Papangelou conditional intensity given by (2.6) which is purely inhibitory, i.e. g(u) ≤ 1 for all u ∈ R d and with finite range R.Then, λ dpp exists uniquely, is an increasing function of β and is such that λ dpp ≤ λ ps .

Numerical study
In this section, we focus on the planar case to investigate the performances of the DPP approximation and compare it with the initial one proposed by Baddeley and Nair (2012).All computations were performed in the R language (R development core team, 2011).The Poisson-saddlepoint approximation as well as the DPP approximation are implemented using root-finding algorithms and in particular we use the R function uniroot for this task.
We have considered 14 different numerical experiments involving Strauss models (S), Strauss hard-core models (SHC), Diggle-Graton models (DG), piecewise Strauss models (PS) and piecewise Strauss hard-core (PSHC) models.The pairwise interaction functions of these models are detailed in Section 2.2.To sum up here are the parameters, that include a continuously varying parameter γ 1 ∈ [0, 1]: For each numerical experiment, we therefore obtain curves of intensity approximation in terms of γ 1 .For γ 1 varying from 0 to 1 by step of 0.05 (the value 0 is exluded for DG models to save time), the true intensity λ is estimated by Monte-Carlo methods.For each set of parameters m realizations of the model are generated on the square [−2R, 1 + 2R] 2 and then clipped to the unit square.That strategy is detailed and justified by Baddeley and Nair (2012).Specifically, the number of points in each realization is averaged to obtain the estimated intensity and its standard error.The simulation results for the Strauss models with β = 50 or 100 were obtained by Baddeley and Nair (2012), where m = 10000 realizations were generated and the exact simulation algorithm was used, implemented in the R function rStrauss of the spatstat package (see Baddeley et al. (2015)).For the Strauss models with β = 200, SHC models, PS and PSHC models, we generate m = 1000 replications and use the rmh function in the spatstat package which implements a Metropolis-Hastings algorithm.Even if we use 10 6 iterations of the algorithm, the results may be slightly biased.For the DG models, the R package spatstat provides an exact simulation algorithm (function rDiggleGraton) and for such models we generate 10000 replications when β = 200 and R = 0.025, 0.05 and when β = 50 and R = 0.15.We used 1000 replications when β = 200 and R = 0.075 to save time.
All results can be found in Figures 3, 4 and 5. Plots provide the same information: we depict intensity approximation λ based on different methods in terms of γ 1 .The dashed curve represents the Poisson-saddlepoint approximation proposed by Baddeley and Nair (2012) and detailed in Section 3.1.The solid curve is the DPP approximation we propose in this paper and is given by (3.8).
Let us first comment Figure 3 dealing with Strauss models.As expected the Poisson-saddlepoint approximation is not efficient when γ 1 is small, i.e. for very repulsive models.This is very significant in particular for the Strauss hard-core model, see Figure 3 (f).The DPP aproximation we propose is more likely able to capture the repulsiveness of the Strauss models.Figure 4 also clearly shows that our approximation is particulalry efficient and outperforms unamibigously the Poisson-saddlepoint approximation.Note that replications for the Diggle-Graton models are generated using an exact algorithm; so the numerical results seem to be exact, except the slight bias induced by clipping the pattern from [−2R, 1 + 2R] 2 to the unit square.
We finally comment Figure 5.When γ 2 = 0.5, i.e.Figures 5 (a)-(b), the results are very satisfactory.Our approximation is able to approximate λ very efficiently for any value of γ 1 .For Figures 5 (c)-(d), γ 2 = 0 which means that points within a distance comprised between 0.05 and 0.1 are forbidden.Such a parameterization tends to create repulsive clusters.When γ 1 = 1 and δ = 0, such a piecewise Strauss model was called annulus model by Stucki and Schuhmacher (2014).This model demonstrates the limitations of our approximation even if when γ 1 is close to zero which means that the model is close to a hard-core process with radius 0.1 our approximation remains satisfactory.

Fig 1 .
Fig 1.Comparison of the exact intensity (small boxplots), the Poisson-saddlepoint approximation (dashed line) and the DPP approximation (solid line) for homogeneous Strauss models with activity parameter β and range of interaction R. Curves and boxplots are reported in terms of the interaction parameter γ 1 ∈ [0, 1].

Fig 4 .
Fig 4. Comparison of the exact intensity (small boxplots obtained by Monte-Carlo method), the Poisson-saddlepoint approximation (dashed line) and the DPP approximation (solid line) for Diggle-Graton models.Curves and boxplots are reported in terms of the interaction parameter γ 1 ∈ [0, 1].
For d ≥ 1, let X be a spatial point process defined on R d , which we see as a random locally finite subset of R d .Local finiteness of X means that X B = X∩B is finite almost surely (a.s.), that is the number of points N (B) of X B is finite R d is bounded.We let N stand for the state space consisting of the locally finite subsets (or point configurations) of R d .Let B(R d ) denote the class of bounded Borel sets in R d .For any B ∈ B(R d ), we denote by |B| its Lebesgue measure.A realization of X B is of the form x = {x 1 , . . ., x m } ⊂ B for some nonnegative finite integer m and we sometimes denotes its cardinal by |x|.For further details about point processes, we refer to Daley and Vere-Jones imsart-ejs ver.2014/10/16 file: appIntDPP.texdate: November 13, 2018 a.s., whenever B ⊂