Bayesian Strategies for Repulsive Spatial Point Processes
Chaoyi Lu Nial Friel





Bayesian Strategies for Repulsive Spatial 
 Point Processes




Chaoyi Lu[Address for correspondence: ], Nial Friel 

School of Mathematics and Statistics, University College Dublin, Ireland

Insight Centre for Data Analytics, University College Dublin, Ireland




April 29, 2024











There is increasing interest to develop Bayesian inferential algorithms for point process models with intractable likelihoods. 
A purpose of this paper is to illustrate the utility of using simulation based strategies, including approximate Bayesian computation (ABC) and Markov chain Monte Carlo (MCMC) methods for this task. 
<cit.> proposed an extended version of an ABC approach for repulsive spatial point processes, including the Strauss point process and the determinantal point process, but their algorithm was not correctly detailed. We explain that is, in general, intractable and therefore impractical to use, except in some restrictive situations.
This motivates us to instead consider an ABC-MCMC algorithm developed by <cit.>. 
We further explore the use of the exchange algorithm, together with the recently proposed noisy Metropolis-Hastings algorithm <cit.>. 
As an extension of the exchange algorithm, which requires a single simulation from the likelihood at each iteration, the noisy Metropolis-Hastings algorithm considers multiple draws from the same likelihood function.
We find that both of these inferential approaches yield good performance for repulsive spatial point processes in both simulated and real data applications and should be considered as viable approaches for the analysis of these models.






  
Keywords and Phrases:  Repulsive spatial point processes; noisy Metropolis-Hastings algorithm; ABC-MCMC algorithm; parallel computation; Strauss point process; Determinantal point process.



§ INTRODUCTION


A spatial point process is a random pattern of points with both the number of points and their locations being random. 
A canonical example is the Poisson point process whose number of points follows Poisson distribution with event locations being independent and identically distributed with density proportional to some intensity function. 
Such a fundamental point process model is a special case of several main classes of spatial point processes. 
These include Cox processes <cit.>, including log Gaussian Cox processes <cit.>, and the shot-noise Cox processes <cit.> for which a Neyman-Scott process <cit.> is a special case. 
Gibbs point processes (GPP) <cit.>, including Markov point processes and pairwise interaction point processes 
where a typical example is the Strauss point process. 
This is the first model which we focus on in this paper. 
The Strauss point process is an example of a doubly-intractable model where both the normalising terms of the posterior and the likelihood functions are unavailable. 
Our second model of interest is the determinantal point process (DPP) <cit.>. 
Although not strictly an intractable likelihood model, it is typically computationally expensive to evaluate. 
A determinantal point process model is usually defined over a Borel set, for example, ℝ^d with d being the event dimension. Simulation from such a process is required to be restricted to a specific area in practice, for example, the unit square [0,1]×[0,1]. 
The restricted version of the DPP model requires a spectral representation of the kernel and such a density function is treated as the true likelihood for the perfect sampler we focus on in our experiments.

Frequentist approaches to infer spatial point process models include those based on maximum likelihood estimation. 
However such approaches are complicated for the Strauss point processes. 
A much quicker alternative inference method is based on the pseudo-likelihood function baddeley2000practicaljensen1991pseudolikelihoodgelfand2010handbook, which is specified in terms of the Papangelou conditional intensity. 
In this paper we focus on Bayesian approaches.
The exchange algorithm, proposed by <cit.> for tackling doubly-intractable problems, is a natural algorithm for practitioners to consider for Gibbs-type likelihood, although it is not well explored in the context of repulsive spatial point processes. It relies on simulating from the intractable likelihood model. However, in practice, it can be difficult to improve the mixing or efficiency of such an algorithm other than the computationally expensive as-long-as-possible implementations which are applied in practice. 
Approximate Bayesian computation (ABC) methods are another popular class of tools to deal with the doubly-intractable problems where a closed-form likelihood function is not required during the implementations.
<cit.> first introduced ABC methods as a rejection technique in population genetics.
A generalized version was produced by <cit.> where a tolerance threshold was introduced when measuring the distance between the summary statistics of observed data and those of the simulated data from the likelihood model. 
However, it instead targets an approximate posterior distribution where a smaller tolerance level leads to better accuracy. 
As such, practitioners of ABC are faced with an efficiency and accuracy trade-off when implementing such a class of methods.

<cit.> adopts an ABC-MCMC approach which can be implemented for doubly-intractable models. 
Following a similar approach, <cit.> propose an ABC approach for spatial point process models where they incorporate an ABC-like rejection sampling step as the proposal step of the MCMC scheme.
A semi-automatic approach proposed by <cit.> is further adopted by <cit.> in order to find an optimal choice of the summary statistic of the observed data for model parameters.
However, we show in Section <ref> that the proposed ABC-MCMC algorithm in <cit.> leads to a stationary distribution with an intractable multiplicative term that depends on model parameters, a fact apparently overlooked by the authors, in effect rendering this algorithm impossible to implement. 

An objective of this paper is to illustrate the utility of using the noisy Metropolis-Hastings (noisy M-H) algorithm <cit.> and a special of it, the exchange algorithm <cit.>. Both algorithms, in common with their ABC counterparts, rely on simulation from the intractable likelihood model. 
In the case of the exchange algorithm, it uses a single draw from the likelihood function at each iteration of the Markov chain, while the noisy M-H algorithm, relies on more than one draw from the likelihood model. In effect, the likelihood draws are used to construct a Monte Carlo estimate of the ratio of the likelihood intractable normalizing constants which is then plugged into the usual Metropolis-Hastings algorithm. The noisy M-H algorithm has not yet received much attention in the context of repulsive spatial point processes, and thus we propose to illustrate its performance by comparing it to some candidate algorithms as outlined above, the exchange and <cit.> ABC-MCMC algorithms. 
Parallel computation can be implemented to carry out multiple draws within the noisy M-H algorithm yielding significant improvement in computational run time and efficiency. 

