This paper studies a cluster point process model defined as follows. Let be a simple locally finite point process defined on the -dimensional Euclidean space ; we can view as a random subset of (for background material on spatial point processes, see Møller and Waagepetersen (2004)). Assume is stationary, that is, its distribution is invariant under translations in . Conditioned on , let be a Poisson process on with intensity function
where and are parameters and
is a probability density function (pdf) on; in our specific models will play the role of a band width (a positive scale parameter). We can identify by a cluster process where conditioned on , the clusters are independent finite Poisson processes on and has intensity function (depending on the ‘offspring’ density relative to the cluster center ).
In the special case where is a stationary Poisson process, is a shot noise Cox process (SNCP), see Møller (2003). Then there may be a large amount of overlap between the clusters unless the intensity of is small as compared to the band width . In this paper, we will instead be interested in repulsive point process models for . This may be an advantage since the repulsiveness of implies less overlap of clusters. Thereby it may be easier to apply statistical methods for cluster detection, and when modelling clustered point pattern data sets a much lower intensity of may be needed as compared to the case of a SNCP. The idea of using a repulsive point process is not new, where Van Lieshout and Baddeley (2002) suggested to use a Markov point process. However, we are in particular interested in the case where is a stationary determinantal point process (DPP) in which case we call a determinantal shot noise Cox process (DSNCP). Briefly, a DPP is a model with repulsion at all scales, cf. Lavancier et al. (2015), Møller and O’Reilly (2021), and the references therein. There are several advantages of using a DPP for : In contrast to a Markov point process, there is no need of MCMC when simulating a DPP, and as we shall see and hence possess nice moment results, which can be used for estimation.
The cluster point process given by (1.1) is a special case of a stationary generalized shot noise Cox process (GSNCP), see Møller and Torrisi (2005), and it may be extended as follows. Suppose conditioned on both
and positive random variablesand is a Poisson process with intensity function
where is a pdf on . In addition, assume that , , and are mutually independent, the are independent identically distributed with mean
and has finite variance, and theare independent identically distributed. Then is still a stationary GSNCP and if also is a stationary DPP we may call a DGSNCP. In fact, all results and statistical methods used in this paper will apply for the DGSNCP when is replaced by in all expressions to follow. The DGSNCP may most naturally be treated in a MCMC Bayesian setting using a similar approach as in Beraha et al. (2022) and the references therein.
In the present paper, we study and exploit for statistical inference the nice moment properties for various DSNCP models as follows. In Section 2 we describe further what it means that is a DPP and present two specific cases of DSNCP models where we let be an isotropic Gaussian density as in the Thomas process (Thomas, 1949), which is the most popular example of a SNCP. Section 3 considers general results for so-called pair correlation and -functions for first the GSNCP model and second the DSNCP model. Section 4 discusses how the results in Section 3 may be used when fitting a parametric DSNCP model in a frequentist setting, and we illustrate this on a real data example.
2 Determinantal shot noise Cox process models
In this section we consider the DSNCP model for and suggest some specific models. In brief, the DPP is specified by a so-called kernel which is usually assumed to be a complex covariance function defined for all ; for details, see Appendix A. We assume is a stationary DPP with intensity , meaning two things: First, if is a bounded Borel set, denotes the cardinality of , and is the Lebesgue measure of , then . Second, for all where denotes the modulus of a complex number . We denote the corresponding complex correlation function by and assume it depends on a correlation/scale parameter so that with
The correlation parameter cannot vary independently of since there is a trade-off between intensity and repulsiveness in order to secure that a DPP model is well defined (Lavancier et al., 2015). For instance, for many DPP models is real, continuous, and stationary, that is, where is a continuous, symmetric, and positive semi-definite function with . Then, if
is square integrable and has Fourier transformwhere is the usual inner product, the DPP is only well-defined for , cf. Lavancier et al. (2015). In case of (2.1), this existence condition of the DPP means that , where for a fixed value of , most repulsiveness is obtained when .
Consider the special case where is a jinc-like DPP, that is, and where is the first order Bessel function of the first kind and denotes usual distance. So, the distribution of depends only on the intensity, and is a most repulsive DPP in the sense of Lavancier et al. (2015), see also Biscio and Lavancier (2016) and Møller and O’Reilly (2021). Christoffersen et al. (2021) used this special case of a DSNCP model in a situation where a realization of but not was observed within a bounded region. They estimated with a minimum contrast procedure based on the pair correlation function (pcf) given by (3.3) in Section 3, where the pcf had to be approximated by numerical methods. Instead we consider more tractable cases, as we shall see in Section 3.
Note that is stationary with intensity . We will consider two specific DSNCP models of where we let be the pdf of , the zero-mean isotropic
-dimensional normal distribution, andis given by one of the following two DPPs.
If is the Gaussian correlation function, then is a Gaussian DPP (Lavancier et al., 2015).
Both these DPP models are well defined if and only if (Lavancier et al., 2015). For a fixed value of , becomes in both cases less and less repulsive as decreases from to , where in the limit is a stationary Poisson process and is a Thomas process. Therefore, we call a Gaussian-DPP-Thomas process in the first case and a Ginibre-DPP-Thomas process in the second case. In both cases, and hence also are stationary and isotropic, although is only stationary when it is the Gaussian correlation function (Appendix B.1 verifies that the distribution of a scaled Ginibre point process is invariant under isometries).
The two first columns of plots in Figure 1 show simulated realizations of a Ginibre-DPP-Thomas process and a Gaussian-DPP-Thomas process within a square region, where , (so we expect to see about 400 points in each simulated point pattern), and in the three rows of plots we have (from bottom to top). In each case, is as large as possible. For comparison, the third column of plots in Figure 1 shows simulations of Thomas processes with the same values of as for the two first columns of plots. So, in each row the three processes have the same expected number of clusters, the same expected cluster sizes, and the same offspring density. For all processes, as increases (that is, decreases and increases) we see that the point patterns look more clustered, since we get less and less cluster centres but larger and larger clusters. We also see that the eye detects less diffuse clusters in the DPP-Thomas processes compared to the Thomas processes, which is in agreement with the fact that cluster centres are repulsive in DPP-Thomas processes whereas they are completely random in Thomas processes. From Figure 1 it can be difficult to make any conclusions about the differences between Gaussian- and Ginibre-DPP-Thomas processes, but we will make further comparisons between these in Sections 3–4.
3 Pair correlation and -functions
3.1 The general setting of stationary GSNCPs
Consider again the general setting in Section 1 of a stationary GSNCP. Henceforth we assume the stationary point process has pair correlation function (pcf) . This means that if are disjoint bounded Borel sets, then
By stationarity, depends only on the lag almost surely (with respect to Lebesgue measure) and for ease of presentation we can assume this is the case for all .
The stationary GSNCP given by (1.1) has intensity and a stationary pcf where
where denotes convolution and is the reflection of , cf. Møller and Torrisi (2005). Thus, decreases as increases; this makes good sense since is the intensity of clusters and the first term in (3.1) corresponds to pairs of points from different clusters whilst the second term is due to pairs of points within a cluster. Furthermore, from (3.1) we obtain Ripley’s -function (Ripley, 1976, 1977)
which will be used in Section 4 for parameter estimation.
In the special case where is a stationary Poisson process (i.e., is a SNCP), we have and (3.1) reduces to . Thus and which is usually interpreted as being a model for clustering. This is of course also the situation if . However, such models for clustering may cause a large amount of overlap between the clusters unless is small as compared to the band width .
3.2 The special setting of DSNCPs
When is a DPP and we let , we have
This is in accordance to intuition: As increases, meaning that decreases and hence that becomes more repulsive, it follows from (3.3) that decreases; and as the band width tends to 0, we see that tends to for every .
When is a Gaussian DPP, we have
and so from (3.3) we obtain that is isotropic with
since . Furthermore, if denotes the volume of the -dimensional unit ball and
is the CDF of a gamma distribution with shape parameterand scale parameter , we obtain
Consider now the case where is a scaled Ginibre point process. The scaled Ginibre point process has some similarity to the Gaussian-DPP, since has the same range in the two processes and
in the case of a scaled Ginibre point process. It thus follows from (3.2), (3.4), and (3.8) that the pcfs of the scaled Ginibre point process and the Gaussian-DPP are of the same form, but in the scaled Ginibre point process corresponds to in the Gaussian-DPP. This shows that the scaled Ginibre point process is more repulsive than the Gaussian-DPP when using the same parameters and , and therefore it will be possible to obtain a larger repulsion between the clusters in a Ginibre-DPP-Thomas process than in a Gaussian-DPP-Thomas process. In fact, if when is a scaled Ginibre point process, then is a most repulsive DPP in the sense of Lavancier et al. (2015). Because in the scaled Ginibre point process corresponds to in the Gaussian-DPP, (3.5)–(3.7) give that when is a scaled Ginibre point process, is isotropic with
Figure 2 shows plots of and for Gaussian- and Ginibre-DPP-Thomas processes for different values of when corresponds to the most repulsive case and without loss of generality we let . Note that in case of a planar stationary Poisson process, and the figure shows that as increases, the processes behave less and less like a planar stationary Poisson process. The pair correlation functions show an increasing degree of clustering at small scales and regularity at larger scales as increases, whereas Ripley’s -function only reveals an increasing degree of clustering. We furthermore see that the Ginibre-DPP-Thomas processes are overall more regular than the corresponding Gaussian-DPP-Thomas processes, especially at larger scales, which again reflects that the cluster centres are more regular in the Ginibre-DPP-Thomas processes.
4 Statistical inference
Suppose , is a bounded observation window, is a point pattern data set, and we want to fit a parametric DSNCP model given by either the Gaussian-DPP-Thomas or Ginibre-DPP-Thomas process. There is a trade-off between and because of the relation . Therefore, when modelling the data, we choose to let , which means that will be as repulsive as possible. That is, and are the unknown parameters.
Likelihood based inference is complicated because of the unobserved process of cluster centres.
Møller and Waagepetersen (2004) showed how a missing data MCMC approach can be used for maximum likelihood estimation in the special case of the Thomas process, and it may be simpler but still rather complicated to use a MCMC Bayesian setting along similar lines as in Beraha et al. (2022). We propose instead to exploit the parametric expressions of the intensity and of the pcf or -function given in Section 3.2 when estimating and .
Specifically, we use a minimum contrast procedure for estimating , where it is preferable to consider Ripley’s -function, since it is easier to estimate than by non-parametric methods, see e.g. Møller and Waagepetersen (2004). Writing to stress the dependence of and for a non-parametric estimate based on , the minimum contrast estimate of is given by
where we use the R-package spatstat (Baddeley et al., 2015) for calculating and the minimum contrast estimate by using default settings for the choice of and . Finally, having estimated , we estimate
from the unbiased estimation equation.
4.2 Model checking
When checking a fitted model, we prefer to use other functional summary statistics than since this was used as part of the estimation procedure. The standard is to consider empirical estimates of theoretical functions known as the empty space function (or spherical contact function) , the nearest-neighbour function , and the -function, which are defined for a stationary point process as follows. Consider any number and an arbitrary point . Then,
Here is only defined for , is the distance from to , and in the definition of when conditioning on it means that follows the reduced Palm distribution of at , see e.g. Møller and Waagepetersen (2004). Since is stationary, the definitions of and do not depend on the choice of .
We have not been able to derive the expressions of , , and for Gaussian- and Ginibre-DPP-Thomas processes; to the best of our knowledge, these expressions are not even known for a Thomas process. We refer to empirical estimates of these theoretical functions as functional summary statistics and use the relevant functions in spatstat to calculate such non-parametric estimates (always using the default settings including settings which account for boundary effects). Figure 3 concerns means of non-parametric estimates , , and calculated from simulations of Ginibre-DPP-Thomas processes and Gaussian-DPP-Thomas processes for the parameters stated in the caption. In agreement with Figure 2 the plots show an increasing degree of clustering as increases. The plots also indicate that the Gaussian-DPP-Thomas processes are more clustered and exhibit more empty space than the corresponding Ginibre-DPP-Thomas processes, and the difference becomes more apparent as decreases.
In order to validate a fitted model, we use 95% global envelopes and global envelope tests based on the extreme rank length as described in Myllymäki et al. (2017), Mrkvička et al. (2018), and Myllymäki and Mrkvička (2019), which is implemented in the R-package GET (Myllymäki and Mrkvička, 2019). These envelopes are based on functional summary statistics calculated from a number of simulations under the fitted model. We use simulations as recommended in the above references.
4.3 An application example
The first point pattern in Figure 4 shows the positions of 448 white oak trees in a square region (scaled to a unit square) of Lansing Woods, Clinton County, Michigan USA, which is part of the lansing data set which is available in spatstat. We will refer to this point pattern as .
By using the method of minimum contrast estimation as described in Section 4.1, we fitted a Gaussian-DPP-Thomas process, Ginibre-DPP-Thomas process, and Thomas process to . The obtained estimates are given in Table 1.
We see that the fitted DPP-Thomas processes expect much fewer clusters than the Thomas process and thus also more points in each cluster. As we expected, the fitted Ginibre-DPP-Thomas process is the one which expects the fewest clusters. Because of its expected 35 clusters with about 12 points on average in each it also seems to be a more sensible cluster process model than the other processes, which have many clusters with only a few points in each cluster. Figure 4 also shows a realization of each fitted model. The behaviour of these realisations are apparently in good agreement with . In order to check whether the models fit to data, we made 95% global envelope tests as described in Section 4.2. Figure 5 shows the results, which indicate that all three models fit well, but the Ginibre-DPP-Thomas process has a much higher -value than the other two processes.
Appendix A Definition of a DPP and some properties
Let be a simple point process defined on and be a complex function defined on so that for every integer and pairwise disjoint bounded Borel sets , we have
where is the determinant of the matrix with ’th entry . Then Macchi (1975) defined to be a DPP with kernel . Note that must be locally finite and the function
is the so-called ’th order intensity function of .
In fact for the DPP , its distribution is unique and completely characterized by the intensity functions of all order, cf. Lemma 4.2.6 in Hough et al. (2009). Thus stationarity of is equivalent to that for all and (Lebesgue almost) all , and isotropy of means that for all rotations matrices and (Lebesgue almost) all .
For later use, consider any numbers and , and the scaled point process . Let be an independent -thinning of (that is, the points in
are independently retained with probabilityand consists of those retained points). It is easily seen that is a DPP with kernel
Appendix B Some results for the scaled Ginibre point processes
In the following assume and is a standard Ginibre point process as defined in Section 2, so we identify with the complex plane . Let be as above and let be its intensity. By (A.2), is the DPP with kernel
and . In Section 2 we used the variation dependent parametrization , which is in one-to-one correspondence to . For the following it is convenient to let and use the variation independent parametrization , which is also in one-to-one correspondence to . Using this parametrization, with a slight abuse of notation we write for the DPP and
for its kernel.
b.1 Invariance under isometries
Below we show that the ’th order intensity function is invariant under translations and rotations, and therefore is stationary and isotropic. In the same way, it can be shown that is invariant under reflections and glide reflections. So the distribution of is invariant under isometries (mappings of the form and where with ; these mappings correspond to translations, rotations, reflections, and glide reflections).
b.2 Spectral decompositions
Spectral representations of the kernel restricted to compact regions are needed for simulation as well as other purposes, cf. Lavancier et al. (2015). The simplest case occurs when we consider restricted to a closed disc around zero. So for , let be the closed disk around zero with radius and the restriction of to . Because is a DPP, is a DPP with kernel
The integral operator corresponding to the kernel
has only one eigenvalue, namely
, and the eigenfunctions are
This follows easily by exploiting the moment properties of two independent zero-mean complex normally distributed random variables and the definition of the complex exponential function ( for ) for the term in (B.1). In other words, we have the spectral representation
Similarly, we see that the integral operator corresponding to has eigenfunctions
with corresponding eigenvalues
and the spectral representation is
Appendix C Simulation procedures
For simulating determinantal point processes, we use the algorithm described in Lavancier et al. (2015) which is a specific case of the simulation algorithm of Hough et al. (2006). We refer to these references for specific details. The algorithm is implemented in spatstat for the models suggested in Lavancier et al. (2015), which include Gaussian DPPs. For these models, it is necessary to approximate the kernel because the spectral representation is unknown. In the case of a scaled Ginibre point process, this approximation is however unnecessary for simulating it on a disc because the spectral representation is known. The simulation is still only approximate because the procedure also involves other approximations including approximating the upper bound for rejection sampling chosen in Lavancier et al. (2015) (an approximation which is in fact not necessary for the models they consider since the expression simplifies in those cases). For simulating a scaled Ginibre point process on a window , we thus use the spectral representation on a disc to simulate the process on and thereafter extract the part which is in .
For simulating DPP-Thomas processes on a window , we first simulate the DPP obtained by restricting to an extended window in order to account for boundary effects. Regarding the extension, we decided to use the default setting from the function rThomas in spatstat which simulates a Thomas process. Given the cluster centers on the extended window, we simulate the clusters of the DPP-Thomas process independently as finite Poisson processes with intensity functions for each . That is, first simulate the number of points in a cluster centered at
from a Poisson distribution with rate. Then sample the independent points in from the -dimensional normal distribution . Finally, the simulation of on is the part of which falls in .
- Spatial point patterns: methodology and applications with r. Chapman & Hall/CRC Press, Boca Raton. Cited by: §4.1.
- MCMC computations for Bayesian mixture models using repulsive point processes. Journal of Computational and Graphical Statistics. Note: To appear. Available at arXiv:2011.06444 Cited by: §1, §4.1.
- Quantifying repulsiveness of determinantal point processes. Bernoulli 22, pp. 2001–2028. Cited by: §2.
- Modelling columnarity of pyramidal cells in the human cerebral cortex. Australian and New Zealand Journal of Statistics 63, pp. 33–54. Cited by: §2.
- The Ginibre point process as a model for wireless networks with repulsion. IEEE Transactions on Wireless Communications 14, pp. 107–121. Cited by: item 2.
- Statistical ensembles of complex, quaternion, and real matrices. Journal of Mathematical Physics 6, pp. 440–449. Cited by: item 2.
- A composite likelihood approach in fitting spatial point process models. Journal of American Statistical Association 101, pp. 1502–1512. Cited by: §4.1.
- Determinantal processes and independence. Probability surveys 3, pp. 206–229. Cited by: Appendix C.
- Zeros of gaussian analytic functions and determinantal point processes. American Mathematical Society, Providence. Cited by: Appendix A.
- Determinantal point process models and statistical inference. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77, pp. 853–877. Cited by: §B.2, Appendix C, §1, item 1, §2, §2, §2, §3.2, §3.2.
- The coincidence approach to stochastic point processes. Advances in Applied Probability 7, pp. 83–122. Cited by: Appendix A.
- Spatial modeling and analysis of cellular networks using the Ginibre point process: A tutorial.. IEICE Transactions on Communic- ations 99, pp. 2247–2255. Cited by: item 2.
- Couplings for determinantal point processes and their reduced Palm distributions with a view to quantifying repulsiveness. Journal of Applied Probability 58, pp. 469–483. Cited by: §1, §2.
- Generalised shot noise Cox processes. Advances in Applied Probability 37, pp. 48–74. Cited by: §1, §3.1.
- Some recent developments in statistics for spatial point patterns. Annual Review of Statistics and Its Applications 4, pp. 317–342. Cited by: §4.1.
- Shot noise Cox processes. Advances in Applied Probability 35, pp. 614–640. Note: Cited by: §1.
- Statistical inference and simulation for spatial point processes. Chapman & Hall/CRC, Boca Raton, Florida. Cited by: §1, §4.1, §4.1, §4.2.
- A one-way ANOVA test for functional data with graphical interpretation. Note: Available at arXiv:1612.03608 External Links: Cited by: §4.2.
- Global envelope tests for spatial processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79, pp. 381–404. Cited by: §4.2.
- GET: global envelopes in R. Note: arXiv preprint arXiv:1911.06583 Cited by: §4.2.
- R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Cited by: §1.
- The second-order analysis of stationary point processes. Journal of Applied Probability 13, pp. 255–266. Cited by: §3.1.
- Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology) 39, pp. 172–212. Cited by: §3.1.
- A generalization of Poisson’s binomial limit in the use in ecology. Biometrika 36, pp. 18–25. Cited by: §1.
Extrapolating and interpolating spatial patterns. In Spatial Cluster Modelling, A. B. Lawson and D. Denison (Eds.), pp. 61–86. Cited by: §1.
- Ggplot2: elegant graphics for data analysis. Springer-Verlag New York. External Links: Cited by: §1.