Second-order variational equations for spatial point processes with a view to pair correlation function estimation

TL;DR

This paper develops second-order variational equations for spatial point processes, enabling closed-form parameter estimation for pair correlation functions, and introduces a method to fit orthogonal series expansions of these functions.

Contribution

It introduces a novel application of second-order variational equations to spatial point processes, facilitating explicit parameter estimation for pair correlation functions.

Findings

Variational equations can be used to derive closed-form estimators.

Method allows fitting orthogonal series expansions of pair correlation functions.

Applicable to log linear parametric models for pair correlation functions.

Abstract

Second-order variational type equations for spatial point processes are established. In case of log linear parametric models for pair correlation functions, it is demonstrated that the variational equations can be applied to construct estimating equations with closed form solutions for the parameter estimates. This result is used to fit orthogonal series expansions of log pair correlation functions of general form.

Figures

Tables

Table 1. Table 1: Square-root of the MISE for different estimates of log ⁡ g 0 subscript 𝑔 0 \log g_{0} , observation windows and models. The figures between brackets correspond to the average of the selected K 𝐾 K s. The NA’s are due to occurrence of non-positive estimates. ( ⋆ : in this setting one replication produced an outlier and is omitted in the root MISE estimation)

	Window	OSE	VSE	KDE
Poisson	${[0,1]}^2$	0.027 (2.1)	0.051 (2.2)	0.093
	${[0,2]}^2$	0.012 (2.0)	0.024 (2.2)	0.037
Thomas	${[0,1]}^2$	0.0995 (3.7)	${0.1418}^{\star }$ (2.7)	0.111
	${[0,2]}^2$	0.044 (4.2)	0.063 (2.9)	0.053
Variance Gamma	${[0,1]}^2$	0.099 (6.5)	0.148 (3.8)	0.110
	${[0,2]}^2$	0.050 (9.6)	0.072 (5.3)	0.057
DPP	${[0,1]}^2$	NA (3)	0.1622 (3.6)	NA
	${[0,2]}^2$	NA (4.1)	0.1582 (5.2)	NA

Equations

$$\mathbb{E}\sum _{u_1,\dots ,u_k\in \mathbf{X}}^{\ne }h(u_1,\dots ,u_k)$$

$$g(u,v)=\frac{{\rho }^{(2)}(u,v)}{\rho (u)\rho (v)}$$

$$\nabla h(w)={\left\{\frac{\partial h}{\partial w_1}(w),\dots ,\frac{\partial h}{\partial w_d}(w)\right\}}^{\top },w\in {\mathbb{R}}^d$$

$$\mathrm{div}h(w)=\sum _{i=1}^d\frac{\partial h_i}{\partial w_i}(w).$$

$$\underset{n\to \infty }{\lim }{\int }_{\partial B_n}g(w)h(w)\cdot \nu (\mathrm{d}w)=0$$

