Orthogonal series estimation of the pair correlation function of a spatial point process

TL;DR

This paper introduces a new orthogonal series estimator for the pair correlation function in spatial point processes, which is consistent, asymptotically normal, and can outperform kernel estimators especially in clustered patterns.

Contribution

The paper proposes a novel orthogonal series estimator for the pair correlation function, addressing bias issues in kernel methods for clustered point patterns.

Findings

Estimator is consistent and asymptotically normal.

Outperforms kernel estimators for Poisson and clustered processes.

Simulation results support theoretical properties.

Abstract

The pair correlation function is a fundamental spatial point process characteristic that, given the intensity function, determines second order moments of the point process. Non-parametric estimation of the pair correlation function is a typical initial step of a statistical analysis of a spatial point pattern. Kernel estimators are popular but especially for clustered point patterns suffer from bias for small spatial lags. In this paper we introduce a new orthogonal series estimator. The new estimator is consistent and asymptotically normal according to our theoretical and simulation results. Our simulations further show that the new estimator can outperform the kernel estimators in particular for Poisson and clustered point processes.

Tables

Table 1. Table 1: Monte Carlo mean, standard error, skewness (S) and kurtosis (K) of g ^ o ​ ( r ) subscript ^ 𝑔 𝑜 𝑟 \hat{g}_{o}(r) using the Bessel basis with the simple smoothing scheme in case of the Thomas process on observation windows W 1 = [ 0 , 1 ] 2 subscript 𝑊 1 superscript 0 1 2 W_{1}=[0,1]^{2} , W 2 = [ 0 , 2 ] 2 subscript 𝑊 2 superscript 0 2 2 W_{2}=[0,2]^{2} and W 3 = [ 0 , 3 ] 3 subscript 𝑊 3 superscript 0 3 3 W_{3}=[0,3]^{3} .

	$r$	$g\left(r\right)$	$\hat{E}\left\{{\hat{g}}_o\left(r\right)\right\}$	${\left[\widehat{\text{var}}\left\{{\hat{g}}_o\left(r\right)\right\}\right]}^{1/2}$	$\hat{\text{S}}\left\{{\hat{g}}_o\left(r\right)\right\}$	$\hat{\text{K}}\left\{{\hat{g}}_o\left(r\right)\right\}$
$W_1$	0.025	3.972	3.961	0.923	1.145	5.240
$W_1$	0.1	1.219	1.152	0.306	0.526	3.516
$W_2$	0.025	3.972	3.959	0.467	0.719	4.220
$W_2$	0.1	1.219	1.187	0.150	0.691	4.582
$W_3$	0.025	3.972	3.949	0.306	0.432	3.225
$W_3$	0.1	1.2187	1.2017	0.0951	0.2913	2.9573

Equations

$$E\{N(B)\}={\int }_B\rho (u)\mathrm{d}u,E\{N(A)N(B)\}={\int }_{A\cap B}\rho (u)\mathrm{d}u+{\int }_A{\int }_B{\rho }^{(2)}(u,v)\mathrm{d}u\mathrm{d}v.$$

$$\text{cov}\{N(A),N(B)\}=\int_{A\cap B}\rho({u}){\mathrm{d}}{u}+\int_{A}\int_{B}\rho({u})\rho({v})\big{\{}g({v},{u})-1\big{\}}{\mathrm{d}}{u}{\mathrm{d}}{v}$$

$${\hat{g}}_k(r;b)=\frac{1}{{\text{sa}}_d{r}^{d-1}}\sum _{u,v\in X_W}^{\ne }\frac{k_b(r-\Vert v-u\Vert )}{\rho (u)\rho (v)|W\cap W_{v-u}|},r\ge 0,$$

$${\hat{g}}_d(r;b)=\frac{1}{{\text{sa}}_d}\sum _{u,v\in X_W}^{\ne }\frac{k_b(r-\Vert v-u\Vert )}{\Vert v-u{\Vert }^{d-1}\rho (u)\rho (v)|W\cap W_{v-u}|},r\ge 0$$

$${\hat{g}}_c(r;b)={\hat{g}}_d(r;b)/c(r;b),c(r;b)={\int }_{-b}^{\min \{r,b\}}{k}_b(t)\mathrm{d}t,$$

$$\textsc{mise}(\hat{g}_{m},w)={\text{sa}_{d}}\int_{I}E\big{\{}\hat{g}_{m}(r;b)-g(r)\big{\}}^{2}w(r-{r_{\min}}){\mathrm{d}}r,$$

$$M(b)={\text{sa}_{d}}\int_{0}^{R}\big{\{}\hat{g}_{m}(r;b)\big{\}}^{2}r^{d-1}{\mathrm{d}}r-2\sum_{\begin{subarray}{c}{u},{v}\in X_{W}\\ \|{v}-{u}\|\leq R\end{subarray}}^{\neq}\frac{\hat{g}_{m}^{-\{{u},{v}\}}(\|{v}-{u}\|;b)}{\rho({u})\rho({v})|W\cap W_{{v}-{u}}|},$$

$$g(r)=\sum _{k=1}^{\infty }{\theta }_k{\varphi }_k(r-r_{\min }),$$

$${\hat{\theta }}_k=\frac{1}{{\text{sa}}_d}\sum _{\begin{array}{c}u,v\in X_W\\ r_{\min }<\Vert u-v\Vert <r_{\min }+R\end{array}}^{\ne }\frac{{\varphi }_k(\Vert v-u\Vert -r_{\min })w(\Vert v-u\Vert -r_{\min })}{\rho (u)\rho (v)\Vert v-u{\Vert }^{d-1}|W\cap W_{v-u}|},$$

$${\hat{g}}_o(r;b)=\sum _{k=1}^{\infty }{b}_k{\hat{\theta }}_k{\varphi }_k(r-r_{\min }),$$

$$\displaystyle\textsc{mise}\big{(}\hat{g}_{o},w\big{)}$$