The structure of this paper is as follows. 
A basic introduction of the Strauss point process model is shown in Section <ref>. 
We provide a description of Markov chain Monte Carlo methods including the noisy Metropolis-Hastings algorithm, and approximate Bayesian computation methods in Section <ref>. 
An overview of the ABC algorithm proposed by <cit.> as well as the corresponding corrected ABC-MCMC algorithm are included in Section <ref>. 
While Section <ref> introduces the determinantal point processes. 
Simulation studies and real data applications are demonstrated, respectively, in Section <ref> and <ref>. 
Conclusions and discussions follow in Section <ref>.



§ STRAUSS POINT PROCESSES



A finite point process X with realization x:= {x_1,x_2,…,x_N} from a finite point configuration space N_f:={x⊂ S:N<∞} on a bounded spatial domain S ⊂ℝ^d is said to be a Gibbs point process (GPP) if it admits a density f(x)∝exp(∑_∅≠y⊆xQ(y)) with respect to a homogeneous Poisson point process with unit intensity. 
Setting exp(-∞)=0, the potential Q(x)∈[-∞,∞) is defined for all non-empty finite point patterns x⊂ S. Conditional on N = n, the joint probability density of x is of the form

    p_n(x_1,x_2,…,x_n|N=n)
       ∝exp(-μ(S))/n!μ(S)^n f({x_1,…,x_n}),

where μ(S):=∫_Sdu is the corresponding intensity measure on S, and the normalising term forms the distribution of N, that is,

    P(N=n)
       = exp(-μ(S))/n!∫_S…∫_Sf({x_1,…,x_n})dμ(x_1)…dμ(x_n).

In most cases, f is specified up to a proportionality constant and we denote f∝ h where h:N_f→ [0,∞) is a known function. For any x∈N_f and u ∈ S, u ∉x, the Papangelou conditional intensity for a point process X with density f is defined as

    λ(u,x)
       = f(x∪{u})/f(x),

which is a fundamental characteristic of a point process and we can interpret λ(u, x)du as the conditional probability of observing a point of the process inside an infinitesimal neighborhood of the location u provided that the rest of X is x. 
We say that a point process X is repulsive if λ(u, x)≥λ(u, y) whenever x⊂y.

A homogeneous pairwise interaction point process is a GPP with density

    f(x)
       ∝ h(x) := β^n(x)∏_{u,v}⊆xψ ({u, v}),

where β > 0 is a constant and n(x) is the number of events in x. Moreover, ψ : S × S → [0,∞) is an interaction function defined as ψ({u, v}):=ψ_0(||u-v||) where ψ_0: (0, ∞)→ [0,∞) is invariant under reflections, rotations and translations. The Strauss point process (SPP) is the simplest non-trivial homogeneous pairwise interaction process with

    ψ_0(||u-v||)
       = γ^1(||u-v|| ≤ R).

The density of the Strauss point process is of the form 

    f(x) 
       ∝ h(x) := β^n(x)γ^s_R(x),

where β>0, 0≤γ≤ 1 and 

    s_R(x)
       := ∑_{ u, v }⊆x1(||u - v|| ≤ R),   u, v ∈ S,

which is the number of R-close pairs of points in x. Further, 1(·) is the indicator function, β is known as the rate and γ is the interaction parameter where smaller values of γ implies stronger repulsion. 
The process is reduced to a homogeneous Poisson point process if we set γ = 1, and the case with γ = 0 is called the hard-core process with the hard-core radius R/2 where each pair of points are not allowed to be closer than R units apart.




§ MARKOV CHAIN MONTE CARLO AND APPROXIMATE BAYESIAN COMPUTATION



In this paper, we focus on the 2-dimensional case, that is, d=2 and, without loss of generality, we restrict our attention to the area S = [0,1]×[0,1].
We consider a posterior distribution π(y|θ) ∝ℒ(y|θ) π(θ) with a likelihood function of the form ℒ(y|θ)= q(y|θ)/z(θ), where y and θ denote the observed data and model parameter(s), respectively. The terms z(θ) and q(·|θ) denote the normalizing constant and the unnormalized likelihood function, respectively. The prior distribution for θ is given by π(θ).
In this situation, the Metropolis-Hastings (M-H) algorithm, outlined in Algorithm <ref> does not apply since the terms z(θ') and z(θ^(t-1)) are both intractable.

In the case of a SPP, the posterior distribution is of the form

    π(β,γ|y)
       ∝β^n(y)γ^s_R(y)/z(β, γ)π(β, γ) = β^n(y)γ^s_R(y)π(β, γ)/z(β, γ),

where z(β, γ) denotes the normalizing constant of the SPP density (<ref>), and π(β,γ) is the joint prior distribution. Once again, it is clear from (<ref>) that the normalising constant in the denominator is intractable and this makes the acceptance ratio α_.85MH in the M-H algorithm <ref> no longer available to us.


In order to deal with doubly-intractable problems, ABC methods
can instead be considered. An important ingredient of ABC algorithms is the requirement to obtain a realization x' ∼ℒ(·|θ'), for model parameter(s) θ'. Additionally, one requires a summary statistic T(·) and a distance function Ψ(T(x'), T(y)), which describes the discrepancy, in some sense, between the observed data y and the realization x'.
One of the first ABC algorithms was proposed by <cit.> and is shown as Algorithm <ref>. In practice, the acceptance threshold ϵ will necessarily not be set to zero. Therefore this algorithm will not target the true posterior distribution and will instead sample from an approximate posterior distribution π_ϵ(θ|y), obtained by marginalising over x, the joint distribution π_ϵ(θ,x|y) ∝ℒ(x|θ)π(θ)1(Ψ(T(x),T(y))≤ϵ). 
The need for a small ϵ to ensure that accuracy of the resulting ABC target distribution to the posterior distribution is balanced with the requirement ensure the algorithm mixes sufficiently.
Thus the balance of accuracy and implementation efficiency is a key aspect practitioners need to consider. 
Implementing ABC methods requires an exact draw from the intractable likelihood function. 
Perfect simulation can be accomplished for the SPP by applying the “dominated coupling from the past” algorithm kendall2000perfect which can be implemented in the R package  baddeley2005spatstat.
An extended version of the ABC adopting with MCMC scheme for spatial point processes is outlined next in Section <ref>.


