Log Gaussian Cox processes on the sphere

by   Jesper Møller, et al.
Aalborg University

We define and study the existence of log Gaussian Cox processes (LGCPs) for the description of inhomogeneous and aggregated/clustered point patterns on the d-dimensional sphere, with d = 2 of primary interest. Useful theoretical properties of LGCPs are studied and applied for the description of sky positions of galaxies, in comparison with previous analysis using a Thomas process. We focus on simple estimation procedures and model checking based on functional summary statistics and the global envelope test.


Structured space-sphere point processes and K-functions

This paper concerns space-sphere point processes, that is, point process...

Testing for complete spatial randomness on three dimensional bounded convex shapes

There is currently a gap in theory for point patterns that lie on the su...

SPHARMA approximations for stationary functional time series on the sphere

In this paper, we focus on isotropic and stationary sphere-cross-time ra...

An anisotropic model for global climate data

We present a new, elementary way to obtain axially symmetric Gaussian pr...

Local inhomogeneous weighted summary statistics for marked point processes

We introduce a family of local inhomogeneous mark-weighted summary stati...

Spatial log-Gaussian Cox processes in Hilbert spaces

A new class of spatial log-Gaussian Cox processes in function spaces is ...

Should we condition on the number of points when modelling spatial point patterns?

We discuss the practice of directly or indirectly assuming a model for t...

1 Introduction

Statistical analysis of point patterns on the sphere have been of interest for a long time, see Lawrence et al. (2016), Møller and Rubak (2016), and the references therein. Although models and methods developed for planar and spatial point processes may be adapted, statistical methodology for point processes on the sphere is still not so developed as discussed in the mentioned references and in the following.

The focus has been on developing functional summary statistics, parametric models, and inference procedures. For homogeneous point patterns,

Robeson et al. (2014) studied Ripley’s -function on the sphere (Ripley, 1976, 1977), and Lawrence et al. (2016) provided a careful presentation of Ripley’s -function and other functional summary statistics such as the empty space function , the nearest-neighbour distance function , and the function, including how to account for edge effects (i.e. when the data is a point pattern observed within a window strictly included in the sphere and unobserved points outside this window may effect the data). See also Møller and Rubak (2016) for details and the connection to reduced Palm distributions. For inhomogeneous point patterns, Lawrence et al. (2016) and Møller and Rubak (2016) studied the pair correlation function and the inhomogeneous -function. The models which have been detailed are rather scarse: Homogeneous Poisson point process models (Raskin, 1994; Robeson et al., 2014); determinantal point processes (Møller and Rubak, 2016; Møller et al., 2018) for regular/repulsive point patterns; inhomogeneous Poisson point process models and Thomas point process models for aggregated/clustered point patterns (Lawrence et al., 2016; Section 2.2-2.3 in the present paper); and inhomogeneous log Gaussian Cox processes (LGCPs) for aggregated/clustered point patterns using the R-INLA approach (Simpson et al., 2016; Section 4.2 in the present paper).

This paper treats inhomogeneous aggregated/clustered point patterns on the sphere and studies how the theory of Cox processes and in particular LGCPs (Møller et al., 1998) can be adapted to analysing such data, using other LGCPs models than in Simpson et al. (2016)

and simpler inference procedures than the R-INLA approach, where we exploit the nice moment properties of LGCPs. In particular, we demonstrate that an inhomogeneous LGCP provides a better description of the sky positions of galaxies than analysed in

Lawrence et al. (2016) by using an inhomogeneous Thomas process. For comparison, as in Lawrence et al. (2016), we use a minimum contrast procedure for parameter estimation, where we discuss the sensibility of the choice of user-specified parameters. No model checking was done for the fitted inhomogeneous Thomas process in Lawrence et al. (2016). Moreover, we show how a thinning procedure can be applied to generate homogeneous point patterns so that the -functions can be used for model checking, and we discuss the sensitivity of this thinning procedure.