$$\displaystyle={\text{sa}_{d}}\sum_{k=1}^{\infty}E(b_{k}\hat{\theta}_{k}-\theta_{k})^{2}={\text{sa}_{d}}\sum_{k=1}^{\infty}\big{[}b_{k}^{2}E\{(\hat{\theta}_{k})^{2}\}-2b_{k}\theta_{k}^{2}+\theta_{k}^{2}\big{]}.$$

$$b_k^{\ast }=\frac{{\theta }_k^2}{E\{({\hat{\theta }}_k{)}^2\}}=\frac{{\theta }_k^2}{{\theta }_k^2+\text{var}({\hat{\theta }}_k)},k\ge 0,$$

$${\hat{\theta }}_{k,n}=\frac{1}{{\text{sa}}_d|W_n|}\sum _{\begin{array}{c}u,v\in X_{W_n}\\ v-u\in B_{r_{\min }}^R\end{array}}^{\ne }\frac{{\varphi }_k(\Vert v-u\Vert -r_{\min })w(\Vert v-u\Vert -r_{\min })}{\rho (u)\rho (v)\Vert v-u{\Vert }^{d-1}{e}_n(v-u)}),$$

$$B_{r_{\min }}^R=\{h\in {\mathbb{R}}^d\mid r_{\min }<\Vert h\Vert <r_{\min }+R\}\text{and}{e}_n(h)=|W_n\cap (W_n{)}_h|/|W_n|.$$

$${\hat{g}}_{o,n}(r;b)=\sum _{k=1}^{K_n}{b}_k({\psi }_n){\hat{\theta }}_{k,n}{\varphi }_k(r-r_{\min })=\frac{1}{{\text{sa}}_d|W_n|}\sum _{\begin{array}{c}u,v\in X_{W_n}\\ v-u\in B_{r_{\min }}^R\end{array}}^{\ne }\frac{w(\Vert v-u\Vert ){\phi }_n(v-u,r)}{\rho (u)\rho (v)\Vert v-u{\Vert }^{d-1}{e}_n(v-u)|},$$

$${\phi }_n(h,r)=\sum _{k=1}^{K_n}{b}_k({\psi }_n){\varphi }_k(\Vert h\Vert -r_{\min }){\varphi }_k(r-r_{\min }).$$

$$E\left\{\sum _{u_1,\dots ,u_k\in X}^{\ne }\text{mathbbm}1(u_1\in A_1,\dots ,u_k\in A_k)\right\}={\int }_{A_1\times \cdots \times A_k}{\rho }^{(k)}(v_1,\dots ,v_k)\mathrm{d}{v}_1\cdots \mathrm{d}{v}_k$$

$$BC_{1}|W_{n}|^{-1}\sum_{k=1}^{\infty}\big{|}b_{k}(\psi_{n})\big{|}+\max_{1\leq k\leq m}\theta_{k}^{2}\sum_{k=1}^{m}(b_{k}(\psi_{n})-1)^{2}+(B^{2}+1)\sum_{k=m+1}^{\infty}\theta_{k}^{2}.$$

$$S_n=\frac{1}{{\text{sa}}_d|W_n|}\sum _{\begin{array}{c}u,v\in X_{W_n}\\ v-u\in B_{r_{\min }}^R\end{array}}^{\ne }\frac{f_n(v-u)}{\rho (u)\rho (v)e_n(v-u)}$$

$$\displaystyle\alpha_{\mathbf{X}}(m;a_{1},a_{2})=\sup\big{\{}$$

$$\displaystyle|B_{1}|\leq a_{1},|B_{2}|\leq a_{2},\mathcal{D}(B_{1},B_{2})\geq m,B_{1},B_{2}\in\mathcal{B}({\mathbb{R}}^{d})\big{\}},$$

$$\mathcal{D}(B_{1},B_{2})=\inf\big{\{}\max_{1\leq i\leq d}|u_{i}-v_{i}|:{u}=(u_{1},\ldots,u_{d})\in B_{1},{v}=(v_{1},\ldots,v_{d})\in B_{2}\big{\}}.$$

$$\{\text{var}({\hat{\theta }}_{k,n}){\}}^{-1/2}({\hat{\theta }}_{k,n}-{\theta }_k)\stackrel{D}{⟶}N(0,1).$$

$$f_n(h)=\frac{{\phi }_n(h,r)w(\Vert h\Vert -r_{\min })}{\Vert h{\Vert }^{d-1}}=\frac{w(\Vert h\Vert -r_{\min })}{\Vert h{\Vert }^{d-1}}\sum _{k=1}^{K_n}{b}_k({\psi }_n){\varphi }_k(\Vert h\Vert -r_{\min }){\varphi }_k(r-r_{\min }),$$

$$K_{n}^{-\omega}\sum_{k=1}^{K_{n}}b_{k}(\psi_{n})\big{|}\phi_{k}(r-{r_{\min}})\phi_{k}(\|{h}\|-{r_{\min}})\big{|}\leq M_{\omega}$$

$$\tau_{n}^{-1}|W_{n}|^{1/2}\big{[}\hat{g}_{o,n}(r;b)-E\{\hat{g}_{o,n}(r;b)\}\big{]}\stackrel{{\scriptstyle D}}{{\longrightarrow}}N(0,1).$$

$${\text{sa}_{d}}I(\psi)=\textsc{mise}(\hat{g}_{o},w)\!-\!\int_{{r_{\min}}}^{{r_{\min}}+R}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\big{\{}g(r)-1\big{\}}^{2}w(r){\mathrm{d}}r=\sum_{k=1}^{K}\big{[}b_{k}(\psi)^{2}E\{(\hat{\theta}_{k})^{2}\}-2b_{k}(\psi)\theta_{k}^{2}\big{]}.$$

$$\sum _{\begin{array}{c}u,v,u^{\prime },v^{\prime }\in X_W\\ v-u,v^{\prime }-u^{\prime }\in B_{r_{\min }}^R\end{array}}^{\ne }\frac{{\varphi }_k(\Vert v-u\Vert -r_{\min }){\varphi }_k(\Vert v^{\prime }-u^{\prime }\Vert -r_{\min })w(\Vert v-u\Vert -r_{\min })w(\Vert v^{\prime }-u^{\prime }\Vert -r_{\min })}{{{\text{sa}}_d}^2\rho (u)\rho (v)\rho (u^{\prime })\rho (v^{\prime })\Vert v-u{\Vert }^{d-1}\Vert v^{\prime }-u^{\prime }{\Vert }^{d-1}|W\cap W_{v-u}||W\cap W_{v^{\prime }-u^{\prime }}|}.$$