Though ABC methods are popular for tackling doubly-intractable problems, it is of interest to compare their performances to both the exchange algorithm and also the noisy Metropolis-Hastings (noisy M-H) algorithm recently proposed by <cit.>. We begin by illustrating the well-known exchange algorithm <cit.> as shown in Algorithm <ref>.
This algorithm samples from an augmented distribution π(θ', x', θ|y) ∝ℒ(y|θ)π(θ)p(θ'|θ)ℒ(x'|θ') whose marginal distribution for θ is our target distribution and the auxiliary function ℒ(x'|θ')=q(x'|θ')/z(θ') is the same likelihood function as ℒ(y|θ)=q(y|θ)/z(θ). Such an augmented distribution leads to a cancellation of the normalizing constants z(θ^(t-1)), z(θ') in the acceptance ratio α_.85Ex of the exchange Algorithm <ref>.

Notice that the auxiliary unnormalized likelihood ratio q(x'|θ^(t-1))/q(x'|θ') in (<ref>) of Algorithm <ref> is actually an unbiased importance sampling estimator of the normalising constant ratio z(θ^(t-1))/z(θ') in (<ref>) of Algorithm <ref>, that is,

    𝔼_x' ∼ℒ(·|θ')(q(x'|θ^(t-1))/q(x'|θ')) = z(θ^(t-1))/z(θ').

By applying a simple Monte Carlo estimator of this expectation, we can obtain an improved unbiased estimator of the normalizing term ratio which can be instead applied at Step 3 of the exchange Algorithm <ref>, as follows. An i.i.d. sample (x_1', x_2', …, x_K') from ℒ(·|θ') can be used to approximate the ratio of the normalizing constant as

    z(θ^(t-1))/z(θ') = 𝔼_x' ∼ℒ(·|θ')(q(x'|θ^(t-1))/q(x'|θ')) ≈1/K∑_k=1^K q(x_k'|θ^(t-1))/q(x_k'|θ').

In turn, this estimator can be plugged into (<ref>) yielding the noisy Metropolis-Hastings (noisy M-H) algorithm <cit.> as shown in Algorithm <ref>. In order to make the algorithm more efficient, parallel computation can be applied to sample K auxiliary chains via the  R package.

Note that for K = 1, the noisy M-H algorithm is equivalent to the exchange algorithm and, as K →∞, the algorithm will become the standard Metropolis-Hastings algorithm. This indicates that both K = 1 and K →∞ leave the target distribution invariant, but this is not guaranteed for 1 < K < ∞. <cit.> proved that a Markov chain resulting from the noisy M-H algorithm will converge to the target posterior density as K →∞ under certain assumptions:

We refer to <cit.> and <cit.> for more details.


 
Note that <cit.> also provide results for the case where the assumption of uniform ergodicity is replaced by the less restrictive assumption of geometric ergodicity and we refer the reader to that paper for more details. In what follows we assume that assumption (A1) holds, although this may be difficult to prove in practice. 
Note that both (A2) and (A3) are satisfied in the case where one uses bounded proposal and prior distributions. Furthermore, <cit.> also showed that (A4) is satisfied for Gibbs random fields for which the SPP is a special case by considering h(x) = β^n(x)γ^s_R(x)=exp(n(x)log(β)+s_R(x)log(γ)). Thus, assuming that all four assumptions above hold, then Theorem 3.1 of <cit.> provides a theoretical guarantee that the Markov chain resulting from the noisy M-H algorithm for a SPP will 
converge to the target posterior density as K →∞. 

Of course, in practice one applies the noisy M-H algorithm with a finite number, K, of auxiliary chains. The objective in this case is to explore the potential improvement in mixing performance compared to the exchange algorithm. This possible accuracy-efficiency trade-off motivates us to also compare both accuracy and efficiency performances of the small-K noisy M-H algorithm to that of the exchange algorithm in this paper.




§ ABC ALGORITHMS FOR REPULSIVE SPATIAL POINT PROCESSES



The ABC algorithm that serves as a candidate algorithm for performance comparisons in this paper is highly related to <cit.> within which an ABC-MCMC-like algorithm was proposed and applied for repulsive spatial point processes.
The method in <cit.> is based on the likelihood-free MCMC idea explored in the ABC-MCMC algorithm proposed by <cit.> as well as a semi-automatic approach proposed by <cit.>. 
Such a semi-automatic approach argues that the optimal choice of the summary statistic T(y) in ABC methods is 𝔼(θ|y) and a linear regression scheme is considered to construct this summary statistic. 
Taking the SPP model as an example here, {(β_l,γ_l)}^L_l=1 are generated from prior distributions through a pilot run with each of the corresponding realizations {x_l}^L_l=1 generated from the SPP likelihood function ℒ(·|β_l,γ_l) (<ref>). 
The radius R = R̂ is estimated by profile pseudo-likelihood method <cit.>. 
Note that taking a log transform of the parameters can facilitate the regression, and thus the pilot draws {θ_l}^L_l=1 in the regression are stored as {(log(β_l),log(γ_l))}^L_l=1. 
The linear regression is implemented for 𝔼(θ_l|y) = a + bη(x_l, y) where η(x_l, y) is a vector of summary statistics

    η(x,y) = (η_1(x,y), η_2(x,y)),
 
with

    η_1(x,y) = log(n(x)) - log(n(y)),  η_2(x,y) = √(K̂_R̂(x)) - √(K̂_R̂(y))^2.
 
Here, K̂_R̂(x) is the empirical estimator of Ripley's K-function ripley1976secondripley1977modelling defined as

    K̂_r(x) = S∑_{ u, v }⊆x1[0 < ||u - v|| ≤ r]/n(n-1)e(u,v),
 
where e(u,v) is an edge correction factor introduced in <cit.>. 
Here we propose to leverage the same “isotropic” edge correction which was also used by <cit.>.

The generalized linear regression under a multi-response Gaussian family is fit with lasso regression, and cross-validation is applied to determine the penalty parameter for the lasso. 
After obtaining â and b̂ using L pilot draws, the distance measures {Ψ(θ̂_l, θ̂_obs) }^L_l = 1 can be calculated in order for setting the acceptance threshold ϵ in <cit.> ABC-MCMC algorithm as shown in Algorithm <ref> by taking p percent estimated percentile of those measures. 
Here, the value of p is set by practitioners prior to implementations, the measure Ψ(θ̂_l, θ̂_obs) := ∑_i (θ̂_l,i - θ̂_obs,i)^2/v̂âr̂(θ̂_i) is the component-wise sum of quadratic loss for the log parameter vector with θ̂_obs = â, and v̂âr̂(θ̂_i) is the sample variance of the ith component of θ̂.




 §.§ Correcting the <cit.> ABC-MCMC Algorithm


Here we explain that the <cit.> ABC-MCMC Algorithm <ref> does not follow exactly the ABC-MCMC scheme proposed in <cit.> or the semi-automatic ABC-MCMC algorithm proposed by <cit.> as shown in Algorithm <ref>.



Instead of an accept-or-stay step shown as Step 1 of Algorithm <ref>, <cit.> propose to leverage an ABC-like rejection sampling step as a proposal step within the Metropolis-Hastings algorithm, that is, Step 1 of Algorithm <ref>. 
The proposal density resulting from Step 1 of this algorithm is

    p_ϵ(θ',x'|θ)
       = p(θ'|θ)ℒ(x'|θ')1_A_ϵ,y(x')/∫_A_ϵ,yθ'p(θ'|θ)ℒ(x'|θ')dx'dθ',

where the set A_ϵ,y:={x'∈N_f: Ψ(θ̂', θ̂_obs) ≤ϵ}, and 1_A(x) is an indicator function giving 1 if x∈ A and giving 0 otherwise.
Note that the proposal density (<ref>) involves a normalizing term which we denote as 
ζ(θ)=∫_A_ϵ,yθ'p(θ'|θ)ℒ(x'|θ')dx'dθ' and which is intractable in general, a point apparently missed by <cit.>.
Note also that <cit.> do not explicitly detail the ergodic distribution resulting from their algorithm. 
However, one can see from the acceptance probability in Step 2 of Algorithm <ref>, that the target distribution must be written as

    π̂_ϵ(θ,x|y)
       ∝π(θ)ℒ(x|θ)1_A_ϵ,y(x)ζ(θ).