$$\displaystyle{\mathbb{E}}\Bigg{\{}\sum_{u,v\in\mathbf{X}\cap W}^{\neq}e(u,v)\nabla\log g(v-u)\,$$

$$-\mathbb{E}\left\{\sum _{u,v\in \mathbf{X}\cap W}^{\ne }e(u,v)\mathrm{div}\mathrm{h}(v-u)\right\},$$

$$t\mapsto h(t)g_0^{\prime }(t)\in L^1({\mathbb{R}}^{+})\text{ and }\underset{n\to \infty }{\lim }\{g_0(n)h(n)-g_0(0)h(0)\}=0$$

$$t\mapsto t^{d-1}h(t)g_0^{\prime }(t)\in L^1({\mathbb{R}}^{+})\text{ and }\underset{n\to \infty }{\lim }\{n^{d-1}{g}_0(n)h(n)-g_0(0)h(0)1(d=1)\}=0.$$

$$\displaystyle{\mathbb{E}}\bigg{\{}\sum_{u,v\in\mathbf{X}\cap W}^{\neq}\frac{e(u,v)}{\|v-u\|^{d-1}}h(\|v-u\|)$$

$$-\mathbb{E}\left\{\sum _{u,v\in \mathbf{X}\cap W}^{\ne }\frac{e(u,v)}{\Vert v-u{\Vert }^{d-1}}{h}^{\prime }(\Vert v-u\Vert )\right\},$$

$$\displaystyle{\mathbb{E}}\bigg{\{}$$

$$-\mathbb{E}\left[\sum _{u,v\in \mathbf{X}\cap W}^{\ne }e(u,v)\left\{(d-1)\frac{h(\Vert v-u\Vert )}{\Vert v-u\Vert }+h^{\prime }(\Vert v-u\Vert )\right\}\right].$$

$$\sum _{u,v\in \mathbf{X}\cap W}^{\ne }\frac{e(u,v)}{\Vert v-u{\Vert }^{d-1}}\mathbf{h}(\Vert v-u\Vert )(\log g_0{)}^{\prime }(\Vert v-u\Vert )+\sum _{u,v\in \mathbf{X}\cap W}^{\ne }\frac{e(u,v)}{\Vert v-u{\Vert }^{d-1}}{\mathbf{h}}^{\prime }(\Vert v-u\Vert )$$

$$S(\mathbf{\beta })=-\mathbb{E}\sum _{u,v\in \mathbf{X}\cap W}^{\ne }\frac{e(u,v)}{\Vert v-u{\Vert }^{d-1}}\mathbf{h}(\Vert v-u\Vert )\frac{\mathrm{d}}{\mathrm{d}{\mathbf{\beta }}^{\top }}(\log g_0{)}^{\prime }(\Vert v-u\Vert ).$$

$$S(\mathbf{\beta })=-{\varsigma }_d{\int }_0^{\infty }\mathbf{h}(t)\left[\frac{\mathrm{d}}{\mathrm{d}{\mathbf{\beta }}^{\top }}(\log g_0{)}^{\prime }(t)\right]g_0(t)\mathrm{d}t,$$

$$S(\mathbf{\beta })=-{\varsigma }_d{\int }_0^{\infty }\mathbf{h}(t)\left[\frac{\mathrm{d}}{\mathrm{d}{\mathbf{\beta }}^{\top }}(\log g_0{)}^{\prime }(t)\right]g_0(t)t^{d-1}\mathrm{d}t.$$

$$\log g_0(t)={\mathbf{\beta }}^{\top }\mathbf{r}(t)={\mathbf{\beta }}^{\top }{\left\{r_1(t),\dots ,r_K(t)\right\}}^{\top }$$

$$\hat{\mathbf{\beta }}=-{\mathbf{A}}^{-1}\mathbf{b}.$$

$$t\mapsto \Vert {\mathbf{r}}^{\prime }(t){\Vert }^2\psi (t)\in L^1({\mathbb{R}}^d)\text{ and }\underset{n\to \infty }{\lim }\psi (n)\mathbf{r}(n{)}^{\top }{\mathbf{r}}^{\prime }(n)-\psi (0)\mathbf{r}(0{)}^{\top }{\mathbf{r}}^{\prime }(0)=0$$

$$t\mapsto t^{d-1}\Vert {\mathbf{r}}^{\prime }(t){\Vert }^2\psi (t)\in L^1({\mathbb{R}}^d)$$

$$\text{ and }\underset{n\to \infty }{\lim }{n}^{d-1}\psi (n)\mathbf{r}(n{)}^{\top }{\mathbf{r}}^{\prime }(n)-\psi (0)\mathbf{r}(0{)}^{\top }{\mathbf{r}}^{\prime }(0)1(d=1)=0.$$

$$\mathbf{A}$$

$$\mathbf{b}$$

$$\mathbf{A}$$

$$\mathbf{b}$$

$$\displaystyle\qquad\qquad+\psi^{\prime}(\|v-u\|)\mathbf{r}^{\prime}(\|v-u\|)+\psi(\|v-u\|)\mathbf{r}^{\prime\prime}(\|v-u\|)\bigg{\}}$$

$$\mathbf{A}\mathbf{\beta }+\mathbf{b}=0$$

$$\log g(w)=\sum _{k=1}^K{\mathbf{\beta }}_k{r}_k(w)={\mathbf{\beta }}^{\top }\mathbf{r}(w)$$

$$\log g_0(t)=\sum _{k=1}^{\infty }{\beta }_k{\varphi }_k(t-r_{\min })$$

$$\widehat{\log g_{0,K}}(t)=\sum _{k=1}^K{\hat{\beta }}_k{\varphi }_k(t-r_{\min }).$$

$${\varphi }_k(r)=\frac{2^{1/2}}{R{J}_{\nu +1}({\alpha }_{\nu ,k})}{J}_{\nu }\left(r{\alpha }_{\nu ,k}/R\right)r^{-\nu },k\ge 1,$$

$$\displaystyle\textsc{mise}\big{(}\widehat{\log g_{0,K}}\big{)}$$

$$\displaystyle={\varsigma_{d}}\sum_{k=1}^{\infty}{\mathbb{E}}(b_{k}\hat{\beta}_{k}-\beta_{k})^{2}={\varsigma_{d}}\sum_{k=1}^{\infty}\big{[}b_{k}^{2}{\mathbb{E}}\{\hat{\beta}_{k}^{2}\}-2b_{k}\beta_{k}{\mathbb{E}}\hat{\beta}_{k}+\beta_{k}^{2}\big{]}.$$

Full text

1 Introduction

Spatial point processes are models for sets of random locations of possibly interacting objects. Background on spatial point processes can be found in Møller and Waagepetersen (2004), Illian et al. (2008) or Baddeley et al. (2015) which gives both an accessible introduction as well as details on implementation in the R package spatstat. Moments of counts of objects for spatial point processes are typically expressed in terms of so-called joint intensity functions or Papangelou conditional intensity functions which are defined via the Campbell or Georgii-Nguyen-Zessin equations (see the aforementioned references or the concise review of intensity functions and Campbell formulae in Section 2). In this paper we consider a third type of equation called variational equations.

A key feature of variational equations compared to Campbell and Georgii-Nguyen-Zessin equations is that they are formulated in terms of the gradient of the log intensity or conditional intensity function rather than the (conditional) intensity itself. Variational equations were introduced for parameter estimation in Markov random fields by Almeida et al. (1993). The authors suggested the terminology ‘variational’ due to the analogy between the derivation of their estimating equation and the variational Euler-Lagrange equations in partial differential equations. The resulting equation consisted in an equilibrium equation involving the gradient of the log conditional probability of the Markov random field. Later, Baddeley and Dereudre (2013) obtained variational equations for Gibbs point processes and exploited them to infer a log-linear parametric model of the conditional intensity function. Coeurjolly and Møller (2014) established a first-order variational equation for general spatial point processes and used it to estimate parameters in a log-linear parametric model for the intensity function.

The first contribution of this paper is to establish second-order variational equations. The second-order properties of a spatial point process are characterized by the so-called pair correlation function which is a normalized version of the second-order joint intensity function. We assume that the pair correlation function is translation invariant and also consider the case when it is isotropic. Since the new variational equations are based on the gradient of the log pair correlation function, they take a particularly simple form for pair correlation functions of log-linear form.

Our second contribution is to propose a new non-parametric estimator of the pair correlation function. The classical approach is to use a kernel estimator, see e.g. Møller and Waagepetersen (2004). More recently, Jalilian et al. (2019) investigated the estimation of the pair correlation function using an orthogonal series expansion. In the setting of their simulation studies, the orthogonal series estimator was shown to be more efficient than the standard kernel estimator. One drawback, however, is that the orthogonal series estimator is not guaranteed to be non-negative. We therefore propose to use our second-order variational equation to estimate coefficients in an orthogonal series expansion of the log pair correlation function. This ensures that the resulting pair correlation function estimator is non-negative. We compare our new estimator with the previous ones in a simulation study and also illustrate its use on real datasets.

2 Background and main results

2.1 Spatial point processes

Throughout this paper we let $\mathbf{X}$ be a spatial point process defined on $\mathbb{R}^{d}$. That is, $\mathbf{X}$ is a random subset of $\mathbb{R}^{d}$ with the property that the intersection of $\mathbf{X}$ with any bounded subset of $\mathbb{R}^{d}$ is of finite cardinality. The joint intensity functions $\rho^{(k)}$, $k\geq 1$, are characterized (when they exist) by the Campbell formulae (equations) (see e.g.  Møller and Waagepetersen, 2004): for any $h:(\mathbb{R}^{d})^{k}\to\mathbb{R}^{+}$ (with $\mathbb{R}^{+}$ the non-negative real numbers)

[TABLE]

More intuitively, for any pairwise distinct points $u_{1},\dots,u_{k}\in\mathbb{R}^{d}$, $\rho^{(k)}(u_{1},\dots,u_{k})\mathrm{d}u_{1}\cdots\mathrm{d}u_{k}$ is the probability that for each $i=1,\ldots,k$, $\mathbf{X}$ has a point in an infinitesimally small region around $u_{i}$ with volume $\mathrm{d}u_{i}$. The intensity function $\rho$ corresponds to the case $k=1$, i.e. $\rho=\rho^{(1)}$. The pair correlation function is obtained by normalizing the second-order joint intensity $\rho^{(2)}$:

[TABLE]

for pairwise distinct $u,v$ and where $g(u,v)$ is set to 0 if $\rho(u)$ or $\rho(v)$ is zero. Intuitively, $g(u,v)>1$ [$g(u,v)<1$] means that presence of a point at $u$ increases [decreases] the probability of observing a further point at $v$ and vice versa. We assume that $\mathbf{X}$ is observed on some bounded domain $W\subset\mathbb{R}^{d}$ with volume $|W|>0$ and without loss of generality we assume that $\rho(u)>0$ for all $u\in W$ (otherwise we just replace $W$ by $\{u\in W|\rho(u)>0\}$ provided the latter set has positive volume).

We will always assume that $\mathbf{X}$ is second-order intensity reweighted stationary (Baddeley et al., 2000), meaning that its pair correlation function $g$ is invariant by translations. We then, with an abuse of notation, write $g(v-u)$ for $g(u,v)$ for any $u,v\in\mathbb{R}^{d}$. We will also consider the case of an isotropic pair correlation function in which case $g(v-u)$ depends only on the distance $\|v-u\|$.

For the presentation of the second-order variational type equation in the next section some additional notation is needed. For a function $h:\mathbb{R}^{d}\to\mathbb{R}$ which is differentiable on $\mathbb{R}^{d}$, we denote by

[TABLE]

the gradient vector with respect to the $d$ coordinates. The inner product is denoted by a ‘$\cdot$’ and for $h:\mathbb{R}^{d}\to\mathbb{R}^{d}$, a multivariate function such that each component is differentiable on $\mathbb{R}^{d}$, we define the divergence operator by

[TABLE]

2.2 Second-order variational equations

In this section, we present our new second-order variational equations. The proofs of the results are given in the Appendices.

Theorem 1**.**

Assume $\mathbf{X}$ is second-order intensity reweighted stationary. Let $h:\mathbb{R}^{d}\to\mathbb{R}^{d}$ be a componentwise continuously differentiable function on $\mathbb{R}^{d}$. Assume that $g$ is continuously differentiable on $\mathbb{R}^{d}$, that $\|h\|\|\nabla g\|\in L^{1}(\mathbb{R}^{d})$, and that there exists a sequence of increasing bounded domains $(B_{n})_{n\geq 1}$ such that $B_{n}\to\mathbb{R}^{d}$ as $n\to\infty$, with piecewise smooth boundary $\partial B_{n}$ and such that

[TABLE]

where $\nu$ stands for the outer normal measure to $\partial B_{n}$. Then

[TABLE]

where $e:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{+}$ denotes the function $e(u,v)=\{\rho(u)\rho(v)|W\cap W_{v-u}|\}^{-1}$ for any $u,v\in\mathbb{R}^{d}$ and where $W_{w}$ denotes the domain $W$ translated by $w\in\mathbb{R}^{d}$.

We note that condition (3) is in particular satisfied if the function $h$ is compactly supported.

We next consider the case where the pair correlation function is isotropic, i.e. for any $u,v\in\mathbb{R}^{d}$ there exists $g_{0}:\mathbb{R}^{+}\to\mathbb{R}^{+}$ such that $g(u,v)=g(v-u)=g_{0}(\|v-u\|)$.

Theorem 2**.**

Assume $\mathbf{X}$ is second-order intensity reweighted stationary with isotropic pair correlation function $g_{0}$. Let $h:\mathbb{R}^{+}\to\mathbb{R}$ be continuously differentiable on $\mathbb{R}^{+}$. Assume that $g_{0}$ is continuously differentiable on $\mathbb{R}^{+}$ and that either

[TABLE]

or

[TABLE]

Then we have the two following cases. If (5) is assumed,

[TABLE]

where $e(u,v)=\{\rho(u)\rho(v)|W\cap W_{v-u}|\}^{-1}$ for any $u,v\in\mathbb{R}^{d}$. Instead, if (6) is assumed,

[TABLE]

We stress that the derivatives involved in Theorem 2 are derivatives with respect to $t\geq 0$. Like for Theorem 1, conditions (5) and (6) are in particular satisfied if $h$ is compactly supported in $(0,\infty)$.

2.3 Sensitivity matrix

In the next section we use empirical versions of (7) and (8) to construct estimating functions for a parametric model of an isotropic pair correlation function $g_{0}$ depending on a $K$-dimensional parameter $\boldsymbol{\beta}$, $K\geq 1$. We here investigate the expression for the associated sensitivity matrices.

Consider functions $h_{1},\ldots,h_{K}$ all fulfilling (5) and possibly depending on $\boldsymbol{\beta}$. By stacking the $K$ equations obtained by applying these functions for $h_{1},\dots,h_{K}$ in (7) we obtain the estimating function

[TABLE]

where $\mathbf{h}$ and $\mathbf{h}^{\prime}$ are vector functions with components $h_{i}$ and $h_{i}^{\prime}$. The sensitivity matrix is obtained as the expectation of the negated derivative (with respect to $\boldsymbol{\beta}$) of (9). After applying (7) once again after differentiation we obtain the sensitivity matrix

[TABLE]

Applying the Campbell theorem and converting to polar coordinates, we obtain

[TABLE]

where ${\varsigma_{d}}$ is the surface area of the $d$-dimensional unit ball. In case of (8) we obtain a similar expression,

[TABLE]

By choosing $\mathbf{h}(t)=-\psi(t)\frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\beta}}(\log g_{0})^{\prime}(t)$ for some real function $\psi$, $S(\boldsymbol{\beta})$ becomes at least positive semi-definite.