The paper is organised as follows. Section 2 provides the setting and needed background material on point processes, particularly on Poisson and Cox processes, the data example of sky positions of galaxies, and the inhomogeneous Thomas process introduced in Lawrence et al. (2016). Section 3 contains the definition and existence conditions for LGCPs on the sphere, studies their useful properties, and compares the fitted Thomas processes and LGCPs for the data example. Finally, Section 4 summarizes our results, establishes some further results, and discusses future directions for research.

2 Background

2.1 Setting

Let denote the -dimensional unit sphere included in the -dimensional Euclidean space , equipped with the usual inner product for points and the usual length . We are mainly interested in the case of . By a region we mean that

is a Borel set. For unit vectors

, let be the geodesic distance on the sphere.

By a point process on , we understand a random finite subset X of . We say that X is isotropic if is distributed as X for any rotation matrix . We assume that X has an intensity function, , and a pair correlation function, , meaning that if are disjoint regions and denotes the cardinality of , then

where is the Lebesgue/surface measure on . We say that X is (first order) homogeneous if is constant, and second order intensity reweighted homogeneous if only depends on . Note that these properties are implied by isotropy of X, and second order intensity reweighted homogeneity allows to define the (inhomogeneous) -function (Baddeley et al., 2000) by

for an arbitrary chosen point , where denotes Lebesgue/surface measure on . For instance,


2.2 Poisson and Cox processes

For a Poisson process on with intensity function ,

is Poisson distributed with mean

, and conditional on , the points in are independent identically distributed with a density proportional to . The process is second order intensity reweighted homogeneous, with -function


Let be a non-negative random field so that almost surely is well-defined and finite, and

has finite variance for all

. We say that is a Cox process driven by if conditional on is a Poisson process on with intensity function . Then has intensity function


and defining the residual driving random field by , where (setting ), has pair correlation function


Thus is second order intensity reweighted homogeneous if is isotropic, that is, is distributed as for any rotation matrix . We shall naturally specify such Cox processes in terms of and .

2.3 Data example and inhomogeneous Thomas process

Figure 0(a) shows the sky positions of 10,546 galaxies from the Revised New General Catalogue and Index Catalogue (RNGC/IC) (Steinicke, 2015). Here, we are following Lawrence et al. (2016) in making a rotation of the original data so that the two circles limit a band around the equator, which is an approximation of the part of the sky obscured by the Milky Way, and the observation window is the complement of the band; 64 galaxies contained in the band are omitted.

(a) Sky positions of galaxies before thinning.
(b) Sky positions of galaxies after thinning.
Figure 1: The sky positions of (a) 10,546 galaxies and (b) 3,285 galaxies obtained after a thinning procedure so that a homogeneous point pattern is expected to be obtained. The circles limit the part of the sky obscured by the Milky Way.

Different plots and tests in the accompanying supporting information to Lawrence et al. (2016) clearly show that the galaxies are aggregated and not well-described by an inhomogeneous Poisson process model. Lawrence et al. (2016) fitted the intensity function


where , is the colatitude, and is the longitude. The term is the inner product of with , and the term allows a gradient in intensity from the poles to the equator, cf. the discussion in Lawrence et al. (2016) and their plot Movie 2 (where we recommend using the electronic version of the paper and the link in the caption to this plot). Throughout this paper, we also use (5).

Lawrence et al. (2016) proposed an inhomogeneous Thomas process, that is, a Cox process with intensity function (5) and isotropic residual driving random field given by


where is a homogeneous Poisson process on with intensity , and where

is the density for a von Mises-Fisher distribution on with mean direction and concentration parameter . They estimated and by a minimum contrast procedure (Diggle and Gratton, 1984; Diggle, 2013), where a non-parametric estimate is compared to

Specifically, the minimum contrast estimate is minimising the contrast


where are user-specified parameters. Lawrence et al. (2016) used the integration interval , where 1.396 is the maximal value of the smallest distance from the north pole respective south pole to the boundary of . They obtained and .

3 Log Gaussian Cox processes

3.1 Definition and existence

Let be a Gaussian random field (GRF), that is,

is normally distributed for any integer

, numbers , and . Let be its mean function and its covariance function. Assuming that almost surely is integrable, the Cox process driven by is called a log Gaussian Cox process (LGCP; Møller et al., 1998; Møller and Waagepetersen, 2004).