This is to ensure that the ratio in α_.85S&G in Algorithm <ref> (comprised of the usual target ratio multiplied by proposal ratio) appears as:

    p(θ|θ')ℒ(x|θ)1_A_ϵ,y(x)/ζ(θ')/p(θ'|θ)ℒ(x'|θ')1_A_ϵ,y(x')/ζ(θ)×π(θ')ℒ(x'|θ')1_A_ϵ,y(x')ζ(θ')/π(θ)ℒ(x|θ)1_A_ϵ,y(x)ζ(θ)   = π(θ')p(θ|θ')/π(θ)p(θ'|θ)
    .

As a consequence, the target distribution (<ref>) resulting from Algorithm <ref> differs from the target distribution, π_ϵ(θ,x|y)∝π(θ)ℒ(x|θ)1_A_ϵ,y(x), resulting from Algorithm <ref>, 
 by a multiplicative term, ζ(θ), which itself depends on the model parameter(s) θ. 
Further, it is not clear what effect the term ζ(θ) has on the target distribution.

In order to target the distribution π_ϵ(θ,x|y), the ratio in α_.85S&G should instead be modified as

    π(θ')ℒ(x'|θ')1_A_ϵ,y(x')/π(θ)ℒ(x|θ)1_A_ϵ,y(x)×p(θ|θ')ℒ(x|θ)1_A_ϵ,y(x)/ζ(θ')/p(θ'|θ)ℒ(x'|θ')1_A_ϵ,y(x')/ζ(θ)    = π(θ')p(θ|θ')ζ(θ)/π(θ)p(θ'|θ)ζ(θ'),

which leads to what we term the corrected <cit.> ABC-MCMC algorithm shown as Algorithm <ref>.

However, as before, the acceptance probability (<ref>) is generally intractable due to the intractability of the terms ζ(θ) and ζ(θ').
Thus, in general, Algorithm <ref> is impossible to implement in practice and we cannot see an immediate way to fix this issue.


Note that it is possible to approximate the intractable term ζ(θ) by considering that 

    ζ(θ)=∫_A_ϵ,yθ'p(θ'|θ)ℒ(x'|θ')dx'dθ'=E_p(θ'|θ)[E_ℒ(x'|θ')(1_A_ϵ,y(x'))],

where Monte Carlo approximations can estimate the double expectation (<ref>) based on multiple auxiliary draws of θ' and further auxiliary draws of x'|θ'.
While this is a compuationally intensive computation, in general, we remark that this might be implemented in parallel, providing some potential efficiency. 
Thus it would be interesting to explore whether Algorithm <ref> implemented with such Monte Carlo approximations is able to bring any improvement in efficiency or accuracy performance compared to a typical ABC or the <cit.> ABC-MCMC Algorithm <ref>.
We leave this as future work.





 §.§ A Specific Case of the Corrected <cit.> ABC-MCMC Algorithm