$$\underset{n\to \infty }{\lim }E(\widehat{{\theta }_{k,n}^2})\to {\left\{{\int }_0^{R}g(r+r_{\min }){\varphi }_k(r)w(r)\mathrm{d}r\right\}}^2={\theta }_k^2$$

$$\hat{I}(\psi)=\sum_{k=1}^{K}\big{\{}b_{k}(\psi)^{2}(\hat{\theta}_{k})^{2}-2b_{k}(\psi)\widehat{\theta_{k}^{2}}\big{\}}$$

Full text

Keywords: Asymptotic normality; Consistency; Kernel estimator; Orthogonal series estimator; Pair correlation function; Spatial point process.

1 Introduction

The pair correlation function is commonly considered the most informative second-order summary statistic of a spatial point process (Stoyan and Stoyan,, 1994; Møller and Waagepetersen,, 2003; Illian et al.,, 2008). Non-parametric estimates of the pair correlation function are useful for assessing regularity or clustering of a spatial point pattern and can moreover be used for inferring parametric models for spatial point processes via minimum contrast estimation (Stoyan and Stoyan,, 1996; Illian et al.,, 2008). Although alternatives exist (Yue and Loh,, 2013), kernel estimation is the by far most popular approach (Stoyan and Stoyan,, 1994; Møller and Waagepetersen,, 2003; Illian et al.,, 2008) which is closely related to kernel estimation of probability densities.

Kernel estimation is computationally fast and works well except at small spatial lags. For spatial lags close to zero, kernel estimators suffer from strong bias, see e.g. the discussion at page 186 in Stoyan and Stoyan, (1994), Example 4.7 in Møller and Waagepetersen, (2003) and Section 7.6.2 in Baddeley et al., (2015). The bias is a major drawback if one attempts to infer a parametric model from the non-parametric estimate since the behavior near zero is important for determining the right parametric model (Jalilian et al.,, 2013).

In this paper we adapt orthogonal series density estimators (see e.g. the reviews in Hall,, 1987; Efromovich,, 2010) to the estimation of the pair correlation function. We derive unbiased estimators of the coefficients in an orthogonal series expansion of the pair correlation function and propose a criterion for choosing a certain optimal smoothing scheme. In the literature on orthogonal series estimation of probability densities, the data are usually assumed to consist of indendent observations from the unknown target density. In our case the situation is more complicated as the data used for estimation consist of spatial lags between observed pairs of points. These lags are neither independent nor identically distributed and the sample of lags is biased due to edge effects. We establish consistency and asymptotic normality of our new orthogonal series estimator and study its performance in a simulation study and an application to a tropical rain forest data set.

2 Background

2.1 Spatial point processes

We denote by $X$ a point process on ${\mathbb{R}}^{d}$, $d\geq 1$, that is, $X$ is a locally finite random subset of ${\mathbb{R}}^{d}$. For $B\subseteq{\mathbb{R}}^{d}$, we let $N(B)$ denote the random number of points in $X\cap B$. That $X$ is locally finite means that $N(B)$ is finite almost surely whenever $B$ is bounded. We assume that $X$ has an intensity function $\rho$ and a second-order joint intensity $\rho^{(2)}$ so that for bounded $A,B\subset{\mathbb{R}}^{d}$,

[TABLE]

The pair correlation function $g$ is defined as $g({u},{v})=\rho^{(2)}({u},{v})/\{\rho({u})\rho({v})\}$ whenever $\rho({u})\rho({v})>0$ (otherwise we define $g({u},{v})=0$). By (1),

[TABLE]

for bounded $A,B\subset{\mathbb{R}}^{d}$. Hence, given the intensity function, $g$ determines the covariances of count variables $N(A)$ and $N(B)$. Further, for locations $u,v\in{\mathbb{R}}^{d}$, $g({u},{v})>1$ ($<1$) implies that the presence of a point at ${v}$ yields an elevated (decreased) probability of observing yet another point in a small neighbourhood of ${u}$ (e.g.  Coeurjolly et al.,, 2016). In this paper we assume that $g$ is isotropic, i.e. with an abuse of notation, $g({u},{v})=g(\|{v}-{u}\|)$. Examples of pair correlation functions are shown in Figure 1.

2.2 Kernel estimation of the pair correlation function

Suppose ${X}$ is observed within a bounded observation window $W\subset{\mathbb{R}}^{d}$ and let ${X}_{W}={X}\cap W$. Let $k_{b}(\cdot)$ be a kernel of the form $k_{b}(r)=k(r/b)/b$, where $k$ is a probability density and $b>0$ is the bandwidth. Then a kernel density estimator (Stoyan and Stoyan,, 1994; Baddeley et al.,, 2000) of $g$ is

[TABLE]

where ${\text{sa}_{d}}$ is the surface area of the unit sphere in ${\mathbb{R}}^{d}$, $\sum^{\neq}$ denotes sum over all distinct points, $1/|W\cap W_{{h}}|$, ${h}\in{\mathbb{R}}^{d}$, is the translation edge correction factor with $W_{{h}}=\{{u}-{h}:{u}\in W\}$, and $|A|$ is the volume (Lebesgue measure) of $A\subset{\mathbb{R}}^{d}$. Variations of this include (Guan, 2007a, )

[TABLE]

and the bias corrected estimator (Guan, 2007a, )

[TABLE]