Define the corresponding mean-zero GRF, , which has also covariance function . Assuming is a Borel function with an upper bounded, then almost sure continuity of implies almost sure integrability of , and so is well-defined. In turn, almost sure continuity of is implied if almost surely is locally sample Hölder continuous of some order

, which means the following: With probability one, for every

, there is a neighbourhood of and a constant so that


The following proposition is a special case of Lang et al. (2016, Corollary 4.5), and it provides a simple condition ensuring (8).

Proposition 0.


be the variogram of a mean-zero GRF . Suppose there exist numbers and such that


Then, for any , is almost surely locally sample Hölder continuous of order .

3.2 Properties

The following proposition is easily verified along similar lines as Proposition 5.4 in Møller and Waagepetersen (2004) using (3)-(4

) and the expression for the Laplace transform of a normally distributed random variable.

Proposition 0.

A LGCP has intensity function and pair correlation functions given by


where and are the mean and covariance functions of the underlying GRF.

In other words, the distribution of is specified by because and specify the distribution of the GRF.

By (10), second order intensity reweighted homogeneity of the LGCP is equivalent to isotropy of the covariance function, that is, depends only on . Parametric classes of isotropic covariance functions on are provided by Gneiting (2013), see Table 2 in Section 4.1, and by Møller et al. (2018). In Section 3.3, we consider the so-called multiquadric covariance function which is isotropic and given by


where .

Proposition 0.

A mean-zero GRF with multiquadric covariance function is locally sample Hölder continuous of any order .


By (11), for ,


Here, letting ,


where the first inequality follows from because , and the second inequality follows from whenever . If , then

where the first inequality follows from (12) and the second from (13). If , then for any ,

where the first inequality follows from (12)-(13) and the second from whenever . Hence, for any and any , (11) satisfies (9) with , , and , where . ∎

In fact locally sample Hölder continuity is satisfied when considering other commonly used parametric classes of covariance functions, but the proof will be case specific as shown in the proof of Proposition 4 in Section 4.1.

3.3 Comparison of fitted Thomas processes and LGCPs for the data example

In this section, to see how well the estimated Thomas process, , obtained in Lawrence et al. (2016) fits the sky positions of galaxies discussed in Section 2.3, we use methods not involving the -function (or the related pair correlation function, cf. (1)) because it was used in the estimation procedure.

The point pattern in Figure 0(b) was obtained by an independent thinning procedure of the point pattern in Figure 0(a), with retention probability at location , where is given by (5) and . If we imagine a realization of on could had been observed so that the independent thinning procedure could had taken place on the whole sphere, then the corresponding thinned Thomas process, , is isotropic. Thus the commonly used -functions would apply: Briefly, for a given isotropic point process on , an arbitrary chosen location , and , by definition is the probability that has a point in a cap with center and geodesic distance from to the boundary of the cap; is the conditional probability that has a point in given that ; and when . For the empirical estimates of the -functions, as we only observe points within the subset specified in Section 2.3, edge correction factors will be needed, where we choose to use the border-corrected estimates from Lawrence et al. (2016).

Empirical estimates can then be used as test functions for the global rank envelope test, which is supplied with graphical representations of global envelopes for as described in Myllymäki et al. (2017). As we combined all three test functions as discussed in Mrkvic̆ka et al. (2016) and Mrkvic̆ka et al. (2017), we followed their recommendation of using simulations of for the calculation of -values and envelopes. Briefly, under the claimed model, with probability 95% we expect (the solid lines in Figure 1(a)1(c)) to be within the envelope (either the dotted or dashed lines in Figure 1(a)1(c) depending on the choice of as detailed below), and the envelope corresponds to a so-called global rank envelope test at level 5%. In Table 1, the limits of the -intervals correspond to liberal and conservative versions of the global rank envelope test (Myllymäki et al., 2017). Table 1 and Figure 1(a)1(c) clearly show that the estimated Thomas process is not providing a satisfactory fit no matter if in the contrast (7) the long integration interval from Lawrence et al. (2016) or the shorter interval (corresponding to 0-10 degrees) is used for parameter estimation. When , larger estimates and are obtained as compared to and from Lawrence et al. (2016). Note that the -values and envelopes are not much affected by the choice of integration interval in the minimum contrast estimation procedure, cf. Table 1 and Figure 1(a)1(c). The envelopes indicate that at short inter-point distances , there is more aggregation in the data than expected under the two fitted Thomas processes.

 Thomas process 0.01% - 1.05% 0.01% - 1.28%
 LGCP 24.02% - 24.09% 0.01% - 1.23%