In restrictive situations, (<ref>) can become tractable.
Here, if we propose an independent proposal distribution such that p(·|θ)=p(·|θ')=p(·) 
for each iteration t of Algorithm <ref>, the intractable terms ζ(θ) and ζ(θ') coincide since,

    ζ(θ)=∫_A_ϵ,yθ'p(θ')ℒ(x'|θ')dx'dθ'≡∫_A_ϵ,yθp(θ)ℒ(x|θ)dxdθ=ζ(θ').

Thus both terms cancel in the ABC-MCMC Algorithm <ref>, which is a special case of the corrected <cit.> ABC-MCMC Algorithm <ref>.

In the case that the proposal distribution p(·) is set as the prior distribution π(θ), Algorithm <ref> becomes the <cit.> ABC algorithm shown as Algorithm <ref>.





§ DETERMINANTAL POINT PROCESSES



We now illustrate another popular and complex family of spatial point process models, the determinantal point process (DPP). After describing this model, we illustrate some variants of this model and related MCMC algorithms, building on those developed in Sections <ref> and <ref>.

A simple locally finite spatial point process X on ℝ^2 is called a determinantal point process with kernel C lavancier2015determinantal if it has a product density function

    ρ^(n)(x_1,…,x_n)
       = det[C](x_1,…,x_n),

where (x_1,…,x_n) ∈ (ℝ^2)^n for n = 1,2, …, and [C](x_1,…,x_n) is the n × n matrix with (i,j)th entry being C(x_i,x_j). Here C is a covariance kernel defined on ℝ^2×ℝ^2 and det(·) denotes determinant of the matrix. We write X ∼ DPP_ℝ^2(C) and, for any compact subset S ⊆ℝ^2, we denote by DPP_S(C), the distribution of the DPP on S with kernel given by the restriction of C to S × S. Note that DPP_S(C) is equivalent to the distribution of X_S=X∩ S. The Poisson point process is a special case of the DPP where C(x,y)=0 whenever x ≠ y. The intensity function is the first order density function: ρ^(1)(x) = C(x,x) for x ∈ℝ^2 and where the pairwise correlation function of X is

    g(x,y) = ρ^(2)(x,y)/ρ^(1)(x)ρ^(1)(y) = 1 - C(x,y)C(y,x)/C(x,x)C(y,y),

if C(x,x) > 0 and C(y,y) > 0. Otherwise it is equal to zero. Due to the fact that C is a real covariance kernel in our setting, the repulsiveness of the DPP is reflected by g ≤ 1 and  

    ρ^(n)(x_1,…,x_n)
       ≤ρ^(1)(x_1)⋯ρ^(1)(x_n),

for any n = 2,3,…. The inequality follows from the fact that the determinant of a covariance matrix never exceeds the product of its diagonal elements. We refer to <cit.> and <cit.> for more details and properties of DPPs.

Since simulation from a DPP over ℝ^2 is practically impossible, <cit.> instead proposed to simulate X_S within a restricted region S ⊆ℝ^2, for example, S = [0,1] × [0,1]. 
Thus the kernel C restricted to S × S has a spectral representation: C(x,y) = ∑^∞_k=1λ_kϕ_k(x)ϕ_k(y) where (x,y) ∈ S × S, and (·) denotes the complex conjugate.
While {ϕ_k}^∞_k=1 are eigenfunctions and {λ_k}^∞_k=1 are the set of eigenvalues which are required to be less than or equal to 1 in order to ensure the existence of the DPP. 
The simulation process is based on the fact that the distribution of the DPP_S(C) is the same as that of the so-called determinantal projection point process DPP_S(χ) with χ(x,y) = ∑^∞_k=1B_kϕ_k(x)ϕ_k(y) := v^*(y)v(x) where (·)^* denotes the conjugate transpose of the vector and {B_k}^∞_k=1 are independent Bernoulli variables with mean λ_k for k = 1,2,…. The corresponding simulation algorithm is shown as Algorithm <ref>. Note that the simulation process can be implemented via the  R package. The sampling process in Steps 2,3 are based on rejection sampling with a uniform density where the acceptance rate is p_i(x)/U with U=sup_y∈ S||v(y)||^2/i being an upper bound applied by <cit.> on p_i(x) for x ∈ S.


In this paper we focus on the case where the covariance kernel is stationary and isotropic. The density of X_S, which has realization (x_1,…,x_n)∈ S^n, is of the form

    f(x_1,…,x_n)
       = exp(|S|-D)det[C̃](x_1,…,x_n),

for n = 0,1,2,… where D=-logP(n=0)=-∑_k=1^∞log(1-λ_k) and C̃(x,y)=∑^∞_k=1λ_k/1-λ_kϕ_k(x)ϕ_k(y) with the setting that det[C̃](x_1,…,x_n)=1 if n =0. Note that both the simulation algorithm and density require the spectral representation of the covariance kernel. However, it is not always the case that such a representation is explicitly known. Thus <cit.> proposed to apply a Fourier series approximation of the kernel and instead consider the approximate DPP: X̂_S∼DPP_S(Ĉ) where

    Ĉ(x,y)
       = ∑_k∈ℤ^2φ(k)exp(2π ik·(x-y)), x,y ∈ S.

Defining C(x,y)=C_0(x-y), then φ(k), which is called spectral density, is the Fourier transform of C_0. The corresponding approximate density is of the form

    f̂(x_1,…,x_n)
       =  exp(|S|-D̂)det[Ĉ̃](x_1,…,x_n),

where Ĉ̃(x,y)=∑_k∈ℤ^2φ̃(k)exp(2π ik·(x-y)), D̂=∑_k∈ℤ^2log(1+φ̃(k)) and φ̃(k) = φ(k)/(1-φ(k)). However, summation over the whole ℤ^2 above is impossible in practice, thus a truncation M is applied by <cit.> and by default in the  package, so that ∑_k∈ℤ_M^2φ(k)>0.99n/|S| where ℤ_M = {-M,-M+1,…,M-1,M}. We refer to <cit.> for more details of this truncation approximation.

Similar to <cit.> and instead of the DPP_ℝ^2(C) over the whole ℝ^2, we also focus on the DPP_S(Ĉ) within the specific region S= [0,1] × [0,1] in our experiments. Thus Algorithm <ref> which is implemented in the  R package can therefore be treated as a perfect sampler from the likelihood function f̂ (<ref>). The covariance kernel we focus on in this paper is the Gaussian kernel defined as

    C(x,y)
       =  τexp(-||x-y||^2/σ^2),

with a spectral density of the form

    φ(k)=τ(√(π)σ)^2exp(-||πσ k||^2).

The existence of such a so-called determinantal point process with a Gaussian kernel (DPPG) is guaranteed by 0≤σ≤σ_max=1/√(πτ) or, equivalently, 0≤τ≤τ_max=1/(πσ^2).
It was shown in <cit.> that realizations simulated from f̂ provides the empirical means of L(r)-r being close to the corresponding theoretical L(r)-r function, where the L-function is the variance stabilizing transformation of the K-function which is defined as K(r):= π r^2-(1-exp(-2r^2/σ^2))πσ^2/2 for the Gaussian model. 
This indicates that the simulations from f̂ are appropriate approximations of the DPP_S(C) with a Gaussian kernel.

In contrast to the SPP, the normalizing term of f̂ is analytically available when the truncation is applied, so it may seem unnecessary to implement the exchange algorithm or the noisy M-H algorithm, since the Metropolis-Hasting algorithm is available to us. However, the normalizing term as well as the covariance kernel Ĉ between each pair of events in the likelihood function can be computationally expensive when M is large. Thus, in order to reduce the computational burden, we propose in this paper to also apply the exchange and noisy M-H algorithms by leveraging (<ref>) as an approximate density of X_S. 

More specifically, by treating f̂ as the likelihood function from which we apply a perfect sampler <ref> and denoting 

    q̂(x_1,…,x_n)
       =  det[Ĉ̃](x_1,…,x_n),

as the corresponding unnormalized likelihood function, we apply the noisy M-H Algorithm <ref> to infer the parameters θ = (τ,σ) of the DPP_S(Ĉ) with a Gaussain kernel and the exchange algorithm is the K=1 special case of such an algorithm. The target posterior distribution is thus of the form

    π(τ,σ|y)
       ∝q̂(y|τ,σ)π(τ,σ).


Define ρ(x):=ρ^(n)(x_1,…,x_n) for n=1,2,…, and then for any x∈ X_S, we consider the truncated density of the form

    f(x|x∈ X_S)
       =  1(x∈ X_S)·ρ(x)/P(x∈ X_S)∝ρ(x),

where 1(x∈ X_S) is the indicator function returning 1 if x∈ X_S and 0 otherwise. 
Here, P(x∈ X_S) is the normalising constant over all the possible events of X_S. 
By leveraging (<ref>) and (<ref>), we substitute all the q̂(·) in the acceptance ratio α_.85DPPG of Algorithm <ref> with the ρ(·),
and obtain the approximate noisy M-H Algorithm <ref>. The approximate exchange algorithm is the K=1 case of such an algorithm.




Note that this approximation does not require knowledge of the specific form of the normalizing constant P(x∈ X_S) in (<ref>) due to the cancellation of the normalizing constants in the acceptance ratio (<ref>). 
In this way, we no longer need to compute the computationally expensive summations shown in (<ref>) for each pair of (x_i,x_j) where i,j=1,…,n;  i≤ j.
Instead directly computing the Gaussian kernel (<ref>) significantly reduces the computation and at the same time does not lose much accuracy if Ĉ is an appropriate approximation of C as discussed by <cit.>. Similar MCMC approximation schemes can also be found, for example, in <cit.>, but what is different to our paper is that, instead of approximating the unnormalized likelihood, they proposed a variety of approximations for the likelihood normalizing constant or constant ratio. The motivations behind these approximations are similar to the interpretation of the unbiased importance sampling estimator of the normalizing constant ratio we discussed in Section <ref> for the exchange algorithm. 

The convergence assumptions (A2) and (A3) shown in Section <ref> can also be satisfied for the noisy M-H Algorithm <ref> by proposing bounded priors and proposals for the model parameters of DPPG. 
While (A4) holds by noticing that exp(|S|-D̂)≤f̂≤exp(|S|-D̂)(∑_k∈ℤ^2φ̃(k))^n for any n=0,1,2,…, where ∑_k∈ℤ_M^2φ(k) tends to τ from below as truncation M→∞. The comparison performances of the M-H, approximate exchange and approximate noisy M-H algorithms 
are shown in Section <ref>. 
The summary statistic 𝔼(θ|y) used in the <cit.> ABC-MCMC Algorithm <ref> for DPPG model is an extension of the SPP case.
Similar to <cit.>, we also leverage a set of {K̂_r(x)} evaluated at 10 equally spaced values over [0.01,0.1] for η_2(x,y) in our experiments, that is, η_2(x,y) = ( η_2,r_1(x,y),η_2,r_2(x,y), …, η_2,r_10(x,y) ) over the set (r_1 = 0.01, r_2 = 0.02,…, r_10 = 0.1), where η_2,r_i(x,y) = √(K̂_r_i(x))-√(K̂_r_i(y))^2 is defined in the same way as the SPP case. 




§ SIMULATION STUDIES



In this section, we compare the performances of the exchange algorithm, the noisy M-H algorithm and the <cit.> ABC-MCMC algorithm (Algorithms <ref>, <ref> and <ref>, respectively). We apply each algorithm to datasets randomly simulated from a SPP and a DPPG.
We focus mainly on noisy M-H algorithms with K≤ J, where J is the number of available processing cores, where in our experiments, J=7. This is for two reasons. 
Firstly, if K>J, then more computation time is needed to synchronize the output of the J likelihood draws, before starting another set of J auxiliary chains. 
Secondly, in practice, we find that not much additional statistical efficiency results if K>J.

Several summary statistics are chosen to assess the performances of the various algorithms in our experiments. 
The effective sample size (ESS) kass1998markov due to autocorrelation is usually defined as ESS(θ) = T/[1+2∑^∞_i=1ν_i(θ)], where T is the posterior sample size, ν_i is the autocorrelation at lag i, and the infinity sum is often truncated at lag i when ν_i(θ)<0.05. 
ESS is used to check the dependence and the autocorrelation of posterior samples, that is, the mixing performance, where larger values of ESS imply better mixing. 
Due to the fact that the computational runtime of the different algorithms varies considerably, we propose to monitor the ESS per second (ESS/sec), to assess the efficiency of the algorithms. 
All experiments are based on a CPU with a 1.80GHz processor and 7 cores. 
We present some basic statistics to compare the accuracy of the output from each algorithm including the posterior mean and standard deviation, posterior box plots or density plots, and the absolute value of the error between the posterior mean and the corresponding ground truth of the model parameters.



 §.§ Strauss Point Process



The SPP dataset y_1 was randomly simulated with parameters β = 200, γ = 0.1 and R = 0.05 on S = [0,1]×[0,1] and y_1 contains 83 point locations as shown in the left plot of Figure <ref>. 
We consider three different values of p for setting the acceptance threshold ϵ in the <cit.> ABC-MCMC algorithm, namely, 0.5, 1 and 2.5. 
The initial states of the β and γ parameters for all the algorithms are set to be β_0 = 190 and γ_0 = 0.2. 
The prior distributions are specified to be uniform distributions: β∼U(50, 400), γ∼U(0,1) and the interaction radius R is estimated by the same profile pseudo-likelihood method as used in <cit.> with the estimated value being R̂ = 0.0508 (right plot of Figure <ref>). 
Bounded uniform proposals are applied for both parameters conditional on current state at each iteration t, 
that is, the proposed state follows:

    β'    ∼U(max(50,β^(t-1)-ϵ_β ),min(400,β^(t-1)+ϵ_β)),
    γ'    ∼U(max(0,γ^(t-1)-ϵ_γ ),min(1,γ^(t-1)+ϵ_γ )),

with ϵ_β and ϵ_γ tuned to be 65 and 0.16, respectively, so that the acceptance rate of the exchange algorithm and the noisy M-H algorithm was around 0.25. 
The acceptance rate of the ABC-MCMC algorithm was around 0.01 to 0.1 depending on the values of p.
All three algorithms were implemented for 0.12 million iterations where the first 0.02 million are burn-in iterations. The Markov chain obtained by the exchange algorithm is known to target the true posterior distribution and thus an extra 1.2-million-iteration implementation of this algorithm was treated as a ground truth with 0.2-million iterations as burn-in. 
 






Table <ref> displays a table of output from different candidate algorithms. 
We denote the ABC-MCMC, exchange and noisy M-H algorithms as “ABC”, “Ex” and “NMH”, respectively. 
While the long run of the exchange algorithm which serves as the ground truth is denoted by “GT”.
Note that E(·), sd(·) and |Bias(·)|, respectively, represents the posterior mean, posterior standard deviation and the absolute value of the error between the posterior mean and the corresponding ground truth. ESS(Ave)/s records the average posterior ESS(·)/s statistics for the two parameters β and γ.

Overall, we can make the following observations, based on the results presented in Table <ref> and Figure <ref>. The ABC algorithms generally yield posterior point estimates that are more biased than the exchange and noisy M-H counterparts. As expected, as p decreases, the ABC algorithm yields improved point estimates, however, this is balanced with the fact that decreasing p corresponds to smaller ESS values and ESS(Ave)/s values. By contrast, the performance of the noisy M-H algorithms in terms of mixing is much improved compared to the ABC algorithms, as illustrated by improved ESS values by comparison. For example, the K=2 noisy M-H case yields a more than 10-fold increase in ESS(Ave)/s compared to the ABC p=0.5 case. 

We also note that the implementation time of the noisy M-H algorithm increases with K, the number of auxiliary chains, as one expects leading to smaller ESS(Ave)/s values. This is especially apparent when K=8. The reason for this, as before, 
is due to the fact that the experiments were run on a laptop with 7 cores and so the K=8-auxiliary-chain case requires another round of parallel computation which leads to significantly more implementation time.











 §.§ Determinantal Point Process with a Gaussian Kernel



Our second simulation study concerns a determinantal point process. 
Figure <ref> shows the observed n=99 locations randomly generated from a DPP with a Gaussian kernel (DPPG) with parameters τ = 100, σ = 0.05. 
The initial states of the two parameters in the Gaussian model are set to be τ_0 = 125 and σ_0 = 0.04 and, similar to Section <ref>, uniform priors and bounded uniform proposals are proposed for the corresponding parameters, that is, τ∼U(50,200) which includes the estimated τ̂ = n/|S|, and σ∼U(0.001, 1/√(50π)) which allows our proposed τ' to lie within our prior range. 
The proposed state in each iteration t follows: τ' ∼U(max(50,τ^(t-1)-ϵ_τ ),min(200,τ^(t-1)+ϵ_τ)) and σ' ∼U(max(0.001,σ^(t-1)-ϵ_σ ),min(1/√(πτ'),σ^(t-1)+ϵ_σ )) with ϵ_τ and ϵ_σ tuned to be 32 and 0.015, respectively.






Both the noisy M-H Algorithm <ref> and the approximate noisy M-H Algorithm <ref>, of which the exchange algorithm and the approximate exchange algorithm are, respectively, the specific K=1 case, are implemented for performance comparisons. 
The ABC-MCMC algorithm steps are similar to those of the SPP cases due to the fact that this algorithm does not require the evaluation of likelihood functions and this is one of the key differences compared to other candidate algorithms. 
Due to the tractability of the likelihood normalizing constant, the Metropolis-Hastings (M-H) algorithm is also available. 
However, considering the increased computational burden required by DPPG implementations compared to the SPP simulation study, 12000-iteration implementations with 2000-iteration burn-in are instead applied for all the algorithms outlined here. 
Indeed Figure <ref> indicates that the mixing is sufficiently adequate to sample from the stationary distribution, though it suggests that the ABC-MCMC algorithm mixes significantly worse than other candidates.
Here the M-H algorithm is additionally implemented for 0.12 million iterations and the corresponding posterior samples with 20000-iteration burn-in are treated as the ground truth. 






Again, we denote “MH”, “Ex”, “NMH” and “ABC”, respectively, for the M-H, exchange, noisy M-H and ABC-MCMC algorithms. 
While “Ex^app” and “NMH^app”, respectively, denote the approximate exchange and the approximate noisy M-H algorithms. 
Table <ref> illustrates that the mixing and posterior point estimates corresponding to the exchange, noisy M-H and ABC-MCMC algorithms are similar to Section <ref>.
However, the ESS(Ave)/s values of the “Ex” and “NMH” algorithms are generally worse than those of the ABC-MCMC cases.
This is due to the increased computational burden required by the auxiliary-draws and the likelihood function evaluation in the “NMH” (or “Ex”) algorithm leading to an increased computational runtime. 
This computational burden is even greater than that of the evaluation of likelihood normalizing constants in the “MH” algorithm.
Thus we can also observe in Table <ref> that the “MH” case provides even better efficiency, since the mixing of this algorithm is shown to be the best among all the candidates.
This is our motivation to further apply the approximate acceptance ratio (<ref>) in Algorithm <ref> 
to reduce the computational runtime of the noisy M-H algorithm, in order that the efficiency, as measured by the average ESS per second,  is significantly improved for the “Ex^app” and “NMH^app” cases in Table <ref>.


The significantly improved efficiency above 
is traded-off against the increase in the bias of the posterior point estimates of parameters resulting from the “Ex” and “NMH” algorithms, both of which are able to provide excellent agreement with the ground truth as shown in both Table <ref> and Figure <ref>. 
Though the “NMH^app” algorithm is able to bring better mixing performance for larger K implementations, the biases are shown to not significantly reduce and are larger than those of the “ABC” cases.
While 
the posterior density of the “Ex^app” and “NMH^app” cases is shown to be comparable with that of the “ABC” p=2.5 and p=1 cases in Figure <ref>.
Thus we do not lose much accuracy when applying the approximate algorithms 
while the efficiency is even better than that of the “MH” algorithm. 
However, following the red and blue entries in Table <ref> as well as the posterior density plots in Figure <ref>, 
one can argue that in general the “MH” algorithm provides the best balance of the accuracy, mixing and efficiency in this DPPG experiment.











§ REAL DATA APPLICATION



In this section, we apply the same experiments as our first simulation study shown in section <ref> for the purpose of comparison of different algorithms to fit the SPP model to the real dataset which was also applied by <cit.>.
The original data contains 13,655 tree locations with 68 species in the Blackwood region of Duke Forest. 
<cit.> processed this dataset by aggregating the species and removing trees which are under 40 dbh (diameter at breast height), considering the fact that the repulsion or inhibition can only be discovered by older trees. 
The left plot of Figure <ref> illustrates the processed dataset which contains 89 tree locations. 








The interaction radius R within the SPP model is again estimated by the profile pseudo-likelihood method and the estimated value of 0.053 is shown in the right plot of the Figure <ref>. 
The initial state of all the algorithms is the same as previously, that is, β_0 = 190 and γ_0 = 0.2. 
The same prior settings in <cit.> are also used here: β∼U(50, 350), γ∼U(0,1) with the bounded proposals being modified to have bounds agreeing with the prior settings. 
After the tuning process, ϵ_β and ϵ_γ are tuned to be 50 and 0.23, respectively. 
Each algorithm was run for 0.12 million iterations. 
Additionally a very long run of the exchange algorithm for 1.2-million iterations serves as our ground truth.











Once again, the labels “ABC”, “Ex” and “NMH” in Table <ref> represents the ABC-MCMC, exchange and noisy M-H algorithm cases, respectively. 
Our overall conclusions based on this experiment are broadly similar to the simulation study presented in Section <ref>. Generally, the performance of the various ABC algorithms is not as competitive as the exchange and noisy M-H algorithms, in terms of the posterior density point estimates and also in terms of the ESS and the ESS(Ave)/s summaries. While Figures <ref> and <ref> illustrate that the posterior density estimates based on the exchange and noisy M-H algorithms are generally closer to the ground truth than those of the ABC algorithms. 















§ CONCLUSION AND DISCUSSION



In this paper, we point out that <cit.> ABC-MCMC algorithm is shown to target an intractable posterior distribution and therefore impractical to implement generally. 
Thus we instead consider the <cit.> ABC-MCMC algorithm for the purpose of algorithm comparisons on repulsive point process models.
We compare this ABC-MCMC algorithm to the exchange algorithm and the recently proposed noisy Metropolis-Hastings algorithm applied on two models: the Strauss point process and the determinantal point process with a Gaussian kernel. 
The SPP is a doubly-intractable model for which perfect sampling from this process is possible. 
While perfect simulation of the DPPG is required to be restricted to a specific area and a Fourier series approximation is further needed in practice. 
The corresponding approximate restricted likelihood function is instead treated as the true likelihood function in our simulation studies. 

For the SPP model, our simulation study reveals that the noisy M-H algorithm is able to obtain almost the same accuracy as the exchange algorithm and at the same time provide better mixing performance even for the two-auxiliary-chain setting. 
In comparison, the ABC-MCMC cases are shown to provide worse performance in terms of both accuracy and mixing, and furthermore no apparent advantage can be observed when considering the resulting increased computational run time.
This indicates improved performance for the exchange and noisy M-H algorithms compared to the ABC-MCMC algorithm.
Moreover, the noisy M-H algorithm is shown to provide better mixing at the expense of increased run-time as the number of auxiliary chains increases. 
All of these characteristics can also be observed in the real data application. 

Regarding the DPPG model, the performance of the exchange, noisy M-H and ABC-MCMC algorithms is similar to that observed for the SPP model except that the efficiency of the exchange and noisy M-H algorithms become worse compared to the ABC-MCMC algorithm.
This is due to the huge computational burden required in the simulation process and in calculating the true likelihood function.
This motivated us to further apply the approximate exchange and approximate noisy M-H algorithms which resulted in improved efficiency compared to their ABC counterparts, without much deterioration in accuracy.
However, due to the tractability of the likelihood normalising constant with truncation approximation applied in practice, the Metropolis-Hastings algorithm is available 
and, in general, this algorithm is shown to obtain the best balance of accuracy, mixing and efficiency among all the candidate algorithms for the DPPG model.

Note that while we do not have theoretical guarantees for the approximate exchange and approximate noisy M-H algorithms, the empirical performance of these algorithms suggest that this is worthy of further study. 
Similar to the SPP cases, the mixing-efficiency trade-off can also be observed for the DPPG model by noticing that the efficiency of the (approximate) noisy M-H algorithm becomes worse as the number of auxiliary chains increases, but such an algorithm is able to provide similar accuracy performance and at the same time significantly better mixing performance for a small number of auxiliary chains compared to the (approximate) exchange algorithm.

Future work comparing the algorithms presented in this paper on other spatial point processes can be considered, for example, Diggle-Gratton point processes <cit.>, Diggle-Gates-Stibbard processes <cit.>, Penttinen processes <cit.> for which perfect simulations are all known to be available.  
Moreover, the approximate acceptance ratio for the DPPG model does not guarantee the convergence of the posterior samples, although our simulation study suggests that the results appeared to converge without much loss of accuracy.




§ ACKNOWLEDGEMENTS

The Insight Centre for Data Analytics is supported by Science Foundation Ireland under Grant Number 12/RC/2289_P2.
We would like to thank Andrew Golightly for helpful discussions.



§ CODE AND DATA

The application code and data for the simulation studies and the real data application is available at <https://github.com/Chaoyi-Lu/Bayesian_Strategies_for_RSPP>.