assuming $k$ has bounded support $[-1,1]$. Regarding the choice of kernel, Illian et al., (2008), p. 230, recommend to use the uniform kernel $k(r)=\mathbbm{1}(|r|\leq 1)/2$, where $\mathbbm{1}(\,\cdot\,)$ denotes the indicator function, but the Epanechnikov kernel $k(r)=(3/4)(1-r^{2})\mathbbm{1}(|r|\leq 1)$ is another common choice. The choice of the bandwidth $b$ highly affects the bias and variance of the kernel estimator. In the planar ($d=2$) stationary case, Illian et al., (2008), p. 236, recommend $b=0.10/\surd{\hat{\rho}}$ based on practical experience where $\hat{\rho}$ is an estimate of the constant intensity. The default in spatstat (Baddeley et al.,, 2015), following Stoyan and Stoyan, (1994), is to use the Epanechnikov kernel with $b=0.15/\surd{\hat{\rho}}$.

Guan, 2007b and Guan, 2007a suggest to choose $b$ by composite likelihood cross validation or by minimizing an estimate of the mean integrated squared error defined over some interval $I$ as

[TABLE]

where $\hat{g}_{m}$, $m=k,d,c$, is one of the aforementioned kernel estimators, $w\geq 0$ is a weight function and ${r_{\min}}\geq 0$. With $I=(0,R)$, $w(r)=r^{d-1}$ and $r_{\min}=0$, Guan, 2007a suggests to estimate the mean integrated squared error by

[TABLE]

where $\hat{g}_{m}^{-\{{u},{v}\}}$, $m=k,d,c$, is defined as $\hat{g}_{m}$ but based on the reduced data $({X}\setminus\{{u},{v}\})\cap W$. Loh and Jang, (2010) instead use a spatial bootstrap for estimating (2). We return to (3) in Section 5.

3 Orthogonal series estimation

3.1 The new estimator

For an $R>0$, the new orthogonal series estimator of $g(r)$, $0\leq{r_{\min}}<r<{r_{\min}}+R$, is based on an orthogonal series expansion of $g(r)$ on $({r_{\min}},{r_{\min}}+R)$ :

[TABLE]

where $\{\phi_{k}\}_{k\geq 1}$ is an orthonormal basis of functions on $(0,R)$ with respect to some weight function $w(r)\geq 0$, $r\in(0,R)$. That is, $\int_{0}^{R}\phi_{k}(r)\phi_{l}(r)w(r){\mathrm{d}}r=\mathbbm{1}(k=l)$ and the coefficients in the expansion are given by $\theta_{k}=\int_{0}^{R}g(r+{r_{\min}})\phi_{k}(r)w(r){\mathrm{d}}r$.

For the cosine basis, $w(r)=1$ and $\phi_{1}(r)=1/\surd{R}$, $\phi_{k}(r)=(2/R)^{1/2}\cos\{(k-1)\pi r/R\}$, $k\geq 2$. Another example is the Fourier-Bessel basis with $w(r)=r^{d-1}$ and $\phi_{k}(r)=2^{1/2}J_{\nu}\left(r\alpha_{\nu,k}/R\right)r^{-\nu}/\{RJ_{\nu+1}(\alpha_{\nu,k})\}$, $k\geq 1$, where $\nu=(d-2)/2$, $J_{\nu}$ is the Bessel function of the first kind of order $\nu$, and $\{\alpha_{\nu,k}\}_{k=1}^{\infty}$ is the sequence of successive positive roots of $J_{\nu}(r)$.

An estimator of $g$ is obtained by replacing the $\theta_{k}$ in (4) by unbiased estimators and truncating or smoothing the infinite sum. A similar approach has a long history in the context of non-parametric estimation of probability densities, see e.g. the review in Efromovich, (2010). For $\theta_{k}$ we propose the estimator

[TABLE]

which is unbiased by the second order Campbell formula, see Section S2 of the supplementary material. This type of estimator has some similarity to the coefficient estimators used for probability density estimation but is based on spatial lags ${v}-{u}$ which are not independent nor identically distributed. Moreover the estimator is adjusted for the possibly inhomogeneous intensity $\rho$ and corrected for edge effects.

The orthogonal series estimator is finally of the form

[TABLE]

where $b=\{b_{k}\}_{k=1}^{\infty}$ is a smoothing/truncation scheme. The simplest smoothing scheme is $b_{k}=\mathbbm{1}[k\leq K]$ for some cut-off $K\geq 1$. Section 3.3 considers several other smoothing schemes.

3.2 Variance of $\hat{\theta}_{k}$

The factor $\|{v}-{u}\|^{d-1}$ in (5) may cause problems when $d>1$ where the presence of two very close points in $X_{W}$ could imply division by a quantity close to zero. The expression for the variance of $\hat{\theta}_{k}$ given in Section S2 of the supplementary material indeed shows that the variance is not finite unless $g(r)w(r-{r_{\min}})/r^{d-1}$ is bounded for ${r_{\min}}<r<{r_{\min}}+R$. If ${r_{\min}}>0$ this is always satisfied for bounded $g$. If ${r_{\min}}=0$ the condition is still satisfied in case of the Fourier-Bessel basis and bounded $g$.

For the cosine basis $w(r)=1$ so if ${r_{\min}}=0$ we need the boundedness of $g(r)/r^{d-1}$. If ${X}$ satisfies a hard core condition (i.e. two points in ${X}$ cannot be closer than some $\delta>0$), this is trivially satisfied. Another example is a determinantal point process (Lavancier et al.,, 2015) for which $g(r)=1-c(r)^{2}$ for a correlation function $c$. The boundedness is then e.g. satisfied if $c(\cdot)$ is the Gaussian ($d\leq 3$) or exponential ($d\leq 2$) correlation function. In practice, when using the cosine basis, we take ${r_{\min}}$ to be a small positive number to avoid issues with infinite variances.

3.3 Mean integrated squared error and smoothing schemes

The orthogonal series estimator (6) has the mean integrated squared error

[TABLE]

Each term in (7) is minimized with $b_{k}$ equal to (cf. Hall,, 1987)

[TABLE]

leading to the minimal value ${\text{sa}_{d}}\sum_{k=1}^{\infty}b_{k}^{*}\text{var}(\hat{\theta}_{k})$ of the mean integrated square error. Unfortunately, the $b_{k}^{*}$ are unknown.