Table 1: Intervals for -values obtained from the global envelope test based on combining the -functions and using either a short or long integration interval when calculating the contrast used for parameter estimation in the Thomas process or the LGCP.
(a) -function (Thomas).
(b) -function (Thomas).
(c) -function (Thomas).
(d) -function (LGCP).
(e) -function (LGCP).
(f) -function (LGCP).
Figure 2: Empirical functional summary statistics and 95% global envelopes under (a)–(c) fitted Thomas process and (d)–(f) fitted LGCPs for the sky positions of galaxies after thinning. The solid lines show versus distance (measured in radians). The dotted and dashed lines limit the global envelopes and correspond to the short and long integration intervals and , respectively, used in the minimum contrast estimation procedure.

We also fitted an inhomogeneous LGCP, , still with intensity function given by (5) and with multiquadric covariance function as in (11). The same minimum contrast procedure as above was used except of course that in the contrast given by (7), the theoretical -function was obtained by combining (1), (10), and (11), where numerical calculation of the integral in (1) was used. The estimates are if is the integration interval, and if . For each choice of integration interval, Figure 3 shows the estimated log pair correlation function, that is, the estimated covariance function of the underlying GRF, cf. (10). Comparing the two pair correlation functions, the one based on the short integration interval is much larger for very short inter-point distances , then rather similar to the other at a short interval of -values, and afterwards again larger, so the fitted LGCP using the short integration interval is more aggregated than the other case. Further, Figure 4 shows the empirical -function and the fitted -functions for both the Thomas process and the LGCP when using the different integration intervals (for ease of comparison, we have subtracted , the theoretical -function for a Poisson process, cf. (2)). The fitted -functions are far away from the empirical -function for large values of . However, having a good agreement for small and modest values of is more important, because the variance of the empirical -function seems to be an increasing function of and the interpretation of the -function becomes harder for large -values. For small and modest -values, using the LGCP model and the short integration interval provides the best agreement between the empirical -function and the theoretical -function, even when is outside the short integration interval.

Figure 3: Plot of estimated pair correlation functions (on a logarithmic scale) versus for the LGCP when different integration intervals were used in the minumum contrast estimation procedure (dotted line: ; dashed line: ).
(a) Integration interval [0,0.175].
(b) Integration interval [0,1.396].
Figure 4: Fitted -functions minus the theoretical Poisson -function versus distance for the sky positions of galaxies when using different integration intervals in the minumum contrast estimation procedure. The solid lines correspond to the empirical functional summary statistics, and the dashed and dotted lines correspond to the theoretical functional summary statistics under the fitted Thomas processes and LGCPs, respectively.

Table 1 and Figure 1(d)1(f) show the results for the fitted LGCPs when making a model checking along similar lines as for the fitted Thomas processes (i.e., based on a thinned LGCP and considering the empirical -functions together with global envelopes). Table 1 shows that the fitted LGCP based on the long integration interval gives an interval of similar low -values as for the fitted Thomas process in Lawrence et al. (2016), and Figure 1(e)1(f) indicate that the data is more aggregated than expected under this fitted LGCP. Finally, the fitted LGCP based on the short integration interval gives -values of about 24%, and the empirical curves of the functional summary statistics appear in the center of the envelopes, cf. Figure 1(d)1(f).