3 Estimation of log linear pair correlation function

We now consider the estimation of an isotropic pair correlation function of the form

[TABLE]

where the functions $r_{k}:\mathbb{R}^{+}\to\mathbb{R}$, $k=1,\dots,K$ are known. Following Section 2.3, the idea is to apply Theorem 2 $K$ times to functions $h_{i}$, $i=1,\ldots,K$, of the form $h_{i}(t)={-\psi(t)}\frac{\partial}{\partial\boldsymbol{\beta}_{i}}(\log g_{o})^{\prime}(t)=-\psi(t)r^{\prime}_{i}(t)$ where the function $\psi:\mathbb{R}^{+}\to\mathbb{R}$ will be justified and specified later. It is then remarkable that we obtain a simple estimating equation of the form $\mathbf{A}\boldsymbol{\beta}+\mathbf{b}=0$. The sensitivity matrix discussed in Section 2.3 is $S(\boldsymbol{\beta})=-{\mathbb{E}}\mathbf{A}$. Provided $\mathbf{A}$ is invertible we obtain the explicit solution

[TABLE]

The matrix $\mathbf{A}$ and the vector $\mathbf{b}$ are specified in the following corollary.

Corollary 1**.**

Let $\psi:\mathbb{R}^{+}\to\mathbb{R}$. Assume that $\psi$ and $r_{k}$ ($k=1,\dots,K$) are respectively continuously differentiable and twice continuously differentiable on $\mathbb{R}^{+}$. Assume either that

