suggests a data-driven lower bound for $k$ in the MCCV case: sequentially
increase $k$ at least until
$\mathcal{R}_{p}(\\{({\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})\\}_{i=1}^{k})=(0,\bar{R})$.

By imposing the restriction that
$\theta=(\theta^{\prime},R)\in\Theta^{\prime}\times\mathcal{R}$, the
innovations (4.15) reduce to

$\displaystyle\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle\mathcal{I}_{p(1-p)^{-1}\lambda_{\theta}}^{f(\lambda_{\theta}(\cdot))}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
(5.4) $\displaystyle=$ $\displaystyle\sum_{x\in{\mathbf{x}}_{i}^{V}\cap
W}f\left(\frac{p\beta_{\theta^{\prime}}(x)}{1-p}\right)-\int_{W\setminus\bigcup_{x\in{\mathbf{x}}_{i}^{T}}b(x,R)}f\left(\frac{p\beta_{\theta^{\prime}}(u)}{1-p}\right)\frac{p\beta_{\theta^{\prime}}(u)}{1-p}\mathrm{d}u,$

i.e. the loss function for estimating $\theta_{0}^{\prime}$ is given by a
combination of
$\mathcal{I}_{p(1-p)^{-1}\lambda_{\theta}}^{f(\lambda_{\theta}(\cdot))}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$,
$\theta=(\theta^{\prime},R)\in\Theta^{\prime}\times\mathcal{R}$,
$i\in\mathcal{T}_{k}$. Since
$\mathcal{I}_{\xi_{\theta_{1}}^{1}}^{h_{\theta_{1}}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=\mathcal{I}_{\xi_{\theta_{2}}^{1}}^{h_{\theta_{2}}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
does not imply that $\theta_{1}=\theta_{2}$ (if
$\theta_{1}=(\theta^{\prime},R_{1})$ and $\theta_{2}=(\theta^{\prime},R_{2})$
these two innovations are the same for any $R_{1},R_{2}\in\mathcal{R}$), the
loss function
$\theta\mapsto\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
is not identifiable for a fixed $i\in\mathcal{T}_{k}$. One would typically
deal with this by fixing a point estimate $\widehat{R}$ of $R_{0}$, most
naturally
$\widehat{R}=\widehat{R}_{p}(\\{({\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})\\}_{i=1}^{k})=\sup\mathcal{R}$,
and then proceed by exploiting the innovations in (5.4) for the estimation of
$\theta^{\prime}$. However, we have seen that, numerically, this is not
necessary when employing any of the loss functions $\mathcal{L}_{1}$,
$\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ in (4.7), (4.8) and (4.9), i.e. we may
let both $R\in\mathcal{R}$ and $\theta^{\prime}\in\Theta^{\prime}$ be free
parameters to be estimated. In other words, the component-wise
unidentifiability seems to not spill over on the loss functions.

######  Remark 5.1.

Regarding the test function family considered here, note in contrast e.g. that
a test function for which $\lim_{x\to 0}f(x)=0$ will instead minimise the
squared/absolute innovations when the hard-core constraint is violated, which
is clearly not what we want here.

We next turn to the special case where
$\beta_{\theta^{\prime}}(\cdot)\equiv\beta\in\Theta^{\prime}=(0,\infty)$ and
$\beta_{\theta_{0}^{\prime}}(\cdot)=\beta_{0}\in\Theta^{\prime}$. Here the
loss function terms become

$\displaystyle\widetilde{\mathcal{I}}_{\xi_{\theta}^{1}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle\widetilde{\mathcal{I}}_{p(1-p)^{-1}\lambda_{\theta}}^{f(\lambda_{\theta}(\cdot))}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
(5.5) $\displaystyle=$ $\displaystyle
I_{i}f\left(\frac{p\beta}{1-p}\right)\left(\\#{\mathbf{x}}_{i}^{V}-\frac{p\beta}{1-p}\left|W\setminus\bigcup_{x\in{\mathbf{x}}_{i}^{T}}b(x,R)\right|\right),\quad
R\in\mathcal{R}.$

If we impose that $|f(x)|>0$, $x>0$, which e.g. holds for $f(x)=x^{-\gamma}$,
$\gamma>0$, then this is 0 if either
$I_{i}={\mathbf{1}}\\{1\leq\\#{\mathbf{x}}_{i}^{T}\leq\\#{\mathbf{x}}-1\\}=0$
or if $\beta$ is given by

$\widehat{\beta}_{i}|R=\frac{(1-p)\\#{\mathbf{x}}_{i}^{V}}{p|W\setminus\bigcup_{x\in{\mathbf{x}}_{i}^{T}}b(x,R)|}=\frac{|W|}{|W\setminus\bigcup_{x\in{\mathbf{x}}_{i}^{T}}b(x,R)|/(1-p)}\overbrace{\frac{\widetilde{\theta}({\mathbf{x}}_{i}^{V},W)}{p}}^{\approx\widetilde{\theta}({\mathbf{x}},W)},\quad
R\in\mathcal{R},$

which essentially is equivalent to a CV-based version of
$\widehat{\beta}_{PL}(R)$, $R\in\mathcal{R}$. It should be emphasised that
$|W\setminus\bigcup_{x\in{\mathbf{x}}_{i}^{T}}b(x,R)|/(1-p)$ is not linear in
$p$ (from a distributional point of view) so we expect the choice of $p$ to be
of significance here. Recalling Theorem 3, we see that the estimate obtained
by minimising either $\mathcal{L}_{1}$ in (4.7) or $\mathcal{L}_{2}$ in (4.8)
here tries to find a pair $(\beta,R)$ such that, on average (in a median sense
in the former case and in a mean sense in the latter case), we estimate
$\\#{\mathbf{x}}_{i}^{V}$ as well as possible by means of
$\widehat{\beta}_{i}|R$, $i=1,\ldots,k$.

######  Remark 5.2.

As an alternative to optimising with respect to $\beta$ and $R$ jointly, one
could consider the profile alternative where one fixes $R\in\mathcal{R}$, e.g.
$R=\sup\mathcal{R}$. When this is the case, minimising (4.7) yields the
estimate
$\widehat{\beta}=\operatorname{med}\\{\widehat{\beta}_{i}|R:i\in\mathcal{T}_{k}\\};$
cf. the estimate in expression (4.11). Similarly, the estimates obtained using
$\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ in (4.8) and (4.9) are here given by
$\widehat{\beta}=\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\widehat{\beta}_{i}|R;$
cf. the estimate in (4.12).

####  5.2.3 Numerical evaluations

We next evaluate our approach numerically in the case where
$\beta_{\theta^{\prime}}(\cdot)\equiv\beta\in\Theta^{\prime}=(0,\infty)$ and
$\beta_{\theta_{0}^{\prime}}(\cdot)=\beta_{0}\in\Theta^{\prime}$. More
specifically, we consider 100 realisations of a hard core model on
$W=[0,1]^{2}$ with parameters $R_{0}=0.05$ and $\beta_{0}=100$; this
particular choice of parameters, which give rise to an average point count of
58.51, was made completely arbitrarily.

We here consider the loss functions $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and
$\mathcal{L}_{3}$ in (4.7), (4.8) and (4.9), in combination with (5.5),
$I_{i}={\mathbf{1}}\\{1\leq\\#{\mathbf{x}}_{i}^{T}\leq\\#{\mathbf{x}}-1\\}$
and MCCV, where $k=400$ and
$p\in\\{0.05,0.07,0.1,0.15,0.2,0.3,0.5,0.7,0.9\\}$. We further focus on the
test functions $f(x)=1/x$ and $f(x)=1/\sqrt{x}$, $x\in\mathbb{R}$. To compare
our approaches to the state of the art, we additionally carry out
pseudolikelihood estimation.

In Figure 3 we report the absolute bias, the variance and the mean squared
error ($\mathrm{MSE}$) for each estimator; as anticipated, the choice of $p$
plays a significant role here. All estimators, irrespective of the choice of
$p$, yield estimated biases of $R_{0}$ which are close to 0; this is also the
the case for the pseudolikelihood estimator. More interestingly, our point
process learning approach does not require that a fixed estimate of $R_{0}$ is
plugged into (5.5) in order to obtain a good estimate of $\beta_{0}$. In the
case of $\beta_{0}$, compared to the pseudolikelihood estimator
($\mathrm{MSE}=410.26$), we can find choices for $p$ such that we always
achieve either a smaller bias ($p=0.07$ with $f(x)=1/\sqrt{x}$, using either
$\mathcal{L}_{1}$ or $\mathcal{L}_{2}$, and $p=0.1$ with $f(x)=1/x$, using
$\mathcal{L}_{2}$), variance (essentially any $p$; note that it decreases with
$p$) or $\mathrm{MSE}$. The superior performance in $\mathrm{MSE}$ with
respect to pseudolikelihood estimation holds for $f(x)=1/\sqrt{x}$ with
$p\in\\{0.05,0.07,0.1\\}$ and either of $\mathcal{L}_{1}$ or
$\mathcal{L}_{2}$; it is minimised with $p=0.07$ and $\mathcal{L}_{2}$, giving
$\mathrm{MSE}=225.14$. Similarly, when $f(x)=1/x$ this is true for
$p\in\\{0.07,0.1,0.15\\}$ and either of $\mathcal{L}_{1}$ or
$\mathcal{L}_{2}$; it is minimised with $p=0.1$ and $\mathcal{L}_{2}$,
yielding $\mathrm{MSE}=202.23$. Hence, it seems that $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$, which perform next to identically, are the preferred loss
functions here. Finally, we have observed that increasing $\beta_{0}$, while
keeping $R_{0}=0.05$ fixed, results in the $\mathrm{MSE}$-optimal $p$ being
slightly smaller; $p=0.07$ in combination with $\mathcal{L}_{2}$ still
performs better than pseudolikelihood estimation in terms of $\mathrm{MSE}$.

Figure 3:  Parameter estimation results based on 100 realisations of a hard
core model on $W=[0,1]^{2}$ with parameters $R_{0}=0.05$ and $\beta_{0}=100$.
Loss functions: $\mathcal{L}_{1}$ (solid curves, –), $\mathcal{L}_{2}$ (dashed
curves, $---$) and $\mathcal{L}_{3}$ (dotted curves, $\cdots$), in combination
with MCCV, where $p\in\\{0.05,0.07,0.1,0.15,0.2,0.3,0.5,0.7,0.9\\}$ and
$k=400$, and the test functions $f(x)=1/x$ (top row) and $f(x)=1/\sqrt{x}$
(bottom row). The red lines represent pseudolikelihood estimation. From left
to right, the plots correspond to absolute bias ($\mathrm{AB}$), variance
($\mathrm{Var}$), and mean squared error ($\mathrm{MSE}$), respectively.

###  5.3 Non-parametric kernel intensity estimation

We next turn to non-parametric intensity estimation. Recalling Section A.1, we
here consider a non-parametric intensity estimator
$\widehat{\rho}_{\theta}(u,{\mathbf{y}})$, $u\in W$,
${\mathbf{y}}\in{\mathcal{X}}$, $\theta\in\Theta$, and our aim is to choose
the tuning parameter $\theta$ in some optimal way, using an observed point
pattern ${\mathbf{x}}\subseteq W$.

From Section 4.3.3 we have that the bivariate innovations here are given by

$\displaystyle\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle\mathcal{I}_{p(1-p)^{-1}\widehat{\rho}_{\theta}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
(5.6) $\displaystyle=$ $\displaystyle\sum_{x\in{\mathbf{x}}_{i}^{V}\cap
W}h_{\theta}(x;{\mathbf{x}}_{i}^{T})-\frac{p}{1-p}\int_{W}h_{\theta}(u;{\mathbf{x}}_{i}^{T})\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})\mathrm{d}u,\quad
i=1,\ldots,k,$

and the tuning parameter selection/estimation may be carried out using these
innovations in e.g. any of the loss functions (4.7), (4.8) and (4.9). Aside
from choosing some suitable estimator $\widehat{\rho}_{\theta}$ and CV
parameters $p$ and $k$, the crucial choice to be made here is clearly the test
function. Motivated by Cronie and van Lieshout (2018), we here consider the
form
$h_{\theta}(u,{\mathbf{x}}_{i}^{T})=f(p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p))$,
for some test function $f\geq 0$. This results in

$\displaystyle\widetilde{\mathcal{I}}_{\xi_{\theta}^{1}}^{h_{\theta}}(W;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle
I_{i}\left(\sum_{x\in{\mathbf{x}}_{i}^{V}}f\left(\frac{\widehat{\rho}_{\theta}(x,{\mathbf{x}}_{i}^{T})p}{1-p}\right)-\int_{W}f\left(\frac{\widehat{\rho}_{\theta}(u,{\mathbf{x}}_{i}^{T})p}{1-p}\right)\frac{\widehat{\rho}_{\theta}(u,{\mathbf{x}}_{i}^{T})p}{1-p}\mathrm{d}u\right),$
(5.7)

where
$I_{i}={\mathbf{1}}\\{1\leq\\#{\mathbf{x}}_{i}^{V}\leq\\#{\mathbf{x}}-1\\}$,
$i=1,\ldots,k$.

####  5.3.1 Kernel intensity estimation

In what follows, we will focus on kernel estimation for point processes in
$\mathbb{R}^{d}$; see Section A.1 for details. More specifically, we will
focus on optimal selection of the bandwidth $\theta\in\Theta=(0,\infty)$,
given some point pattern ${\mathbf{x}}\subseteq W\subseteq\mathbb{R}^{d}$. Our
main objective here is to study what effect our point process learning
approach has on bandwidth selection and to do so we will use the approach of
Cronie and van Lieshout (2018) as benchmark since, generally speaking, it
outperformed its predecessors/competitors.

######  Remark 5.3.

It should be emphasised that we may normalise intensity estimators to become
density estimators, so the conclusions below apply equally well to density
estimation (in the context of generalised random sampling).

Recall from Section A.1 that Cronie and van Lieshout (2018) conjectured that
(A.3) with $f(x)=1/x$ results in (A.3) being a convex function, when using a
Gaussian kernel and no edge correction, i.e. $w_{\theta}(u,x)\equiv 1$.
Assuming that their conjecture is true, also $\mathcal{L}_{2}(\theta)$ in
combination with (5.7) would be convex under these conditions and consequently
there would be a global minimiser. Hence, we here consider a Gaussian kernel
estimator, using no edge correction when we select the bandwidth, and let the
test function be given by

$h_{\theta}(u,{\mathbf{x}}_{i}^{T})=f(p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p))=1/(p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p)).$

This puts us in a setting which is equivalent to that in Cronie and van
Lieshout (2018) and, consequently, we may explicitly analyse how our point
process learning framework improves the performance of the approach of Cronie
and van Lieshout (2018) by simply comparing the performances of the two
approaches. We further note that the above choice of test function is a
natural and convenient one for our purposes but many other choices may perform
equally well/better. Taking the observations in Section 5.1 into account, we
will also look closer at the choice

$h_{\theta}(u,{\mathbf{x}}_{i}^{T})=f(p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p))=1/\sqrt{p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p)}.$

######  Remark 5.4.

As an aside, it should be noted that the Poisson process likelihood leave-one-
out CV approach (see Section A.1 for details) can be straightforwardly altered
such that any of our CV-based approaches would replace the leave-one-out part
in it. One would maximise
$\theta\mapsto\frac{1}{k}\sum_{i=1}^{k}I_{i}(\sum_{x\in{\mathbf{x}}_{i}^{V}}\log(\widetilde{\rho}_{\theta}(x,{\mathbf{x}}_{i}^{T}))-\int_{W}\widetilde{\rho}_{\theta}(u,{\mathbf{x}}_{i}^{T})\mathrm{d}u)$
to select the bandwidth.

To evaluate the proposed setup numerically, we let the study region be given
by $W=[0,1]^{2}$ and we consider 100 realisations of each of three different
models, which constitute a subset of the models evaluated in Cronie and van
Lieshout (2018). The models considered, which represent aggregation, complete
randomness and inhibition are the following.

  * •

Recalling Section 2.4.3, we here consider a log-Gaussian Cox process (LGCP)
with random intensity $\Lambda(u)=\exp\\{Z(u)\\}$, $u=(u_{1},u_{2})\in
W=[0,1]^{2}$, where $Z$ is a Gaussian random field with mean function
$u=(u_{1},u_{2})\mapsto 10+80u_{1}$, $u\in W$, and exponential covariance
function $(u,v)\mapsto\sigma^{2}\exp\\{-r\|u-v\|_{2}\\}$, $u,v\in W$, with
$(\sigma^{2},r)=(2\log 5,50)$. Consequently, the intensity function is given
by $\rho(u)=(10+80u_{1})\operatorname{e}^{\sigma^{2}/2}$, $u\in W$, and the
expected number of points in $W$ is 250.

  * •

We consider a linear trend Poisson process (recall Section 2.4.2) with
intensity function $\rho(u)=\rho(u_{1},u_{2})=10+au_{1}$, $u\in W=[0,1]^{2}$,
where $a=480$, so that the expected number of points is given by 250.

  * •

We consider a homogeneous determinantal point process (DPP) with kernel
$(u,v)\mapsto\sigma^{2}\exp\\{-r\|u-v\|_{2}\\}$, $u,v\in W=[0,1]^{2}$, where
$(\sigma^{2},\beta)=(250,50)$; recall Section 2.4.5. By applying independent
thinning to it, using the retention probability function
$u=(u_{1},u_{2})\mapsto(10+80u_{1})/90$, $u\in W$, we obtain an inhomogeneous
determinantal point process with intensity function
$\rho(u)=\rho(u_{1},u_{2})=\sigma^{2}(10+80u_{1})/90$, $u\in W$, and expected
total point count given by 138.9.

As noted above, we let the kernel $\kappa$ be a Gaussian kernel and we let
$w_{\theta}(u,x)\equiv 1$, i.e. we use no edge correction, when we carry out
the bandwidth selection. Then, when we generate the final intensity estimates
based on the selected bandwidths, $\widehat{\theta}$, we use the local edge
correction
$w_{\widehat{\theta}}(u,x)=\int_{W}\kappa_{\widehat{\theta}}(u-x)\mathrm{d}x$.
Moreover, to measure the performance, for each model and bandwidth selection
approach we report estimates of the integrated absolute bias ($\mathrm{IAB}$),
the integrated squared bias ($\mathrm{ISB}$), the integrated variance
($\mathrm{IV}$) and the mean integrated squared error ($\mathrm{MISE}$):

$\displaystyle\mathrm{IAB}=$
$\displaystyle\int_{W}|\widehat{\mathbb{E}}[\widehat{\rho}_{\widehat{\theta}}(u,X)]-\rho(u)|\mathrm{d}u,$
$\displaystyle\mathrm{ISB}=$
$\displaystyle\int_{W}(\widehat{\mathbb{E}}[\widehat{\rho}_{\widehat{\theta}}(u,X)]-\rho(u))^{2}\mathrm{d}u,$
$\displaystyle\mathrm{IV}=$
$\displaystyle\int_{W}\widehat{\operatorname{Var}}(\widehat{\rho}_{\widehat{\theta}}(u,X))\mathrm{d}u,$
$\displaystyle\mathrm{MISE}=$ $\displaystyle\mathrm{ISB}+\mathrm{IV}.$

For a given model, these estimates are obtained by averaging with respect to
the outcomes for the 100 simulated realisations.

What we will show here is that, in terms of $\mathrm{MISE}$, our new approach
quite substantially outperforms the state of the art, which here is
represented by the approach of Cronie and van Lieshout (2018). What we
specifically do is to considerate the following settings:

  * •

For each model, we use each of the loss functions (4.7), (4.8) and (4.9) with
$n=1$, the innovations (5.6) and the test functions
$h_{\theta}(u,{\mathbf{x}}_{i}^{T})=f(p\widehat{\rho}_{\theta}(u;{\mathbf{x}}_{i}^{T})/(1-p))$,
where $f(x)=1/x$; recall (5.7). In the case of MCCV, we consider the sequence
$p=0.1,0.3,0.5,0.7,0.9$ and $k=400$, whereas in the case of multinomial CV we
evaluate the performance for the sequence $k=2,3,\ldots,10$. The numerical
output for the MCCV case can be found in Figure 4 and the output for the
multinomial CV case can be found in Figure 5.

  * •

For each model, we also consider the approach of Cronie and van Lieshout
(2018), i.e. the loss function (A.3) with $f(x)=1/x$, which has been
implemented in the function `bw.CvL` in the R package spatstat (Baddeley et
al., 2015). The numerical results can be found in Table 1.

  * •

Finally, to shed some light on the choice of test function, we evaluate the
loss function $\mathcal{L}_{2}(\theta)$ in combination with (5.7), i.e. (4.8)
with $n=1$, using the test function $f(x)=1/\sqrt{x}$ as well as MCCV with
$p=0.1,0.3,0.5,0.7,0.9$ and $k=400$. The corresponding numerical results can
be found in Figure 6.

First of all, regarding the choice $k=400$ in the MCCV case, it should be
stated here that we have observed that increasing $k$ beyond 100 essentially
has little/no effect on the chosen performance measures so our general
suggestion is to fix $k\geq 100$. Choosing $k$ is clearly a question of
computational cost and as a computationally cheaper alternative in the MCCV
case, one could instead sequentially increase $k$ and stop once the loss
function shows signs of converges (theoretically, the convergence is
guaranteed by the law of large numbers).

Most importantly, comparing Figure 4, Figure 6 and Figure 5 with Table 1, the
first thing we note is that regardless of the choice of $p$, $k$ and model,
all of the point process learning approaches outperform the approach of Cronie
and van Lieshout (2018) in terms of $\mathrm{MISE}$. Although the approach of
Cronie and van Lieshout (2018) performs slightly better in terms of bias, it
performs comparatively poorly in terms of variance, which is consequently the
reason for its higher $\mathrm{MISE}$; it is worth emphasising that it is
precisely the lower variance which ensures that the Cronie and van Lieshout
(2018) approach outperforms its predecessors (Cronie and van Lieshout, 2018,
Moradi et al., 2019), e.g. the bandwidth selection approach in (A.2). We do
however hypothesise that if $p\to 0$ in the MCCV case (possibly in combination
with $k\to\infty$), or $k\to\infty$ in the multinomial CV case (e.g. in
combination with $\mathcal{L}_{3}$), we would reach the same bias level as the
Cronie and van Lieshout (2018) approach, but still with a significantly lower
$\mathrm{MISE}$. In conclusion, the approach of Cronie and van Lieshout (2018)
is improved substantially by framing it within our point process learning
approach.

Having compared the performances of the new approaches to the state of the
art, we next look closer at the different choices we have made. First of all,
it seems that $\mathcal{L}_{3}$ favours a lower bias over a lower
variance/$\mathrm{MISE}$, whereas $\mathcal{L}_{2}$ favours the opposite;
$\mathcal{L}_{1}$ seems to offer some middle-ground between the two. We
further see in Figure 4 and Figure 6 that, in the case of MCCV,
$p\in[0.5,0.7]$ tends to be a safe/good choice, which balances the trade-off
between bias and variance, irrespective of the degree of clustering/inhibition
of the underlying model. Comparing the two associated test functions, i.e.
$f(x)=1/x$ (Figure 4) and $f(x)=1/\sqrt{x}$ (Figure 6), we see that the latter
reduces the bias with respect to the former, but at the cost of increasing the
variance and, consequently, $\mathrm{MISE}$. Moreover, it seems that the
performance of multinomial CV in terms of $\mathrm{MISE}$ is the best when
$k=2$ (see Figure 5), which is equivalent to considering MCCV with $p=0.5$ and
$k=1$. Since multinomial CV overall performs only slightly poorer than the
MCCV approach with $p=0.5$ and $k=400$, irrespective of the choice of $k$, we
draw the conclusion that the most relevant aspect of the choice of CV setting
in the bandwidth selection context is that we split the data roughly in half.

Turning to the computational aspects, we always have to compute $k$
innovations, irrespective of whether we use $k$-fold multinomial CV or MCCV.
This e.g. means that 2-fold MCCV is roughly 400 times faster than MCCV with
$k=400$ and $p=0.5$. It should be stated, however, that parallelisation may
reduce this speed difference by roughly as many times as there are cores to be
accessed. Using a standard laptop, we have experienced that for a single
pattern of moderate size the speed difference between the two approaches is
barely noticeable. However, since multinomial CV is faster than MCCV, it would
be one’s go-to method if computational aspects are the main priority, whereas
MCCV is the go-to method if precision is what one is going for.

The fact that in the case of clustered models the performance in terms of MISE
is improved when the training sets have few points, i.e. when $p$ is large,
corresponds precisely with observations made in Moradi et al. (2019). We
hypothesise that the reason for clustered processes requiring more thinning is
that here the points need to be sufficiently sparse for the approach to
capture the general trends of the underlying intensity.

Figure 4:  Performance of the loss functions $\mathcal{L}_{1}$ (solid curves,
–), $\mathcal{L}_{2}$ (dashed curves, $---$) and $\mathcal{L}_{3}$ (dotted
curves, $\cdots$), using MCCV with $p=0.1,0.3,0.4,0.7,0.9$ and $k=400$
together with the test function $f(x)=1/x$. Columns: LGCP (left), Poisson
(middle) and DPP (right). Top row: $\mathrm{IAB}$ (grey curve, right axis) and
$\mathrm{IV}$ (black curve, left axis). Bottom row: $\mathrm{ISB}$ (grey
curve, right axis) and $\mathrm{MISE}$ (black curve, left axis).

Figure 5:  Performance of the loss functions $\mathcal{L}_{1}$ (solid curves,
–), $\mathcal{L}_{2}$ (dashed curves, $---$) and $\mathcal{L}_{3}$ (dotted
curves, $\cdots$), using multinomial CV with $k=2,3,\ldots,10$ and the test
function $f(x)=1/x$. Columns: LGCP (left), Poisson (middle) and DPP (right).
Top row: $\mathrm{IAB}$ (grey curve, right axis) and $\mathrm{IV}$ (black
curve, left axis). Bottom row: $\mathrm{ISB}$ (grey curve, right axis) and
$\mathrm{MISE}$ (black curve, left axis).

Figure 6:  Performance of the loss function $\mathcal{L}_{2}$, using MCCV with $p=0.1,0.3,0.4,0.7,0.9$ and $k=400$ together with the test function $f(x)=1/\sqrt{x}$. Columns: LGCP (left), Poisson (middle) and DPP (right). Top row: $\mathrm{IAB}$ (grey curve, right axis) and $\mathrm{IV}$ (black curve, left axis). Bottom row: $\mathrm{ISB}$ (grey curve, right axis) and $\mathrm{MISE}$ (black curve, left axis).  Table 1: The numerical performance of the method of Cronie and van Lieshout (2018), i.e. (A.3) in combination with $f(x)=1/x$, for each of the models above. | IAB | ISB | IV | MISE
---|---|---|---|---
LGCP | 19.48 | 963.47 | 17597.99 | 18561.47
Poisson | 15.80 | 921.82 | 4408.21 | 5330.04
DPP | 9.14 | 276.75 | 2002.55 | 2279.31

####  5.3.2 Real data example

We next consider an application of our point process learning-based kernel
intensity estimation approach, in the context of a point pattern
${\mathbf{x}}\subseteq W\subseteq\mathbb{R}^{2}$. More precisely, we consider
MCCV with $(k,p)=(400,0.7)$, the test function $f(x)=1/x$ and the loss
function $\mathcal{L}_{2}$. Figure 7 shows the estimated intensity for the
dataset, which consists of the locations of 3605 tropical rain forest trees of
the species Beilschmiedia pendula (Lauraceae), sampled on Barro Colorado
Island, Panama (Hubbell and Foster, 1983). This dataset, which is accessible
through the R package spatstat (Baddeley et al., 2015), has previously been
analysed in various papers (see e.g. Møller and Waagepetersen, 2007). We
obtain a bandwidth of $56.65$ (meters) and from Figure 7 we can see that the
selected bandwidth leads to an estimated intensity which adapts well to the
inhomogeneity of the spatial locations of the trees.

Figure 7: Kernel intensity estimate for a dataset consisting of the locations
of tropical rain forest trees on Barro Colorado Island, Panama.

###  5.4 Fixed size samples

For completeness, we here look closer at the setting in Section 2.4.1, i.e.
the case where we condition on the total point count $X(S)=N\geq 1$.

We here have the parametric intensity function family
$\rho_{\theta}(u)=Nf_{1}(u;\theta)$, $u\in S$, where $f_{1}(\cdot;\theta)$,
$\theta\in\Theta$, is a parametric family of density functions and
$f_{1}(\cdot;\theta_{0})$ is the true common marginal density function of
$x_{1},\ldots,x_{N}$. From (4.13) we obtain the innovations

$\displaystyle\mathcal{I}_{(1-p)\rho_{\theta}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{T})=$
$\displaystyle\sum_{x\in{\mathbf{x}}_{i}^{T}}h_{\theta}(x)-(1-p)\int_{S}h_{\theta}(u)\rho_{\theta}(u)\mathrm{d}u$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}{\mathbf{1}}\\{x_{i}\in{\mathbf{x}}_{i}^{T}\\}h_{\theta}(x)-(1-p)N\int_{S}h_{\theta}(u)f_{1}(u;\theta)\mathrm{d}u;$

we have to require that the test functions are such that the expectation
$\int_{S}h_{\theta}(u)f_{1}(u;\theta)\mathrm{d}u$ is finite. In particular, by
setting
$h_{\theta}(\cdot)=((1-p)Nf_{1}(\cdot;\theta))^{-1}\nabla_{\theta}f_{1}(\cdot;\theta)$
in this unbiased estimating equation (provided that the gradient
$\nabla_{\theta}f_{1}(\cdot;\theta)$ exists), since $f_{1}(\cdot;\theta)$ is a
density function, we obtain

$\displaystyle\mathcal{I}_{(1-p)\rho_{\theta}}^{\nabla_{\theta}f_{1}(\cdot;\theta)/((1-p)Nf_{1}(\cdot;\theta))}(S;{\mathbf{x}}_{i}^{T})=$
$\displaystyle\frac{1}{(1-p)N}\sum_{i=1}^{N}{\mathbf{1}}\\{x_{i}\in{\mathbf{x}}_{i}^{T}\\}\frac{\nabla_{\theta}f_{1}(x_{i};\theta)}{f_{1}(x_{i};\theta)}-\nabla_{\theta}\int_{S}f_{1}(u;\theta)\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\frac{{\mathbf{1}}\\{x_{i}\in{\mathbf{x}}_{i}^{T}\\}}{1-p}\frac{\nabla_{\theta}f_{1}(x_{i};\theta)}{f_{1}(x_{i};\theta)}-1,$

assuming that integration and differentiation may be interchanged. This is
connected to the likelihood score function based on the sample
${\mathbf{x}}_{i}^{T}$. Hence, we obtain a form of subsampled maximum
likelihood estimation approach. Note that we here have not imposed the
assumption that $x_{1},\ldots,x_{N}$ are independent, i.e. that we are
considering a Binomial point process.

We next consider the case where we want to fit a multivariate density
$f_{N}(u_{1},\ldots,u_{N};\theta)$, $u_{1},\ldots,u_{N}\in S$,
$\theta\in\Theta$, to $x_{1},\ldots,x_{N}$, where we assume that the true
joint density of $x_{1},\ldots,x_{N}$ is given by
$f_{N}(u_{1},\ldots,u_{N};\theta_{0})$, $u_{1},\ldots,u_{N}\in S$. Recall from
Section 2.4.1 the associated Papangelou conditional intensity
$\lambda_{\theta}(u,{\mathbf{x}})=nf_{1}(u|{\mathbf{x}};\theta)$, $u\in S$,
${\mathbf{x}}=\\{x_{1},\ldots,x_{n-1}\\}\subseteq S^{n-1}$, $1\leq n\leq N-1$,
where $n=1$ yields
$f_{1}(u|{\mathbf{x}};\theta)=f_{1}(u|\emptyset;\theta)=f_{1}(u;\theta)$. Here
the innovations in (4.15) become

$\displaystyle\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle\mathcal{I}_{p(1-p)^{-1}\lambda_{\theta}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$
$\displaystyle=$
$\displaystyle\sum_{x\in{\mathbf{x}}_{i}^{V}}h_{\theta}(x;{\mathbf{x}}_{i}^{T})-\frac{p}{1-p}\int_{S}h_{\theta}(u;{\mathbf{x}}_{i}^{T})\lambda_{\theta}(u;{\mathbf{x}}_{i}^{T})\mathrm{d}u$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{N}{\mathbf{1}}\\{x_{j}\in{\mathbf{x}}_{i}^{V}\\}h_{\theta}(x_{j};{\mathbf{x}}_{i}^{T})-\frac{p}{1-p}(\\#{\mathbf{x}}_{i}^{T}+1)\int_{S}h_{\theta}(u;{\mathbf{x}}_{i}^{T})f_{1}(u|{\mathbf{x}}_{i}^{T};\theta)\mathrm{d}u,$

and we use
$\widetilde{\mathcal{I}}_{\xi_{\theta}^{1}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=I_{i}\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$,
where
$I_{i}={\mathbf{1}}\\{1\leq\\#{\mathbf{x}}_{i}^{T}\leq\\#{\mathbf{x}}-1\\}={\mathbf{1}}\\{0<\\#{\mathbf{x}}_{i}^{T}<N\\}$,
to build the loss function. Arguing as in the marginal case, if we set
$h_{\theta}(\cdot;{\mathbf{x}}_{i}^{T})=\nabla_{\theta}f_{1}(\cdot|{\mathbf{x}}_{i}^{T};\theta)/f_{1}(\cdot|{\mathbf{x}}_{i}^{T};\theta)p$,
we obtain a subsampled conditional likelihood estimation approach;

$\displaystyle\mathcal{I}_{\xi_{\theta}^{1}}^{h_{\theta}}(S;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})=$
$\displaystyle\sum_{j=1}^{N}\frac{{\mathbf{1}}\\{x_{j}\in{\mathbf{x}}_{i}^{V}\\}}{p}\frac{\nabla_{\theta}f_{1}(x_{j}|{\mathbf{x}}_{i}^{T};\theta)}{f_{1}(x_{j}|{\mathbf{x}}_{i}^{T};\theta)}-\frac{(\\#{\mathbf{x}}_{i}^{T}+1)}{1-p}.$

Note that, as we have seen in the previous sections, by changing the chosen
test function we also change the emphasis put on either the bias or the
variance of the resulting estimators. Hence, we may obtain estimators which
behave quite differently than likelihood-type estimators.

##  6 Discussion

In this paper we present a general theory for statistical learning for point
processes, point process learning, which essentially consists of two
components. The first is what we refer to as bivariate innovations, which
essentially are measures of discrepancy between two point processes. These are
used to “predict” one point process from another point process and they are
motivated by integral formulas such as Campbell, Campbell-Mecke and Georgii-
Nguyen-Zessin (GNZ) formulas, which relate expectations of sums over the
points of a point process to integrals with respect to various distributional
characteristics of the point process, including factorial moment measures and
Papangelou conditional intensity functions. A particular instance of our
bivariate innovations reduce to the “classical” innovations of Baddeley et al.
(2005, 2008), which were introduced to define point process residuals. The
second component of our statistical learning approach is thinning-based cross-
validation for point processes. Due to the many appealing distributional
properties of independent thinnings, which are obtained by independently
retaining the points of a point process according to some given probability
function, we study the case of independent thinning-based cross-validation in
detail. By using the bivariate innovations to predict cross-validation-
generated validation sets, using the corresponding training sets, we connect
the two components to form our (supervised) statistical learning approach.
Through a range of applications, including kernel intensity estimation,
parametric intensity estimation and Papangelou conditional intensity fitting,
we show that our learning approach generally performs better than the
classical statistical state of the art for point processes.

Our hope is that this paper becomes the starting point for a research
direction which deals with various aspects of statistical learning for point
processes. Aside from the range of future/ongoing work topics indicated here
and there in this paper, we are currently also looking into further topics. To
begin with, we are studying additional theoretical properties of bivariate
innovations, in particular asymptotic properties, where there are clear
connections to recent work in stochastic geometry (see e.g. Biscio et al.,
2018, Biscio and Waagepetersen, 2019, Błaszczyszyn et al., 2019, Biscio et
al., 2020). Moreover, a natural extension of our work in this paper is to
apply our point process learning framework to marked point processes, since
these enable regression-type estimation in the context of dependent samples.
Further ongoing work on extensions is mentioned in Section 6.2 below; the
sections below discuss various alternatives and extensions for the framework
developed in this paper.

###  6.1 Alternative cross-validation approaches

As a sort of mix between the cross-validation (CV) methods in Definition 4.2
and Definition 4.3, one could consider letting ${\mathbf{x}}_{i}^{V}$ be a
$p$-thinning of ${\mathbf{x}}$ with retention probability $p_{i}=i/k$ and
letting ${\mathbf{x}}_{i}^{T}={\mathbf{x}}\setminus{\mathbf{x}}_{i}^{V}$,
$i=1,\ldots,k\geq 1$. On the one hand, the possible advantage over MCCV is
that we only have to choose the parameter $k$. On the other hand, it is not
likely that it would perform better than an "optimally" chosen pair $(k,p)$
for MCCV. A further variant of the CV method in Definition 4.2 is to consider
an "empirical Bayes" type CV approach, where the retention probability used
would be estimated non-parametrically, e.g. by means of a scaled intensity
estimate, using (a part of) the data. The issue here is that we do not
actually employ independent thinning to generate the training and validation
sets. It is not clear whether there are any actual benefits of doing this, but
this may be worth exploring further.

Another natural and intuitive approach is what we call domain partitioning CV.
The idea is to split the study region $W$ into $k\geq 1$ pieces $W_{i}$,
$i=1,\ldots,k$, and let the $i$th validation set be
${\mathbf{x}}_{i}^{V}={\mathbf{x}}\cap W_{i}$; see the left panel of Figure 8
for an illustration. This may be obtained through a sequence of $k$
independent thinnings, where in the $i$th run we use the retention probability
$p_{i}(u)={\mathbf{1}}\\{u\in W_{i}\\}$, $u\in W$, $i=1,\ldots,k$. Note the
philosophical difference between domain partitioning CV and independent
thinning based CV strategies. In the former we "block", or sample "from left
to right/top to bottom", since we first split $W$ and then generate the
validation sets. In the latter, however, we sample "from above"; see the right
panel of Figure 8.

Figure 8:  A point pattern ${\mathbf{x}}$ in a study region
$W\subseteq\mathbb{R}^{2}$. Left: A partition $\\{W_{i}\\}_{i=1}^{k}$ of $W$,
which in turn yields the validation sets
${\mathbf{x}}_{i}^{V}={\mathbf{x}}\cap W_{i}$, $i=1,\ldots,k\geq 2$. Right: A
validation set ${\mathbf{x}}_{i}^{V}$ (illustrated by $\Delta$) which is given
by a $p$-thinning of ${\mathbf{x}}$.

However natural and intuitive domain partitioning CV may seem, there are
problems related to our point process learning approach. The main problem is
related to the retention probability for ${\mathbf{x}}_{i}^{V}$ being 1 and
the retention probability for ${\mathbf{x}}_{i}^{T}$ being 0. If
$\Xi_{\Theta}^{n}=\\{\xi_{\theta}^{n}:\theta\in\Theta\\}$ and
$\mathcal{H}_{\Theta}=\\{h_{\theta}:\theta\in\Theta\\}$ are of the form (3.1),
then (3.9) tells us that we should have
$\xi(u_{1},\ldots,u_{n};X_{i}^{T})=\prod_{j=1}^{n}p(u_{j})(\prod_{j=1}^{n}(1-p(u_{j})))^{-1}\lambda_{Y}^{(n)}(u_{1},\ldots,u_{n};X_{i}^{T})={\mathbf{1}}\\{u_{1},\ldots,u_{n}\in
W\\}{\mathbf{1}}\\{u_{1},\ldots,u_{n}\notin
W\\}^{-1}\lambda_{Y}^{(n)}(u_{1},\ldots,u_{n};X_{i}^{T})$ a.e., which is
infinite if any $u_{j}$ lies in $W_{i}$ and 0 if any $u_{j}$ lies outside
$W_{i}$; this would mean the our loss functions would be infinite. Moreover,
what we here are dealing with is "extrapolation/interpolation", i.e.
prediction outside the observation domain which is generally a hard task. One
may also speculate about whether there emerge other problems here; the folds
may have highly varying point counts and such an approach may introduce
edge/boundary effects which one may have to correct for (Cronie and Särkkä,
2011, Chiu et al., 2013, Diggle, 2014, Baddeley et al., 2015). Finally, it is
not clear how exactly the partitioning of $W$ should be chosen; how many folds
there should be, whether they should be equally sized etc.

###  6.2 Optimal loss functions, test functions and hyperparameters

A burning question is how one should choose the test function family and the
CV parameter $p$ here. The current paper has taken a proof-of-concept
approach, with the aim of showing that there exist test functions with
associated values for $p$ which result in (substantially) improved fits in
different contexts. As previously indicated, these may be viewed as
hyperparameters and, quite naturally, one may ask how a hyperparameter may be
chosen "optimally".

Turning to the literature, in the context of univariate innovation-based
parametric intensity estimation (see Section A.2), more precisely composite
likelihood estimation, Guan et al. (2015) exploited estimating equation theory
to find an optimal test function in the form of a Fredholm integral equation.
Further, Coeurjolly et al. (2016) gave an extensive account on test function
selection in the context of using classical innovations for the estimation
(see Section A.3), so-called Takacs-Fiksel estimation. In addition, motivated
by Guan et al. (2015), they found a similar Fredholm integral equation
approach to finding a (semi-)optimal test function. For our purposes, however,
since our point process learning framework (in general) cannot be expressed
through estimating equations, we believe these findings to be of limited use
here.

A further idea for finding an optimal test function is to apply calculus of
variations to find a minimiser of the variances in Theorem 2. Aside from the
possible associated mathematical intractability of such an approach, the
optimality may be model specific and the solution may not solve the issue with
the remaining hyperparameters. Arguing from the point of view of Taylor
expansion, another potential idea would be to let the test function
$h_{\theta}$ be given by a combination of derivatives of $\xi_{\theta}^{n}$.

####  6.2.1 Data-driven hyperparameter selection

A further approach, which we are currently exploring in a parallel paper,
embraces the statistical learning paradigm to a larger extent than the
approaches indicated above. It is (arguably) also more general in the sense
that it deals with hyperparameters in general, and thereby test functions and
CV parameters in particular. Its algorithm, which is laid out below, is
motivated by an algorithm commonly encountered in the statistical learning
literature (James et al., 2013). Our aim here is to have a completely data-
driven approach with a minimal amount of hard-coded parameters. Roughly
speaking, it states that out of a class of potential loss functions, the one
which (in some suitable sense) best predicts the validation sets
${\mathbf{x}}_{i}^{V}$, $i=1,\ldots,k$, will be labelled the optimal one;
clearly, in some cases this is equivalent to finding an optimal family of test
functions $\mathcal{H}_{\Theta}$.

Consider a collection of loss functions of either of the forms

$\displaystyle\\{\mathcal{L}_{\gamma}(\theta;{\mathbf{z}},{\mathbf{y}}):{\mathbf{z}},{\mathbf{y}}\in{\mathcal{X}},\theta\in\Theta\\},$
$\displaystyle\\{\mathcal{L}_{\gamma}(\theta;{\mathbf{y}}):{\mathbf{y}}\in{\mathcal{X}},\theta\in\Theta\\},$

where $\gamma\in\Theta_{\gamma}\subseteq\mathbb{R}^{l_{\gamma}}$,
$l_{\gamma}\geq 1$, is the associated hyper-parameter. In addition, consider
some measure $\mathcal{G}:\Theta\times{\mathcal{X}}\to[0,\infty)$ of goodness
of fit, where small means a good fit. Most naturally here,
$\mathcal{L}_{\gamma}(\theta;{\mathbf{z}},{\mathbf{y}})=|\widetilde{\mathcal{I}}_{\xi_{\theta}^{n}}^{h_{\theta}^{\gamma}}(W^{n};{\mathbf{z}},{\mathbf{y}})|$
if $\xi_{\theta}^{n}$ and $h_{\theta}$ are of the form (3.1) or
$\mathcal{L}_{\gamma}(\theta;{\mathbf{y}})=|\widetilde{\mathcal{I}}_{\xi_{\theta}^{n}}^{h_{\theta}^{\gamma}}(W^{n};{\mathbf{y}})|$
if $\xi_{\theta}^{n}$ and $h_{\theta}$ are of the form (3.2), where
$\mathcal{H}_{\Theta}^{\gamma}=\\{h_{\theta}^{\gamma}:\theta\in\Theta\\}$,
$\gamma\in\Theta_{\gamma}$, are different families of test functions, say
$h_{\theta}^{\gamma}(\cdot)=\xi_{\theta}^{n}(\cdot)^{\gamma}$,
$\gamma\in\Theta_{\gamma}\subseteq\mathbb{R}$. However, $\gamma$ can also e.g.
incorporate the CV parameter $p$, if the algorithm below is run for different
CV settings. Similarly, as goodness of fit criterion we may e.g. consider a
univariate innovation-based one,
$\mathcal{G}(\theta;{\mathbf{y}})=|\mathcal{I}_{\xi_{\theta}^{n}}^{h_{E}}(W^{n};{\mathbf{y}})|^{i}$,
$i\in\\{1,2\\}$, for some fixed test function $h_{E}(\cdot)$.

######  Remark 6.1.

In a parallel paper we exploit the current data-driven approach in the context
of penalisation, most notably regularisation (James et al., 2013, Yue and Loh,
2015, Rajala et al., 2018, Choiruddin et al., 2018). Here we let
$\mathcal{L}_{\gamma}(\theta;{\mathbf{z}},{\mathbf{y}})=|\widetilde{\mathcal{I}}_{\xi_{\theta}^{n}}^{h_{\theta}}(W;{\mathbf{z}},{\mathbf{y}})|+\gamma
R(\theta)$ or
$\mathcal{L}_{\gamma}(\theta;{\mathbf{y}})=|\widetilde{\mathcal{I}}_{\xi_{\theta}^{n}}^{h_{\theta}}(W;{\mathbf{y}})|+\gamma
R(\theta)$, where $\mathcal{H}_{\Theta}=\\{h_{\theta}:\theta\in\Theta\\}$ is a
fixed family of test functions, $\gamma\geq 0$ and $R(\cdot)$ is a suitable
penalty, e.g. a regularisation penalty such as the elastic-net penalty (Zou
and Hastie, 2005). Another penalisation setting which is parallely being
looked into is smoothness penalisation in the context of non-parametric
intensity estimation; specifically, the case of soap film smoothing (Wood et
al., 2008) is currently being studied.

The general approach described in the algorithm below illustrates how an
“optimal” $\mathcal{L}_{\gamma}$ may be found in a data-driven way.

  1. 1.

When the loss functions are of the form
$\mathcal{L}_{\gamma}(\theta;{\mathbf{z}},{\mathbf{y}})$, as in Section 4.1.2,
consider triples
$({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{E})$,
$i=1,\ldots,k$, which have been obtained through two layers of CV-partitioning
of ${\mathbf{x}}$; here ${\mathbf{x}}_{i}^{E}$, $i=1,\ldots,k$, are generated
by a retention probability $p_{E}\in(0,1)$ and
$({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})$,
${\mathbf{x}}_{i}^{T}\cup{\mathbf{x}}_{i}^{V}={\mathbf{x}}\setminus{\mathbf{x}}_{i}^{E}$,
$i=1,\ldots,k$, are generated by a retention probability $p\in(0,1)$, both
using some of the CV-approaches in Definition 4.2 and Definition 4.3.

When the loss functions are of the form
$\mathcal{L}_{\gamma}(\theta;{\mathbf{y}})$, we instead consider CV-splittings
$({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})$,
${\mathbf{x}}_{i}^{V}={\mathbf{x}}\setminus{\mathbf{x}}_{i}^{T}$,
$i=1,\ldots,k$, using one of the CV-approaches in Definition 4.2 and 4.3 with
some retention probability $p\in(0,1)$.

  2. 2.

For each $i=1,\ldots,k$:

    1. (a)

For each $\gamma\in\Theta_{\gamma}$, find a minimiser
$\widehat{\theta}_{i}(\gamma)\in\Theta$ of
$\mathcal{L}_{\gamma}(\theta;({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V}))$ or
$\mathcal{L}_{\gamma}(\theta;{\mathbf{x}}_{i}^{T})$.

    2. (b)

Consider the $i$th "prediction error path"
$\mathcal{G}_{i}(\gamma)=\mathcal{G}(\widehat{\theta}_{i}(\gamma);{\mathbf{x}}_{i}^{E})$,
$\gamma\in\Theta_{\gamma}$, or
$\mathcal{G}_{i}(\gamma)=\mathcal{G}(\widehat{\theta}_{i}(\gamma);{\mathbf{x}}_{i}^{V})$,
$\gamma\in\Theta_{\gamma}$ (when $\mathcal{G}$ is innovation-based, in keeping
with Baddeley et al. (2005), these should be called residuals).

  3. 3.

Let $\widehat{\gamma}$ be a minimiser of
$\gamma\mapsto\frac{1}{k}\sum_{i=1}^{k}\mathcal{G}_{i}(\gamma)$,
$\gamma\in\Theta_{\gamma}$.

  4. 4.

The optimal loss function is given by
$\mathcal{L}_{\widehat{\gamma}}(\theta;{\mathbf{z}},{\mathbf{y}})$,
${\mathbf{z}},{\mathbf{y}}\in{\mathcal{X}}$, $\theta\in\Theta$, or
$\mathcal{L}_{\widehat{\gamma}}(\theta;{\mathbf{y}})$,
${\mathbf{y}}\in{\mathcal{X}}$, $\theta\in\Theta$.

  5. 5.

Generate a new collection of training-validation splittings
$({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})$, $i=1,\ldots,k$, of the full
dataset ${\mathbf{x}}$ (possibly using a different $k\geq 1$ and a different
retention probability $p\in(0,1)$) and consider some test function generated
by a combination of
$\mathcal{L}_{\widehat{\gamma}}(\theta;{\mathbf{x}}_{i}^{T})$, $i=1,\ldots,k$,
or
$\mathcal{L}_{\widehat{\gamma}}(\theta;{\mathbf{x}}_{i}^{V},{\mathbf{x}}_{i}^{T})$,
$i=1,\ldots,k$, e.g. as was done in Section 4.2.1, or, in the former case, the
($p\to 0$) limiting case,
$\mathcal{L}_{\widehat{\gamma}}(\theta;{\mathbf{x}})$. We obtain the final
estimate $\widehat{\theta}$ by minimising this loss function.

A further option is to (additionally) let the choice of $\xi_{\theta}^{n}$
vary with $\gamma$ as well, e.g. if we want to compare different models which
are represented by different families $\Xi_{\Theta}^{n}$ – in some sense we
here would jointly carry out fitting, model selection and model validation. In
the context of large and complex datasets, one could e.g. think of hybrid
Gibbs processes (Baddeley et al., 2013) or pairwise interaction processes (van
Lieshout, 2000) with an increasing number of components.

##  7 Acknowledgements

We are grateful to Achmad Choiruddin (Institut Teknologi Sepuluh Nopember,
Indonesia) for feedback on an early draft and to Mathew Penrose (University of
Bath, UK) for supplying us with details on the proof of Last and Penrose
(2017, Exercise 5.9).

## References

  * (1)
  * Anderes et al. (2020) Anderes, E., Møller, J., Rasmussen, J. G. et al. (2020).  Isotropic covariance functions on graphs and their edges, Annals of Statistics 48(4): 2478–2503.
  * Ang et al. (2012) Ang, Q. W., Baddeley, A. and Nair, G. (2012).  Geometrically corrected second order analysis of events on a linear network, with applications to ecology and criminology, Scandinavian Journal of Statistics 39(4): 591–617.
  * Arlot and Celisse (2010) Arlot, S. and Celisse, A. (2010).  A survey of cross-validation procedures for model selection, Statistics surveys 4: 40–79.
  * Babu and Feigelson (1996) Babu, G. J. and Feigelson, E. D. (1996).  Spatial point processes in astronomy, Journal of Statistical Planning and Inference 50(3): 311–326.
  * Baccelli et al. (2020) Baccelli, F., Blaszczyszyn, B. and Karray, M. (2020).  Random Measures, Point Processes, and Stochastic Geometry, Inria.
  * Baccelli and Brémaud (2013) Baccelli, F. and Brémaud, P. (2013).  Elements of queueing theory: Palm Martingale calculus and stochastic recurrences, Vol. 26, Springer Science & Business Media.
  * Baddeley et al. (2008) Baddeley, A., Møller, J. and Pakes, A. G. (2008).  Properties of residuals for spatial point processes, Annals of the Institute of Statistical Mathematics 60(3): 627–649.
  * Baddeley et al. (2000) Baddeley, A., Møller, J. and Waagepetersen, R. (2000).  Non-and semi-parametric estimation of interaction in inhomogeneous point patterns, Statistica Neerlandica 54(3): 329–350.
  * Baddeley et al. (2015) Baddeley, A., Rubak, E. and Turner, R. (2015).  Spatial Point Patterns: Methodology and Applications with R, CRC Press.
  * Baddeley et al. (2013) Baddeley, A., Turner, R., Mateu, J. and Bevan, A. (2013).  Hybrids of gibbs point process models and their implementation, Journal of Statistical Software 55(11): 1–43.
  * Baddeley et al. (2005) Baddeley, A., Turner, R., Møller, J. and Hazelton, M. (2005).  Residual analysis for spatial point processes (with discussion), Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67(5): 617–666.
  * Barr and Schoenberg (2010) Barr, C. D. and Schoenberg, F. P. (2010).  On the Voronoi estimator for the intensity of an inhomogeneous planar poisson process, Biometrika 97(4): 977–984.
  * Bayisa et al. (2020) Bayisa, F. L., Ådahl, M., Rydén, P. and Cronie, O. (2020).  Large-scale modelling and forecasting of ambulance calls in northern sweden using spatio-temporal log-gaussian cox processes, Spatial Statistics 39: 100471.
  * Beneš and Rataj (2004) Beneš, V. and Rataj, J. (2004).  Stochastic Geometry: Selected Topics, Kluwer Academic Publishers.
  * Besag (1974) Besag, J. (1974).  Spatial interaction and the statistical analysis of lattice systems, Journal of the Royal Statistical Society: Series B (Methodological) 36(2): 192–225.
  * Biscio et al. (2020) Biscio, C. A. N., Chenavier, N., Hirsch, C. and Svane, A. M. (2020).  Testing goodness of fit for point processes via topological data analysis, Electronic Journal of Statistics 14(1): 1024 – 1074.
  * Biscio et al. (2018) Biscio, C. A. N., Poinas, A. and Waagepetersen, R. (2018).  A note on gaps in proofs of central limit theorems, Statistics & Probability Letters 135: 7–10.
  * Biscio and Waagepetersen (2019) Biscio, C. A. N. and Waagepetersen, R. (2019).  A general central limit theorem and a subsampling variance estimator for $\alpha$-mixing point processes, Scandinavian Journal of Statistics 46(4): 1168–1190.
  * Breiman (2001) Breiman, L. (2001).  Statistical modeling: The two cultures (with comments and a rejoinder by the author), Statistical science 16(3): 199–231.
  * Brémaud (1981) Brémaud, P. (1981).  Point processes and queues: martingale dynamics, Vol. 50, Springer.
  * Błaszczyszyn et al. (2019) Błaszczyszyn, B., Yogeshwaran, D. and Yukich, J. E. (2019).  Limit theory for geometric statistics of point processes having fast decay of correlations, The Annals of Probability 47(2): 835 – 895\.
  * Chaiban et al. (2019) Chaiban, C., Biscio, C., Thanapongtharm, W., Tildesley, M., Xiao, X., Robinson, T. P., Vanwambeke, S. O. and Gilbert, M. (2019).  Point pattern simulation modelling of extensive and intensive chicken farming in thailand: Accounting for clustering and landscape characteristics, Agricultural Systems 173: 335–344.
  * Chaudhuri et al. (2021) Chaudhuri, S., Moradi, M. and Mateu, J. (2021).  On the trend detection of time-ordered intensity images of point processes on linear networks, Communication in Statistics- Simulation and Computation .
  * Chiu et al. (2013) Chiu, S. N., Stoyan, D., Kendall, W. S. and Mecke, J. (2013).  Stochastic Geometry and its Applications, John Wiley & Sons.
  * Choiruddin et al. (2018) Choiruddin, A., Coeurjolly, J.-F. and Letué, F. (2018).  Convex and non-convex regularization methods for spatial point processes intensity estimation, Electronic Journal of Statistics 12(1): 1210–1255.
  * Christensen (2019) Christensen, R. (2019).  Advanced linear modeling: statistical learning and dependent data, Springer Nature.
  * Coeurjolly et al. (2016) Coeurjolly, J.-F., Guan, Y., Khanmohammadi, M. and Waagepetersen, R. (2016).  Towards optimal takacs–fiksel estimation, Spatial Statistics 18: 396–411.
  * Coeurjolly and Lavancier (2019) Coeurjolly, J.-F. and Lavancier, F. (2019).  Understanding spatial point patterns through intensity and conditional intensities, in D. Coupier (ed.), Stochastic Geometry, Lecture Notes in Mathematics, vol 2237, Springer, pp. 45–85.
  * Coeurjolly et al. (2017) Coeurjolly, J.-F., Møller, J. and Waagepetersen, R. (2017).  A tutorial on palm distributions for spatial point processes, International Statistical Review 85(3): 404–420.
  * Coeurjolly and Rubak (2013) Coeurjolly, J.-F. and Rubak, E. (2013).  Fast covariance estimation for innovations computed from a spatial gibbs point process, Scandinavian Journal of Statistics 40(4): 669–684.
  * Cronie et al. (2020) Cronie, O., Moradi, M. and Mateu, J. (2020).  Inhomogeneous higher-order summary statistics for point processes on linear networks, Statistics and computing 30(5): 1221–1239.
  * Cronie et al. (2013) Cronie, O., Nyström, K. and Yu, J. (2013).  Spatiotemporal modeling of swedish scots pine stands, Forest Science 59(5): 505–516.
  * Cronie and Särkkä (2011) Cronie, O. and Särkkä, A. (2011).  Some edge correction methods for marked spatio-temporal point process models, Computational Statistics & Data Analysis 55(7): 2209–2220.
  * Cronie and van Lieshout (2015) Cronie, O. and van Lieshout, M. N. M. (2015).  A J-function for inhomogeneous spatio-temporal point processes, Scandinavian Journal of Statistics 42: 562–579.
  * Cronie and van Lieshout (2016a) Cronie, O. and van Lieshout, M. N. M. (2016a).  Bandwidth selection for kernel estimators of the spatial intensity function, arXiv preprint arXiv:1611.10221 .
  * Cronie and van Lieshout (2016b) Cronie, O. and van Lieshout, M. N. M. (2016b).  Summary statistics for inhomogeneous marked point processes, Annals of the Institute of Statistical Mathematics 68: 905–928.
  * Cronie and van Lieshout (2018) Cronie, O. and van Lieshout, M. N. M. (2018).  A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions, Biometrika 105(2): 455–462.
  * Daley and Vere-Jones (2003) Daley, D. J. and Vere-Jones, D. (2003).  An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods, second edn, Springer-Verlag New York.
  * Daley and Vere-Jones (2008) Daley, D. J. and Vere-Jones, D. (2008).  An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure, second edn, Springer-Verlag New York.
  * Decreusefond and Vasseur (2018) Decreusefond, L. and Vasseur, A. (2018).  Stein’s method and papangelou intensity for Poisson or Cox process approximation, arXiv preprint arXiv:1807.02453 .
  * Diggle (2014) Diggle, P. (2014).  Statistical Analysis of Spatial and Spatio-Temporal Point Patterns, 3rd edn, Taylor & Francis/CRC Press, Boca Raton.
  * Flint et al. (2019) Flint, I., Privault, N. and Torrisi, G. L. (2019).  Functional inequalities for marked point processes, Electronic Journal of Probability 24(116): 1–40.
  * Ghorbani et al. (2020) Ghorbani, M., Cronie, O., Mateu, J. and Yu, J. (2020).  Functional marked point processes: a natural structure to unify spatio-temporal frameworks and to analyse dependent functional data, TEST pp. 1–40.
  * González et al. (2016) González, J. A., Rodriguez-Cortes, F. J., Cronie, O. and Mateu, J. (2016).  Spatio-temporal point process statistics: a review, Spatial Statistics 18: 505–544.
  * Guan et al. (2015) Guan, Y., Jalilian, A. and Waagepetersen, R. (2015).  Quasi-likelihood for spatial point processes, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(3): 677–697.
  * Guan and Shen (2010) Guan, Y. and Shen, Y. (2010).  A weighted estimating equation approach for inhomogeneous spatial point processes, Biometrika 97(4): 867–880.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R. and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer Science & Business Media.
  * Hough et al. (2009) Hough, J. B., Krishnapur, M., Peres, Y. and Virág, B. (2009).  Zeros of Gaussian analytic functions and determinantal point processes, Vol. 51 of University Lecture Series, American Mathematical Society, Providence, RI.
  * Hubbell and Foster (1983) Hubbell, S. P. and Foster, R. B. (1983).  Diversity of Canopy Trees in a Neotropical Forest and Implications for Conservation, in S. Sutton, T. Whitmore and A. Chadwick (eds), Tropical rain forest: ecology and management, Blackwell Scientific Publications, pp. 25–41.
  * Iftimi et al. (2019) Iftimi, A., Cronie, O. and Montes, F. (2019).  Second-order analysis of marked inhomogeneous spatiotemporal point processes: Applications to earthquake data, Scandinavian Journal of Statistics 46(3): 661–685.
  * Illian et al. (2008) Illian, J., Penttinen, A., Stoyan, H. and Stoyan, D. (2008).  Statistical Analysis and Modelling of Spatial Point Patterns, John Wiley & Sons.
  * James et al. (2013) James, G., Witten, D., Hastie, T. and Tibshirani, R. (2013).  An Introduction to Statistical Learning, Vol. 112, Springer.
  * Kallenberg (2017) Kallenberg, O. (2017).  Random Measures, Theory and Applications, Springer.
  * Kerscher (2000) Kerscher, M. (2000).  Statistical analysis of large-scale structure in the universe, in K. Mecke and D. Stoyan (eds), Statistical physics and spatial statistics, Vol. 554, Springer, pp. 36–71.
  * Kulesza and Taskar (2012) Kulesza, A. and Taskar, B. (2012).  Determinantal point processes for machine learning, Foundations and Trends® in Machine Learning 5(2–3): 123–286.
  * Last and Penrose (2017) Last, G. and Penrose, M. (2017).  Lectures on the Poisson process, Vol. 7, Cambridge University Press.
  * Lavancier et al. (2015) Lavancier, F., Møller, J. and Rubak, E. (2015).  Determinantal point process models and statistical inference, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(4): 853–877.
  * Lawrence et al. (2016) Lawrence, T., Baddeley, A., Milne, R. K. and Nair, G. (2016).  Point pattern analysis on a region of a sphere, Stat 5(1): 144–157.
  * Loader (1999) Loader, C. (1999).  Local Regression and Likelihood, Springer, New York.
  * Macchi (1975) Macchi, O. (1975).  The coincidence approach to stochastic point processes, Advances in Applied Probability 7(1): 83–122.
  * Marsan and Lengline (2008) Marsan, D. and Lengline, O. (2008).  Extending earthquakes’ reach through cascading, Science 319(5866): 1076–1079.
  * Mateu et al. (2020) Mateu, J., Moradi, M. and Cronie, O. (2020).  Spatio-temporal point patterns on linear networks: Pseudo-separable intensity estimation, Spatial Statistics 37: 100400.
  * McSwiggan et al. (2017) McSwiggan, G., Baddeley, A. and Nair, G. (2017).  Kernel density estimation on a linear network, Scandinavian Journal of Statistics 44(2): 324–345.
  * Meyer et al. (2012) Meyer, S., Elias, J. and Höhle, M. (2012).  A space–time conditional intensity model for invasive meningococcal disease occurrence, Biometrics 68(2): 607–616.
  * Møller and Rubak (2016) Møller, J. and Rubak, E. (2016).  Functional summary statistics for point processes on the sphere with an application to determinantal point processes, Spatial Statistics 18: 4–23.
  * Møller et al. (1998) Møller, J., Syversveen, A. R. and Waagepetersen, R. (1998).  Log Gaussian Cox processes, Scandinavian Journal of Statistics 25(3): 451–482.
  * Møller and Waagepetersen (2004) Møller, J. and Waagepetersen, R. (2004).  Statistical Inference and Simulation for Spatial Point Processes, CRC Press.
  * Møller and Waagepetersen (2007) Møller, J. and Waagepetersen, R. (2007).  Modern statistics for spatial point processes, Scandinavian Journal of Statistics 34(4): 643–684.
  * Møller and Waagepetersen (2017) Møller, J. and Waagepetersen, R. (2017).  Some recent developments in statistics for spatial point patterns, Annual Review of Statistics and Its Application 4: 317–342.
  * Moradi (2018) Moradi, M. (2018).  Spatial and Spatio-Temporal Point Patterns on Linear Networks, PhD dissertation, University Jaume I.
  * Moradi et al. (2019) Moradi, M., Cronie, O., Rubak, E., Lachieze-Rey, R., Mateu, J. and Baddeley, A. (2019).  Resample-smoothing of Voronoi intensity estimators, Statistics and Computing 29(5): 995––1010.
  * Moradi and Mateu (2020) Moradi, M. and Mateu, J. (2020).  First-and second-order characteristics of spatio-temporal point processes on linear networks, Journal of Computational and Graphical Statistics 29(3): 432–443.
  * Moradi et al. (2020) Moradi, M., Mateu, J. and Comas, C. (2020).  Directional analysis for point patterns on linear networks, Stat p. e323.
  * Moradi et al. (2018) Moradi, M., Rodriguez-Cortes, F. and Mateu, J. (2018).  On kernel-based intensity estimation of spatial point patterns on linear networks, Journal of Computational and Graphical Statistics 27(2): 302–311.
  * Ogata (1998) Ogata, Y. (1998).  Space-time point-process models for earthquake occurrences, Annals of the Institute of Statistical Mathematics 50(2): 379–402.
  * Ogata and Katsura (1991) Ogata, Y. and Katsura, K. (1991).  Maximum likelihood estimates of the fractal dimension for random spatial patterns, Biometrika 78(3): 463–474.
  * Okabe and Sugihara (2012) Okabe, A. and Sugihara, K. (2012).  Spatial Analysis along Networks: Statistical and Computational Methods, John Wiley & Sons.
  * Pelletier (2005) Pelletier, B. (2005).  Kernel density estimation on riemannian manifolds, Statistics & probability letters 73(3): 297–304.
  * Poinas et al. (2019) Poinas, A., Delyon, B., Lavancier, F. et al. (2019).  Mixing properties and central limit theorem for associated point processes, Bernoulli 25(3): 1724–1754.
  * Prokešová et al. (2017) Prokešová, M., Dvořák, J. and Jensen, E. B. V. (2017).  Two-step estimation procedures for inhomogeneous shot-noise cox processes, Annals of the Institute of Statistical Mathematics 69(3): 513–542.
  * Rajala et al. (2018) Rajala, T., Murrell, D. J. and Olhede, S. C. (2018).  Detecting multivariate interactions in spatial point patterns with gibbs models and variable selection, Journal of the Royal Statistical Society: Series C (Applied Statistics) 67(5): 1237–1273.
  * Rakshit et al. (2019) Rakshit, S., Davies, T. M., Moradi, M., McSwiggan, G., Nair, G., Mateu, J. and Baddeley, A. (2019).  Fast kernel smoothing of point patterns on a large network using two-dimensional convolution, International Statistical Review 87(3): 531–556.
  * Rakshit et al. (2017) Rakshit, S., Nair, G. and Baddeley, A. (2017).  Second-order analysis of point patterns on a network using any distance metric, Spatial Statistics 22: 129–154.
  * Ripley (1991) Ripley, B. D. (1991).  Statistical Inference for Spatial Processes, Cambridge university press.
  * Robeson et al. (2014) Robeson, S. M., Li, A. and Huang, C. (2014).  Point-pattern analysis on the sphere, Spatial Statistics 10: 76–86.
  * Schneider and Weil (2008) Schneider, R. and Weil, W. (2008).  Stochastic and Integral Geometry, Probability and Its Applications, Springer, Dordrecht.
  * Shirai and Takahashi (2003) Shirai, T. and Takahashi, Y. (2003).  Random point fields associated with certain fredholm determinants i: fermion, poisson and boson point processes, Journal of Functional Analysis 205(2): 414–463.
  * Silverman (1986) Silverman, B. W. (1986).  Density estimation for statistics and data analysis, Vol. 26, Chapman & Hall/CRC press.
  * Stoyan and Penttinen (2000) Stoyan, D. and Penttinen, A. (2000).  Recent applications of point process methods in forestry statistics, Statistical science 15(1): 61–78.
  * Toreti et al. (2019) Toreti, A., Cronie, O. and Zampieri, M. (2019).  Concurrent climate extremes in the key wheat producing regions of the world, Scientific reports 9(1): 1–8.
  * van Lieshout (2000) van Lieshout, M. N. M. (2000).  Markov Point Processes and Their Applications, Imperial College Press/World Scientific.
  * van Lieshout (2006) van Lieshout, M. N. M. (2006).  Campbell and moment measures for finite sequential spatial processes, in M. Hušková and M. Janžur (eds), Proceedings Prague Stochastics, Matfyzpress, Prague, p. 215–224.
  * van Lieshout (2011) van Lieshout, M. N. M. (2011).  A J-function for inhomogeneous point processes, Statistica Neerlandica 65(2): 183–201.
  * van Lieshout (2012) van Lieshout, M. N. M. (2012).  On estimation of the intensity function of a point process, Methodology and Computing in Applied Probability 14: 567–578.
  * van Lieshout (2020) van Lieshout, M. N. M. (2020).  Infill asymptotics and bandwidth selection for kernel estimators of spatial intensity functions, Methodology and computing in applied probability 22(3): 995–1008.
  * Vapnik (2013) Vapnik, V. (2013).  The Nature of Statistical Learning Theory, Springer science & business media.
  * Waagepetersen (2007) Waagepetersen, R. (2007).  An estimating function approach to inference for inhomogeneous Neyman–Scott processes, Biometrics 63(1): 252–258.
  * Wood et al. (2008) Wood, S. N., Bravington, M. V. and Hedley, S. L. (2008).  Soap film smoothing, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(5): 931–955.
  * Yang et al. (2019) Yang, J., Rao, V. and Neville, J. (2019).  A stein–papangelou goodness-of-fit test for point processes, The 22nd International Conference on Artificial Intelligence and Statistics, pp. 226–235.
  * Yue and Loh (2015) Yue, Y. and Loh, J. M. (2015).  Variable selection for inhomogeneous spatial point process models, Canadian Journal of Statistics 43(2): 288–305.
  * Zou and Hastie (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67(2): 301–320.

##  Appendix A Univariate innovations

In the subsequent subsections we will consider different settings where the
statistical inference fits into the univariate innovations framework. We
stress that there may be existing approaches which do not fit into this
framework and, also, that there may be approaches which do fit in but which
(accidentally) have not been included; we have focused on approaches which
deal with inhomogeneous point processes.

###  A.1 Non-parametric intensity estimation

As the intensity function $\rho(u)$, $u\in W$, governs the univariate
distributions of the points of $X\cap W$, it is not hard to imagine that its
estimation forms the very foundation of subsequent statistical analyses. When
we do not have access to any covariates (on $W$), we have to resort to non-
parametric estimation. Although a considerable amount of progress has been
made regarding non-parametric estimation of intensity functions (van Lieshout,
2012, Baddeley et al., 2015, Cronie and van Lieshout, 2018, Moradi, 2018,
Moradi et al., 2019, van Lieshout, 2020), there is still significant room for
improvement.

We may represent non-parametric intensity estimation using the general
estimator structure in (3.1):

$\widehat{\rho}_{\theta}(u,{\mathbf{y}})=\xi_{\theta}^{1}(u;{\mathbf{y}}),\quad
u\in W,\quad{\mathbf{y}}\in{\mathcal{X}},$

where $\theta\in\Theta\subseteq\mathbb{R}^{l}$, $l\geq 1$, is the
tuning/smoothing parameter associated to the estimator. Given an observed
point pattern ${\mathbf{x}}$, to obtain an estimate/selection
$\widehat{\theta}_{W}({\mathbf{x}})$ we consider
$\widehat{\rho}_{\theta}(u,{\mathbf{x}})=\xi_{\theta}^{1}(u;{\mathbf{x}})$ and
minimise
$\mathcal{L}(\theta;{\mathbf{x}})=\mathcal{I}_{\widehat{\rho}_{\theta}}^{h_{\theta}}(W;{\mathbf{x}})^{2}$,
$\theta\in\Theta$, where

$\mathcal{I}_{\widehat{\rho}_{\theta}}^{h_{\theta}}(W;{\mathbf{x}})=\sum_{x\in{\mathbf{x}}\cap
W}h_{\theta}(x;{\mathbf{x}}\setminus\\{x\\})-\int_{A}h_{\theta}(u;{\mathbf{x}})\widehat{\rho}_{\theta}(u;{\mathbf{x}})\mathrm{d}u,$

using some family of test functions $h_{\theta}$, $\theta\in\Theta$. We then
plug $\widehat{\theta}_{W}({\mathbf{x}})$ into the intensity estimator, i.e.
$\widehat{\rho}_{\widehat{\theta}_{W}({\mathbf{x}})}(u,{\mathbf{x}})$, $u\in
W$, to obtain the final estimate of the intensity function.

######  Remark A.1.

If the estimator has no tuning parameter $\theta$, which e.g. is the case for
the Voronoi intensity estimator (Barr and Schoenberg, 2010, Moradi et al.,
2019)

$\displaystyle\widehat{\rho}_{\theta}(u,{\mathbf{y}})$
$\displaystyle=\sum_{x\in{\mathbf{y}}}\frac{{\mathbf{1}}\\{u\in\mathcal{V}_{x}\\}}{|\mathcal{V}_{x}|},\quad
u\in W,$ $\displaystyle\mathcal{V}_{x}$ $\displaystyle=\\{v\in W:d(v,x)\leq
d(v,y)\text{ for all }y\in{\mathbf{y}}\setminus\\{x\\}\\},\quad
x\in{\mathbf{y}},$

we can simply set $\widehat{\rho}_{\theta_{1}}=\widehat{\rho}_{\theta_{2}}$
for any $\theta_{1},\theta_{2}\in\Theta$ or $\Theta=\emptyset$. In the case of
the so-called Resample-smoothed Voronoi intensity estimator of Moradi et al.
(2019) on the other hand, there are two tuning parameters to be selected.

####  A.1.1 Kernel estimation

The arguably most prominent intensity estimation technique is kernel
estimation (Illian et al., 2008, van Lieshout, 2012, Diggle, 2014, Baddeley et
al., 2015). E.g., when $W\subseteq S=\mathbb{R}^{d}$, $d\geq 1$, in its most
typical form a kernel estimate is given by

$\displaystyle\widehat{\rho}_{\theta}(u,{\mathbf{x}})=\sum_{x\in{\mathbf{x}}}\frac{\kappa_{\theta}(u-x)}{w_{\theta}(u,x)}=\sum_{x\in{\mathbf{x}}}\frac{\theta^{-d}\kappa((u-x)/\theta)}{w_{\theta}(u,x)},\quad
u\in W,$ (A.1)

where the kernel $\kappa$ is a symmetric density function on $\mathbb{R}^{d}$,
$w_{\theta}(u,x)$ is an edge correction term which compensates for possible
interaction between points inside and outside $W$, and the smoothing parameter
$\theta\in\Theta=(0,\infty)$ is the so-called bandwidth. A common edge
correction factor is the local factor
$w_{\theta}(u,x)=\int_{W}\kappa_{\theta}(v-x)\mathrm{d}v$, which also ensures
that
$\int_{W}\widehat{\rho}_{\theta}(u,{\mathbf{x}})\mathrm{d}u=\\#{\mathbf{x}}$.
It should be emphasised that for point processes in other (non-Euclidean)
domains $S$, kernel functions $\kappa$, and thereby kernel estimators, may
generally take other forms and can be quite abstract entities (Pelletier,
2005, McSwiggan et al., 2017, Moradi, 2018, Mateu et al., 2020). In certain
cases, however, there are straightforward extensions of (A.1); for a point
pattern on a linear network $S=L$, we may e.g. consider the kernel estimate
obtained by leaving the kernel unchanged in (A.1) and simply replacing the
edge correction $w_{\theta}(u,x)$ by
$w_{\theta}(u,x)=\int_{L}\kappa_{\theta}(v-x)\mathrm{d}v$, where the
integration is with respect to arc length, i.e. 1-dimensional Hausdorff
measure in $\mathbb{R}^{2}$ (Rakshit et al., 2019).

In kernel estimation, the main challenge is optimal bandwidth selection; often
the choice of kernel is of much less importance than the choice of bandwidth
(Silverman, 1986). To the best of our knowledge, the two best-performing
approaches currently available for point processes in $\mathbb{R}^{d}$ are the
Poisson process likelihood leave-one-out cross-validation approach (Loader,
1999, Baddeley et al., 2015) and the Campbell formula-based approach of Cronie
and van Lieshout (2016a, 2018). Given a point pattern ${\mathbf{x}}\subseteq
W$, the Poisson process likelihood leave-one-out cross-validation approach
maximises

$\displaystyle\theta\mapsto\sum_{x\in{\mathbf{x}}\cap
W}\log\widehat{\rho}_{\theta}(x,{\mathbf{x}}\setminus\\{x\\})-\int_{W}\widehat{\rho}_{\theta}(u,{\mathbf{x}})\mathrm{d}u$
(A.2)

in order to obtain an optimal bandwidth; the Poisson process log-likelihood
function is a direct consequence of (2.15). To express this as a loss
function, we may either multiply (A.2) by $-1$ or consider its derivative with
respect to $\theta$, assuming sufficient differentiability; the latter results
in a univariate innovation of the kind given in Definition (3.1). This
approach is particularly suited when the data come from a Poisson process
(Cronie and van Lieshout, 2018). Moreover, Cronie and van Lieshout (2018)
proposed to select the bandwidth by minimising the squared innovation

$\mathcal{L}(\theta)=\mathcal{I}_{\widehat{\rho}_{\theta}}^{f(\widehat{\rho}_{\theta}(\cdot))}(W;{\mathbf{x}})^{2}=\left(\sum_{x\in{\mathbf{x}}\cap
W}f(\widehat{\rho}_{\theta}(x,{\mathbf{x}}\setminus\\{x\\}))-\int_{W}f(\widehat{\rho}_{\theta}(u,{\mathbf{x}}))\widehat{\rho}_{\theta}(u,{\mathbf{x}})\mathrm{d}u\right)^{2},$

where $f:\mathbb{R}\to\mathbb{R}$. They studied in detail the choice
$f(x)=1/x$, $x\in\mathbb{R}$, which results in the so-called inverse/Stoyan-
Grabarnik statistic/diagnostic, and found that

$\displaystyle\mathcal{L}(\theta)=\left(\sum_{x\in{\mathbf{x}}\cap
W}f(\widehat{\rho}_{\theta}(x,{\mathbf{x}}))-\int_{W}f(\widehat{\rho}_{\theta}(u,{\mathbf{x}}))\widehat{\rho}_{\theta}(u,{\mathbf{x}})\mathrm{d}u\right)^{2}$
(A.3)

with the choice $f(x)=1/x$, $x\in\mathbb{R}$, gave rise to (A.3) being the
square of a (conjectured) monotonic function of $\theta\geq 0$ when using
$w_{\theta}(u,x)\equiv 1$ (no edge correction) and a Gaussian kernel; the
unmodified version is not monotonic. Note that $f(x)=1/x$, $x\in\mathbb{R}$,
sets the integral in (A.3) to $|W|$ when
$\widehat{\rho}_{\theta}(u,{\mathbf{x}})>0$, $u\in W$, so we choose the
bandwidth by estimating the (known) size of the study region by the sum of
reciprocal intensity estimates. Moradi et al. (2019) indicated that this
particular choice of test function favours a low variance with respect to the
variance-bias trade off, which makes it particularly suited for aggregated
point processes.

###  A.2 Parametric product density and Palm-likelihood estimation

The general idea in parametric product density estimation is to fit some
parametric family of functions
$\rho_{\theta}^{(n)}(u_{1},\ldots,u_{n})=\xi_{\theta}^{1}(u_{1},\ldots,u_{n})\geq
0$, $u_{1},\ldots,u_{n}\in W$, $\theta\in\Theta$, to the observed point
pattern ${\mathbf{x}}\subseteq W$. The corresponding innovations
$\mathcal{I}_{\rho_{\theta}^{(n)}}^{h_{\theta}}(\cdot;{\mathbf{y}})$,
$\theta\in\Theta$, ${\mathbf{y}}\in{\mathcal{X}}$, and loss function are given
by

$\displaystyle\mathcal{I}_{\rho_{\theta}^{(n)}}^{h_{\theta}}(A;{\mathbf{y}})=$
$\displaystyle\sum_{(x_{1},\ldots,x_{n})\in{\mathbf{y}}_{\neq}^{n}\cap
A}h_{\theta}(x_{1},\ldots,x_{n})$ (A.4)
$\displaystyle-\int_{A}h_{\theta}(u_{1},\ldots,u_{n})\rho_{\theta}^{(n)}(u_{1},\ldots,u_{n})\mathrm{d}u_{1}\cdots\mathrm{d}u_{n},\quad
A\subseteq S^{n},$ $\displaystyle\mathcal{L}(\theta)=$
$\displaystyle\mathcal{I}_{\rho_{\theta}^{(n)}}^{h_{\theta}}(W^{n};{\mathbf{x}})^{2},$

respectively, and by minimising $\mathcal{L}(\theta)$, $\theta\in\Theta$, we
consequently obtain an estimate
$\widehat{\theta}=\widehat{\theta}_{W}({\mathbf{x}})\in\Theta$ of
$\theta_{0}$, which yields the final intensity estimate
$\rho_{\widehat{\theta}_{W}({\mathbf{x}})}^{(n)}(u_{1},\ldots,u_{n})$,
$u_{1},\ldots,u_{n}\in W$.

####  A.2.1 Intensity estimation

In particular, when $n=1$, i.e. when we are carrying out parametric intensity
estimation by means of some parametric family of intensity functions
$\rho_{\theta}(u)=\xi_{\theta}^{1}(u)\geq 0$, $u\in W$, $\theta\in\Theta$, we
obtain the innovations

$\displaystyle\mathcal{I}_{\rho_{\theta}}^{h_{\theta}}(A;{\mathbf{y}})=$
$\displaystyle\sum_{x\in{\mathbf{y}}\cap
A}h_{\theta}(x)-\int_{A}h_{\theta}(u)\rho_{\theta}(u)\mathrm{d}u,\quad
A\subseteq S,\quad\theta\in\Theta,\quad{\mathbf{y}}\in{\mathcal{X}}.$ (A.5)

Viewed as an unbiased estimating equation approach (recall Lemma 3.1), the
corresponding loss function here essentially corresponds to the quasi-
likelihood approach of Guan et al. (2015). Moreover, when one has the access
to a vector of $l\geq 1$ covariates $z(u)=(z_{1}(u),\ldots,z_{l}(u))^{\top}$
measured at location $u\in W$, it is often natural to model the intensity
function parametrically through these covariates. The arguably most common and
popular model choice for the intensity function is the log-linear form

$\displaystyle\rho_{\theta}(u)=\exp\\{\theta^{\top}z(u)\\},\quad u\in W,$
(A.6)

where
$\theta=(\theta_{1},\ldots,\theta_{l})^{\top}\in\Theta\subseteq\mathbb{R}^{l}$
is a real-valued $l$-dimensional parameter (Coeurjolly and Lavancier, 2019).
Note that an intercept may be included by setting $z_{1}(\cdot)\equiv 1$ and,
in addition, (projections of) the location $u\in W$ itself may also be
included in $z(u)$ if explicit spatial dependence is warranted. When we assume
that $X$ is a Poisson process, the standard procedure is to maximise the log-
likelihood function (see (2.15)) with respect to $\theta\in\Theta$, which for
the intensity model (A.6) is given by

$\displaystyle\ell(\theta;{\mathbf{x}})=\sum_{x\in{\mathbf{x}}\cap
W}\log\rho_{\theta}(u)-\int_{W}\rho_{\theta}(u)\mathrm{d}u=\sum_{x\in{\mathbf{x}}\cap
W}\theta^{\top}z(u)-\int_{W}\exp\\{\theta^{\top}z(u)\\}\mathrm{d}u.$ (A.7)

This has gradient

$\nabla\ell(\theta;{\mathbf{x}})=\sum_{u\in{\mathbf{x}}\cap
W}z(u)-\int_{W}z(u)\rho_{\theta}(u)\mathrm{d}u\in\mathbb{R}^{l},$ (A.8)

and solving $\nabla\ell(\theta;{\mathbf{x}})=0\in\mathbb{R}^{l}$ is equivalent
to minimising an $l$-dimensional vector of squared innovations of the form
(A.5), where entry $i=1\ldots,l$ is based on the test function
$h_{i}(u)=z_{i}(u)$, $u\in W$. Even when $X$ is not a Poisson process, we
obtain that (A.8) remains an unbiased estimating equation (recall Lemma 3.1).
Hence, the maximum of (A.7) still makes sense for non-Poissonian models, an
the obtained estimator may be interpreted as a composite likelihood estimator
(Waagepetersen, 2007, Guan and Shen, 2010, Guan et al., 2015).

####  A.2.2 Palm-likelihood estimation

In certain cases, the go-to approach is Palm-likelihood estimation (Baddeley
et al., 2015, Coeurjolly and Lavancier, 2019). Following Prokešová et al.
(2017), for a SOIRS point process $X$, one instance of log-Palm-likelihood
estimation is obtained by maximising

$\displaystyle\ell_{P}(\theta;{\mathbf{x}})=$
$\displaystyle\mathop{\sum\nolimits^{\neq}}_{x_{1},x_{2}\in{\mathbf{x}}\cap
W}{\mathbf{1}}\\{x_{1}\in W_{\ominus r}\\}{\mathbf{1}}\\{d(x_{1},x_{2})\leq
r\\}\log\rho_{\theta}(x_{1})\rho_{\theta}(x_{2})g_{\theta}^{(2)}(x_{1}-x_{2})$
$\displaystyle-\int_{W}\int_{W}{\mathbf{1}}\\{u\in W_{\ominus
r}\\}{\mathbf{1}}\\{d(u_{1},u_{2})\leq
r\\}\rho_{\theta}(u_{1})\rho_{\theta}(u_{2})g_{\theta}^{(2)}(u_{1}-u_{2})\mathrm{d}u$

for some $r>0$, where $W_{\ominus r}$ is an $r$-trimming of $W$, and
$\rho_{\theta}(\cdot)$ and $g_{\theta}^{(2)}(\cdot)$ are the intensity
function and the pair correlation function of the model parametrised by
$\theta\in\Theta$. Prokešová et al. (2017) considered a two-step procedure
where the intensity was estimated externally and then plugged into the above
log-Palm-likelihood function, thus replacing $\rho_{\theta}(\cdot)$. The name
Palm-likelihood stems from the original approach (Ogata and Katsura, 1991) for
stationary point processes in $\mathbb{R}^{d}$, where the idea is to apply
Poisson process likelihood heuristics to the intensity functions of the first
order reduced Palm processes $X_{u}^{!}$; recall expression (2.8). Assuming
sufficient differentiability, if we let

$\displaystyle\xi_{\theta}^{2}(u_{1},u_{2})=$
$\displaystyle\rho_{\theta}(u_{1})\rho_{\theta}(u_{2})g_{\theta}^{(2)}(u_{1}-u_{2}),$
$\displaystyle h_{\theta}(u_{1},u_{2})=$ $\displaystyle{\mathbf{1}}\\{u\in
W_{\ominus r}\\}{\mathbf{1}}\\{d(u_{1},u_{2})\leq
r\\}\frac{\nabla(\rho_{\theta}(u_{1})\rho_{\theta}(u_{2})g_{\theta}^{(2)}(u_{1}-u_{2}))}{\rho_{\theta}(u_{1})\rho_{\theta}(u_{2})g_{\theta}^{(2)}(u_{1}-u_{2})},$

then the corresponding score function is given by the $l$-dimensional vector
of innovations

$\displaystyle\nabla\ell_{P}(\theta;{\mathbf{x}})=\mathcal{I}_{\xi_{\theta}^{2}}^{h_{\theta}}(W^{2};{\mathbf{x}})=$
$\displaystyle=$
$\displaystyle\mathop{\sum\nolimits^{\neq}}_{x_{1},x_{2}\in{\mathbf{x}}\cap
W}{\mathbf{1}}\\{x_{1}\in W_{\ominus r}\\}{\mathbf{1}}\\{d(x_{1},x_{2})\leq
r\\}\frac{\nabla(\rho_{\theta}(x_{1})\rho_{\theta}(x_{2})g_{\theta}^{(2)}(x_{1}-x_{2}))}{\rho_{\theta}(x_{1})\rho_{\theta}(x_{2})g_{\theta}^{(2)}(x_{1}-x_{2})}$
$\displaystyle-\int_{W}\int_{W}{\mathbf{1}}\\{u\in W_{\ominus
r}\\}{\mathbf{1}}\\{d(u_{1},u_{2})\leq
r\\}\nabla(\rho_{\theta}(u_{1})\rho_{\theta}(u_{2})g_{\theta}^{(2)}(u_{1}-u_{2}))\mathrm{d}u\in\mathbb{R}^{l},$

provided that we can interchange differentiation and integration. Hence,
carrying out Palm-likelihood is equivalent to minimising a vector of squared
innovations of the form (A.4).

######  Remark A.2.

One could, potentially, consider a similar approach for other (higher-order)
reduced Palm product densities/correlation functions, where other (higher-
order) innovations would be exploited.

###  A.3 Papangelou conditional intensity estimation

The idea here is to assume that the observed point pattern
${\mathbf{x}}\subseteq W$ is a realisation of a point process $X$ with
Papangelou conditional intensity $\lambda_{\theta_{0}}$,
$\theta_{0}\in\Theta$, which is a member of some parametric family of
Papangelou conditional intensity functions
$\Xi_{\Theta}^{1}=\\{\xi_{\theta}^{1}:\theta\in\Theta\\}=\\{\lambda_{\theta}:\theta\in\Theta\\}$,
$\Theta\subseteq\mathbb{R}^{l}$, $l\geq 1$. Hence, we obtain classical
innovations and the innovation-based loss function

$\displaystyle\mathcal{L}(\theta)=\mathcal{I}_{\lambda_{\theta}}^{h_{\theta}}(W;{\mathbf{x}})^{2}=$
$\displaystyle\left(\sum_{x\in{\mathbf{x}}\cap
W}h_{\theta}(x;{\mathbf{x}}\setminus\\{x\\})-\int_{W}h_{\theta}(u;{\mathbf{x}})\lambda_{\theta}(u;{\mathbf{x}})\mathrm{d}u\right)^{2},\quad\theta\in\Theta,$

for some suitable family of test functions $\mathcal{H}_{\Theta}$. In the
literature, this is sometimes referred to as Takacs–Fiksel estimation (see
e.g. Møller and Waagepetersen, 2017, Coeurjolly and Lavancier, 2019, and the
references therein). The most prominent special case hereof, namely pseudo-
likelihood estimation (Besag, 1974, Ripley, 1991, van Lieshout, 2000, Møller
and Waagepetersen, 2004, Baddeley et al., 2015), is obtained by setting
$h_{\theta}(\cdot)=\nabla\lambda_{\theta}(\cdot)/\lambda_{\theta}(\cdot)\in\mathbb{R}^{l}$,
where $\nabla\lambda_{\theta}(\cdot)$ is the (well-defined) gradient of
$\lambda_{\theta}$. Note that this yields the score function of the
(log-)pseudo-likelihood function, which in turn may be expressed as the loss
function

$\displaystyle\mathcal{L}(\theta)=$
$\displaystyle\int_{W}\lambda_{\theta}(u;{\mathbf{x}})\mathrm{d}u-\sum_{x\in{\mathbf{x}}\cap
W}\log\lambda_{\theta}(x;{\mathbf{x}}\setminus\\{x\\}),\quad\theta\in\Theta.$

###  A.4 Minimum contrast estimation of second-order summary statistics

We next turn to the case where the general parametrised estimator family is of
the form (3.2) and is given by $\xi_{\theta}^{2}=\rho_{\theta}^{(2)}$,
$\theta\in\Theta$, for some family of second-order product densities. Assuming
that the point process $X$ in $S=\mathbb{R}^{d}$ has the second-order product
density $\rho_{\theta_{0}}^{(2)}$ for some $\theta_{0}\in\Theta$, we consider
the test function

$h(u_{1},u_{2};r)={\mathbf{1}}\\{d(u_{1},u_{2})\leq
r\\}/(|W|\rho(u_{1})\rho(u_{2})),\quad u_{1},u_{2}\in W,$

where $r\geq 0$, $\rho(u)$, $u\in W$, is the true (unknown) intensity of $X$
and $d(u_{1},u_{2})=\|u_{1}-u_{2}\|_{2}$ is the Euclidean metric on
$\mathbb{R}^{d}$. This yields the innovations

$\displaystyle\mathcal{I}_{\rho_{\theta}^{(2)}}^{h(\cdot;r)}(W\times\mathbb{R}^{d};X)=$
$\displaystyle\frac{1}{|W|}\sum_{(x_{1},x_{2})\in
X_{\neq}^{2}\cap(W\times\mathbb{R}^{d})}\frac{{\mathbf{1}}\\{d(x_{1},x_{2})\leq
r\\}}{\rho(x_{1})\rho(x_{2})}$ (A.9)
$\displaystyle-\frac{1}{|W|}\int_{W}\int_{b(u_{1},r)}g_{\theta}^{(2)}(u_{1},u_{2})\mathrm{d}u_{1}\mathrm{d}u_{2}\quad\theta\in\Theta,$

where $b(u,r)$ denotes the closed Euclidean $r$-ball around
$u\in\mathbb{R}^{d}$ and
$g_{\theta}^{(2)}(u_{1},u_{2})=\rho_{\theta}^{(2)}(u_{1},u_{2})/(\rho(u_{1})\rho(u_{2}))$
is the pair correlation function corresponding to $\theta$. If the true
second-order product density of $X$ is given by $\rho_{\theta_{0}}^{(2)}$ for
some $\theta_{0}\in\Theta$, by Lemma 3.1 we obtain that
$\mathbb{E}[\mathcal{I}_{\rho_{\theta_{0}}^{(2)}}^{h}(W;X)]=0$. Under the
assumption of second-order intensity reweighted stationarity (recall Section
2.3.1), i.e. when the intensity is positive and the pair correlation functions
satisfy $g_{\theta}^{(2)}(u_{1},u_{2})=g_{\theta}^{(2)}(u_{1}-u_{2})$,
$u_{1},u_{2}\in\mathbb{R}^{d}$, $\theta\in\Theta$, we have that the integral
term in (A.9) becomes

$K_{\rm inhom}(r;\theta)=\int_{b(v,r)}g_{\theta}^{(2)}(u)\mathrm{d}u$

for any $W\subseteq\mathbb{R}^{d}$, $|W|>0$, and any $\theta\in\Theta$; these
functions are constant as functions of $v\in\mathbb{R}^{d}$. Note that the
function $K_{\rm inhom}(r)=K_{\rm inhom}(r;\theta_{0})$, $r\geq 0$, is the
inhomogeneous $K$-function of $X$ (Baddeley et al., 2000), which a measure of
interaction. For a Poisson process $X$ we have that $K_{\rm
inhom}(r)=|b(v,r)|$, $v\in\mathbb{R}^{d}$, so $K_{\rm inhom}(r)>|b(v,r)|$
indicates clustering/aggregation between points with inter-point distance at
most $r$ whereas $K_{\rm inhom}(r)<|b(v,r)|$ indicates inhibition/regularity.

In practice one considers
$\mathcal{I}_{\rho_{\theta}^{(2)}}^{h(\cdot;r)}(W\times W;X)$ because we only
observe events within $W$. Consequently, one needs to include an edge-
correction/normalisation factor $e(\cdot,\cdot)$ in the test function. This is
done by simply changing the test function to

$\displaystyle\widetilde{h}_{e}(u_{1},u_{2};r)=e(u_{1},u_{2})h(u_{1},u_{2};r).$

An edge-correction factor adjusts for (unobserved) interactions between events
in $X\cap W$ and those in $X\cap\mathbb{R}^{d}\setminus W$. Ideally, including
$e(\cdot,\cdot)$ will ensure that
$\mathcal{I}_{\rho_{\theta}^{(2)}}^{\widetilde{h}(\cdot;r)}(W^{2};X)$ has
expectation 0. For further details, see e.g. (Baddeley et al., 2000, Møller
and Waagepetersen, 2004, Baddeley et al., 2015). When employing
$\mathcal{I}_{\rho_{\theta}^{(2)}}^{\widetilde{h}_{e}(\cdot;r)}(W^{2};X)$,
i.e. using $\widetilde{h}_{e}$ as test function, we refer to the sum-term in
(A.9) as a non-parametric estimator of $K_{\rm inhom}(r)$.

Assuming second-order intensity reweighted stationarity, by minimising the
loss function

$\mathcal{L}(\theta)=\int_{0}^{r_{max}}\mathcal{I}_{\rho_{\theta}^{(2)}}^{\widetilde{h}_{e}(\cdot;r)}(W^{2};X)^{2}\mathrm{d}r,\quad\theta\in\Theta,$

for some $r_{max}>0$, we obtain a minimum contrast estimator
$\widehat{\theta}_{W}(X)$ of $\theta_{0}$. There are also further developments
of this concept, where metrics other than the $L_{2}$-distance are used to
measure discrepancies between $K_{\rm inhom}(r;\theta)$, $r\geq 0$, and its
estimator (Baddeley et al., 2015, Diggle, 2014).

##  Appendix B Proofs

###  B.1 Proof of Theorem 1

###### Proof of Theorem 1.

Starting with expression (2.22), the form of the Papangelou conditional
intensity is a direct consequence of Decreusefond and Vasseur (2018, Theorem
4.7) and (2.9). The result on the product densities follows from e.g. Baccelli
et al. (2020, Proposition 2.3.24).

The structure of the proof of the prediction formula in (1) follows the lines
of the proof of Last and Penrose (2017, Exercise 5.9). Consider the random
measure representation of $X$, where there are random variables
$N=X(S)\in\\{0,\ldots,\infty\\}$ and $X_{1},\ldots X_{N}\in S$ such that
$X(A)=\sum_{i=1}^{N}\delta_{X_{i}}(A)=\sum_{i=1}^{N}{\mathbf{1}}\\{X_{i}\in
A\\}$, $A\subseteq S$. An independent thinning of $X$ has the same
distribution as $Z(\cdot)=\sum_{i=1}^{N}B_{i}\delta_{X_{i}}(\cdot)$, where,
conditionally on $N$ and $X_{1},\ldots,X_{N}$, the random variables
$B_{1},\ldots,B_{N}$ are mutually independent and, for any $i=1,\ldots,N$,
conditionally on $X_{i}$, the random variable $B_{i}$ is Bernoulli distributed
with parameter $p(X_{i})$. Similarly, $Y=X\setminus Z$ has the random measure
representation
$Y(\cdot)=X(\cdot)-Z(\cdot)=\sum_{i=1}^{N}(1-B_{i})\delta_{X_{i}}(\cdot)$. For
any $m\geq n$, let $\mathcal{A}_{m}$ be the set of all $n$-tuples of distinct
integers $i_{1},\ldots,i_{n}\in\\{1,\ldots,m\\}$; if $m$ is infinite, we let
$i_{1},\ldots,i_{n}$ be finite. It now follows that

$\displaystyle\mathbb{E}\left[\mathop{\sum\nolimits^{\neq}}_{x_{1},\ldots,x_{n}\in
Z}h(x_{1},\ldots,x_{n},Y)\prod_{i=1}^{n}(1-p(x_{i}))\right]=$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y)\prod_{j=1}^{n}B_{i_{j}}(1-p(X_{i_{j}}))\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y)\prod_{j=1}^{n}\mathbb{E}[1-B_{i_{j}}|X]\prod_{j=1}^{n}B_{i_{j}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\})\prod_{j=1}^{n}\mathbb{E}[1-B_{i_{j}}|X]\prod_{j=1}^{n}B_{i_{j}}\right],$

where we have used that $p(X_{i_{j}})=\mathbb{E}[B_{i_{j}}|X]$ and that
$Y\cap\\{X_{i_{1}},\ldots,X_{i_{n}}\\}=\emptyset$, i.e.
$Y=Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\}$, when $B_{i_{j}}=1$ for all
$j=1,\ldots,n$. By the conditional independence of the $B_{i_{j}}$’s, we have
that
$\prod_{j=1}^{n}\mathbb{E}[1-B_{i_{j}}|X]=\mathbb{E}[\prod_{j=1}^{n}(1-B_{i_{j}})|X]$,
and writing
$\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)=h(X_{i_{1}},\ldots,X_{i_{n}},Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\})$,
we obtain

$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\})\prod_{j=1}^{n}\mathbb{E}[1-B_{i_{j}}|X]\prod_{j=1}^{n}B_{i_{j}}\right]=$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\left.\mathbb{E}\left[\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}B_{i_{j}}\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\left.\mathbb{E}\left[\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}B_{i_{j}}\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}B_{i_{j}}\right],$