(a) -function (Thomas).
(b) -function (Thomas).
(c) -function (Thomas).
(d) -function (LGCP).
(e) -function (LGCP).
(f) -function (LGCP).
Figure 5: Empirical functional summary statistics and 95% global envelopes under (a)–(c) fitted Thomas process and (d)–(f) fitted LGCPs for the sky positions of galaxies after thinning. The 1000 solid lines show versus distance (measured in radians) when the independent thinning procedure is repeated 1000 times. The dotted and dashed lines limit the global envelopes and correspond to the short and long integration intervals and , respectively, used in the minimum contrast estimation procedure.
Figure 6: Intervals for -values obtained from the global envelope test based on combining the -functions when repeating the independent thinning procedure 1000 times. The three lower solid lines are very close and therefore appear as one thick solid line in the plot. Each pair of dashed and solid lines from below to the top corresponds to conservative and liberal -values for the fitted Thomas process, using first the long and second the short integration interval, and for the fitted LGCP, using first the long and second the short integration interval.

Finally, we have studied how sensible the results will be when using the independent thinning procedure. Figure 5 is similar to Figure 2 except that the 1000 solid lines show the empirical estimates when we repeat the thinning procedure 1000 times. The empirical estimates do not vary much in the 1000 cases. Moreover, Figure 6 is similar to Table 1 and shows the intervals for -values obtained in the 1000 cases. The intervals for the fitted Thomas process, using the short or long integration interval, and for the fitted LGCP, using the long integration interval, are effectively not varying in the 1000 cases (the three lower solid curves correspond to liberal -values and visually they appear as one curve which is close to 0: the curve in the LGCP case varies around , and the two curves in the Thomas process case vary around ). The intervals for the fitted LGCP, using the short integration interval, are much shorter than for the three other fitted models, and they vary from 0.75% - 0.77% to 41.44% - 41.48%, with only 12 out of the 1000 intervals below the 5% level. Thus, the results are not so sensitive to the independent thinning procedure and the conclusion, namely that the fitted LGCP, using the short integration interval, is fitting well whilst the other fitted models do not, is almost the same in the 1000 cases.

4 Further results and concluding remarks

4.1 Existence

As noticed in Section 3.1, under mild conditions for the mean function , almost sure locally sample Hölder continuity of the mean-zero GRF is a sufficient condition for establishing the existence of a LGCP driven by . In turn, almost sure locally sample Hölder continuity of is implied by the condition (9) in Proposition 1, and Proposition 3 states that the multiquadric covariance function is satisfying this condition. As shown in the following proposition, (9) is satisfied for any of the parametric classes of isotropic covariance functions on given by Gneiting (2013), see Table 2. In brief, these are the commonly used isotropic covariance functions which are expressible on closed form.

Model Correlation function Parameter range
Powered exponential ,
Matérn ,
Generalized Cauchy , ,
Dagum ,
multiquadric ,
Sine power
Askey ,
-Wendland ,
-Wendland ,
Table 2: Parametric models for an isotropic correlation function on , where , is the gamma function, is the modified Bessel function of the second kind, and for . For the powered exponential, Matérn, generalized Cauchy, Dagum, multiquadric, and sine power models, , whilst for the spherical, Askey, -Wendland, and -Wendland models, . For each model, the specified parameter range ensures that is well-defined, cf. Gneiting (2013), and hence for any , is an isotropic covariance function.
Proposition 0.

Any isotropic covariance function as given by Table 2 is satisfying (9) for some and .


We have already verified this for the multiquadric model, cf. Proposition 3. For any of models in Table 2, the variance is strictly positive, so dividing the left and right side in (9) by , we can without loss of generality assume the covariance function to be a correlation function, i.e. , cf. the caption to Table 2. Hence we need only to show the existence of numbers , and so that the condition


is satisfied, where .

Consider the powered exponential model. If and , then , and so we obtain (14) because

where the first inequality follows from . Further, the case is the special case of the Matérn model with shape parameter , which is considered below. For the remaining eight models in Table 2, we establish a limit


for some and , where depends on the specific model. The limit (15) implies that (14) holds for some , with , and . Note that .

For example, consider the Matérn model. By DLMF (2018, Equation 10.27.4), when ,


Thus, defining , the expression of the Matérn correlation function in Table 2 can be rewritten as


and so (15) holds, with and .