In practice we consider a parametric class of smoothing schemes $b(\psi)$. For practical reasons we need a finite sum in (6) so one component in $\psi$ will be a cut-off index $K$ so that $b_{k}(\psi)=0$ when $k>K$. The simplest smoothing scheme is $b_{k}(\psi)=\mathbbm{1}(k\leq K)$. A more refined scheme is $b_{k}(\psi)=\mathbbm{1}(k\leq K)\hat{b}_{k}^{*}$ where $\hat{b}_{k}^{*}=\widehat{\theta_{k}^{2}}/(\hat{\theta}_{k})^{2}$ is an estimate of the optimal smoothing coefficient $b_{k}^{*}$ given in (8). Here $\widehat{\theta_{k}^{2}}$ is an asymptotically unbiased estimator of $\theta_{k}^{2}$ derived in Section 5. For these two smoothing schemes $\psi=K$. Adapting the scheme suggested by Wahba, (1981), we also consider $\psi=(K,c_{1},c_{2})$, $c_{1}>0,c_{2}>1$, and $b_{k}(\psi)=\mathbbm{1}(k\leq K)/(1+c_{1}k^{c_{2}})$. In practice we choose the smoothing parameter $\psi$ by minimizing an estimate of the mean integrated squared error, see Section 5.

3.4 Expansion of $g(\cdot)-1$

For large $R$, $g({r_{\min}}+R)$ is typically close to one. However, for the Fourier-Bessel basis, $\phi_{k}(R)=0$ for all $k\geq 1$ which implies $\hat{g}_{o}({r_{\min}}+R)=0$. Hence the estimator cannot be consistent for $r={r_{\min}}+R$ and the convergence of the estimator for $r\in({r_{\min}},{r_{\min}}+R)$ can be quite slow as the number of terms $K$ in the estimator increases. In practice we obtain quicker convergence by applying the Fourier-Bessel expansion to $g(r)-1=\sum_{k\geq 1}\vartheta_{k}\phi_{k}(r-{r_{\min}})$ so that the estimator becomes $\tilde{g}_{o}(r;b)=1+\sum_{k=1}^{\infty}b_{k}\hat{\vartheta}_{k}\phi_{k}(r-{r_{\min}})$ where $\hat{\vartheta}_{k}=\hat{\theta}_{k}-\int_{0}^{{r_{\min}}+R}\phi_{k}(r)w(r){\mathrm{d}}r$ is an estimator of $\vartheta_{k}=\int_{0}^{R}\{g(r+{r_{\min}})-1\}\phi_{k}(r)w(r){\mathrm{d}}r$. Note that $\text{var}(\hat{\vartheta}_{k})=\text{var}(\hat{\theta}_{k})$ and $\tilde{g}_{o}(r;b)-E\{\tilde{g}_{o}(r;b)\}=\hat{g}_{o}(r;b)-E\{\hat{g}_{o}(r;b)\}$. These identities imply that the results regarding consistency and asymptotic normality established for $\hat{g}_{o}(r;b)$ in Section 4 are also valid for $\tilde{g}_{o}(r;b)$.

4 Consistency and asymptotic normality

4.1 Setting

To obtain asymptotic results we assume that ${X}$ is observed through an increasing sequence of observation windows $W_{n}$. For ease of presentation we assume square observation windows $W_{n}=\times_{i=1}^{d}[-na_{i},na_{i}]$ for some $a_{i}>0$, $i=1,\ldots,d$. More general sequences of windows can be used at the expense of more notation and assumptions. We also consider an associated sequence $\psi_{n}$, $n\geq 1$, of smoothing parameters satisfying conditions to be detailed in the following. We let $\hat{\theta}_{k,n}$ and $\hat{g}_{o,n}$ denote the estimators of $\theta_{k}$ and $g$ obtained from ${X}$ observed on $W_{n}$. Thus

[TABLE]

where

[TABLE]

Further,

[TABLE]

where

[TABLE]

In the results below we refer to higher order normalized joint intensities $g^{(k)}$ of ${X}$. Define the $k$’th order joint intensity of $X$ by the identity

[TABLE]

for bounded subsets $A_{i}\subset{\mathbb{R}}^{d}$, $i=1,\ldots,k$, where the sum is over distinct ${u}_{1},\ldots,{u}_{k}$. We then let $g^{(k)}({v}_{1},\ldots,{v}_{k})=\rho^{(k)}({v}_{1},\ldots,{v}_{k})/\{\rho({v}_{1})\cdots\rho({v}_{k})\}$ and assume with an abuse of notation that the $g^{(k)}$ are translation invariant for $k=3,4$, i.e. $g^{(k)}({v}_{1},\ldots,{v}_{k})=g^{(k)}({v}_{2}-{v}_{1},\ldots,{v}_{k}-{v}_{1})$.

4.2 Consistency of orthogonal series estimator

Consistency of the orthogonal series estimator can be established under fairly mild conditions following the approach in Hall, (1987). We first state some conditions that ensure (see Section S2 of the supplementary material) that $\text{var}(\hat{\theta}_{k,n})\leq C_{1}/|W_{n}|$ for some $0<C_{1}<\infty$:

1. V1

There exists $0<\rho_{\min}<\rho_{\max}<\infty$ such that for all ${u}\in{\mathbb{R}}^{d}$, $\rho_{\min}\leq\rho({u})\leq\rho_{\max}$. 2. V2

For any ${h},{h}_{1},{h}_{2}\in{B_{r_{\min}}^{R}}$, $g({h})w(\|{h}\|-{r_{\min}})\leq C_{2}\|{h}\|^{d-1}$ and $g^{(3)}({h}_{1},{h}_{2})\leq C_{3}$ for constants $C_{2},C_{3}<\infty$. 3. V3

A constant $C_{4}<\infty$ can be found such that \sup_{{h}_{1},{h}_{2}\in{B_{r_{\min}}^{R}}}\int_{{\mathbb{R}}^{d}}\Big{|}g^{(4)}({h}_{1},{h}_{3},{h}_{2}+{h}_{3})-g({h}_{1})g({h}_{2})\Big{|}{\mathrm{d}}{h}_{3}\leq C_{4}.