where the last equality follows from the "pulling out known factors" property
of conditional expectations; $N$ and $\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)$ are
measurable with respect to the $\sigma$-algebra generated by $X$. By the law
of total expectation it follows that

$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}B_{i_{j}}\right]=$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}B_{i_{j}}\right|X\right]\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\left.\mathbb{E}\left[\prod_{j=1}^{n}B_{i_{j}}\right|X\right]\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\right|X\right]\prod_{j=1}^{n}p(X_{i_{j}})\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[\left.\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\prod_{j=1}^{n}p(X_{i_{j}})\right|X\right]\right]$
$\displaystyle=$
$\displaystyle\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{\infty}}\mathbb{E}\left[{\mathbf{1}}\\{N\geq\max\\{i_{1},\ldots,i_{n}\\}\\}\tilde{h}_{i_{1},\ldots,i_{n}}(X,Y)\prod_{j=1}^{n}(1-B_{i_{j}})\prod_{j=1}^{n}p(X_{i_{j}})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\})\prod_{j=1}^{n}(1-B_{i_{j}})\prod_{j=1}^{n}p(X_{i_{j}})\right],$

where we have used the fact that
$\mathbb{E}[\prod_{j=1}^{n}B_{i_{j}}|X]=\prod_{j=1}^{n}\mathbb{E}[B_{i_{j}}|X]=\prod_{j=1}^{n}p(X_{i_{j}})$
by the conditional independence of the $B_{i_{j}}$’s, as well as the above-
mentioned property of conditional expectations for
$\prod_{j=1}^{n}p(X_{i_{j}})$ and the $\sigma$-algebra generated by $X$.
Exploiting the representation
$Y(\cdot)=X(\cdot)-Z(\cdot)=\sum_{i=1}^{N}(1-B_{i})\delta_{X_{i}}(\cdot)$, we
finally obtain that