For the remaining seven models, the values of and in (15) are straightforwardly derived, using L’Hospital’s rule when considering the generalized Cauchy, Askey, -Wendland, and -Wendland models. The values are given in Table 2.

Generalized Cauchy
Sine power
Spherical 1
Askey 1
-Wendland 2
-Wendland 2
Table 3: List of and values for the last seven models considered in the proof of Proposition 4.

4.2 Moment properties, Palm distribution, and statistical inference

As demonstrated, a LGCP on the sphere possesses useful theoretical properties, in particular moment properties as provided by Proposition 2. We exploited these expressions for the intensity and pair correlation function when dealing with the inhomogeneous -function, which concerns the second moment properties of a second order intensity reweighted homogeneous point process.

Proposition 2 extends as follows. For a general point process on , the

-th order pair correlation function

is defined for integers and multiple disjoint regions by

provided this multiple integral is well-defined and finite. The following proposition follows immediately as in Møller and Waagepetersen (2003, Theorem 1).

Proposition 0.

For any integer , a LGCP on has -th order pair correlation function given by


where is the covariance function of the underlying GRF.

Higher-order pair correlation functions as given by (16) may be used for constructing further functional summaries e.g. along similar lines as the third-order characteristic studied in Møller et al. (1998) and the inhomogeneous -function studied in Cronie and van Lieshout (2015, Section 5.3).

It would also be interesting to exploit the following result for the reduced Palm distribution of a LGCP on the sphere as discussed in Coeurjolly et al. (2017) in the case of LGCPs on . Intuitively, the reduced Palm distribution of a point process on at a given point corresponds to the distribution of conditional on that ; see Lawrence et al. (2016) and Møller and Rubak (2016). Along similar lines as in Coeurjolly et al. (2017, Theorem 1), we obtain immediately the following result.

Proposition 0.

Consider a LGCP whose underlying GRF has mean and covariance functions and . For any , the reduced Palm distribution of at is a LGCP, where the underlying GRF has unchanged covariance function but mean function for .

Bayesian analysis of inhomogeneous LGCPs on the sphere has previously been considered in Simpson et al. (2016)

using the R-INLA approach: Their covariance model for the underlying GRF is an extension of the Matérn covariance function defined as the covariance function at any fixed time to the solution to a stochastic partial differential equation which is stationary in time and isotropic on the sphere. Using the R-INLA approach, a finite approximation of their LGCP based on a triangulation of the sphere is used, which for computational efficiency requires

to be an integer, where is the shape parameter (the expression in their paper is a typo). Simpson et al. (2016) considered an inhomogeneous LGCP defined only on the world’s oceans, with and , so . In the original Matérn model (Gneiting, 2013, Table 1), and the covariance function has a nice expression in term of a Bessel function, but otherwise the covariance function can only be expressed by an infinite series in terms of spherical harmonics (Simpson, 2009, Section 6.4). In contrast, we have focused on covariance functions which are expressible on closed form so that we can easily work with the pair correlation function . Indeed it could be interesting to apply the R-INLA approach as well as other likelihood based methods of inference (Møller and Waagepetersen, 2017), though they are certainly much more complicated than the simple inference procedure considered in the present paper.

For comparison with the analysis of the sky positions of galaxies in Lawrence et al. (2016), we used a minimum contrast estimation procedure, but other simple and fast methods such as composite likelihood (Guan, 2006; Møller and Waagepetersen, 2017, and the references therein)

could have been used as well. It is well-known that such estimation procedures can be sensitive to the choice of user-specified parameters. We demonstrated this for the choice of integration interval in the contrast, where using a short interval and an inhomogeneous LGCP provided a satisfactory fit, in contrast to using an inhomogeneous Thomas process or a long interval. It could be interesting to use more advanced estimation procedures such as maximum likelihood and Bayesian inference. This will involve a time-consuming missing data simulation-based approach

(Møller and Waagepetersen, 2004, 2017, and the references therein).

4.3 Modelling the sky positions of galaxies