[TABLE]

or

[TABLE]

If (12) is assumed, we define the $(K,K)$ matrix $\mathbf{A}$ and the vector $\mathbf{b}\in\mathbb{R}^{K}$ by

[TABLE]

where again the edge effect factor is $e(u,v)=\{\rho(u)\rho(v)|W\cap W_{v-u}|\}^{-1}$ for any $u,v\in\mathbb{R}^{d}$. Instead, in case of (13), we define

[TABLE]

Then, the equation

[TABLE]

is an unbiased estimating equation.

Proof.

The proof consists in applying Theorem 2 with $h(t)=-\psi(t)r_{k}^{\prime}(t)$ for $k=1,\dots,K$ and in noticing that $(\log g_{0})^{\prime}(t)=\boldsymbol{\beta}^{\top}\mathbf{r}^{\prime}(t)=\mathbf{r}^{\prime}(t)^{\top}\boldsymbol{\beta}$. ∎

We note that if $\psi$ is compactly supported in $[0,\infty)$, then (12) or (13) are always valid assumptions. Another special case is also interesting: let $d>1$ and $\psi=1$, then (13) is true if for any $k,l=1,\dots,K$, $t\mapsto t^{d-1}r^{\prime}_{k}(t)^{2}\in L^{1}(\mathbb{R}^{d})$ and $\lim_{n\to\infty}n^{d-1}r_{k}(n)r_{l}^{\prime}(n)=0$. This simple condition is for instance satisfied if the $r_{k}$’s’ are exponential covariance functions.