$\displaystyle\mathbb{E}\left[\sum_{i_{1},\ldots,i_{n}\in\mathcal{A}_{N}}h(X_{i_{1}},\ldots,X_{i_{n}},Y\setminus\\{X_{i_{1}},\ldots,X_{i_{n}}\\})\prod_{j=1}^{n}(1-B_{i_{j}})\prod_{j=1}^{n}p(X_{i_{j}})\right]=$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathop{\sum\nolimits^{\neq}}_{x_{1},\ldots,x_{n}\in
Y}h(x_{1},\ldots,x_{n},Y\setminus\\{x_{1},\ldots,x_{n}\\})\prod_{i=1}^{n}p(x_{i})\right],$

which proves (1).

Next, let
$\psi_{0}(\breve{{\mathbf{x}}})=\\{x:(x,m)\in\breve{{\mathbf{x}}}\cap
S\times\\{0\\}\\}$, $\breve{{\mathbf{x}}}\in\breve{{\mathcal{X}}}$, and
consider any non-negative or integrable
$h:S^{n}\times{\mathcal{X}}\to\mathbb{R}$. Applying the GNZ formula to the
left-hand side of (1) yields

$\displaystyle\mathbb{E}\left[\mathop{\sum\nolimits^{\neq}}_{x_{1},\ldots,x_{n}\in
Z}h(x_{1},\ldots,x_{n},Y)\prod_{i=1}^{n}(1-p(x_{i}))\right]=$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathop{\sum\nolimits^{\neq}}_{(x_{1},m_{1}),\ldots,(x_{n},m_{n})\in\breve{X}}h(x_{1},\ldots,x_{n},\psi_{0}(\breve{X}))\prod_{i=1}^{n}m_{i}(1-p(x_{i}))\right]$
$\displaystyle=$
$\displaystyle\int_{S^{n}}\mathbb{E}\left[\prod_{i=1}^{n}(1-p(u_{i}))h(u_{1},\ldots,u_{n};\psi_{0}(\breve{X}))\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\right]\mathrm{d}u_{1}\cdots\mathrm{d}u_{n}$
$\displaystyle=$
$\displaystyle\int_{S^{n}}\mathbb{E}\left[\prod_{i=1}^{n}(1-p(u_{i}))h(u_{1},\ldots,u_{n};Y)\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\right]\mathrm{d}u_{1}\cdots\mathrm{d}u_{n},$