The first part of V2 is needed to ensure finite variances of the $\hat{\theta}_{k,n}$ and is discussed in detail in Section 3.2. The second part simply requires that $g^{(3)}$ is bounded. The condition V3 is a weak dependence condition which is also used for asymptotic normality in Section 4.3 and for estimation of $\theta_{k}^{2}$ in Section 5.

Regarding the smoothing scheme, we assume

1. S1

B=\sup_{k,\psi}\big{|}b_{k}(\psi)\big{|}<\infty and for all $\psi$, \sum_{k=1}^{\infty}\big{|}b_{k}(\psi)\big{|}<\infty. 2. S2

$\psi_{n}\rightarrow\psi^{*}$ for some $\psi^{*}$, and \lim_{\psi\rightarrow\psi^{*}}\max_{1\leq k\leq m}\big{|}b_{k}(\psi)-1\big{|}=0 for all $m\geq 1$. 3. S3

|W_{n}|^{-1}\sum_{k=1}^{\infty}\big{|}b_{k}(\psi_{n})\big{|}\rightarrow 0.

E.g. for the simplest smoothing scheme, $\psi_{n}=K_{n}$, $\psi^{*}=\infty$ and we assume $K_{n}/|W_{n}|\rightarrow 0$.

Assuming the above conditions we now verify that the mean integrated squared error of $\hat{g}_{o,n}$ tends to zero as $n\rightarrow\infty$. By (7), \textsc{mise}\big{(}\hat{g}_{o,n},w\big{)}/{\text{sa}_{d}}=\sum_{k=1}^{\infty}\big{[}b_{k}(\psi_{n})^{2}\text{var}(\hat{\theta}_{k})+\theta_{k}^{2}\{b_{k}(\psi_{n})-1\}^{2}\big{]}. By V1-V3 and S1 the right hand side is bounded by

[TABLE]

By Parseval’s identity, $\sum_{k=1}^{\infty}\theta_{k}^{2}<\infty$. The last term can thus be made arbitrarily small by choosing $m$ large enough. It also follows that $\theta_{k}^{2}$ tends to zero as $k\rightarrow\infty$. Hence, by S2, the middle term can be made arbitrarily small by choosing $n$ large enough for any choice of $m$. Finally, the first term can be made arbitrarily small by S3 and choosing $n$ large enough.

4.3 Asymptotic normality

The estimators $\hat{\theta}_{k,n}$ as well as the estimator $\hat{g}_{o,n}(r;b)$ are of the form

[TABLE]

for a sequence of even functions $f_{n}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$. We let $\tau^{2}_{n}=|W_{n}|\text{var}(S_{n})$.

To establish asymptotic normality of estimators of the form (11) we need certain mixing properties for $X$ as in Waagepetersen and Guan, (2009). The strong mixing coefficient for the point process $X$ on ${\mathbb{R}}^{d}$ is given by (Ivanoff,, 1982; Politis et al.,, 1998)

[TABLE]

where $\mathcal{B}({\mathbb{R}}^{d})$ denotes the Borel $\sigma$-field on ${\mathbb{R}}^{d}$, $\mathcal{F}_{X}(B_{i})$ is the $\sigma$-field generated by ${X}\cap B_{i}$ and

[TABLE]

To verify asymptotic normality we need the following assumptions as well as V1 (the conditions V2 and V3 are not needed due to conditions N2 and N4 below):

1. N1

The mixing coefficient satisfies $\alpha_{{X}}(m;(s+2R)^{d},\infty)=O(m^{-d-\varepsilon})$ for some $s,\varepsilon>0$. 2. N2

There exists a $\eta>0$ and $L_{1}<\infty$ such that $g^{(k)}({h}_{1},\ldots,{h}_{k-1})\leq L_{1}$ for $k=2,\ldots,2(2+\lceil\eta\rceil)$ and all ${h}_{1},\ldots,{h}_{k-1}\in{\mathbb{R}}^{d}$. 3. N3

$\liminf_{n\rightarrow\infty}\tau^{2}_{n}>0$. 4. N4

There exists $L_{2}<\infty$ so that $|f_{n}({h})|\leq L_{2}$ for all $n\geq 1$ and ${h}\in{B_{r_{\min}}^{R}}$.

The conditions N1-N3 are standard in the point process literature, see e.g. the discussions in Waagepetersen and Guan, (2009) and Coeurjolly and Møller, (2014). The condition N3 is difficult to verify and is usually left as an assumption, see Waagepetersen and Guan, (2009), Coeurjolly and Møller, (2014) and Dvořák and Prokešová, (2016). However, at least in the stationary case, and in case of estimation of $\hat{\theta}_{k,n}$, the expression for $\text{var}(\hat{\theta}_{k,n})$ in Section S2 of the supplementary material shows that $\tau_{n}^{2}=|W_{n}|\text{var}(\hat{\theta}_{k,n})$ converges to a constant which supports the plausibility of condition N3. We discuss N4 in further detail below when applying the general framework to $\hat{\theta}_{k,n}$ and $\hat{g}_{o,n}$. The following theorem is proved in Section S3 of the supplementary material.

Theorem 1

*Under conditions V1, N1-N4, \tau_{n}^{-1}|W_{n}|^{1/2}\big{\{}S_{n}-E(S_{n})\big{\}}\stackrel{{\scriptstyle D}}{{\longrightarrow}}N(0,1). *

4.4 Application to $\hat{\theta}_{k,n}$ and $\hat{g}_{o,n}$

In case of estimation of $\theta_{k}$, $\hat{\theta}_{k,n}=S_{n}$ with $f_{n}({h})=\phi_{k}(\|{h}\|-{r_{\min}})w(\|{h}\|-{r_{\min}})/\|{h}\|^{d-1}$. The assumption N4 is then straightforwardly seen to hold in the case of the Fourier-Bessel basis where $|\phi_{k}(r)|\leq|\phi_{k}(0)|$ and $w(r)=r^{d-1}$. For the cosine basis, N4 does not hold in general and further assumptions are needed, cf. the discussion in Section 3.2. For simplicity we here just assume ${r_{\min}}>0$. Thus we state the following

Corollary 1

Assume V1, N1-N4, and, in case of the cosine basis, that ${r_{\min}}>0$. Then