The Thomas process is a mechanistic model since it has an interpretation as a cluster point process (Lawrence et al., 2016). The original Thomas process in (i.e., using a 3-dimensional isotropic zero-mean normal distribution as the density for a point relative to its cluster centre) may perhaps appear natural for positions of galaxies in the 3-dimensional space – but we question if the Thomas process from Lawrence et al. (2016) is a natural model for the sky positions because these points are obtained by projecting clusters of galaxies in space to a sphere which may not produce a clear clustering because of overlap. Rather, we think both the inhomogeneous Thomas and the inhomogeneous LGCP should be viewed as empirical models for the data. Moreover, it could be investigated if the Thomas process replaced by another type of Neyman-Scott process (Neyman and Scott, 1958, 1972) or a (generalised) shot noise Cox process (Møller, 2003; Møller and Torrisi, 2005) would provide an adequate fit.

As a specific model, we only considered the multiquadric covariance function for the underlying GRF of the LGCP. A more flexible model could be the spectral model studied in Møller et al. (2018)

, where a parametric model for the eigenvalues of the spectral representation of an isotropic covariance function

() in terms of spherical harmonics is used (incidentally, the eigenvalues are also known for the multiquadric covariance function if ). The spectral representation allows a Karhunen-Loève representation which could be used for simulation. However, for the data analysed in this paper, we found it easier and faster just to approximate the GRF , using a finite grid so that each is replaced by if is the nearest grid point to

, and using the singular value decomposition when simulating the finite random field

. When using spherical angles as in (5), a regular grid over can not be recommended, since the density of grid points will large close to the poles and small close to equator. Using a regular grid avoids this problem, nonetheless, there are only five regular grids on the sphere (Coxeter, 1973). We used a nearly-regular grid consisting of 4098 points on the sphere (Szalay et al., 2007, and the references therein) and computed using the R package mvmesh (Nolan, 2016).