since the reference measure on the mark space is the counting measure on
${\mathcal{M}}=\\{0,1\\}$. On the other hand, applying the GNZ formula to the
right hand side of (1) yields

$\displaystyle\mathbb{E}\left[\mathop{\sum\nolimits^{\neq}}_{x_{1},\ldots,x_{n}\in
Y}h(x_{1},\ldots,x_{n},Y\setminus\\{x_{1},\ldots,x_{n}\\})\prod_{i=1}^{n}p(x_{i})\right]=$
$\displaystyle=$
$\displaystyle\int_{S^{n}}\mathbb{E}\left[\prod_{i=1}^{n}p(u_{i})h(u_{1},\ldots,u_{n};Y)\lambda_{Y}^{(n)}(u_{1},\ldots,u_{n};Y)\right]\mathrm{d}u_{1}\cdots\mathrm{d}u_{n}.$

The equality of these two expressions for arbitrary
$h:S^{n}\times{\mathcal{X}}\to\mathbb{R}$ yields that for almost every
$u_{1},\ldots,u_{n}\in S^{n}$ and every non-negative or integrable
$h^{*}:{\mathcal{X}}\to\mathbb{R}$ we have that

$\mathbb{E}\left[h^{*}(Y)\left(\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})-\frac{\prod_{i=1}^{n}p(u_{i})}{\prod_{i=1}^{n}(1-p(u_{i}))}\lambda_{Y}^{(n)}(u_{1},\ldots,u_{n};Y)\right)\right]=0,$