[TABLE]

For $\hat{g}_{o,n}(r;b)=S_{n}$,

[TABLE]

where $\varphi_{n}$ is defined in (10). In this case, $f_{n}$ is typically not uniformly bounded since the number of not necessarily decreasing terms in the sum defining $\varphi_{n}$ in (10) grows with $n$. We therefore introduce one more condition:

1. N5

There exist an $\omega>0$ and $M_{\omega}<\infty$ so that

[TABLE]

for all ${h}\in{B_{r_{\min}}^{R}}$.

Given N5, we can simply rescale: $\tilde{S}_{n}:=K_{n}^{-\omega}S_{n}$ and $\tilde{\tau}^{2}_{n}:=K_{n}^{-2\omega}\tau^{2}_{n}$. Then, assuming $\liminf_{n\rightarrow\infty}\tilde{\tau}_{n}^{2}>0$, Theorem 1 gives the asymptotic normality of $\tilde{\tau}_{n}^{-1}|W_{n}|^{1/2}\{\tilde{S}_{n}-E(\tilde{S}_{n})\}$ which is equal to $\tau_{n}^{-1}|W_{n}|^{1/2}\{S_{n}-E(S_{n})\}$. Hence we obtain

Corollary 2

Assume V1, N1-N2, N5 and $\liminf_{n\rightarrow\infty}K_{n}^{-2\omega}\tau_{n}^{2}>0$. In case of the cosine basis, assume further ${r_{\min}}>0$. Then for $r\in({r_{\min}},{r_{\min}}+R)$,

[TABLE]

In case of the simple smoothing scheme $b_{k}(\psi_{n})=\mathbbm{1}(k\leq K_{n})$, we take $\omega=1$ for the cosine basis. For the Fourier-Bessel basis we take $\omega=4/3$ when $d=1$ and $\omega=d/2+2/3$ when $d>1$ (see the derivations in Section S6 of the supplementary material).

5 Tuning the smoothing scheme

In practice we choose $K$, and other parameters in the smoothing scheme $b(\psi)$, by minimizing an estimate of the mean integrated squared error. This is equivalent to minimizing

[TABLE]

In practice we must replace (12) by an estimate. Define $\widehat{\theta^{2}_{k}}$ as

[TABLE]

Then, referring to the set-up in Section 4 and assuming V3,

[TABLE]

(see Section S4 of the supplementary material) and hence $\widehat{\theta^{2}_{k,n}}$ is an asymptotically unbiased estimator of $\theta_{k}^{2}$. The estimator is obtained from $(\hat{\theta}_{k})^{2}$ by retaining only terms where all four points ${u},{v},{u}^{\prime},{v}^{\prime}$ involved are distinct. In simulation studies, $\widehat{\theta_{k}^{2}}$ had a smaller root mean squared error than $(\hat{\theta}_{k})^{2}$ for estimation of $\theta_{k}^{2}$.

Thus

[TABLE]

is an asymptotically unbiased estimator of (12). Moreover, (13) is equivalent to the following slight modification of Guan, 2007a ’s criterion (3):

[TABLE]

For the simple smoothing scheme $b_{k}(K)=\mathbbm{1}(k\leq K)$, (13) reduces to

[TABLE]

where $\hat{b}^{*}_{k}=\widehat{\theta_{k}^{2}}/(\hat{\theta}_{k})^{2}$ is an estimator of $b^{*}_{k}$ in (8).

In practice, uncertainties of $\hat{\theta}_{k}$ and $\widehat{\theta_{k}^{2}}$ lead to numerical instabilities in the minimization of (13) with respect to $\psi$. To obtain a numerically stable procedure we first determine $K$ as

[TABLE]

That is, $\hat{K}$ is the first local minimum of (14) larger than 1 and smaller than an upper limit $K_{\max}$ which we chose to be 49 in the applications. This choice of $K$ is also used for the refined and the Wahba smoothing schemes. For the refined smoothing scheme we thus let $b_{k}=\mathbbm{1}(k\leq\hat{K})\hat{b}_{k}^{*}$. For the Wahba smoothing scheme $b_{k}=\mathbbm{1}(k\leq\hat{K})/(1+\hat{c}_{1}k^{\hat{c}_{2}})$, where $\hat{c}_{1}$ and $\hat{c}_{2}$ minimize $\sum_{k=1}^{\hat{K}}\left\{(\hat{\theta}_{k})^{2}/(1+c_{1}k^{c_{2}})^{2}-2\widehat{\theta_{k}^{2}}/(1+c_{1}k^{c_{2}})\right\}$ over $c_{1}>0$ and $c_{2}>1$.

6 Simulation study

We compare the performance of the orthogonal series estimators and the kernel estimators for data simulated on $W=[0,1]^{2}$ or $W=[0,2]^{2}$ from four point processes with constant intensity $\rho=100$. More specifically, we consider $n_{\text{sim}}=1000$ realizations from a Poisson process, a Thomas process (parent intensity $\kappa=25$, dispersion standard deviation $\omega=0.0198$), a Variance Gamma cluster process (parent intensity $\kappa=25$, shape parameter $\nu=-1/4$, dispersion parameter $\omega=0.01845$, Jalilian et al.,, 2013), and a determinantal point process with pair correlation function $g(r)=1-\exp\{-2(r/\alpha)^{2}\}$ and $\alpha=0.056$. The pair correlation functions of these point processes are shown in Figure 1.

For each realization, $g(r)$ is estimated for $r$ in $({r_{\min}},{r_{\min}}+R)$, with ${r_{\min}}=10^{-3}$ and $R=0.06,0.085,0.125$, using the kernel estimators $\hat{g}_{k}(r;b)$, $\hat{g}_{d}(r;b)$ and $\hat{g}_{c}(r;b)$ or the orthogonal series estimator $\hat{g}_{o}(r;b)$. The Epanechnikov kernel with bandwidth $b=0.15/\surd{\hat{\rho}}$ is used for $\hat{g}_{k}(r;b)$ and $\hat{g}_{d}(r;b)$ while the bandwidth of $\hat{g}_{c}(r;b)$ is chosen by minimizing Guan, 2007a ’s estimate (3) of the mean integrated squared error. For the orthogonal series estimator, we consider both the cosine and the Fourier-Bessel bases with simple, refined or Wahba smoothing schemes. For the Fourier-Bessel basis we use the modified orthogonal series estimator described in Section 3.4. The parameters for the smoothing scheme are chosen according to Section 5.