The results above are e.g. applicable to the case of a pair correlation function for a log Gaussian Cox process with covariance function given by a sum of known correlation functions scaled by unknown variance parameters. Assuming a known correlation function is on the other hand quite restrictive. However, any log pair correlation function can be approximated well on a finite interval using a suitable basis function expansion so that we can effectively represent it as a log linear model. We exploit this in Section 4 where we consider the case where the functions $r_{k}$ are basis functions on a bounded real interval.

Remark 1**.**

In the more general case of a translation invariant pair correlation function of log linear form

[TABLE]

where $r_{k}(w)$, $k=1,\dots,K$, are $K$ known functions assumed to be continuously differentiable on $\mathbb{R}^{d}$, we can also obtain an estimating equation of the form (18) using Theorem 1 instead of Theorem 2. We omit the details.

Remark 2**.**

In applications of (14)-(15) for $d=2$ or (16)-(17) for $d\geq 1$ the division by $\|v-u\|^{d-1}$ or $\|v-u\|$ may lead to numerical instability for pairs of close points $u$ and $v$. This can be mitigated by a proper choice of the function $\psi$. In the spatial case of $d=2$ we propose to define $\psi(t)=(t/b)^{2}(1-(t/b))^{2}\mathbf{1}(t\in[0,b])$ for some $b>0$. With this choice of $\psi$ the divisors $\|v-u\|^{d-1}=\|v-u\|$ cancel out preventing very large or infinite variances of (14)-(17).