which concludes the proof.

∎

###  B.2 Proof of Lemma 3.1

###### Proof of Lemma 3.1.

This is an immediate consequence of the GNZ formula (2.7) and the Campbell
formula (2.3). ∎

###  B.3 Proof of Theorem 2

###### Proof of Theorem 2.

For ease of notation, we sometimes write $\mathrm{d}u$ for
$\mathrm{d}u_{1}\cdots\mathrm{d}u_{n}$.

#### When $h$ and $\xi$ are of the form (3.2)

When $h$ and $\xi$ are of the form (3.2), by the Campbell formula we have that

$\displaystyle\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)]=$
$\displaystyle\mathbb{E}\left[\sum_{(x_{1},\ldots,x_{n})\in Z_{\neq}^{n}\cap
A}h(x_{1},\ldots,x_{n})\right]-\int_{A}h(u_{1},\ldots,u_{n})\xi(u_{1},\ldots,u_{n})\mathrm{d}u$
$\displaystyle=$
$\displaystyle\int_{A}h(u_{1},\ldots,u_{n})\left(\rho_{Z}^{(n)}(u_{1},\ldots,u_{n})\mathrm{d}u-\xi(u_{1},\ldots,u_{n})\right)\mathrm{d}u,\qquad
A\subseteq S^{n}.$

