Jean-François Coeurjolly
email: jean-francois.coeurjolly@upmf-grenoble.fr.

Jesper Møller

Rasmus Waagepetersen

Palm distributions for log Gaussian Cox processes

Keywords: J-function, joint intensities, Laplace approximation, nearestneighbour distribution function, spatial point process

This paper establishes a remarkable result regarding Palm distributions for a log Gaussian Cox process: the reduced Palm distribution for a log Gaussian Cox process is itself a log Gaussian Cox process which only differs from the original log Gaussian Cox process in the intensity function. This new result is used to study functional summaries for log Gaussian Cox processes.

Introduction

Palm distributions (see e.g. [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF][START_REF] Daley | An Introduction to the Theory of Point Processes, Volume II: General Theory and Structure[END_REF] are important in the theory and application of spatial point processes. Intuitively speaking, for a prespecified location in space, the Palm distribution of a point process, with respect to this location, plays the role of the conditional distribution of the point process given that the aforementioned location is occupied by a point of the point process.

The present paper focuses on log Gaussian Cox processes [START_REF] Møller | Log Gaussian Cox processes[END_REF] which provide very flexible, useful, and popular models for modeling spatial patterns in e.g. biology and spatial epidemiology. The paper establishes a surprisingly simple characterization of Palm distributions for such a process: The reduced n-point Palm distribution is for any n ≥ 1 itself a log Gaussian Cox process that only differs from the original log Gaussian Cox process in its intensity function (not to be confused with the random intensity function generating this kind of Cox process). This result can be exploited for functional summaries as discussed later. The simplicity and completeness of this result is remarkable when compared with Palm distributions for other common classes of spatial point processes. Reduced Palm distributions for Gibbs point processes are also themselves Gibbs point processes but with densities only known up to a normalizing constant. For shot-noise Cox processes [START_REF] Møller | Shot noise Cox processes[END_REF] one-point reduced Palm distributions have a simple characterization as cluster processes similar to shot-noise Cox processes but this is not the case for n-point Palm distributions when n > 1.

The paper is organized as follows. Section 2 reviews the general definition of reduced Palm distributions of any order and relates this to Cox processes. Section 3 establishes our characterization result for log Gaussian Cox processes. Section 4 applies this result to functional summaries for stationary log Gaussian Cox processes, in particular the so-called F , G, and J-functions, where we establish some new theoretical results, consider how to calculate F , G, and J using Laplace approximations, and discuss an application. Section 5 concludes the paper.

Palm distributions

Our general setting is as follows. For ease of exposition we view a point process as a random locally finite subset X of a Borel set S ⊆ R d , d ≥ 1; for measure theoretical details, see e.g. [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] or [START_REF] Daley | An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods[END_REF]. Denoting X B = X ∩ B the restriction of X to a set B ⊆ S, local finiteness of X means that X B is finite almost surely (a.s.) whenever B is bounded. We denote N the state space consisting of the locally finite subsets (or point configurations) of S. We use the generic notation h for an arbitrary non-negative measurable function defined on N , S n , or S n × N for n = 1, 2, . . .. Furthermore, B 0 is the family of all bounded Borel subsets of S. Finally, recall that the void probabilities P(X K = ∅), K ⊆ S compact, uniquely determine the distribution of X.

Factorial moment measures and Palm distributions

This section provides the general definition of reduced Palm distributions of any order. For finite point processes specified by a density, a simpler and more explicit definition is available as reviewed in [START_REF] Coeurjolly | Conditioning in spatial point processes[END_REF].

For n = 1, 2, . . . and B i ∈ B 0 , the n-th order factorial moment measure α (n) is defined by

α (n) (B 1 × B 2 × • • • × B n ) = E = x 1 ,...,xn∈X 1(x 1 ∈ B 1 , . . . , x n ∈ B n ),
where 1(•) denotes the indicator function and = over the summation sign means that x 1 , . . . , x n are pairwise distinct. If α (n) has a density ρ (n) with respect to Lebesgue measure, ρ (n) is called the n-th order joint intensity function and is determined up to a Lebesgue nullset. Therefore, we can assume that ρ (n) (x 1 , . . . , x n ) is invariant under permutations of x 1 , . . . , x n , and we need only to consider the case where x 1 , . . . , x n ∈ S are pairwise distinct. Then ρ (n) (x 1 , . . . , x n ) dx 1 • • • dx n can be interpreted as the approximate probability for X having a point in each of infinitesimally small regions around x 1 , . . . , x n of volumes dx 1 , . . . dx n , respectively. We also write ρ(u) for the intensity function ρ (1) (u).

Moreover, for any measurable F ⊆ N , define the n-th order reduced Campbell measure C (n)! as the measure on S n × N given by

C (n)! (B 1 × B 2 × • • • × B n × F ) = E = x 1 ,...,xn∈X 1(x 1 ∈ B 1 , . . . , x n ∈ B n , X \ {x 1 , . . . , x n } ∈ F ).
Note that C (n)! (• × F ), as a measure on S n , is absolutely continuous with respect to α (n) , with a density P !

x 1 ,...,xn (F ) which is determined up to an α (n) nullset, and

α (n) (N ) = C (n)! (B 1 × B 2 × • • • × B n × N )
. By the so-called Campbell-Mecke formula/theorem, we can assume that P !

x 1 ,...,xn (•) is a point process distribution on N , called the n-th order reduced Palm distribution given x 1 , . . . , x n (see e.g. [START_REF] Daley | An Introduction to the Theory of Point Processes, Volume II: General Theory and Structure[END_REF]. We denote by X ! x 1 ,...,xn a point process distributed according to P ! x 1 ,...,xn . Again we need only to consider the case where x 1 , . . . , x n are pairwise distinct. Then P ! x 1 ,...,xn can be interpreted as the conditional distribution of X \ {x 1 , . . . , x n } given that x 1 , . . . , x n ∈ X.

If ρ (n) exists, then by standard measure theoretical arguments we obtain the extended Campbell-Mecke formula

E = x 1 ,...,xn∈X h(x 1 , . . . , x n , X \ {x 1 , . . . , x n }) = S • • • S Eh(x 1 , . . . , x n , X ! x 1 ,...,xn )ρ (n) (x 1 , . . . , x n ) dx 1 • • • dx n . (1)
Suppose ρ (m+n) exists for an m ≥ 1 and n ≥ 1. Then, for pairwise distinct u 1 , . . . , u m , x 1 , . . . , x n ∈ S, it follows easily by expressing α (m+n) as an expectation of the form (1) that X ! x 1 ,...,xn has m-th order joint intensity function

ρ (m) x 1 ,...,xn (u 1 , . . . , u m ) = ρ (m+n) (u 1 ,...,um,x 1 ,...,xn) ρ (n) (x 1 ,...,xn) if ρ (n) (x 1 , . . . , x n ) > 0, 0 otherwise. (2)
We also write ρ x 1 ...,xn for the intensity function ρ

(1)

x 1 ...,xn .

Cox processes

Let Λ = {Λ(x)} x∈S be a nonnegative random field such that Λ is locally integrable a.s., that is, for any B ∈ B 0 , the integral B Λ(x) dx exists and is finite a.s. In the sequel, X conditional on Λ is assumed to be a Poisson process with intensity function Λ; we say that X is a Cox process driven by Λ. We also assume that Λ has moments of any order n = 1, 2, . . .. Then the joint intensities of X exist: For any n = 1, 2, . . . and pairwise distinct x 1 , . . . , x n ∈ S,

ρ (n) (x 1 , . . . , x n ) = E n i=1 Λ(x i ) . (3)
The following lemma, which is verified in Appendix A, gives a characterization of the reduced Palm distributions and their void probabilities.

Lemma 1. Let X be a Cox process satisfying the conditions above. Then, for any n = 1, 2, . . ., pairwise distinct x 1 , . . . , x n ∈ S, and compact K ⊆ S,

ρ (n) (x 1 , . . . , x n )E h x 1 , . . . , x n , X ! x 1 ,...,xn = E h(x 1 , . . . , x n , X) n i=1 Λ(x i ) (4) and ρ (n) (x 1 , . . . , x n )P(X ! x 1 ,...,xn ∩ K = ∅) = E exp - K Λ(u) du n i=1 Λ(x i ) .
(5)

3 Reduced Palm distributions for log Gaussian Cox processes

For the remainder of this paper, let X be a Cox process driven by Λ = {Λ(x)} x∈S , where Λ(x) = exp{Y (x)} and Y = {Y (x)} x∈S is a Gaussian random field with mean function µ and covariance function c so that Λ is locally integrable a.s. (simple conditions ensuring this are given in [START_REF] Møller | Log Gaussian Cox processes[END_REF]. Then X is a log Gaussian Cox process (LGCP) as introduced by [START_REF] Coles | A lognormal model for the cosmological mass distribution[END_REF] in astronomy and independently by [START_REF] Møller | Log Gaussian Cox processes[END_REF] in statistics.

For distinct x, y ∈ S, define the so-called pair correlation function g(x, y) = ρ (2) (x, y)/{ρ(x)ρ(y)} (the following result shows that ρ > 0 in the present case). By Møller et al. (1998, Theorem 1),

ρ(x) = exp{µ(x)+c(x, x)/2}, g(x, y) = exp{c(x, y)}, x, y ∈ S, x = y, (6) and for pairwise distinct x 1 , . . . , x n ∈ S, ρ (n) (x 1 , . . . , x n ) = n i=1 ρ(x i ) 1≤i<j≤n g(x i , x j ) (7) is strictly positive. For u, x 1 , . . . , x n ∈ S, define µ x 1 ,...,xn (u) = µ(u) + n i=1 c(u, x i ).
Combining ( 2) and ( 6)-( 7), we obtain for any pairwise distinct u 1 , . . . , u m , x 1 , . . . , x n ∈ S with m > 0 and n > 0,

ρ (m) x 1 ,...,xn (u 1 , . . . , u m ) = m i=1 ρ x 1 ,...,xn (u i ) 1≤i<j≤m g(u i , u j ) , (8)
where

ρ x 1 ,...,xn (u) = exp {µ x 1 ,...,xn (u) + c(u, u)/2} .
Thereby the following proposition follows.

Proposition 1. For the LGCP X and any pairwise distinct x 1 , . . . , x n ∈ S, X ! x 1 ,...,xn has m-th order joint intensity (8) which agrees with the m-th order joint intensity function for an LGCP with mean function µ x 1 ,...,xn and covariance function c for the underlying Gaussian random field.

Proposition 1 indicates that also X !

x 1 ,...,xn could be an LGCP. A sufficient condition, considered by [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], is the existence of a number a = a

(B) > 1 for each set B ∈ B 0 such that E exp a B Λ(u) du < ∞. (9)
However, we have not been successful in verifying this condition which seems too strong to hold for any of the covariance function models we have considered, including when c is constant (then X is a mixed Poisson process) or weaker cases of correlation, e.g. if c is a stationary exponential covariance function. The case where c is constant is closely related to the log normal distribution which is not uniquely determined by its moments [START_REF] Heyde | On a property of the lognormal distribution[END_REF]. Accordingly we use instead Lemma 1 when establishing the following theorem, which implies that the LGCPs X and X !

x 1 ,...,xn share the same pair correlation function and differ only in their intensity functions.

Theorem 1. For pairwise distinct x 1 , . . . , x n ∈ S, X !

x 1 ,...,xn is an LGCP with underlying Gaussian random field Y x 1 ,...,xn , where Y x 1 ,...,xn has mean function µ x 1 ,...,xn and covariance function c.

Let Ỹ = Y -µ be the centered Gaussian random field with covariance function c. Theorem 1 is a consequence of the fact that the probability measure of Y x 1 ,...,xn is absolutely continuous with respect to the one of Ỹ, with density exp n i=1 ỹ(x i ) -n i,j=1 c(x i , x j )/2 when ỹ is a realization of Ỹ. This result is related to the Cameron-Martin-Girsanov formula for onedimensional Gaussian processes. A short selfcontained proof covering our spatial setting is given in Appendix A.

Often we consider a non-negative covariance function c or equivalently g ≥ 1, which is interpreted as 'attractiveness of the LGCP at all ranges', but even more can be said: A coupling between X and X !

x 1 ,...,xn is obtained by taking Y x 1 ,...,xn (x) = Y (x) + n i=1 c(x, x i ). Thus, if c ≥ 0 and we are given pairwise distinct points x 1 , . . . , x n ∈ S, we can consider X as being included in X !

x 1 ,...,xn , since X can be obtained by an independent thinning of X ! x 1 ,...,xn , with inclusion probabilities exp{-n i=1 c(x, x i )}, x ∈ X \ {x 1 , . . . , x n }. This property clearly shows the attractiveness of the LGCP if c ≥ 0 (equivalently g ≥ 1).

Functional summaries for stationary log Gaussian Cox processes

Throughout this section, let S = R d and assume that the LGCP X is stationary, i.e., its distribution is invariant under translations in R d . By ( 6

K(r) = 1 ρ E # X ! o ∩ B(o, r) = B(o,r) g(x).
Thus, ρK(r) is the expected number of further points in X within distance r of a typical point in X. If g(x) depends only on the distance x then g and K are in one-to-one correspondence;

(ii) the empty space function given by

F (r) = P {X ∩ B(o, r) = ∅} ,
which is the probability that X has a point within distance r of an arbitrary fixed location;

(iii) the nearest-neighbour distribution function given by

G(r) = P X ! o ∩ B(o, r) = ∅ ,
which is the probability that X has a further point within distance r of a typical point in X;

(iv) the J-function given by

J(r) = 1 -G(r) 1 -F (r) ,
with the convention a/0 = 0 for any a ≥ 0.

Section 4.1 establishes some new results for these theoretical functions and Section 4.2 discusses how they can be calculated using a Laplace approximation. Section 4.3 illustrates this calculation and Section 4.4 discusses an application for a real dataset.

New formulae for G and J

By conditioning on Y, we see that

1 -F (r) = E exp - B(o,r) exp{Y (x)} dx . ( 10
)
Using the Slivnyak-Mecke formula, [START_REF] Møller | Log Gaussian Cox processes[END_REF] showed that

1 -G(r) = 1 ρ E exp Y (o) - B(o,r) exp{Y (x)} dx . ( 11
)
Since the nearest-neighbour distribution function for X is the same as the empty space function for X ! o , which is an LGCP with underlying Gaussian random field Y o (x) = Y (x) + c(x), and since g(x) = exp{c(x)}, we obtain an alternative expression

1 -G(r) = E exp - B(o,r) g(x) exp{Y (x)} dx . ( 12
)
Therefore, we also obtain a new expression for the J-function,

J(r) = E exp -B(o,r) g(x) exp{Y (x)} dx E exp -B(o,r) exp{Y (x)} dx . ( 13
)
Van Lieshout ( 2011) established for a general stationary point process the approximation J(r) -1 ≈ -ρ{K(r) -ω d r d }, where ω d = |B(o, 1)| and r → ω d r d is the K-function for a stationary Poisson process. It is therefore not so surprising that often empirical J and K-functions lead to the same practical interpretations. In particular, if for our LGCP c ≥ 0, i.e., g ≥ 1, then we have K(r) -ω d r d ≥ 0, and so we expect that J(r) ≤ 1. Indeed Van Lieshout (2011) verified this in the case of an LGCP with g ≥ 1. This result immediately follows by the new expression (13).

Laplace approximation

Since Laplace's pioneering work (see e.g. [START_REF] Stigler | Memoir on the probability of the causes of events, Laplace, Pierre Simon[END_REF], Laplace approximations of complex integrals have gained much attention in probability and statistics, in particular when considering integrals involving Gaussian random fields (see e.g. [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF]. This section discusses a Laplace approximation of 1 -G(r); a Laplace approximation of 1 -F (r) can be obtained along similar lines.

For ∆ > 0, consider a grid of quadrature points,

G(∆) = {(∆i 1 , . . . , ∆i d ) | i 1 , . . . , i d ∈ Z},
and for v ∈ G(∆), let A ∆ v = [v 1 -∆/2, v 1 + ∆/2[× • • • × [v d -∆/2, v d + ∆/2
[ be the grid cell associated with v. Then for any non-negative Borel function : R d → R, we use the numerical quadrature approximation

B(o,r) exp{Y (x)} (x) dx ≈ v∈G(∆)∩B(o,r) w v (v) exp{Y (v)}, (14)
where the quadrature weight

w v = |A ∆ v ∩ B(o, r)|.
Denote by f , M and Σ the density, the mean vector and the covariance matrix of the normally distributed vector {Y (v)} v∈G(∆)∩B(o,r) . Then ( 12) and ( 14) give

1 -G(r) ≈ R m exp    - v∈G(∆)∩B(o,r) w v g(v) exp(y v )    f (y) dy = R m exp{h(y)} dy ( 15
)
where y is the vector (y v ) v∈G(∆)∩B(o,r) of dimension m = #{G(∆) ∩ B(o, r)} and

h(y) = - v∈G(∆)∩B(o,r) w v g(v) exp(y v ) - 1 2 (y -M ) Σ -1 (y -M ) - 1 2 log{(2π) m |Σ|}.
The gradient vector for h is

∇h(y) = -d(y) -Σ -1 (y -M ), (16)
where d(y) = {w v g(v) exp(y v )} v∈G(∆)∩B(o,r) , and minus one times the Hessian matrix for h is

H(y) = D(y) + Σ -1 ,
where D(y) is the diagonal matrix with entries {d(y)} v , vG(∆) ∩ B(o, r)}.

Since H(y) is a positive definite matrix, h has a unique maximum at a point ŷ, which can be found using Newton-Raphson iterations

y (l+1) = y (l) + H -1 {y (l) }∇h{y (l) }. (17)
Therefore, the logarithm of the Laplace approximation of the right hand side in (15) (see e.g. [START_REF] Stigler | Memoir on the probability of the causes of events, Laplace, Pierre Simon[END_REF] gives

log{1 -G(r)} ≈ - v∈G(∆)∩B(o,r) w v g(v) exp(ŷ v ) + 1 2 (ŷ -M ) d(ŷ) - 1 2 log |D(ŷ)Σ + I|. ( 18
)
where I is the m × m identity matrix. For the computation of Σ -1 (y -M ) in ( 16) we solve LL z = y -M where L is the Cholesky factor of Σ. In the same way, considering the QR decomposition, Q(y)R(y) for y ∈ R m , of the matrix D(y)Σ + I, the computation of H -( 1) {y (l) }∇h{y (l) } in ( 17) is done by first solving Q{y (l) }R{y (l) }z = ∇h{y (l) } and second by evaluating Σz. Finally, in (18), |D(ŷ)Σ + I| = |R(ŷ)|.

Numerical illustration

To illustrate the Laplace approximations of the G and J-functions (Section 4.2) we consider three planar stationary LGCPs with intensity ρ = 50 and spherical covariance function

c(x) =      σ 2 1 -2 π x α 1 -x α 2 + sin -1 x α if x ≤ α, 0 otherwise,
with variance σ 2 = 4 and scale parameters α = 0.1, 0.2, 0.3, respectively. We evaluate the approximations of G(r) and J(r) at r ∈ R, where R is the set of 50 equispaced values between 0.01 and 0.25. For r ∈ R, we define the grid G(∆ r ) with ∆ r = 2r/q, where q is a fixed integer. Such a choice implies that #{G(∆ r ) ∩ [-r, r] 2 } = q 2 , and so we have at least q 2 π/4 quadrature points in B(o, r). For a given q, denote by G q , F q , and J q the corresponding Laplace approximations of G, F , and J, respectively. Figure 1 shows the resulting curves with q = 16. To see how far these Cox processes deviate from the Poisson case (which would correspond to σ 2 = 0), we also plot the G-function in the Poisson case, namely 1 -G(r) = exp(-ρπr 2 ). To study the role of q, we report in Table 1 the maximal differences max r∈R |G 16 (r) -G q (r)| and max r∈R |J 16 (r)-J q (r)| for q = 4, 8, 12. As expected, each difference decreases as q increases and is already very small when q = 12 (less than 4×10 -3 except for the J-function and α = 1). This justifies our choice q = 16 in Figure 1. 1.0 r J function q q q q q q α=0.3 α=0.2 α=0.1

Poisson

Figure 1: three planar stationary LGCPs with intensity ρ = 50 and c given by a spherical covariance function, with variance σ 2 = 4 and scale parameters α = 0.1, 0.2, 0.3, respectively, Laplace approximations of the Gfunction (left) and the J-function (right). 1: For the same three LGCPs as in Figure 1, maximal differences between the Laplace approximations G q and G 16 , and between the Laplace approximations J q and J 16 , with q = 4, 8, 12, that is, the evaluation of max r∈R |H 16 (r) -H q (r)| for H = G, J. Results are multiplied by 10 3 .

q = 4 q = 8 q = 12 α = 0.1, G 59
The Laplace approximation of the G-function could also be derived using (11). To check the agreement of the numerical approximations based on ( 11) and ( 12), respectively, Table 2 shows the maximal difference between the two approximations of first the G-function and second the J-function. In agreement with the theoretical developments, in both cases, the difference does not exceed 4 × 10 -4 when q = 16. q = 4 q = 8 q = 12 q = 16 α = 0.1, G 3.8 8.4 4.4 3.2 J 15 8.5 6.1 3.9 α = 0.2, G 4.7 3.1 3.9 1.5

J 4.7 3.1 2.1 1.7 α = 0.3, G 0.1 1.9 1.4 1.1 J 0.2 2.0 1.5 1.3
Table 2: For the same three LGCPs as in Figure 1, maximal differences between the Laplace approximations of the G-function based on (11) respective ( 12), with q = 4, 8, 12, 16, and similarly for the J-function. Results are multiplied by 10 4 .

Scots pine saplings dataset

The left panel in Figure 2 shows the locations of 126 Scots pine saplings in a 10 by 10 metre square. The dataset is included in the R package spatstat as finpines, and it has previously been analyzed by [START_REF] Penttinen | Marked point processes in forest statistics[END_REF], Stoyan and[START_REF] Stoyan | Fractals, Random Shapes, and Point Fields: Methods of Geometrical Statistics[END_REF][START_REF] Møller | Log Gaussian Cox processes[END_REF]. In the first two papers a Matérn cluster process is fitted, using the K-function (or its equivalent Lfunction) and its nonparametric estimate both for parameter estimation and model checking, while the third paper considered an LGCP with exponential covariance function and used the pair correlation function for parameter estimation and the F and G-functions for model checking. [START_REF] Møller | Log Gaussian Cox processes[END_REF] concluded that both models provide a reasonable fit although when also including a third-order functional summary (i.e., one based on X ! o,x ) the LGCP model showed a better fit. Below we extend this analysis by using the J-function and the approximation established in Section 4.2.

We fitted both models by minimum contrast estimation (method kppm in spatstat) which compares a non-parametric estimate of the K-function with its theoretical value. When approximating the J-function for the LGCP, we used q = 12; no improvements were noticed with higher values of q. The right panel in Figure 2 shows the theoretical J-functions for the two fitted models together with a non-parametric estimate of the J-function obtained from data, considering 50 equispaced distances (r-values) between 0 and 0.9 meter (for the exact expression of the J-function for the Matérn cluster process, see e.g. [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]. Clearly, the fitted LGCP provides a better fit than the fitted Matérn cluster process. Indeed, the maximal difference between the non-parametric estimate and the theoretical J-function equals 0.43 for the Matérn cluster model and 0.20 for the LGCP model. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.

Concluding remarks

We expect that our results for the reduced Palm distributions for an LGCP can be exploited further regarding third-order and higher order functional summaries (one such characteristic was briefly mentioned in Section 4.4), parameter estimation procedures, model checking, etc. This is discussed in some detail below. For likelihood based inference, suppressing in the notation any dependence on a parametric model and assuming a realization x of an LGCP is observed within a region W ⊂ R d of Lebesgue measure |W | ∈ (0, ∞), the likelihood function is given by the density

f (x) = E exp |W | - W exp{Y (u)} du u∈x exp{Y (u)}
with respect to the unit rate Poisson process on W . This density expression has no explicit form even for simple covariance function models for the underlying Gaussian random field. For maximum likelihood estimation and prediction of the Gaussian random field, rather elaborate and time-consuming Markov chain Monte Carlo (MCMC) methods [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]) may be used; alternatively, a Bayesian approach based on integrated nested Laplace approximations or MCMC methods may be used [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF][START_REF] Taylor | INLA or MCMC? a tutorial and comparative evaluation for spatial prediction in log-Gaussian Cox processes[END_REF]. Using the results of this paper, we have the following new expression for the density:

f (x) = ρ (n) (x 1 , . . . , x n )E exp |W | - W exp{Y x 1 ,...,xn (u)} du .
Here the expression for the n-th order intensity ρ (n) is explicit, but it remains to investigate if the expectation, which apparently has a simpler expression, may be easier to approximate.

For model checking, when considering non-parametric estimates of functional summaries together with simulated confidence bounds under a fitted LGCP model (such as extreme rank envelopes, see [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF][START_REF] Myllymäki | Global envelope tests for spatial processes[END_REF], it could be pertinent to include the theoretical expressions of the functional summaries for LGCPs obtained in the present paper.

Finally, recall that for any point process, the pair correlation function (when it exists) is invariant under independent thinning. Could this property be exploited in connection to LGCPs where we know how the pair correlation function is related to those of the reduced Palm distributions?

A Proofs

Proof of Lemma 1: By conditioning on Λ, (1) becomes

n i=1 Λ(x i ) dx 1 • • • dx n Λ (19) = S • • • S E h(x 1 , . . . , x n , X) n i=1 Λ(x i ) dx 1 • • • dx n . (20)
Here, in (19) we use that X given Λ is a Poisson process and apply the extended Slivnyak-Mecke theorem [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF], and in (20) we use Fubini's theorem. Combining ( 1) and ( 20), we deduce (4). Finally, (5) follows from (4) with h(x 1 , . . . , x n , x) = 1(x ∩ K = ∅).

Proof of Theorem 1: By ( 5) and ( 6)-( 7), we just have to show that for any compact K ⊆ S and pairwise distinct points x 1 , . . . , x n ∈ S, This will follow by verifying that the distribution of { Ỹ (u)+ n i=1 c(u, x i )} u∈S is absolutely continuous with respect to the distribution of Ỹ = { Ỹ (u)} u∈S , with density exp n i=1 ỹ(x i ) -n i,j=1 c(x i , x j )/2 when ỹ is a realization of Ỹ. Since the distribution of a random field is determined by its finite dimensional distributions, we just need to verify the agreement of the characteristic functions of the probability measures Q 1 and Q 2 given by = exp e T Σe/2 E exp i{ Ỹ (u)} T u∈I (t -ie) = exp e T Σe/2 + ie T Σt -t T Σt/2 -e T Σe/2 = exp ie T Σt -t T Σt/2 . The last expression is the characteristic function of Q 1 which concludes the proof. For the second last equality in the above derivation we considered { Ỹ (u)} u∈U as a complex Gaussian vector and used the expression for its characteristic function.

  )-(7), this is equivalent to stationarity of the underlying Gaussian random field Y, that is, the intensity ρ is constant and the pair correlation function g(x, y) = g(x -y) is translation invariant, where x, y ∈ R d are distinct, and g(x) = exp{c(x)} and c(x) = c(o, x) for x ∈ R d , where o denotes the origin in R d . It is custom to call P ! o the reduced Palm distribution at a typical point, noticing that for any x ∈ R d , X ! o and X ! x -x = {y -x : y ∈ X ! x } are identically distributed.Denote B(o, r) the ball in R d of radius r > 0 and centered at o. Popular tools for exploratory purposes as well as model fitting and model checking are based on the following functional summaries where r > 0 (see e.g.[START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF]):(i) the pair correlation function g and the related Ripley's K-functiongiven by

Figure 2 :

 2 Figure 2: Left panel: Locations of 126 Scots pine saplings in a 10 by 10 metre square. Right panel: Non-parametric estimate of the J-function (solid curve) and fitted J-functions for the Matérn cluster process (dashed curve) and the LGCP with exponential covariance function (dotted curve).

  1 , . . . , x n , X \ {x 1 , . . . , x n }) 1 , . . . , x n , X)

  i , x j )/2 -K exp µ(u) + Ỹ (u) du .

  i , x j )/2 , for any Borel set B ⊆ R n+m , any pairwise distinct locations u 1 , . . . , u m ∈ R d \ {x 1 , . . . , x n }, with m > 0 and U = {x 1 , . . . , x n , u 1 , . . . , u m }. Let Σ = {c(u, v)} u,v∈U denote the covariance matrix of { Ỹ (u)} u∈U , and let c x 1 ,...,xn = { n i=1 c(u, x i )} u∈U = Σe, where e consists of n 1's followed by m 0's. For t ∈ R n+m , the characteristic function of Q 2 is (with i 2 = -1) E exp i{ Ỹ (u)} T u∈I t + { Ỹ (u)} T u∈I e + e T Σe/2

Acknowledgments

We thank the two reviewers, the associate editor and the editor for their careful reading of our paper and their useful comments and suggestions which helped us to improve the paper.

J. Møller and R. Waagepetersen were supported by the Danish Council for Independent Research -Natural Sciences, grant 12-124675, "Mathematical and Statistical Analysis of Spatial Data", and by the "Centre for Stochastic Geometry and Advanced Bioimaging", funded by grant 8721 from the Villum Foundation. J.-F. Coeurjolly was supported by ANR-11-LABX-0025 PERSYVAL-Lab (2011, project OculoNimbus).