Remark 3**.**

The quantities (14)-(17) depend on the unknown intensity function. If the intensity function is constant equal to $\rho>0$ we can multiply (18) by $\rho^{2}$ whereby the resulting estimating equation no longer depends on $\rho$. Thus $g_{0}$ can be estimated without estimating $\rho$. Otherwise, the intensity function has to be estimated first, for instance in a parametric way, see Guan et al. (2015), and plugged into (14)-(17).

4 Variational orthogonal series estimation of the pair correlation function

In this section we consider the estimation of an isotropic pair correlation function $g_{0}$ on a bounded interval $[r_{\min},r_{\min}+R]$, $0\leq r_{\min}<\infty$ and $0<R<\infty$, using a series expansion of $\log g_{0}$. Let $\{\phi_{k}\}_{k\geq 1}$ denote an orthonormal basis of functions on $[0,R]$ with respect to some weight function $w(\cdot)\geq 0$, i.e. $\int_{0}^{R}\phi_{k}(t)\phi_{l}(t)w(t)\mathrm{d}t=\delta_{kl}$. Provided $\log g_{0}$ is square integrable (with respect to $w(\cdot)$) on $[r_{\min},r_{\min}+R]$, we have the expansion

[TABLE]

where the coefficients $\beta_{k}$ are defined by $\beta_{k}=\int_{0}^{R}g_{0}(t+r_{\min})\phi_{k}(t)w(t)\mathrm{d}t$.

We propose to approximate $\log g_{0}$ by truncating the infinite sum up to some $K\geq 1$ and obtain estimates $\hat{\beta}_{1},\ldots,\hat{\beta}_{K}$ using (18). The resulting estimate thus becomes

[TABLE]

In the sequel this estimator is referred to as the variational (orthogonal series) estimator (VSE for short). The approach is related to Zhao (2018) who also considers an estimating equation approach to estimate a pair correlation function of the form (20) but for a number $m>1$ of independent point processes on $\mathbb{R}$. The approach in Zhao (2018) further does not yield closed form expressions for the estimates of the coefficients.

Orthogonal series estimators have already been considered by Jalilian et al. (2019) who expand $g_{0}-1$ instead of $\log g_{0}$. They propose very simple unbiased estimators of the coefficients but the resulting estimator of $g_{0}$, referred to as the OSE in the sequel, is not guaranteed to be non-negative.

4.1 Implementation of the VSE

Examples of orthogonal bases include the cosine basis with $w(r)=1$, $\phi_{1}(r)=1/\surd{R}$ and $\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

[TABLE]

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)$. In the context of the variational equation (18) we need that the basis functions $\phi_{k}$ have non-zero derivatives in order to estimate $\beta_{k}$. This is not the case for $\phi_{1}$ of the cosine basis. We therefore consider in the following the Fourier-Bessel basis.

Let $b_{k}=1[k\leq K]$, $k\geq 1$. The mean integrated squared error (MISE) for $\log g_{0}$ of the VSE over the interval $[{r_{\min}},R+{r_{\min}}]$ is

[TABLE]