Hence, $\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)]=0$ for any (bounded)
$A\subseteq S^{n}$ and function $h$ if and only if

$\xi(u_{1},\ldots,u_{n})=\rho_{Z}^{(n)}(u_{1},\ldots,u_{n}),$

for $|\cdot|^{n}$-almost every $(u_{1},\ldots,u_{n})\in S^{n}$; see e.g.
Møller and Waagepetersen (2004, Section 2.3.3).

We further have that

$\displaystyle\operatorname{Var}(\mathcal{I}_{\xi}^{h}(A;Z,Y))=$
$\displaystyle\operatorname{Var}\left(\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap A}h(x_{1},\ldots,x_{n})\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap A}h(x_{1},\ldots,x_{n})\right)^{2}\right]$
$\displaystyle-\left(\int_{A}h(u_{1},\ldots,u_{n})\rho_{Z}^{(n)}(u_{1},\ldots,u_{n})\mathrm{d}u\right)^{2},$

where, by Poinas et al. (2019, Equation (B.3)),

$\displaystyle\mathbb{E}\left[\left(\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap A}h(x_{1},\ldots,x_{n})\right)^{2}\right]=$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\mathop{\sum\nolimits^{\neq}}_{x_{1},\ldots,x_{n}\in
Z}{\mathbf{1}}\\{(x_{1},\ldots,x_{n})\in
A\\}h(x_{1},\ldots,x_{n})\right)^{2}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\sum_{{\mathbf{y}}\subseteq
Z}n!{\mathbf{1}}\\{\\#{\mathbf{y}}=n\\}{\mathbf{1}}\\{{\mathbf{y}}\in
A\\}h({\mathbf{y}})\right)^{2}\right]$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{n}\frac{(n!)^{2}}{(2n-j)!}\binom{n}{j}\binom{2n-j}{n}\int_{S^{2n-j}}h(u_{1},\ldots,u_{n})h(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j})\times$
$\displaystyle\times{\mathbf{1}}\\{(u_{1},\ldots,u_{n})\in
A\\}{\mathbf{1}}\\{(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j})\in A\\}\times$
$\displaystyle\times\rho_{Z}^{(2n-j)}(u_{1},\ldots,u_{2n-j})\mathrm{d}u_{1}\cdots\mathrm{d}u_{2n-j}.$