Supported by The Danish Council for Independent Research Natural Sciences, grant DFF – 7014-00074 ”Statistics for point processes in space and beyond”, and by the ”Centre for Stochastic Geometry and Advanced Bioimaging”, funded by grant 8721 from the Villum Foundation.



  • Baddeley et al. [2000] A. J. Baddeley, J. Møller, and R. Waagepetersen. Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54:329–350, 2000.
  • Coeurjolly et al. [2017] J.-F. Coeurjolly, J. Møller, and R. Waagepetersen. Palm distributions for log Gaussian Cox processes. Scandinavian Journal of Statistics, 44:192–203, 2017.
  • Coxeter [1973] H. S. M. Coxeter. Regular Polytopes. Methuen, London, 1973.
  • Cronie and van Lieshout [2015] O. Cronie and M. N. M. van Lieshout. A -function for inhomogeneous spatio-temporal point processes. Scandinavian Journal of Statistics, 42:562–579, 2015.
  • Diggle [2013] P. J. Diggle. Statistical Analysis of Spatial and Spatio-temporal Point Patterns. CRC Press, Boca Raton, 2013.
  • Diggle and Gratton [1984] P. J. Diggle and R. J. Gratton. Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 46:193–227, 1984.
  • [7] DLMF. NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.18 of 2018-03-27. URL http://dlmf.nist.gov/. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
  • Gneiting [2013] T. Gneiting. Strictly and non-strictly positive definite functions on spheres. Bernoulli, 19:1327–1349, 2013.
  • Guan [2006] Y. Guan. A composite likelihood approach in fitting spatial point process models. Journal of the American Statistical Association, 101:1502–1512, 2006.
  • Lang et al. [2016] A. Lang, J. Potthoff, M. Schlather, and D. Schwab. Continuity of random fields on Riemannian manifolds. Manuscript available at arXiv: 1607.05859, 2016.
  • Lawrence et al. [2016] T. Lawrence, A. Baddeley, R. K. Milne, and G. Nair. Point pattern analysis on a region of a sphere. Stat, 5:144–157, 2016.
  • Møller [2003] J. Møller. Shot noise Cox processes. Advances in Applied Probability, 35:614–640, 2003.
  • Møller and Rubak [2016] J. Møller and E. Rubak. Functional summary statistics on the sphere with an application to determinantal point processes. Spatial Statistics, 18:4–23, 2016.
  • Møller and Torrisi [2005] J. Møller and G. Torrisi. Generalised shot noise Cox processes. Advances in Applied Probability, 37:48–74, 2005.
  • Møller and Waagepetersen [2017] J. Møller and R. Waagepetersen. Some recent developments in statistics for spatial point patterns. Annual Review of Statistics and Its Applications, 4:317–342, 2017.
  • Møller and Waagepetersen [2003] J. Møller and R. P. Waagepetersen. Spatial Statistics and Computational Methods, chapter An introduction to simulation-based inference for spatial point processes, pages 37–60. Lecture Notes in Statistics 173, Springer-Verlag, New York, 2003.
  • Møller and Waagepetersen [2004] J. Møller and R. P. Waagepetersen. Statistical Inference and Simulation for Spatial Point Processes. Boca Raton, FL, Chapman & Hall/CRC., 2004.
  • Møller et al. [1998] J. Møller, R. Waagepetersen, and A. R. Syversveen. Log Gaussian Cox processes. Scandinavian Journal of Statistics, 25:451–482, 1998.
  • Møller et al. [2018] J. Møller, M. Nielsen, E. Porcu, and E. Rubak. Determinantal point process models on the sphere. Bernoulli, 24:1171–1201, 2018.
  • Mrkvic̆ka et al. [2016] T. Mrkvic̆ka, S. Soubeyrand, M. Myllymäki, P. Grabarnik, and U. Hahn. Monte Carlo testing in spatial statistics, with applications to spatial residuals. Spatial Statistics, 18:40–53, 2016.
  • Mrkvic̆ka et al. [2017] T. Mrkvic̆ka, M. Myllymäki, and U. Hahn. Multiple Monte Carlo testing, with applications in spatial points. Statistics and Computing, 27:1239–1255, 2017.
  • Myllymäki et al. [2017] M. Myllymäki, T. Mrkvic̆ka, P. Grabarnik, H. Seijo, and U. Hahn. Global envelope tests for spatial processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79:381–404, 2017.
  • Neyman and Scott [1958] J. Neyman and E. L. Scott. Statistical approach to problems of cosmology. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 20:1–43, 1958.
  • Neyman and Scott [1972] J. Neyman and E. L. Scott. Processes of Clustering and Applications. John Wiley & Sons, New York, Stochastic Processes in Lewis PAW edition, 1972.
  • Nolan [2016] J. P. Nolan. An R package for modeling and simulating generalized spherical and related distributions. Journal of Statistical Distributions and Applications, 3:3–14, 2016.
  • Raskin [1994] R. G. Raskin. Spatial analysis on the sphere: a review, 1994. NGCIA technical report. (Available from http://eprints.cdlib.org/uc/item/5748n2xz). [Accessed February 1, 2018].
  • Ripley [1976] B. D. Ripley. The second-order analysis of spatial point processes. Journal of Applied Probability, 13:255–266, 1976.
  • Ripley [1977] B. D. Ripley. Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology), 39:172–212, 1977.
  • Robeson et al. [2014] S. M. Robeson, A. Li, and C. Huang. Point-pattern analysis on the sphere. Spatial Statistics, 10:76–86, 2014.
  • Simpson [2009] D. Simpson. Krylov Subspace Methods for Approximating Functions of Symmetric Positive Definite Matrices with Applications to Applied Statistics and Anomalous Diffusion. PhD thesis, Queensland University of Technology, 2009.
  • Simpson et al. [2016] D. Simpson, J. B. Illian, F. Lindgren, S. H. Sørbye, and H. Rue. Going off grid: computationally efficient inference for log-Gaussian Cox processes. Biometrika, 103:49–70, 2016.
  • Steinicke [2015] W. Steinicke. Revised new general catalogue and index catalogue, 2015. (Available from http://www.klima-luft.de/steinicke/index.htm). [Accessed February 1, 2018].
  • Szalay et al. [2007] A. S. Szalay, J. Gray, G. Fekete, P. Z. Kunszt, P. Kukol, and A. Thakar. Indexing the sphere with the hierarchical triangular mesh. Manuscript available at arXiv: 0701164, 2007.