Jalilian et al. (2019) chose $K$ by minimizing an estimate of the MISE for $g_{0}$. We have, however, not been able to construct a useful estimate of (21). Instead we choose $K$ by maximizing a composite likelihood cross-validation criterion

[TABLE]

where $\widehat{\log g_{0,K}}^{-\{u,v\}}$ is the estimate of $\log g_{0}$ obtained using all pairs of points in $\mathbf{X}$ except $(u,v)$ and $(v,u)$. This is a simplified version of the cross-validation criterion introduced by Guan (2007a) in the context of non-parametric kernel estimation of the pair correlation function.

For computational simplicity and to guard against overfitting we choose inspired by Jalilian et al. (2019) the first local maximum of $\mathrm{CV}(K)$ larger than or equal to two rather than looking for a global maximum. Note that when $\mathbf{A}$ and $\mathbf{b}$ in (18) have been obtained for one value of $K$, then we obtain the $\mathbf{A}$ and $\mathbf{b}$ for $K+1$ by just adding one new row/column to the previous $\mathbf{A}$ and one new entry to the previous $\mathbf{b}$.

4.2 Simulation study

We study the performance of our variational estimator using simulations of point processes with constant intensity 200 on $W=[0,1]^{2}$ or $W=[0,2]^{2}$. We consider the case of a Poisson process for which the pair correlation function is constant equal to one, a Thomas process (parent intensity $\kappa=25$, dispersal standard deviation $\omega=0.0198$ and offspring intensity $\mu=8$), a variance Gamma cluster process (parent intensity $\kappa=25$, shape parameter $\nu=-1/4$, dispersion parameter $\omega=0.01845$ and offspring intensity $\mu=8$), and a determinantal point process (DPP) with exponential kernel $K(r)=\exp(-r/\alpha)$ and $\alpha=0.039$. The pair correlation functions for the four point process models are shown in Figures 2 and 3 in the usual scale as well as in the log scale. The Thomas and variance Gamma processes are clustered with pair correlation functions bigger than one while the DPP is repulsive with pair correlation function less than one. In all cases we consider $R=0.125$ and we let $r_{\min}=0$ for Poisson, Thomas, and variance Gamma. For the DPP the log pair correlation function is not well-defined for $r=0$ and we therefore use $r_{\min}=0.01$ in case of the DPP. We use (14) and (15) for computing $\mathbf{A}$ and $\mathbf{b}$ and referring to Remark 2 we let $b=r_{\min}+R$. For each point process we generate 500 simulations.

4.2.1 Estimates of coefficients

Equations (14) and (15) are derived from (7) in which $g_{0}$ is the true pair correlation function. In practice, when considering a truncated version of (20), the estimating equation (18) is not unbiased which results in bias of the coefficient estimates. This is exemplified in case of the Thomas process in the left plot of Figure 1 which shows boxplots of the first two coefficient estimates when (20) is truncated to $K=2$. In the right plot, (20) is truncated to $K=8$ which means that the truncated version of (20) is very close to the Thomas pair correlation function. Accordingly, the bias of the estimates is much reduced. However, the estimation variance increases when $K$ is increased. This emphasizes the importance of selecting an appropriate trade-off between bias an variance. The plots in Figure 1 also show how the variance of the coefficient estimates decreases when the observation window $W$ is increased from $[0,1]^{2}$ to $[0,2]^{2}$.

4.2.2 Comparison of estimators

In addition to our new VSE, we also for each simulation consider the OSE proposed by Jalilian et al. (2019) (using the Fourier-Bessel basis and their so-called simple smoothing scheme) and a standard non-parametric kernel density estimate (KDE) with bandwidth chosen by cross-validation (Guan, 2007b; Jalilian and Waagepetersen, 2018). Figures 2 and 3 depict means of the simulated OSE and VSE estimates of $g_{0}$ and $\log g_{0}$ as well as 95% pointwise envelopes. The variational estimator has larger variability whereas the bias can be smaller or larger than the OSE depending on the model.

Table 1 summarizes the root MISE (square root of (21)) for the three estimators across the four models.

The root MISEs are larger for the variational estimator than for the OSE and the KDE except in the Poisson case where the KDE has larger MISE than the VSE. Table 1 also reports the average of the selected $K$ for the variational estimator and the OSE. The averages of the selected $K$’s are pretty similar for the Poisson and DPP models while the OSE tends to select higher $K$ than the variational method for the Thomas and variance Gamma point processes.