From the simulations we estimate the mean integrated squared error (2) with $w(r)=1$ of each estimator $\hat{g}_{m}$, $m=k,d,c,o$, over the intervals $[{r_{\min}},0.025]$ (small spatial lags) and $[{r_{\min}},{r_{\min}}+R]$ (all lags). We consider the kernel estimator $\hat{g}_{k}$ as the baseline estimator and compare any of the other estimators $\hat{g}$ with $\hat{g}_{k}$ using the log relative efficiency $e_{I}(\hat{g})=\log\{\widehat{\textsc{mise}}_{I}(\hat{g}_{k})/\widehat{\textsc{mise}}_{I}(\hat{g})\}$, where $\widehat{\textsc{mise}}_{I}(\hat{g})$ denotes the estimated mean squared integrated error over the interval $I$ for the estimator $\hat{g}$. Thus $e_{I}(\hat{g})>0$ indicates that $\hat{g}$ outperforms $\hat{g}_{k}$ on the interval $I$. Results for W=$[0,1]^{2}$ are summarized in Figure 2.

For all types of point processes, the orthogonal series estimators outperform or does as well as the kernel estimators both at small lags and over all lags. The detailed conclusions depend on whether the non-repulsive Poisson, Thomas and Var Gamma processes or the repulsive determinantal process are considered. Orthogonal-Bessel with refined or Wahba smoothing is superior for Poisson, Thomas and Var Gamma but only better than $\hat{g}_{c}$ for the determinantal point process. The performance of the orthogonal-cosine estimator is between or better than the performance of the kernel estimators for Poisson, Thomas and Var Gamma and is as good as the best kernel estimator for determinantal. Regarding the kernel estimators, $\hat{g}_{c}$ is better than $\hat{g}_{d}$ for Poisson, Thomas and Var Gamma and worse than $\hat{g}_{d}$ for determinantal. The above conclusions are stable over the three $R$ values considered. For $W=[0,2]^{2}$ (see Figure S1 in the supplementary material) the conclusions are similar but with more clear superiority of the orthogonal series estimators for Poisson and Thomas. For Var Gamma the performance of $\hat{g}_{c}$ is similar to the orthogonal series estimators. For determinantal and $W=[0,2]^{2}$, $\hat{g}_{c}$ is better than orthogonal-Bessel-refined/Wahba but still inferior to orthogonal-Bessel-simple and orthogonal-cosine. Figures S2 and S3 in the supplementary material give a more detailed insight in the bias and variance properties for $\hat{g}_{k}$, $\hat{g}_{c}$, and the orthogonal series estimators with simple smoothing scheme. Table S1 in the supplementary material shows that the selected $K$ in general increases when the observation window is enlargened, as required for the asymptotic results. The general conclusion, taking into account the simulation results for all four types of point processes, is that the best overall performance is obtained with orthogonal-Bessel-simple, orthogonal-cosine-refined or orthogonal-cosine-Wahba.

To supplement our theoretical results in Section 4 we consider the distribution of the simulated $\hat{g}_{o}(r;b)$ for $r=0.025$ and $r=0.1$ in case of the Thomas process and using the Fourier-Bessel basis with the simple smoothing scheme. In addition to $W=[0,1]^{2}$ and $W=[0,2]^{2}$, also $W=[0,3]^{2}$ is considered. The mean, standard error, skewness and kurtosis of $\hat{g}_{o}(r)$ are given in Table 1 while histograms of the estimates are shown in Figure S3. The standard error of $\hat{g}_{o}(r;b)$ scales as $|W|^{1/2}$ in accordance with our theoretical results. Also the bias decreases and the distributions of the estimates become increasingly normal as $|W|$ increases.

7 Application

We consider point patterns of locations of Acalypha diversifolia (528 trees), Lonchocarpus heptaphyllus (836 trees) and Capparis frondosa (3299 trees) species in the 1995 census for the $1000\text{m}\times 500\text{m}$ Barro Colorado Island plot (Hubbell and Foster,, 1983; Condit,, 1998). To estimate the intensity function of each species, we use a log-linear regression model depending on soil condition (contents of copper, mineralized nitrogen, potassium and phosphorus and soil acidity) and topographical (elevation, slope gradient, multiresolution index of valley bottom flatness, ncoming mean solar radiation and the topographic wetness index) variables. The regression parameters are estimated using the quasi-likelihood approach in Guan et al., (2015). The point patterns and fitted intensity functions are shown in Figure S5 in the supplementary material.

The pair correlation function of each species is then estimated using the bias corrected kernel estimator $\hat{g}_{c}(r;b)$ with $b$ determined by minimizing (3) and the orthogonal series estimator $\hat{g}_{o}(r;b)$ with both Fourier-Bessel and cosine basis, refined smoothing scheme and the optimal cut-offs $\hat{K}$ obtained from (15); see Figure 3.

For Lonchocarpus the three estimates are quite similar while for Acalypha and Capparis the estimates deviate markedly for small lags and then become similar for lags greater than respectively 2 and 8 meters. For Capparis and the cosine basis, the number of selected coefficients coincides with the chosen upper limit 49 for the number of coefficients. The cosine estimate displays oscillations which appear to be artefacts of using high frequency components of the cosine basis. The function (14) decreases very slowly after $K=7$ so we also tried the cosine estimate with $K=7$ which gives a more reasonable estimate.

Acknowledgement

Rasmus Waagepetersen is supported by the Danish Council for Independent Research — Natural Sciences, grant ”Mathematical and Statistical Analysis of Spatial Data”, and by the ”Centre for Stochastic Geometry and Advanced Bioimaging”, funded by the Villum Foundation.

Supplementary material

Supplementary material includes proofs of consistency and asymptotic normality results and details of the simulation study and data analysis.