#### When $h$ and $\xi$ are of the form (3.1)

We here let $h$ and $\xi$ be of the form (3.1) and start by defining

$\displaystyle H_{1}(A)=$ $\displaystyle\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap A}h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\}),$
$\displaystyle H_{2}(A)=$
$\displaystyle\int_{A}h(u_{1},\ldots,u_{n};Y)\xi(u_{1},\ldots,u_{n};Y)\mathrm{d}u,$
$\displaystyle\mu_{1}(A)=$ $\displaystyle\mathbb{E}\left[H_{1}(A)\right],$
$\displaystyle\mu_{2}(A)=$ $\displaystyle\mathbb{E}\left[H_{2}(A)\right],\quad
A\subseteq S^{n},$

where we note that

$\displaystyle\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)]=$
$\displaystyle\mu_{1}(A)-\mu_{2}(A),$
$\displaystyle\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)^{2}]=$
$\displaystyle\mathbb{E}[H_{1}(A)^{2}]+\mathbb{E}[H_{2}(A)^{2}]-2\mathbb{E}[H_{1}(A)H_{2}(A)],$
$\displaystyle\operatorname{Var}(\mathcal{I}_{\xi}^{h}(A;Z,Y))=$
$\displaystyle\mathbb{E}[H_{1}(A)^{2}]+\mathbb{E}[H_{2}(A)^{2}]-2\mathbb{E}[H_{1}(A)H_{2}(A)]-(\mu_{1}(A)-\mu_{2}(A))^{2}.$
(B.1)

Next, recall the associated marked point process $\breve{X}$ in (2.20), with
Papangelou conditional intensity $\breve{\lambda}^{(n)}(\cdot)$. Given
$\psi_{0}(\breve{{\mathbf{x}}})=\\{x:(x,m)\in\breve{{\mathbf{x}}}\cap
S\times\\{0\\}\\}$, $\breve{{\mathbf{x}}}\in\breve{{\mathcal{X}}}$, by the GNZ
formula,

$\displaystyle\mu_{1}(A)=$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{(x_{1},\ldots,x_{n})\in Z_{\neq}^{n}\cap
A}h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{((x_{1},m_{1}),\ldots,(x_{n},m_{n}))\in\breve{X}_{\neq}^{n}\cap(A\times{\mathcal{M}}^{n})}\prod_{i=1}^{n}m_{i}h(x_{1},\ldots,x_{n};\psi_{0}(\breve{X}\setminus\\{(x_{1},m_{1}),\ldots,(x_{n},m_{n})\\}))\right]$
$\displaystyle=$
$\displaystyle\int_{A}\sum_{m_{1},\ldots,m_{n}\in\\{0,1\\}}\mathbb{E}\left[\prod_{i=1}^{n}m_{i}h(u_{1},\ldots,u_{n};\psi_{0}(\breve{X}))\breve{\lambda}^{(n)}((u_{1},m_{1}),\ldots,(u_{n},m_{n});\breve{X})\right]\mathrm{d}u$
$\displaystyle=$
$\displaystyle\int_{A}\mathbb{E}\left[h(u_{1},\ldots,u_{n};\psi_{0}(\breve{X}))\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\right]\mathrm{d}u$
$\displaystyle=$
$\displaystyle\int_{A}\mathbb{E}\left[h(u_{1},\ldots,u_{n};Y)\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\right]\mathrm{d}u,$

since the reference measure on the mark space is given by the counting measure
on the mark space ${\mathcal{M}}=\\{0,1\\}$. On the other hand, by the Fubini-
Tonelli theorem,

$\displaystyle\mu_{2}(A)=$
$\displaystyle\mathbb{E}\left[\int_{A}h(u_{1},\ldots,u_{n};Y)\xi(u_{1},\ldots,u_{n};Y)\mathrm{d}u\right]$
$\displaystyle=$
$\displaystyle\int_{A}\mathbb{E}\Big{[}h(u_{1},\ldots,u_{n};Y)\xi(u_{1},\ldots,u_{n};Y)\Big{]}\mathrm{d}u.$

Hence, $\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)]=0$ for any (bounded)
$A\subseteq S^{n}$ if and only if

$\displaystyle\mathbb{E}\left[h(u_{1},\ldots,u_{n};Y)\left(\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})-\xi(u_{1},\ldots,u_{n};Y)\right)\right]=0,$

for $|\cdot|^{n}$-almost every $(u_{1},\ldots,u_{n})\in S^{n}$; see e.g.
Møller and Waagepetersen (2004, Section 2.3.3). Moreover, under the assumption
that
$\mathbb{E}[\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})^{2}]<\infty$
and $\mathbb{E}[h(u_{1},\ldots,u_{n};Y)^{2}]<\infty$, $L_{2}$-projection
yields that

$\xi(u_{1},\ldots,u_{n};Y)=\mathbb{E}[\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})|Y].$

We next turn to the variance. Similarly to Poinas et al. (2019, Equation
(B.3)), we find that

$\displaystyle\mathbb{E}[H_{1}(A)^{2}]=$ $\displaystyle=$
$\displaystyle\mathbb{E}\Bigg{[}\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}}\sum_{(y_{1},\ldots,y_{n})\in
Z_{\neq}^{n}}{\mathbf{1}}\\{(x_{1},\ldots,x_{n}),(y_{1},\ldots,y_{n})\in
A\\}\times$ $\displaystyle\times
h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})h(y_{1},\ldots,y_{n};Y\setminus\\{y_{1},\ldots,y_{n}\\})\Bigg{]}$
$\displaystyle=$ $\displaystyle
n!^{2}\mathbb{E}\left[\sum_{{\mathbf{x}}=\\{x_{1},\ldots,x_{n}\\}\subseteq
Z}\sum_{{\mathbf{y}}=\\{y_{1},\ldots,y_{n}\\}\subseteq
Z}{\mathbf{1}}\\{{\mathbf{x}},{\mathbf{y}}\in
A\\}h({\mathbf{x}};Y\setminus{\mathbf{x}})h({\mathbf{y}};Y\setminus{\mathbf{y}})\right]$
$\displaystyle=$ $\displaystyle
n!^{2}\sum_{j=0}^{n}\mathbb{E}\Bigg{[}\sum_{{\mathbf{x}}=\\{x_{1},\ldots,x_{n}\\}\subseteq
Z}\sum_{{\mathbf{y}}=\\{y_{1},\ldots,y_{n}\\}\subseteq
Z}{\mathbf{1}}\\{\\#({\mathbf{x}}\cap{\mathbf{y}})=j\\}{\mathbf{1}}\\{{\mathbf{x}},{\mathbf{y}}\in
A\\}h({\mathbf{x}};Y\setminus{\mathbf{x}})h({\mathbf{y}};Y\setminus{\mathbf{y}})\Bigg{]},$

where the factor $n!^{2}$ comes from the fact that when we go from $n$-subsets
to $n$-tuples we count the same thing $n!$ times; we can rearrange
$(x_{1},\ldots,x_{n})$ in $n!$ different ways. In the last sum, assuming that
$\\#({\mathbf{x}}\cap{\mathbf{y}})=j$, i.e. that ${\mathbf{x}}$ and
${\mathbf{y}}$ have $j$ elements $x_{i}=y_{i^{\prime}}\in Z$ in common, there
are $\binom{n}{j}$ ways in which the elements in
${\mathbf{x}}\cap{\mathbf{y}}$ can be chosen from ${\mathbf{x}}$ and
${\mathbf{y}}$. The remaining $2n-j$ elements now need to be assigned to
${\mathbf{x}}\setminus({\mathbf{x}}\cap{\mathbf{y}})$ and
${\mathbf{y}}\setminus({\mathbf{x}}\cap{\mathbf{y}})$. There are
$\binom{2n-j}{n-j}=\binom{2n-j}{n}$ ways to assign elements of
$({\mathbf{x}}\cap{\mathbf{y}})^{c}$ to
${\mathbf{x}}\setminus({\mathbf{x}}\cap{\mathbf{y}})$ so that
$\\#{\mathbf{x}}=n$; the remaining elements will automatically be assigned to
${\mathbf{y}}\setminus({\mathbf{x}}\cap{\mathbf{y}})$. In the following, we
let $z_{1},\ldots,z_{j}$ denote the elements in
${\mathbf{x}}\cap{\mathbf{y}}$, $z_{j+1},\ldots,z_{n}$ the elements only in
${\mathbf{x}}$, and $z_{n+1},\ldots,z_{2n-j}$ the ones only in ${\mathbf{y}}$.
Consequently,

$\displaystyle\mathbb{E}[H_{1}(A)^{2}]=$ $\displaystyle=$ $\displaystyle
n!^{2}\sum_{j=0}^{n}\binom{n}{j}\binom{2n-j}{n}\mathbb{E}\Bigg{[}\sum_{\\{z_{1},\ldots,z_{2n-j}\\}\subseteq
Z}{\mathbf{1}}\\{\\{z_{1},\ldots,z_{n}\\}\in A\\}\times$
$\displaystyle\times{\mathbf{1}}\\{\\{z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j}\\}\in
A\\}h(z_{1},\ldots,z_{n};Y\setminus\\{z_{1},\ldots,z_{n}\\})\times$
$\displaystyle\times
h(z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j};Y\setminus\\{z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j}\\})\Bigg{]}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{n}\binom{n}{j}\binom{2n-j}{n}\frac{n!^{2}}{(2n-j)!}\mathbb{E}\Bigg{[}\sum_{(z_{1},\ldots,z_{2n-j})\subseteq
Z_{\neq}^{2n-j}}{\mathbf{1}}\\{(z_{1},\ldots,z_{n})\in A\\}\times$
$\displaystyle\times{\mathbf{1}}\\{(z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j})\in
A\\}h(z_{1},\ldots,z_{n};Y\setminus\\{z_{1},\ldots,z_{n}\\})\times$
$\displaystyle\times
h(z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j};Y\setminus\\{z_{1},\ldots,z_{j},z_{n+1},\ldots,z_{2n-j}\\})\Bigg{]},$

where we note that

$\frac{(n!)^{2}}{(2n-j)!}\binom{n}{j}\binom{2n-j}{n}=\binom{n}{j}\frac{(n!)^{2}}{(2n-j)!}\frac{(2n-j)!}{n!(n-j)!}=j!\binom{n}{j}^{2}.$

By applying the GNZ formula to each term in the sum in the last equation, it
follows that

$\displaystyle\mathbb{E}[H_{1}(A)^{2}]=$ (B.2) $\displaystyle=$
$\displaystyle\sum_{j=0}^{n}j!\binom{n}{j}^{2}\int_{S^{2n-j}}{\mathbf{1}}\\{(u_{1},\ldots,u_{n}),(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j})\in
A\\}$
$\displaystyle\times\mathbb{E}\Big{[}h(u_{1},\ldots,u_{n};Y)h(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j};Y)\breve{\lambda}^{(2n-j)}((u_{1},1),\ldots,(u_{2n-j},1);\breve{X})\Big{]}$
$\displaystyle\times\mathrm{d}u_{1}\cdots\mathrm{d}u_{2n-j}.$

We further have that