We have also compared the computing time to evaluate the OSE and VSE. The OSE is generally cheaper except when the number of points and $R$ are large, see also the case of Capparis in Section 4.3.

4.3 Data example

We apply the three estimators to the data example considered in Jalilian et al. (2019). That is, 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 et al., 1996; Condit, 1998). The intensity functions for the point patterns are estimated as in Jalilian et al. (2019) using log-linear regression models depending on various soil and topographical variables. The estimated pair correlation functions are shown in Figure 4.

In case of Capparis and Lonchocarpus, VSE and OSE are quite similar while the kernel estimate is markedly different from the other estimates for small lags. For Acalypha, all estimates differ for small lags. The three estimates are very similar for large spatial lags for all species. The selected number $K$ for the VSE are 3, 9 and 5 for Acalypha, Capparis, and Lonchocarpus, while OSE selects $K=7$ for all species. In the case of Capparis, the computation time (4200 seconds) is higher for the OSE than for the VSE (1244 seconds) due to the high number of points for this species.

5 Discussion

In this paper we derive variational equations based on second order properties of a spatial point process. It is remarkable that in case of log-linear parametric models for the pair correlation function, it is possible to derive variational estimating equations which have closed form solutions for the unknown parameters. We exploit this to construct new variational orthogonal series type estimators for the pair correlation function. In contrast to previous kernel and orthogonal series estimators, our new estimate is guaranteed to be non-negative. For large data sets, the new estimator is further computationally faster than the previous orthogonal series estimate. However, in terms of accuracy as measured by MISE, the new estimator does not outperform the previous estimators. In the data example, the new estimator and the OSE gave similar results.

We believe there is further scope for exploring variational equations. For instance, one could use non-orthogonal bases for expanding the log pair correlation function instead of the orthogonal Fourier-Bessel basis used in this work. One might e.g. in future work consider so-called frames (Christensen, 2008) or spline bases.

**Acknowledgments

** Rasmus Waagepetersen’s and Francisco Cuevas-Pachecho’s research was 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. The research of J.-F. Coeurjolly is funded by the Natural Sciences and Engineering Research Council of Canada.

The BCI forest dynamics research project was made possible by National Science Foundation grants to Stephen P. Hubbell: DEB-0640386, DEB-0425651, DEB-0346488, DEB-0129874, DEB-00753102, DEB-9909347, DEB-9615226, DEB-9615226, DEB-9405933, DEB-9221033, DEB-9100058, DEB-8906869, DEB-8605042, DEB-8206992, DEB-7922197, support from the Center for Tropical Forest Science, the Smithsonian Tropical Research Institute, the John D. and Catherine T. MacArthur Foundation, the Mellon Foundation, the Celera Foundation, and numerous private individuals, and through the hard work of over 100 people from 10 countries over the past two decades. The plot project is part of the Center for Tropical Forest Science, a global network of large-scale demographic tree plots.

The BCI soils data set were collected and analyzed by J. Dalling, R. John, K. Harms, R. Stallard and J. Yavitt with support from NSF DEB021104, 021115, 0212284, 0212818 and OISE 0314581, STRI and CTFS. Paolo Segre and Juan Di Trani provided assistance in the field. The covariates dem, grad, mrvbf, solar and twi were computed in SAGA GIS by Tomislav Hengl (http://spatial-analyst.net/).

Appendix A Proof of Theorem 1

Proof.

Using the Campbell theorem (1) and since $\nabla\log g=(\nabla g)/g$, we start with

[TABLE]

Using first the invariance by translation of $h$ and $\nabla g$, second Fubini’s theorem, and third a change of variables, this reduces to

[TABLE]

By assumption, we have using the dominated convergence theorem,

[TABLE]

We can now use the standard trace theorem (see e.g. Evans and Gariepy (1992)) and obtain

[TABLE]

From (3), we deduce from the dominated convergence theorem that

[TABLE]

Finally, using successively a change of variable and the Campbell theorem we get

[TABLE]

which proves (4). ∎

Appendix B Proof of Theorem 2

Proof.

Both (7) and (8) are proved similarly. We focus only on (8) and follow the proof of Theorem 1. Using the Campbell theorem (1), the fact $(\log g_{0})^{\prime}=g_{0}^{\prime}/g_{0}$ and finally a change to polar coordinates, we have

[TABLE]

Using the dominated convergence theorem, partial integration and (6) we have

[TABLE]

A change to polar coordinates and the Campbell theorem again lead to

[TABLE]

∎