$\displaystyle\mathbb{E}[H_{2}(A)^{2}]=$ $\displaystyle=$
$\displaystyle\int_{A}\int_{A}\mathbb{E}[h(u_{1},\ldots,u_{n};Y)h(v_{1},\ldots,v_{n};Y)\xi(u_{1},\ldots,u_{n};Y)\xi(v_{1},\ldots,v_{n};Y)]\mathrm{d}u\mathrm{d}v$
(B.3)

and

$\displaystyle\mathbb{E}[H_{1}(A)H_{2}(A)]=$ $\displaystyle=$
$\displaystyle\mathbb{E}\Bigg{[}\sum_{(x_{1},\ldots,x_{n})\in Z_{\neq}^{n}\cap
A}h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})\times$
$\displaystyle\times\int_{A}h(v_{1},\ldots,v_{n};Y)\xi(v_{1},\ldots,v_{n};Y)\mathrm{d}v\Bigg{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}\Bigg{[}\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap
A}h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})\times$
$\displaystyle\times\Bigg{(}\int_{A}h(v_{1},\ldots,v_{n};(Y\setminus\\{x_{1},\ldots,x_{n}\\})\cup\\{x_{1},\ldots,x_{n}\\})\times$
$\displaystyle\times\xi(v_{1},\ldots,v_{n};(Y\setminus\\{x_{1},\ldots,x_{n}\\})\cup\\{x_{1},\ldots,x_{n}\\})\mathrm{d}v\Bigg{)}\Bigg{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}\left[\sum_{(x_{1},\ldots,x_{n})\in
Z_{\neq}^{n}\cap
A}\widetilde{h}(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})\right].$

where

$\displaystyle\widetilde{h}(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})=$
$\displaystyle h(x_{1},\ldots,x_{n};Y\setminus\\{x_{1},\ldots,x_{n}\\})$
$\displaystyle\times\int_{A}h(v;(Y\setminus\\{x_{1},\ldots,x_{n}\\})\cup\\{x_{1},\ldots,x_{n}\\})$
$\displaystyle\times\xi(v;(Y\setminus\\{x_{1},\ldots,x_{n}\\})\cup\\{x_{1},\ldots,x_{n}\\})\mathrm{d}v.$

Hence,

$\displaystyle\mathbb{E}[H_{1}(A)H_{2}(A)]=$ $\displaystyle=$
$\displaystyle\int_{A}\mathbb{E}\Bigg{[}\widetilde{h}(u_{1},\ldots,u_{n};Y)\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\Bigg{]}\mathrm{d}u$
$\displaystyle=$
$\displaystyle\int_{A}\int_{A}\mathbb{E}\Bigg{[}h(u_{1},\ldots,u_{n};Y)h(v_{1},\ldots,v_{n};Y\cup\\{u_{1},\ldots,u_{n}\\})\times$
$\displaystyle\times\xi(v_{1},\ldots,v_{n};Y\cup\\{u_{1},\ldots,u_{n}\\})\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\Bigg{]}\mathrm{d}u\mathrm{d}v$
(B.4)

and, consequently, by combining (B.1) with (B.2), (B.3) and (B.4), the
variance is given by

$\displaystyle\operatorname{Var}(\mathcal{I}_{\xi}^{h}(A;Z,Y))=\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)^{2}]-\mathbb{E}[\mathcal{I}_{\xi}^{h}(A;Z,Y)]^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{n}j!\binom{n}{j}^{2}\int_{S^{2n-j}}{\mathbf{1}}\\{(u_{1},\ldots,u_{n}),(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j})\in
A\\}\mathbb{E}\Big{[}h(u_{1},\ldots,u_{n};Y)$ $\displaystyle\times
h(u_{1},\ldots,u_{j},u_{n+1},\ldots,u_{2n-j};Y)\breve{\lambda}^{(2n-j)}((u_{1},1),\ldots,(u_{2n-j},1);\breve{X})\Big{]}\mathrm{d}u_{1}\cdots\mathrm{d}u_{2n-j}$
$\displaystyle+\int_{A}\int_{A}\mathbb{E}[h(u_{1},\ldots,u_{n};Y)h(v_{1},\ldots,v_{n};Y)\times$
$\displaystyle\times\xi^{n}(u_{1},\ldots,u_{n};Y)\xi^{n}(v_{1},\ldots,v_{n};Y)]\mathrm{d}u\mathrm{d}v$
$\displaystyle-2\Bigg{(}\int_{A}\int_{A}\mathbb{E}\Big{[}h(u_{1},\ldots,u_{n};Y)h(v_{1},\ldots,v_{n};Y\cup\\{u_{1},\ldots,u_{n}\\})\times$
$\displaystyle\times\xi(v_{1},\ldots,v_{n};Y\cup\\{u_{1},\ldots,u_{n}\\})\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})\Big{]}\mathrm{d}u\mathrm{d}v\Bigg{)}$
$\displaystyle-\Bigg{(}\mathbb{E}\left[h(u_{1},\ldots,u_{n};Y)\left(\breve{\lambda}^{(n)}((u_{1},1),\ldots,(u_{n},1);\breve{X})-\xi(u_{1},\ldots,u_{n};Y)\right)\right]\Bigg{)}^{2}$

when $h$ and $\xi$ are of the form (3.1).

∎

###  B.4 Proof of Theorem 3

###### Proof of Theorem 3.

By (4.7) and (5.1), we have that

$\mathcal{L}_{1}(\theta)=\frac{1}{k}\sum_{i=1}^{k}|\mathcal{I}_{\xi_{\theta}^{1}}^{h}(W;{\mathbf{x}}_{i}^{T})|=\frac{1}{k}\sum_{i=1}^{k}\left|\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)-\theta(1-p)\int_{W}h(u)\mathrm{d}u\right|{\mathbf{1}}\\{\\#{\mathbf{x}}_{i}^{T}\geq
1\\},$

which has the same minimum with respect to $\theta$ as

$\frac{k\mathcal{L}_{1}(\theta)}{(1-p)|\int_{W}h(u)\mathrm{d}u|}=\sum_{i=1}^{k}\left|\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}-\theta\right|{\mathbf{1}}\\{\\#{\mathbf{x}}_{i}^{T}\geq
1\\}.$

The derivative of the last function is defined for almost all $\theta$ by

$\sum_{i=1}^{k}\mathrm{sgn}\left(\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}-\theta\right){\mathbf{1}}\\{\\#{\mathbf{x}}_{i}^{T}\geq
1\\},$

where for all $x\in\mathbb{R}$, $\mathrm{sgn}(x)={\mathbf{1}}\\{x\geq
0\\}+{\mathbf{1}}\\{x<0\\}$. This sum is null only if on the set of all
$i\in\mathcal{T}_{k}$, half of the $(\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x))/((1-p)\int_{W}h(u)\mathrm{d}u)$ are greater than $\theta$ and half are
less than $\theta$. Therefore we have that a minimiser is given by
$\widehat{\theta}_{1}(\\{({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})\\}_{i=1}^{k},p,W,h)=\operatorname{med}\\{\widehat{\theta}_{i}^{h}:i\in\mathcal{T}_{k}\\}$.

Elementary calculus shows that the derivatives of $\mathcal{L}_{2}$ and
$\mathcal{L}_{3}$ with respect to $\theta$ are null if $\theta$ satisfies

$\sum_{i=1}^{k}\left(\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)-\theta(1-p)\int_{W}h(u)\mathrm{d}u\right){\mathbf{1}}\\{\\#{\mathbf{x}}_{i}^{T}\geq
1\\}=0$

which is equivalent to finding $\theta$ such that

$\sum_{i=1}^{k}\left(\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}-\theta\right){\mathbf{1}}\\{\\#{\mathbf{x}}_{i}^{T}\geq
1\\}=0.$

Therefore, for $j=2,3$ we have that

$\displaystyle\widehat{\theta}_{j}(\\{({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})\\}_{i=1}^{k},p,W,h)=$
$\displaystyle\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}$ $\displaystyle=$
$\displaystyle\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}.$

Further, when $j=2,3$, for any point configuration ${\mathbf{x}}$ we have that

$\displaystyle\mathbb{E}[\widehat{\theta}_{j}(\\{({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})\\}_{i=1}^{k},p,W,h)]=$
$\displaystyle\mathbb{E}\left[\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\left.\mathbb{E}\left[\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x){\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}\right|\mathcal{T}_{k}\right]\right]$
$\displaystyle=$ $\displaystyle\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}\mathbb{E}\left[\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\mathbb{E}[{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}|\mathcal{T}_{k}]\right]$
$\displaystyle=$ $\displaystyle\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}\mathbb{E}[\mathbb{E}[{\mathbf{1}}\\{x\in{\mathbf{x}}_{1}^{T}\\}|\mathcal{T}_{k}]]=\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)}{\int_{W}h(u)\mathrm{d}u}$

by the law of total expectation, since the validation/training set assignments
are independent and identically distributed with
$\mathbb{P}(x\in{\mathbf{x}}_{i}^{T})=1-p$, $i=1,\ldots,k$,
$x\in{\mathbf{x}}$. It is straightforward to see that the same holds true if
$\mathcal{T}_{k}=\\{1,\ldots,k\\}$. ∎

###  B.5 Proof of Lemma 5.1

We here want to show that

$\operatorname{Var}\left(\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\widehat{\theta}_{i}^{h}\right)\geq\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\widehat{\theta}_{i}^{h}\right)=\frac{p}{k(1-p)}\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)^{2}}{(\int_{W}h(u)\mathrm{d}u)^{2}}.$

Staring with the right hand side, we have that

$\displaystyle\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\widehat{\theta}_{i}^{h}\right)=$
$\displaystyle\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x){\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}\right)=\frac{\operatorname{Var}\left(\sum_{x\in{\mathbf{x}}\cap
W}b_{x}h(x)\right)}{k^{2}(1-p)^{2}(\int_{W}h(u)\mathrm{d}u)^{2}},$

where, for any $x\in X$,
$b_{x}=\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\cap W\\}$ follows
a binomial distribution with parameters $k$ and $(1-p)$, with
$\mathbb{E}[b_{x}]=k(1-p)$ and
$\mathbb{E}[b_{x}^{2}]=kp(1-p)+(k(1-p))^{2}=k(1-p)(p+k(1-p))$. Since the
events $x\in{\mathbf{x}}_{i}^{T}$ and $y\in{\mathbf{x}}_{i^{\prime}}^{T}$ are
independent for all $i,i^{\prime}=1,\ldots,k$ and $x\neq y$, $b_{x}$ and
$b_{y}$ are independent when $x\neq y$, whereby

$\displaystyle\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\widehat{\theta}_{i}^{h}\right)=$
$\displaystyle\frac{\sum_{x\in{\mathbf{x}}\cap
W}\operatorname{Var}(b_{x})h(x)^{2}}{k^{2}(1-p)^{2}(\int_{W}h(u)\mathrm{d}u)^{2}}=\frac{p}{k(1-p)}\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)^{2}}{(\int_{W}h(u)\mathrm{d}u)^{2}}.$

We next turn to the inequality. Since $\sum_{x\in\emptyset}h(x)=0$, we a.s.
have that

$\displaystyle\widehat{\theta}_{j}(\\{({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})\\}_{i=1}^{k},p,W,h)=$
$\displaystyle\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\widehat{\theta}_{i}^{h}=\frac{1}{\\#\mathcal{T}_{k}}\sum_{i=1}^{k}\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}$ $\displaystyle=$
$\displaystyle\frac{k}{\\#\mathcal{T}_{k}}\frac{1}{k}\sum_{i=1}^{k}\frac{\sum_{x\in{\mathbf{x}}_{i}^{T}\cap
W}h(x)}{(1-p)\int_{W}h(u)\mathrm{d}u}$ $\displaystyle=$
$\displaystyle\frac{k}{\\#\mathcal{T}_{k}}\frac{\sum_{x\in{\mathbf{x}}\cap
W}h(x)\frac{1}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u},$

whereby

$\displaystyle\operatorname{Var}(\widehat{\theta}_{j}(\\{({\mathbf{x}}_{i}^{T},{\mathbf{x}}_{i}^{V})\\}_{i=1}^{k},p,W,h))$
$\displaystyle=\frac{\operatorname{Var}\left(\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)}{(1-p)^{2}(\int_{W}h(u)\mathrm{d}u)^{2}}.$

The numerator is given by

$\displaystyle\operatorname{Var}\left(\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)=$
$\displaystyle\mathbb{E}\left[\left(\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}h(x)\frac{1}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)^{2}\right]$
$\displaystyle-\mathbb{E}\left[\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right]^{2},$
(B.5)

where, by Theorem 3,

$\mathbb{E}\left[\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right]=(1-p)\sum_{x\in{\mathbf{x}}\cap
W}h(x).$

We further have that

$\displaystyle\operatorname{Var}\left(\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)=$
$\displaystyle\mathbb{E}\left[\left(\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)^{2}\right]$
$\displaystyle-\mathbb{E}\left[\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right]^{2}$
(B.6)

and

$\mathbb{E}\left[\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right]=\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}\mathbb{E}[{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}]=(1-p)\sum_{x\in{\mathbf{x}}\cap
W}h(x).$

Since $k\geq\\#\mathcal{T}_{k}$, from the two expectation equalities and
(B.5)-(B.6) it follows that

$\operatorname{Var}\left(\frac{k}{\\#\mathcal{T}_{k}}\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right)\geq\operatorname{Var}\left(\sum_{x\in{\mathbf{x}}\cap
W}\frac{h(x)}{k}\sum_{i=1}^{k}{\mathbf{1}}\\{x\in{\mathbf{x}}_{i}^{T}\\}\right),$

which concludes the proof.

###  B.6 Proof of Lemma 5.2

###### Proof of Lemma 5.2.

By Campbell’s theorem, we have

$\displaystyle\mathbb{E}[\widehat{\theta}_{j}(\\{(X_{i}^{T},X_{i}^{V})\\}_{i=1}^{k},p,W,h)]=$
$\displaystyle\mathbb{E}\left[\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\left.\mathbb{E}\left[\frac{\sum_{x\in
X\cap W}h(x){\mathbf{1}}\\{x\in
X_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}\right|\mathcal{T}_{k}\right]\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\frac{1}{\\#\mathcal{T}_{k}}\sum_{i\in\mathcal{T}_{k}}\frac{(1-p)\theta_{0}\int_{W}h(u)\mathrm{d}u}{(1-p)\int_{W}h(u)\mathrm{d}u}\right]=\theta_{0}$

and

$\displaystyle\mathbb{E}\left[\frac{1}{k}\sum_{i=1}^{k}\frac{\sum_{x\in X\cap
W}h(x){\mathbf{1}}\\{x\in X_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{k}\sum_{i=1}^{k}\left.\mathbb{E}\left[\frac{\sum_{x\in
X\cap W}h(x){\mathbf{1}}\\{x\in
X_{i}^{T}\\}}{(1-p)\int_{W}h(u)\mathrm{d}u}\right|\mathcal{T}_{k}\right]\right]=\theta_{0}.$

Hence, by the last two equations, we can follow the same calculus as in the
proof of Lemma 5.1 to prove that

$\operatorname{Var}\widehat{\theta}_{j}((X_{i}^{T},X_{i}^{V}),p,W,h)\geq\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\widehat{\theta}((X_{i}^{T},X_{i}^{V}),p,W,h)\right).$

For all $x\in X$, let $b_{x}=\sum_{i\in\mathcal{T}_{k}}{\mathbf{1}}\\{x\in
X_{i}^{T}\cap W\\}$ which, conditionally on $\mathcal{T}_{k}=k$, follows a
binomial distribution with parameters $k$ and $(1-p)$ so that
$\mathbb{E}[b_{x}]=k(1-p)$ and
$\mathbb{E}[b_{x}^{2}]=kp(1-p)+(k(1-p))^{2}=k(1-p)(p+k(1-p))$. Thus,

$\displaystyle\operatorname{Var}\left(\frac{1}{k}\sum_{i=1}^{k}\widehat{\theta}((X_{i}^{T},X_{i}^{V}),p,W,h)\right)=$
$\displaystyle=$ $\displaystyle\frac{\operatorname{Var}\left(\sum_{x\in X\cap
W}b_{x}h(x)\right)}{(k(1-p)\int_{W}h(u)\mathrm{d}u)^{2}}$ $\displaystyle=$
$\displaystyle\frac{\sum_{x\in X\cap
W}h(x)^{2}\mathbb{E}[b_{x}^{2}]+\sum_{x,y\in X\cap
W}^{\neq}h(x)h(y)\mathbb{E}[b_{x}b_{y}]-\mathbb{E}[\sum_{x\in X\cap
W}b_{x}h(x)]^{2}}{(k(1-p)\int_{W}h(u)\mathrm{d}u)^{2}}$ $\displaystyle=$
$\displaystyle\frac{\frac{p+k(1-p)}{k(1-p)}\theta_{0}\int_{W}h(u)^{2}\mathrm{d}u+\int_{W}\int_{W}h(u_{1})h(u_{2})\rho^{(2)}(u_{1},u_{2})\mathrm{d}u_{1}\mathrm{d}u_{2}-(\theta_{0}\int_{W}h(u)\mathrm{d}u)^{2}}{(\int_{W}h(u)\mathrm{d}u)^{2}}$
$\displaystyle=$
$\displaystyle\left(\frac{p}{(1-p)k}+1\right)\theta_{0}\frac{\int_{W}h(u)^{2}\mathrm{d}u}{(\int_{W}h(u)\mathrm{d}u)^{2}}+\theta_{0}^{2}\left(\frac{\int_{W}\int_{W}h(u_{1})h(u_{2})g_{X}^{(2)}(u_{1},u_{2})\mathrm{d}u_{1}\mathrm{d}u_{2}}{(\int_{W}h(u)\mathrm{d}u)^{2}}-1\right).$

∎