2. Main results

As for the estimators (1.1) and (1.2), throughout this paper we consider a point process X with intensity $1/t$ , and we apply Newton–Cotes quadratures to functions evaluated at the points of X.

Let $n\in \mathbb{N}$ be given. We recall the definition of the nth-order Newton–Cotes estimator $\hat{V}_n({\kern1.3pt}f)$ from [Reference Kiderlen and Dorph-Petersen6] for a fixed realization of X. On the interval from a point $x_0\in X$ to its nth right neighbour in X, say $x_n$ , the function f is approximated by a polynomial of degree at most $n\in\mathbb{N}$ passing through the points $\{x_j,f(x_j)\}_{j=0}^{n}$ , where $x_1<\cdots<x_{n-1}$ are the ordered points in $X\cap (x_0,x_n)$ . $\hat{V}_n({\kern1.3pt}f)$ is then an average of the integral of the concatenation of such approximations (composite rule) with respect to the starting point chosen. The estimator $\hat{V}_n({\kern1.3pt}f)$ turns out to be a weighted average of f over all points in X,

(2.1) \begin{equation} \hat{V}_n({\kern1.3pt}f)=\sum_{x\in X} \alpha(x;\,X)\, f(x),\end{equation}

where the weights satisfy $\alpha(x;\,X)=t\alpha(x/t;\,X/t)$ for all $x\in X$ , where (when considered random) $X/t$ is of unit intensity; see (3.4) for details. We will see in Remark 5.1 that $\alpha(x;\,X)=t$ when $X=t(U+\mathbb{Z})$ is an equidistant process, and therefore, Newton–Cotes estimators of any order coincide with (1.2) in the equidistant case.

When applying the estimator on randomized sampling points, we work under the general assumption, stated below, that a typical distance between two consecutive points has finite positive and negative moments of all orders.

Assumption 2.1

(2.2) \begin{equation} {\mathbb{E}}^0 h_1^j <\infty \qquad \textit{for all } j\in\mathbb{Z}.\end{equation}

Here ${\mathbb{E}}^0$ is the expectation under the Palm distribution of X, that is, the distribution of X given that $0\in X$ (see e.g. [Reference Schneider and Weil10, Section 3.3]), and $h_1$ is the lag between 0 and its right neighbour in X. Note that (2.2) holds for X if and only if it holds for aX for all $a>0$ . Assumption 2.1 is certainly not necessary for the results to hold for a given n, but finding a necessary and sufficient condition appears to be quite technical. Our first result shows the unbiasedness of $\hat{V}_n({\kern1.3pt}f)$ .

Theorem 2.1. Let $n\in \mathbb{N}$ be given and assume that X is a stationary point process such that Assumption 2.1 is satisfied. Then $\hat{V}_n({\kern1.3pt}f)$ is unbiased:

\[ {\mathbb{E}} \hat{V}_n({\kern1.3pt}f)=\int_\mathbb{R} f(x) {\text{d}} x\]

for all integrable and real-valued functions f with compact support.

This will be shown in Section 5, where we also argue that Assumption 4.1 below, which is weaker than Assumption 2.1, is sufficient to ensure unbiasedness. We also remark that unbiasedness is known to hold for $n=1$ without integrability conditions and for $n=2$ under a condition weaker than Assumption 4.1; see Example 1 and Corollary 3 in [Reference Kiderlen and Dorph-Petersen6].

As in the case of classical quadrature, high-order quadrature reduces the discretization error when the measurement function is smooth. We adopt a smoothness condition which is in widespread use in stereological applications. For $m,p\in\mathbb{N}_0 \cup \{\infty\}$ , we say that a measurable function f with compact support is (m,p)-piecewise smooth if it is in $C^{m-1}(\mathbb{R})$ for $m\ge 1$ , and all derivatives up to order $m+p$ exist and are continuous except in at most finitely many points, where they may have finite jumps. Hence, if f is (m, p)-piecewise smooth, m is the smallest order of derivative of f which may have jumps; see e.g. [Reference Kiêu7] for details on such functions. For our results to hold, we require that $p\geq 1$ ; however, the exact value of p is otherwise irrelevant. We therefore state all results for (m, 1)-piecewise smooth functions. We say that a function f is exactly (m, 1)-piecewise smooth if it is (m,1)-piecewise smooth with discontinuous mth derivative. We let $\smash{{D}_{f^{(m)}}}$ denote the finite set of discontinuity points of $f^{(m)}$ , with

\begin{equation*} a\mapsto {J}_{f^{(m)}} (a) = \lim_{x\to a^+}f^{(m)}(x)-\lim_{x\to a^-}f^{(m)}(x)\end{equation*}

denoting the corresponding jump function.

Our second result expresses the discretization error

(2.3) \begin{equation} R^{(n)}({\kern1.3pt}f) = \hat{V}_n({\kern1.3pt}f) - \int_\mathbb{R} f(x) {\text{d}} x\end{equation}

in terms of higher-order derivatives of f. We state it for a realization of X; that is, we consider X as a deterministic, locally finite set of distinct points with convex hull $\text{Conv}(X)=\mathbb{R}$ .

Theorem 2.2. (Refined Peano kernel theorem for Newton–Cotes estimation.)

Let $n\in\mathbb{N}$ be fixed. Given X and $m\leq n$ , there exists a function $K_m$ such that

\begin{equation*} R^{(n)}({\kern1.3pt}f) = \int_\mathbb{R} f^{(m+1)}(r) K_m(r) {\text{d}} r + { \sum_{a\in {D}_{f^{(m)}} }} {J}_{f^{(m)}}(a) K_m(a)\end{equation*}

for all (m, 1)-piecewise smooth functions $f\,:\,\mathbb{R}\to\mathbb{R}$ .

As initially considered and shown in [Reference Matheron9], the variance of (1.2) in the equidistant case depends on jumps of high-order derivatives of the measurement function; see also [Reference Kiêu7, Chapter 5]. This is outlined in the following for comparison with the general case. Let $*$ denote the convolution operator and let the reflection $\check{f}$ of f be defined as $\check{f}(x)=f({-}x)$ . When the measurement function f is (m,1)-piecewise smooth, it can be shown [Reference Kiêu7, Corollary 5.8] that the so-called covariogram $g=f*\check f$ of f is $(2m+1,1)$ -piecewise smooth. When f is exactly (m,1)-piecewise smooth, one usually decomposes the variance of $\hat{V}({\kern1.3pt}f)$ as

(2.4) \begin{equation} \mathbb{V}\text{ar} (\hat{V}({\kern1.3pt}f)) = \mathbb{V}\text{ar}_E (\hat{V}({\kern1.3pt}f)) + Z(t) + o\big(t^{2m+2}\big)\end{equation}

when $t\downarrow0$ . The Zitterbewegung Z(t), which is of order $t^{2m+2}$ , is a finite sum of terms oscillating around 0, $o(t^{2m+2})$ is a low-order remainder, and the extension term

(2.5) \begin{equation} \mathbb{V}\text{ar}_E (\hat{V}({\kern1.3pt}f)) = t^{2m+2} g^{(2m+1)}(0^+) c_m\end{equation}

explains the overall trend of the variance. Here

\begin{equation*}c_m = -\frac{2 B_{2m+2}}{(2m+2)!}\neq 0, \end{equation*}

where $B_k$ is the kth Bernoulli number (see Section 6 below), and as such $c_m$ does not depend on t or the function f, other than through its order of smoothness.

Motivated by stereological applications, and adopting the naming from [Reference Ziegel, Baddeley, Dorph-Petersen and Jensen14], we will mainly work with two classes of point process models. Both models are defined as scalings of unit-intensity processes.

If X is from the perturbed model (with nondegenerate perturbations), the variance of (1.1) satisfies $\mathbb{V}\text{ar} (\hat{V}_0({\kern1.3pt}f)) = t^2 c' +Z_0(t) + o(t^2)$ when $m=0$ and $\mathbb{V}\text{ar} (\hat{V}_0({\kern1.3pt}f)) = t^3 c'' + o(t^3)$ when $m\geq 1$ as $t\downarrow 0$ . This was shown in [Reference Ziegel, Baddeley, Dorph-Petersen and Jensen14, Proposition 1], apart from the missing Zitterbewegung term $Z_0(t)$ of order $t^2$ in the first equation, which was omitted there as it was erroneously claimed that the last term in [Reference Ziegel, Baddeley, Dorph-Petersen and Jensen14, Eq. (A3)] is of order $o(t^{2m+2})$ . Hence, for all $m\geq 1$ , the rate of decrease of $V_0({\kern1.3pt}f)$ in the non-equidistant case is strictly smaller than in the equidistant case; cf. (2.5) for the latter. The behaviour is even worse in the model with cumulative errors, as $\mathbb{V}\text{ar} (\hat{V}_0({\kern1.3pt}f)) = t c''' + o(t)$ for all $m\geq 0$ ; see [Reference Ziegel, Baddeley, Dorph-Petersen and Jensen14, Proposition 2].

In order to formulate corresponding rates of decrease for Newton–Cotes estimators, we need the notion of an admissible point process. The Peano kernel in the definition of an admissible point process is explicitly given in (3.7) with $m=n$ .

From the definition (3.7) of $K_m$ it is easily seen that X is weakly/strongly admissible if and only if aX is weakly/strongly admissible for all constants $a>0$ . Admissibility is closely related to the ergodicity properties of the stationary field $K_n$ , and hence to those of X. In fact, if $K_n$ has an exponentially decaying $\alpha$ -mixing coefficient (see, for instance, [Reference Doukhan5, Subsection 1.3.2] for the definition of this coefficient), then [Reference Doukhan5, Theorem 3(1), p. 9] and the fact that ${\mathbb{E}} K_n(0)^{2+\varepsilon}<\infty$ , $\varepsilon>0$ , imply that $H_n(s)$ is exponentially decaying, and hence X is strongly n-admissible for all $n\in\mathbb{N}$ .

The covariance function $H_n$ need not be decaying for X to be strongly n-admissible. When X is from the perturbed model, the covariance function is closely related to Bernoulli functions (see Section 6), which are 1-periodic functions integrating to 0 on an interval of unit length. This is used in Lemma 7.1 to show that X is strongly n-admissible for all $n\in \mathbb{N}$ . Concerning the model with cumulative errors, we show in Lemma 7.3 that $H_n$ is indeed exponentially decaying when it is assumed that the holding times of the process have finite exponential moments. The proof relies on a result from [Reference Asmussen1] concerning the convergence rate of convolutions of the renewal measure of a pure renewal process.

Theorem 2.3. Let Let $n\in \mathbb{N}$ be given and assume that X is a stationary point process with intensity $1/t>0$ such that Assumption 2.1 holds. If f is (m, 1)-piecewise smooth and ${k=\min\{m,n\}}$ , the variance of the estimator (2.1) obeys

(2.6) \begin{equation} \mathbb{V}\text{ar}(\hat{V}_{n}({\kern1.3pt}f))\leq c t^{2k+2}\end{equation}

for some constant c, which does not depend on t.

If $m>n$ and X is strongly n-admissible, then

(2.7) \begin{equation} \mathbb{V}\text{ar}(\hat{V}_{n}({\kern1.3pt}f)) \leq c' t^{2n+3}\end{equation}

for some constant cʹ, which does not depend on t.

If f is exactly (m,1)-piecewise smooth with $m<n$ , the decrease rate in (2.6) is optimal. This is also true in the case $m=n$ if X is weakly n-admissible; see Remark 6.1.

When using the trapezoidal estimator, that is, $n=1$ , we have exact expressions for the asymptotic behaviour of the variance when X is from the perturbed model and the model with cumulative errors. In the perturbed case, the rate of decrease of the upper bound in (2.7) is optimal if the perturbations $E_i$ are nondegenerate.

Theorem 2.4. Let X be from the perturbed model with intensity 1/t, and let $\mu_k$ be the kth moment of the perturbations $E_i$ . Assume that the measurement function f is exactly (m, 1)-piecewise smooth with covariogram $g=f*\,\check{f}$ . Then, for $t\downarrow 0$ ,

(2.8) \begin{align} \mathbb{V}\textit{ar} (\hat{V}_{1}({\kern1.3pt}f)) = -t^{2} g'(0^+) (\mu_2 + \tfrac{1}{6}) + Z_0(t) + o(t^2) &&\qquad\qquad\qquad\textit{for}\quad m=0,\end{align}

(2.9) \begin{align} \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) = t^{4} g^{(3)}(0^+) \tfrac1{12} (2\mu_2 + 2\mu_4 + \tfrac{1}{30}) + Z_1(t) + o(t^4) &&\quad\textit{for}\quad m=1,\end{align}

(2.10) \begin{align}\mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) = t^{5} g^{(4)}(0)\tfrac1{8} (2\mu_4 + \mu_2\mu_4 - \mu_2^3 - \mu_3^2) + o(t^5) &&\,\,\,\,\quad\quad\textit{for}\quad m\geq 2,\end{align}

where the Zitterbewegung $Z_m(t)$ is given by (7.7). It is of order $t^{2m+2}$ , and it is a finite sum of terms oscillating around 0. Moreover, if $E_i$ has a density with a finite number of finite jumps and $m\geq2$ , the remainder $o(t^5)$ is explicitly given by

(2.11) \begin{equation}\begin{aligned} t^6 g^{(5)}(0^+) \tfrac1{720} \bigl(&{-}34\mu_2 - 90\mu_2^2 + 110\mu_4 + 180\mu_2\mu_4 \\ & - 180\mu_2^3 - 170\mu_3^2 + 8\mu_6 - \tfrac{1}{21}\bigr) + Z_2(t) + o(t^6).\end{aligned}\end{equation}

We compare these findings with the equidistant case. The Zitterbewegung in (2.4) is not present in the decomposition of Theorem 2.4 when $m\geq 2$ , or rather it is of lower order and thus is part of the low-order remainder. As the Bernoulli numbers satisfy $B_2=\tfrac16$ , $B_4=-\tfrac{1}{30}$ , and $B_6 = \tfrac1{42}$ , the extension term (2.5) becomes $\mathbb{V}\text{ar}_E \hat{V}({\kern1.3pt}f) = -t^{2} g'(0^+)\cdot \tfrac16$ , $\mathbb{V}\text{ar}_E \hat{V}({\kern1.3pt}f) = t^{4} g^{(3)}(0^+)\cdot \tfrac{1}{12} \cdot \tfrac{1}{30}$ , and $\mathbb{V}\text{ar}_E \hat{V}({\kern1.3pt}f) = - t^{6} g^{(5)}(0^+) \cdot \tfrac{1}{21} \cdot \tfrac{1}{720}$ for $m=0,1,2$ , respectively. Hence, the extension term of the trapezoidal estimator with perturbed sampling can come arbitrarily close to (2.5) if the errors $E_i$ are sufficiently small. Under the model with cumulative errors, a corresponding statement holds as a consequence of the following result.

Theorem 2.5. Let X be from the model with cumulative errors with intensity $1/t$ , and let the i.i.d. holding times $\omega_i$ (of the unit-intensity process) satisfy ${\mathbb{E}} {\text{e}}^{\eta \omega_1}<\infty$ for some $\eta>0$ . Define $\nu_k$ as the kth moment of $\omega_1$ . Let the measurement function f be exactly (m, 1)-piecewise smooth with covariogram $g=f*\,\check{f}$ . Then, for $t\downarrow 0$ ,

(2.12) \begin{align} \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) & = -t^{2} g'(0^+) \tfrac1{6} \nu_3 + o(t^2) &&\qquad\qquad\!\qquad\textit{for}\quad m=0,\end{align}

(2.13) \begin{align} \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) & = t^{4} g^{(3)}(0^+) \tfrac1{12}\tfrac1{30} (6\nu_5 - 5 \nu_3^2) + o(t^4) && \textit{for}\quad m=1.\end{align}

3. The Peano kernel representation

In this section we consider a locally finite set $X\subset \mathbb{R}$ such that $\text{Conv}(X)=\mathbb{R}$ , and an integrable function $f\,:\,\mathbb{R} \to \mathbb{R}$ with compact support which is known at all points in X. For any $x\in X$ and $j\in\mathbb{Z}$ we define $s_j(x)=s_j(x;\ X)$ as the jth successor (predecessor for $j < 0$ ) of x in X, with $s_0(x)=x$ by definition. Hence, for $j\geq 1$ , $s_j(x)$ and $s_{-j}(x)$ are the unique points in $X\cap(x,\infty)$ and $X\cap({-}\infty,x)$ , respectively, such that $\# (X \cap (x,s_j(x)])=\# (X \cap [s_{-j}(x),x))=j$ . Note that

(3.1) \begin{equation} s_j(x+t;\ X+t)=s_j(x;\ X)+t\end{equation}

for all $t\in \mathbb{R}$ . For all $x\in X$ and $j\in \mathbb{Z}$ we let $h_j(x)=h_j(x;\ X)= s_j(x;\ X)-s_{j-1}(x;\ X)$ be the distance from the jth successor (predecessor) of x to its left neighbour in X. By (3.1),

(3.2) \begin{equation} h_j(x+t;\ X+t) = h_j(x;\ X)\end{equation}

for all $t\in \mathbb{R}$ . We now recall the principle of Newton–Cotes quadrature, adapted to an infinite set of nodes; see [Reference Kiderlen and Dorph-Petersen6] for details. On the interval $[x,s_n(x)]$ , $x\in X$ , the function f is approximated by a polynomial of degree at most $n\in\mathbb{N}$ passing through the points $\{s_j(x),f(s_j(x))\}_{j=0}^{n}$ . The integral of this polynomial on $[x,s_n(x)]$ is

\[ I_x^{(n)}({\kern1.3pt}f)=I_x^{(n)}({\kern1.3pt}f;\ X)=\sum_{j=0}^{n} \beta_j^{(n)}(x) f(s_j(x)),\]

where

(3.3) \begin{equation} \beta_j^{(n)}(x)=\beta_j^{(n)}(x;\ X) = \int_{x}^{s_n(x)} \prod_{\substack{k=0 \\ k\neq j}}^{n}\frac{y-s_k(x)}{s_j(x)-s_k(x)} {\text{d}} y\end{equation}

for $x\in X$ . The approximation

\[ \hat{V}_n({\kern1.3pt}f) =\tfrac1n \sum_{x\in X} I_x^{(n)}({\kern1.3pt}f) = \sum_{x\in X} \alpha(x)f(x)\]

is then an average of the sum of the integral approximations $I_x^{(n)}$ with respect to the starting point chosen. Here

(3.4) \begin{equation} \alpha(x)=\alpha(x;\ X)=\frac{1}{n} \sum_{j=0}^n \beta_j^{(n)}(s_{-j}(x)) .\end{equation}

Remark 3.1. From [Reference Stoer and Bulirsch13, Theorem 2.1.1.1] the integral approximation on an interval $[x,s_n(x)]$ is exact whenever $f=p$ is a polynomial of degree at most n. That is, $R_x^{(n)}(p)=0$ , with the discretization error $R_x^{(n)}$ defined by

(3.5) \begin{equation} R_x^{(n)}({\kern1.3pt}f) = R_x^{(n)}({\kern1.3pt}f;\ X)= I_x^{(n)}({\kern1.3pt}f) - \int_{x}^{s_n(x)}\ f(y) {\text{d}} y ,\end{equation}

$x\in X$ .

As shown in the supplementary material [Reference Stehr and Kiderlen12, Lemma SM1.1], $\beta_j^{(n)}$ is a rational function of point-increments, and (3.2) then implies that

(3.6) \begin{equation} \beta_j^{(n)}(x+t;\ X+t) = \beta_j^{(n)}(x;\ X) \quad\text{and}\quad \alpha(x+t;\ X+t) = \alpha(x;\ X)\end{equation}

for all $t\in\mathbb{R}$ and $x\in X$ .

We are now ready to prove the refined Peano kernel theorem as stated in Theorem 2.2. Given n, X, and $m\in \mathbb{N}_0$ , the mth Peano kernel from Theorem 2.2 is defined as

(3.7) \begin{equation} K_m(r)= K_m(r;\ X)=\frac{1}{m!n} \sum_{x\in X}\boldsymbol{1}_{{(x,s_n(x)]}}(r) R_x^{(n)} \left( (\!\cdot-r)_+^m\right) .\end{equation}

The mapping $x \mapsto (x-r)_+^m$ should be understood as

\[ (x-r)_+^m = \begin{cases} (x-r)^m & \text{for} \quad x>r, \\[5pt] 0 & \text{for} \quad x\leq r . \end{cases}\]

Hence, $K_m$ is piecewise polynomial of degree at most $m+1$ with coefficients determined by X.

Proof of Theorem 2.2. Fix $n\in\mathbb{N}$ and note that $nR^{(n)} ({\kern1.3pt}f)=\sum_{x\in X}R^{(n)}_x ({\kern1.3pt}f)$ by (2.3) and (3.5). For all $x\in X$ and $y\in[x,s_n(x)]$ , an induction argument using the refined partial integration formula [Reference Kiêu7, Lemma 4.1] yields

\begin{equation*}\begin{aligned} f(y^-) = {} & \sum_{k=0}^m \frac{f^{(k)}(x^+)}{k!}(y-x)^k \\ &+ \frac{1}{m!} [r]{\sum_{a\in {D}_{f^{(m)}} \cap(x,y)}}\quad {J}_{f^{(m)}} (a) (y-a)^m + \frac{1}{m!}\int_x^y f^{(m+1)}(t)(y-t)^m {\text{d}} t\end{aligned}\end{equation*}

for all (m,1)-piecewise smooth functions f, with $m\in\mathbb{N}_0$ . We now assume $m\leq n$ . Using the linearity of $R^{(n)}_{x}$ , the fact that all polynomials of order at most n are integrated exactly, and the fact that $R_x^{(n)}$ commutes with integration, we find that (with all expressions considered as functions of y)

\begin{align*} &m! R_x^{(n)}({\kern1.3pt}f) \\ &\quad = R_x^{(n)}\Bigg({ \sum_{a\in {D}_{f^{(m)}} \cap(x,y)}}\quad {J}_{f^{(m)}}(a) (y-a)^m + \int_x^y f^{(m+1)}(t)(y-t)^m {\text{d}} t \Bigg) \\ &\quad = {\sum_{a\in {D}_{f^{(m)}} \cap (x,s_n(x)]}}\quad {J}_{f^{(m)}}(a) R_x^{(n)}\bigl( (y-a)_+^m \bigr) + \int_x^{s_n(x)} f^{(m+1)}(t)R_x^{(n)}\bigl((y-t)_+^m \bigr) {\text{d}} t .\end{align*}

Changing the summation order, (3.7) implies that

\begin{equation*} R^{(n)}({\kern1.3pt}f) = \frac1n \sum_{x\in X} R_x^{(n)}({\kern1.3pt}f) = [r]{\sum_{a\in {D}_{f^{(m)}} }} {J}_{f^{(m)}}(a) K_m(a) + \int_\mathbb{R} f^{(m+1)}(t) K_m(t) {\text{d}} t,\end{equation*}

as claimed.

Before proceeding, we state a useful lemma on continuity properties of the Peano kernel. For $r\in\mathbb{R}$ we have

\[ K_m(r) = \frac{1}{m!n}\sum_{x\in X} \boldsymbol{1}_{{(x,s_1(x)]}}(r) \sum_{i=1-n}^0 R_{s_i(x)}^{(n)}\left( (\!\cdot-r)_+^m\right) .\]

The following result is a simple consequence of this representation and the fact that polynomials of degree at most n are approximated exactly.

Lemma 3.1. Fix $n\in\mathbb{N}$ and a locally finite point-set X with $\text{Conv}(X)=\mathbb{R}$ . Then, for all $x\in X$ and $m\in \mathbb{N}$ , the function $K_m$ is differentiable on $(x,s_1(x))$ with derivative $-K_{m-1}$ and jump

\begin{equation*} {J}_{K_m}(x) =\frac{1}{m!n} R_x^{(n)} \bigl( (\!\cdot-x)^m\bigr).\end{equation*}

In particular, $K_m$ is $(m-1)$ times continuously differentiable for all $1\leq m \leq n$ .

4. Integrability properties

To argue that $\hat{V}_n({\kern1.3pt}f)$ is an unbiased estimator for $\int f(x) {\text{d}} x$ when applied to randomized sampling points, we recall the notion of the Palm distribution of a stationary point process $X\subset \mathbb{R}$ . It can be interpreted as the conditional distribution of X given that $0\in X$ . We denote it by $\mathbb{P}^0$ with the corresponding expectation denoted by ${\mathbb{E}}^0$ . When considering the point process X under its Palm distribution, we will often suppress the dependence on the point $0\in X$ in various expressions; i.e. under $\mathbb{P}^0$ we write, for instance,

\[ s_i=s_i(0), \qquad h_i=h_i(0), \qquad \beta_j^{(n)} = \beta_j^{(n)}(0)\]

for all $i\in\mathbb{Z}$ and $j=0,\ldots,n$ . In addition we write $\mathbf{h}=(h_1,\ldots,h_n)$ and, for $i\in\mathbb{Z}$ , $\mathbf{h}(s_i)=(h_1(s_i),\ldots,h_n(s_i))=(h_{i+1},\ldots,h_{i+n})$ under $\mathbb{P}^0$ . As mentioned in Section 1, a weaker assumption than (2.2) is sufficient to ensure the unbiasedness of the estimator.

Assumption 4.1 For a given $n\in \mathbb{N}$ we assume that

(4.1) \begin{equation} {\mathbb{E}}^0 \left[ \frac{ \mathbf{h}^{\mathbf{m}} }{ \mathbf{h}^{\mathbf{m}'} } \right] <\infty\end{equation}

for all multi-indices $\mathbf{m},\mathbf{m}'\in\mathbb{N}_0^n$ with $|\mathbf{m}|\in \{n+1,n+2\}$ and $|\mathbf{m}'|=n$ , where $|\mathbf{m}|= |(m_1,\ldots,m_n)| =\sum_{k=1}^n m_k$ .

Using Hölder’s inequality and [Reference Kiderlen and Dorph-Petersen6, Eq. (13)], one shows that Assumption 2.1 is stronger than Assumption 4.1. In [Reference Stehr and Kiderlen12, Lemma SM1.1] of the supplementary material it is shown that the weight $\beta_j^{(n)}(x)$ is a rational function of the point-increments $(h_1(x),\ldots,h_n(x))$ , $x\in X$ , where the numerator is a homogeneous polynomial of degree $n+1$ , and the denominator is a nonvanishing homogeneous polynomial of degree n with nonnegative coefficients. From the fact that the Palm distribution is invariant under bijective point shifts [Reference Kiderlen and Dorph-Petersen6, Eq. (13)], it is easily seen that ${\mathbb{E}}^0|{\beta_j^{(n)}(s_{-j})}| < \infty$ for all $j\in\{0,\ldots,n\}$ when Assumption 4.1 is satisfied, and consequently

(4.2) \begin{equation} {\mathbb{E}}^0 |\alpha(0)| < \infty;\end{equation}

see [Reference Stehr and Kiderlen12, Lemma SM1.2] of the supplementary material. We conclude that either of the two assumptions is sufficient to guarantee the Palm-integrability of $\alpha(0)$ , which will be used in the proof of Theorem 2.1.

To argue for the variance bounds presented in Theorem 2.3 we need higher-order moment and translation-invariance properties of the Peano kernel $K_m$ defined in (3.7).

Lemma 4.1. Let $n\in\mathbb{N}$ be given and assume that X is a stationary point process. Then, for all $m\in\mathbb{N}_0$ , $K_m$ is a stationary stochastic process. If Assumption 2.1 holds, $K_m(0)$ has finite absolute moments of all (positive) orders. Moreover, if X has intensity $\gamma$ , $K_m$ satisfies

(4.3) \begin{equation} {\mathbb{E}} K_m(0) = \gamma {\mathbb{E}}^0 J_{K_{m+1}}(0) =\frac{\gamma}{(m+1)!n} {\mathbb{E}}^0 R_0^{(n)} \big( ({\cdot})^{m+1}\big)\end{equation}

for all $m\in \mathbb{N}_0$ . In particular, ${\mathbb{E}} K_m(0)=0$ for all $m<n$ .

Proof. Fix $n\in\mathbb{N}$ . For any $r,s\in\mathbb{R}$ and any locally finite point-set X, the Peano kernel satisfies

(4.4) \begin{equation} K_m(r+s;\ X) = K_m(r;\ X-s) .\end{equation}

This follows from the definition of $K_m$ and

\[ R_{x}^{(n)}\left((\!\cdot-(r+s))_+^m;\ X\right) = R_{x-s}^{(n)}\left((\!\cdot-r)_+^m;\ X-s\right) , \quad \quad x\in X,\]

which in turn is a consequence of (3.1) and (3.6). By (4.4), the stationarity of $K_m$ is inherited from the stationarity of the point process X.

We now prove that $K_m(0)$ has finite absolute moments. Let $k\in \mathbb{N}$ be given. For arbitrary $r\in \mathbb{R}$ put $I_r=\{x\in X\,:\, r\in (x,s_n(x)]\}$ . Using Hölder’s inequality and some rather crude upper bounds we obtain from (3.5) and (3.7) that

\[ \big\lvert K_m^k(0)\big\rvert \leq \sum_{x\in I_0} \big\lvert R_{x}^{(n)}\left(({\cdot})_+^m\right)\big\rvert^k \leq \sum_{x\in I_0} (s_n(x))^{km} \bigg( \sum_{j=0}^n \lvert\beta_j^{(n)}(x)\rvert + s_n(x)\bigg)^k .\]

By the refined Campbell theorem [Reference Schneider and Weil10, Theorem 3.5.3], (3.6), and (3.1), it follows that

\[ {\mathbb{E}} \big\lvert K_m^k(0)\big\rvert \leq \gamma {\mathbb{E}}^0 \int_0^{s_n} x^{km} \bigg( \sum_{j=0}^n \lvert\beta_j^{(n)}\rvert + x\bigg)^k {\text{d}} x \leq \gamma {\mathbb{E}}^0 s_n^{km+1} \bigg( \sum_{j=0}^n \lvert\beta_j^{(n)}\rvert + s_n\bigg)^k ,\]

where $\gamma$ is the intensity of X. By the supplementary material [Reference Stehr and Kiderlen12, Lemma SM1.1], Assumption 2.1, and the fact that $s_n=\sum_{j=1}^n h_j$ under $\mathbb{P}^0$ , the variables $s_n$ and $\beta_j^{(n)}$ have finite absolute moments of all orders under $\mathbb{P}^0$ . This implies that ${\mathbb{E}} \lvert K_m^k(0)\rvert<\infty$ .

Equation (4.3) is a simple consequence of the refined Campbell theorem [Reference Schneider and Weil10, Theorem 3.5.3], Lemma 3.1, and [Reference Kiderlen and Dorph-Petersen6, Eq. (13)].

5. Unbiasedness of Newton–Cotes estimators

Proof of Theorem 2.1. Fix $n\in\mathbb{N}$ and let $X\subseteq \mathbb{R}$ be a stationary point process with finite and positive intensity $\gamma$ . As $\alpha$ satisfies (3.6) and $\alpha(0)$ is Palm-integrable by (4.2), [Reference Kiderlen and Dorph-Petersen6, Theorem 1] can be applied. It states that

(5.1) \begin{equation} {\mathbb{E}} \hat{V}_n({\kern1.3pt}f) = \gamma {\mathbb{E}}^0[\alpha(0)] \int_\mathbb{R} f(x) {\text{d}} x\end{equation}

holds for all integrable functions $f\,:\,\mathbb{R}\to\mathbb{R}$ with compact support. Hence, if we can show that ${\mathbb{E}}^0[\alpha(0)]=\gamma^{-1}$ , we have shown that $\hat{V}_n({\kern1.3pt}f)$ is unbiased.

For $s\in \mathbb{R}$ we reuse the notation $I_s$ from the end of the previous section. When f is an integrable function and $|f|\leq 1$ , (3.5) implies

\[ \sum_{x\in I_s} \lvert R_x^{(n)}({\kern1.3pt}f)\rvert \leq \sum_{x\in I_s} \bigg(\sum_{j=0}^{n} \lvert\beta_{j}^{(n)}(x)\rvert +(s_{n}(x)-x)\bigg).\]

The refined Campbell theorem [Reference Schneider and Weil10, Theorem 3.5.3], (3.6), and (3.1) imply

\[ {\mathbb{E}} \sum_{x\in I_s} \lvert R_x^{(n)}({\kern1.3pt}f)\rvert \leq \gamma {\mathbb{E}}^0 \int_{s-s_n}^{s} \bigg( \sum_{j=0}^n \lvert\beta_{j}^{(n)}\rvert + s_n \bigg) {\text{d}} x =\gamma {\mathbb{E}}^0 \bigg[s_n\sum_{j=0}^n \lvert\beta_j^{(n)}\rvert + s_n^2\bigg]<\infty,\]

where the finiteness follows from the supplementary material [Reference Stehr and Kiderlen12, Lemma SM1.1] and Assumption 4.1, which is weaker than Assumption 2.1. Note that the finite upper bound does not depend on s.

Now let $r>0$ be given and consider the function $f_r=\boldsymbol{1}_{{[0,r]}}$ . Recall that

\[ R^{(n)}({\kern1.3pt}f_r)=\hat{V}_n({\kern1.3pt}f_r)-\int_\mathbb{R} f_r(x) {\text{d}} x =\frac1 n \sum_{x\in X} R_x^{(n)}({\kern1.3pt}f_r)\]

is the error of the nth Newton–Cotes estimator. The Newton–Cotes approximation on an interval $[x,s_n(x)]$ is exact for all polynomials of degree at most n, and in particular, it is exact for constant functions. Hence, $R_x^{(n)}({\kern1.3pt}f_r)=0$ whenever $[x,s_n(x)]\cap \{0,r\}=\emptyset$ . This implies

\[ |{\mathbb{E}} R^{(n)}({\kern1.3pt}f_r)|\leq {\mathbb{E}}\sum_{x\in I_0} |R_x^{(n)}({\kern1.3pt}f)|+{\mathbb{E}}\sum_{x\in I_r} |R_x^{(n)}({\kern1.3pt}f)|\leq 2C,\]

for some finite $C\in\mathbb{R}$ which is independent of r. Equation (5.1) now implies

\[ 0=\lim_{r\to\infty}\tfrac1r {\mathbb{E}} R^{(n)}({\kern1.3pt}f_r)= \gamma{\mathbb{E}}^0\alpha(0)-1,\]

so ${\mathbb{E}}^0\alpha(0)=1/\gamma$ as asserted.

Remark 5.1. If $X=t \bigl(U+\mathbb{Z}\bigr)$ is the equidistant point process, $\alpha(x)=t$ for all $x\in X$ . In fact, since the Palm version of X is the deterministic set $t\mathbb{Z}$ , (3.6) yields

\[ \alpha(x;\ X) = \alpha(x-x;\ X-x) = \alpha(0;\ t\mathbb{Z})\]

for all $x\in X$ . Hence $\alpha = \alpha(x)$ is deterministic, and $\hat{V}_n({\kern1.3pt}f) = \alpha \sum_{x\in X} f(x)$ . Assumption 2.1 is trivially satisfied, so Theorem 2.1 implies the well-known fact that $\hat{V}_n({\kern1.3pt}f)$ is unbiased for $\int f {\text{d}} x$ . This is equivalent to $\alpha=t$ .

6. Asymptotic variance behaviour of Newton–Cotes estimators

In this and the following section we derive variance expressions showing the exact dependence on the mean point distance $t>0$ . To this end, we will consider the Peano kernel and its associated covariance function applied to the unit-intensity scaling of the process X. That is, let the point process of interest, X, have intensity $1/t$ , and define its unit-intensity scaled process by $X^u = X/t$ . Hence the Peano kernel $K_m({\cdot};\ X^u)$ with respect to $X^u$ remains unchanged when X is rescaled. To avoid intricate notation, we will write ${\mathbb{E}}_u$ , $\text{Cov}_u$ , $\mathbb{V}\text{ar}_u$ when expectation, covariance, and variance are understood with respect to $X^u$ ; for instance ${\mathbb{E}}_u K_m(r)= {\mathbb{E}} K_m(r;\ X^u)$ . Similarly, the Palm distribution of $X^u$ and its expectation are denoted $\mathbb{P}^0_u$ and ${\mathbb{E}}^0_u$ , respectively. Lastly, the covariance function of $K_m$ applied to $X^u$ is denoted by $ H_m^u(s)=\text{Cov}_u\bigl(K_m(s + r),K_m(r) \bigr)$ for all $s,r\in\mathbb{R}$ .

Before proving Theorem 2.3, we recall the variance decomposition of the estimator (1.2) in the equidistant case, as it shows a great resemblance to the new non-equidistant set-up. First we introduce the periodic Bernoulli functions $P_m$ , which we define as in [Reference Knopp8, Paragraph 297]. Let $(\tilde P_m)_{m=0}^\infty$ be the sequence of rescaled Bernoulli polynomials defined inductively by $\tilde{P}_0(x)=1$ , $\tilde{P}_1(x)=x-\frac{1}{2}$ , and $\tilde{P}_{\!m+1}^{\prime}=\tilde{P}_{m}$ , with

\[\tilde{P}_{m+1}(0)=\tilde{P}_{m+1}(1)=\frac{1}{(m+1)!}B_{m+1},\]

for $m\in\mathbb{N}$ , where $B_m$ is the mth Bernoulli number. This normalization is chosen as in [Reference Kiêu7] for ease of comparison with the results there. Then $P_m(x)=\tilde{P}_m(x-\lfloor x \rfloor)$ is the mth Bernoulli polynomial, evaluated at the fractional part of $x\in\mathbb{R}$ . Note that $P_m$ is continuous for all $m\neq 1$ . When the measurement function f is (m,1)-piecewise smooth, the variance decomposes as

(6.1) \begin{equation} \mathbb{V}\text{ar} (\hat{V}({\kern1.3pt}f)) = - t^{2m+2} \quad{ \sum_{a\in {D}_{g^{(2m+1)}} } } {J}_{g^{(2m+1)}}(a) P_{2m+2}\big(\tfrac{a}{t}\big) + o\big(t^{2m+2}\big)\end{equation}

as $t\downarrow 0$ (see [Reference Kiêu7, Chapter 5]). Here, $g=f*\check f$ is the covariogram of f, and the term $o(t^{2m+2})$ can explicitly be given as $- t^{2m+2}\int_\mathbb{R} g^{(2m+2)}(s) P_{2m+2}\big(\tfrac{s}{t}\big) {\text{d}} s$ . When the point process X is not equidistant, we find a similar variance representation involving the Peano kernels instead of the periodic Bernoulli functions.

Proposition 6.1. Let $n\in\mathbb{N}$ be given, and assume that X is a stationary point process with intensity $1/t$ such that Assumption 2.1 holds. If f is (m, 1)-piecewise smooth and $k = \min\{m,n\}$ , then

(6.2) \begin{equation}\begin{aligned} ({-}1)^{k+1} \mathbb{V}\text{ar} (\hat{V}_{n}({\kern1.3pt}f)) ={} & t^{2k+2} \quad{\sum_{a\in {D}_{g^{(2k+1)}} } } {J}_{g^{(2k+1)}}(a) H_{k}^u\big(\tfrac{a}{t}\big) \\ &+ t^{2k+2} \int_\mathbb{R} g^{(2k+2)}(s) H_{k}^u\big(\tfrac{s}{t}\big) {\text{d}} s .\end{aligned}\end{equation}

If $k=m<n$ or X is weakly n-admissible, the variance behaviour is determined by the first term, as

\[ \int_\mathbb{R} g^{(2k+2)}(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s = o(1)\]

for $t\downarrow 0$ .

Proof. The definition of $\alpha(x)$ and elementary calculations give

\[ \alpha(x;\ X)=t\alpha(x/t;\ X^u)\]

for $x\in X$ , so putting $f_t(x)=f(tx)$ we see that

(6.3) \begin{equation} \hat{V}_{n}({\kern1.3pt}f)=t\hat{V}_{n}({\kern1.3pt}f_t;\ X^u),\end{equation}

where the latter estimator is given in terms of the unit-intensity process $X^u$ . As $k\leq n$ , Theorem 2.2 implies

\begin{equation*} R^{(n)} ({\kern1.3pt}f_t;\ X^u) = \int_\mathbb{R} f_t^{(k+1)}(s) K_k(s;\ X^u) {\text{d}} s + {\sum_{a\in {D}_{f_t^{(k)}} } } {J}_{f_t^{(k)}}(a) K_k(a;\ X^u).\end{equation*}

Using $\,f_t'(x)=t f'(tx)$ whenever the derivative is defined, we arrive at

\begin{equation*} R^{(n)} ({\kern1.3pt}f_t;\ X^u) = t^{k}\int_\mathbb{R} f^{(k+1)}(s) K_k\big(\tfrac{s}{t};\ X^u\big) {\text{d}} s +t^k {\sum_{a\in {D}_{f^{(k)}} } } {J}_{f^{(k)}}(a) K_k\big(\tfrac{a}{t};\ X^u\big).\end{equation*}

Hence, using (6.3) and the unbiasedness of $\hat{V}_{n}$ , we get

(6.4) \begin{align}\nonumber \mathbb{V}\text{ar} (\hat{V}_{n} ({\kern1.3pt}f))& = t^2\mathbb{V}\text{ar}_u (\hat{V}_{n} ({\kern1.3pt}f_t)) = t^2{\mathbb{E}}_u (R^{(n)} ({\kern1.3pt}f_t))^2 \\ & = t^{2k+2}{\mathbb{E}}_u\Big(\!\int_\mathbb{R} f^{(k+1)}(s) K_k\big(\tfrac{s}{t}\big) {\text{d}} s + {\sum_{a\in {D}_{f^{(k)}} } } {J}_{f^{(k)}}(a) K_k\big(\tfrac{a}{t}\big) \Big)^2.\end{align}

An application of [Reference Kiêu7, Proposition 5.7] yields

(6.5) \begin{equation}\begin{aligned} f^{(k+1)}*\,\check{f^{(k+1)}}(x) = ({-}1)^{k+1}g^{(2k+2)}(x) \\ & - {\sum_{a\in {D}_{f^{(k)}} } } J_{f^{(k)}} (a) f^{(k+1)}(a-x) - {\sum_{a\in D_{f^{(k)}} }} J_{f^{(k)}} (a) f^{(k+1)}(a+x),\end{aligned}\end{equation}

and furthermore the jumps of $g^{(2k+1)}$ are given by

\[ J_{\!g^{(2k+1)}}(a) = {\sum_{b\in {D_{f^{(k)}}}}} J_{\!f^{(k)}} (b) J_{\!f^{(k)}} (b-a);\]

see [Reference Kiêu7, Eq. (5.12)]. The stationarity and square-integrability of $K_k$ from Lemma 4.1 implies that ${\mathbb{E}}_u K_k(r)$ and $H_k^u(s) = \text{Cov}_u(K_k(r),K_k(s+r))$ are both finite and independent of ${r\in\mathbb{R}}$ . Equation (6.2) now follows by expanding (6.4), applying (6.5), and using the structure of $ J_{g^{(2k+1)}}$ together with Fubini’s theorem. The latter may be applied thanks to the square-integrability of $K_k$ and the fact that $f^{(k+1)}$ is bounded with compact support.

We now show that

\[\lim_{t\downarrow 0}\int_\mathbb{R} g^{(2k+2)}(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s = 0\]

if $k=m<n$ or X is weakly n-admissible. The weak admissibility assumption yields

(6.6) \begin{equation} \lim_{t\downarrow 0} \int_0^1 H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s = 0\end{equation}

for $k=n$ . Equation (6.6) also holds for $k=m<n$ without additional assumptions. In fact, for $k<n$ we have $K^{\prime}_{k+1}=-K_k$ by Lemma 3.1, and thus, using Fubini’s theorem,

\[ \Bigl| \int_0^t H_k^u(s) {\text{d}} s \Bigr| = \left| \text{Cov}_u\big(K_k(0), K_{k+1}(0)\big) - \text{Cov}_u\big(K_k(0), K_{k+1}(t) \big) \right| \leq c<\infty ,\]

where Hölder’s inequality and the stationarity of the Peano kernels have been used to show that the constant c is independent of t. A substitution allows us to derive (6.6) from this.

Now fix $k\leq n$ and let $\epsilon > 0$ be given. As $g^{(2k+2)}$ is integrable and bounded, there is a simple function $\phi$ such that $\phi \leq g^{(2k+2)}$ and

\[ 0 \leq \int_\mathbb{R} g^{(2k+2)}(s) {\text{d}} s - \int_\mathbb{R} \phi(s) {\text{d}} s < \frac{\epsilon}{2C},\]

where the finite constant $C>0$ satisfies $\sup_{s\in\mathbb{R}}|H_k^u(s)| \leq C$ . This implies that

\[ \left|\int_\mathbb{R} g^{(2k+2)}(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s - \int_\mathbb{R} \phi(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s \right| < \frac{\epsilon}{2} .\]

As $\phi$ is simple, (6.6) implies that

\[\lim_{t\downarrow 0}\int_\mathbb{R} \phi(s)H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s=0.\]

We conclude that

\[\bigg|\int_\mathbb{R} \phi(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s\bigg| < \frac{\epsilon}{2}\]

for sufficiently small $t>0$ , and hence

\[ \left|\int_\mathbb{R} g^{(2k+2)}(s) H_k^u\big(\tfrac{s}{t}\big) {\text{d}} s \right|<\epsilon\]

for such small $t>0$ .

Proof of Theorem 2.3. Recall that $k=\min\{m,n\}$ . By Lemma 4.1 there exists $C<\infty$ such that $\sup_{s\in\mathbb{R}}|H_k^u(s)|\leq C$ , and we immediately see from (6.2) that

\begin{equation*} \mathbb{V}\text{ar} (\hat{V}_{n} ({\kern1.3pt}f)) \leq t^{2k+2} \Bigg( C \| g^{(2k+2)}\|_{\infty} \lambda(\text{supp}\,g) + C {\sum_{a\in D_{g^{(2k+1)}} }} \lvert J_{g^{(2k+1)}}(a) \rvert \Bigg),\end{equation*}

where $\lambda(\text{supp}\,g)<\infty$ is the Lebesgue measure of the support of g. As g is $(2k+1,1)$ -piecewise smooth by [Reference Kiêu7, Corollary 5.8], the t-independent constant is finite, and (2.6) therefore follows.

For the stronger result (2.7), note that $m>n$ , and hence g is $(2n+3,1)$ -piecewise smooth; in particular, $g^{(2n+2)}$ is continuous. An application of Proposition 6.1 to the (n,1)-piecewise smooth function f and a substitution gives

(6.7) \begin{equation}({-}1)^{n+1} \mathbb{V}\text{ar} (\hat{V}_{n}({\kern1.3pt}f)) = t^{2n+3} \int_\mathbb{R} g^{(2n+2)}(st) H_n^u(s) {\text{d}} s .\end{equation}

Let $b>0$ satisfy $\text{supp}\,g \subset [{-}b,b]$ . As $g^{(2n+3)}$ is bounded and measurable, $g^{(2n+2)}$ is absolutely continuous. As $H_n^u$ is bounded and hence integrable on $[{-}b/t,b/t]$ for any $t>0$ , the function V given by $V(s)= \int_{-b/t}^s H_n^u(y) {\text{d}} y$ is also absolutely continuous on $[{-}b/t,b/t]$ with derivative $H_n^u$ almost everywhere; see e.g. [Reference De Barra4, Section 9.3] for details on absolutely continuous functions. Furthermore, as X and therefore $X^u$ are assumed strongly n-admissible, V is bounded by a t-independent constant Cʹ, say. Partial integration for absolutely continuous functions [Reference De Barra4, Theorem 9] shows that

\[ \int_\mathbb{R} g^{(2n+2)}(st) H_n^u(s) {\text{d}} s = \int_{-b/t}^{b/t} g^{(2n+2)}(st) H_n^u(s) {\text{d}} s = -t \int_{-b/t}^{b/t} g^{(2n+3)}(st) V(s) {\text{d}} s ,\]

where we used that $g^{(2n+2)}$ vanishes at $\pm b$ . Returning to (6.7) we find

\[ \mathbb{V}\text{ar} (\hat{V}_{n}({\kern1.3pt}f)) \leq t^{2n+3} t \int_{-b/t}^{b/t} |g^{(2n+3)}(st)| |V(s)| {\text{d}} s \leq t^{2n+3} 2 b \|g^{(2n+3)}\|_\infty C' . \]

This proves the assertion.

Remark 6.1. If f is exactly (m,1)-piecewise smooth with $m\leq n$ , and (6.6) is satisfied with $k=m$ , the variance

\begin{equation*} \mathbb{V}\text{ar} (\hat{V}_{n}({\kern1.3pt}f)) = ({-}1)^{m+1} t^{2m+2} \enspace{\sum_{a\in D_{g^{(2m+1)}} }} J_{g^{(2m+1)}}(a) H_m^u\big(\tfrac{a}{t}\big) + o\big(t^{2m+2}\big)\end{equation*}

is exactly of order $t^{2m+2}$ . This is easily seen by assuming that

\begin{equation*} {\sum_{a\in D_{g^{(2m+1)}} }} J_{g^{(2m+1)}}(a) H_m^u\big(\tfrac{a}{t}\big) \to 0\end{equation*}

as $t \downarrow 0$ , and using that $g^{(2m+1)}$ has a jump at 0. Applying (6.6) yields a contradiction. In particular, the decrease rate in (2.6) is optimal if $m<n$ or X is weakly n-admissible.

7. Variance behaviour under perturbed and cumulative sampling

In this section the general findings will be exemplified and made more explicit for the perturbed model and the model with cumulative errors introduced in Examples 2.1 and 2.2, respectively.

7.1. Perturbed sampling

To construct the perturbed model, we let U be uniform on (0,1) and independent of the sequence of i.i.d. variables $\{E_i\}_{i\in\mathbb{Z}}$ , where $|E_i|<\frac12$ almost surely and ${\mathbb{E}} E_1 = 0$ . The perturbed model is the stationary point process $X=\{x_i\}_{i\in\mathbb{Z}}$ for which $x_i=t\bigl(U+i+E_i\bigr)$ , for all $i\in\mathbb{Z}$ . Under its Palm distribution, we have

(7.1) \begin{equation} h_k = t\bigl(1+ E_k - E_{k-1}\bigr) \leq 2t ,\end{equation}

$k\in\mathbb{Z}$ , so (2.2) is equivalent to

(7.2) \begin{equation} {\mathbb{E}} (1+E_1-E_0)^{-j}<\infty\end{equation}

for all $j\in\mathbb{N}$ . For instance, (7.2) holds if there is $\epsilon > 0$ such that $|E_0|\leq \tfrac12 - \epsilon$ almost surely. For the scaled perturbed model $X^u$ we define the shifted kernel $K_m^*$ by

\[ K_m^*(r)=K_m(r+U;\ X^u)=K_m(r;\ X^u-U)\]

for $m\in\mathbb{N}_0$ . Note that it only depends on the perturbations $\{E_i\}$ and not on the initial uniform translation, and thus it is not (necessarily) a stationary process. However, by the i.i.d. structure of $\{E_i\}$ and the fact that $\beta_j^{(n)}$ is a rational function of point-increments, we see that

(7.3) \begin{equation} K_{m}^{*}(r) \overset{\mathcal{D}}{=} K_{m}^{*}(r+k)\end{equation}

for all $k\in\mathbb{Z}$ . This can be used to show that $X^u$ (and equivalently X) is strongly admissible.

Lemma 7.1. Let $n\in\mathbb{N}$ be given and assume that X is a stationary point process from the perturbed model such that (7.2) holds. Then, for all $m\in\mathbb{N}_0$ and $r\geq 2n+2$ , $H_m^u(r)=H_m^u(r+1)$ and $\int_r^{r+1} H_m^u(s) {\text{d}} s =0$ . In particular, X is strongly n-admissible.

Proof. Fix $n\in\mathbb{N}$ and $r\geq 2n+2$ . For such r, the independence between U and $\{E_i\}$ yields

(7.4) \begin{equation} {\mathbb{E}}_u[K_m(0)K_m(r)] = \int_0^1 {\mathbb{E}}[K_m^*({-}u)K_m^*(r-u)] {\text{d}} u = \int_0^1 {\mathbb{E}}[K_m^*({-}u)]{\mathbb{E}}[K_m^*(r-u)] {\text{d}} u .\end{equation}

Equation (7.3) and (7.4) imply that ${\mathbb{E}}_u[K_m(0)K_m(r)]={\mathbb{E}}_u[K_m(0)K_m(r+1)]$ , and by stationarity of $K_m$ we conclude that $H_m^u(r)=H_m^u(r+1)$ .

Returning to (7.4), we find by Fubini’s theorem, a substitution, and the stationarity of $K_m$ that

\begin{align*} \int_r^{r+1} {\mathbb{E}}_u[K_m(0)K_m(s)] {\text{d}} s \\ & = \int_0^1 {\mathbb{E}}[K_m^*({-}u)] \int_0^1 {\mathbb{E}}[K_m^*(r+1-u-s)] {\text{d}} s {\text{d}} u \\ & = \int_0^1 {\mathbb{E}}[K_m^*({-}u)] {\mathbb{E}}_u[K_m(r+1-u)] {\text{d}} u = ({\mathbb{E}}_u K_m(0))^2, \end{align*}

which yields the asserted properties of $H_m^u$ . This clearly implies that $X^u$ and equivalently X are strongly n-admissible.

In order to obtain explicit leading terms in Theorem 2.4 we state in the following a connection between the covariance function $H_m^u$ and certain periodic Bernoulli functions. For our purpose, it is enough to consider $m\in\{0,1\}$ , but we also state that the result holds for all m, when no perturbations are present.

Lemma 7.2. Let $n=1$ and let X be from the perturbed model. Then

(7.5) \begin{equation} H_m^u(r) = ({-}1)^m {\mathbb{E}} [P_{2m+2}(r + E_1 - E_0)]\end{equation}

for $m=0,1$ and all $|r| \geq 4$ .

If $X = t \bigl(U + \mathbb{Z}\bigr)$ (hence $X^u=U+\mathbb{Z}$ ) is the equidistant model, then

(7.6) \begin{equation} H_m^u(r) = ({-}1)^m P_{2m+2} (r)\end{equation}

for all $n\in\mathbb{N}$ , $m\leq n$ , and $r\in\mathbb{R}$ .

The proof of Lemma 7.2 can be found in the supplementary material of this paper; see [Reference Stehr and Kiderlen12, Corollaries SM2.1 and SM2.2]. As a consequence of (7.6), for the equidistant model $X=t\bigl(U+\mathbb{Z}\bigr)$ , the variance representation (6.2) found using the Peano kernels coincides with the classical variance representation (6.1) found using Euler–McLaurin formulae.

Before turning to the proof of Theorem 2.4, we emphasize that the integrability condition (2.2), or equivalently the condition (7.2), was omitted in the statement of the theorem as we work with the trapezoidal rule. In fact, the unbiasedness of $\hat{V}_{1}({\kern1.3pt}f)$ for all stationary point processes and integrable, compactly supported functions f was already noted in the paragraph following the statement of Theorem 2.1. By (7.1), the weights satisfy

\[ \beta_0^{(1)}(x) =\beta_1^{(1)}(x) = \tfrac12 h_1(x) \leq t ,\]

$x\in X$ , which replaces the condition (2.2) in all the arguments in Sections 5 and 6. The assumptions of Proposition 6.1 are thus satisfied.

Proof of Theorem 2.4. Let $m\in\{0,1\}$ . The $(2m+1)$ th derivative of the covariogram g is an odd function, implying $ J_{g^{(2m+1)}}(0) = 2g^{(2m+1)}(0^+)$ . As X is strongly admissible, Proposition 6.1 in combination with (7.5) yields the variance decomposition

\[ \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) = ({-}1)^{m+1} t^{2m+2} 2 g^{(2m+1)}(0^+) H_m^u(0) +Z_m(t) + o(t^{2m+2}) ,\]

where the Zitterbewegung $Z_m(t)$ is given by

(7.7) \begin{equation} Z_m(t) = -t^{2m+2} \quad{\sum_{a\in D_{g^{(2m+1)}} \setminus\{0\}}}\quad J_{g^{(2m+1)}}(a) {\mathbb{E}} [P_{2m+2}\big(\tfrac{a}{t} + E_1 - E_0\big)] .\end{equation}

The facts that $Z_m$ is a finite sum of terms each oscillating around 0 and that it is of order $t^{2m+2}$ follow from arguments similar to those of [Reference Kiêu7, Section 5.2] as f is assumed to be exactly (m,1)-piecewise smooth. By the refined Campbell theorem [Reference Schneider and Weil10, Theorem 3.5.3] and the facts that ${\mathbb{E}}_u K_0(0)=0$ and ${\mathbb{E}}_u K_1(0) = \tfrac12{\mathbb{E}}^0_u [R_0^{(1)}(({\cdot})^2)]$ by (4.3), we find that $H_m^u(0) = \mathbb{V}\text{ar}_u(K_m(0))$ satisfies

(7.8) \begin{align} H_0^u(0) & = {\mathbb{E}}^0_u \int_0^{h_1} \big(\tfrac12 h_1 - y\big)^2 {\text{d}} y = \tfrac{1}{12} {\mathbb{E}}^0_u h_1^3\ , \end{align}

(7.9) \begin{align} H_1^u(0) & = {\mathbb{E}}^0_u \int_0^{h_1} \big(\tfrac12 h_1 y - \tfrac12 y^2\big)^2 {\text{d}} y - \big(\tfrac{1}{12} {\mathbb{E}}^0_u h_1^3\big)^2 = \tfrac{1}{120} {\mathbb{E}}^0_u h_1^5 - \tfrac{1}{144} \big({\mathbb{E}}^0_u h_1^3\big)^2 .\end{align}

Using (7.1), it is elementary to conclude (2.8) and (2.9).

Now let $m\geq 2$ be given and define $\tilde H_1^u$ by $\tilde H_1^u(s) = H_1^u(s) + {\mathbb{E}}[P_4(s+E_1 - E_0)]$ . By Lemma 7.2, $\tilde H_1^u(s)$ vanishes for $|r|>4$ . Since $g^{(4)}$ is continuous, an application of Proposition 6.1 to the (1,1)-piecewise smooth function f, Fubini’s theorem, and the refined partial integration formula [Reference Kiêu7, Lemma 4.1] yield

(7.10) \begin{equation}\begin{aligned} \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) ={} & t^5 \int_\mathbb{R} g^{(4)}(st) \tilde H_1^u(s) {\text{d}} s - t^6 \int_\mathbb{R} g^{(6)}(s) {\mathbb{E}}[P_6(\tfrac{s}{t} + E_1 - E_0)] {\text{d}} s \\[6pt] & - t^{6} {\sum_{a\in D_{g^{(5)}} }} J_{g^{(5)}}(a) {\mathbb{E}} [P_{6}(\tfrac{a}{t} + E_1 - E_0)].\end{aligned}\end{equation}

As the last two terms in (7.10) are of order $o(t^5)$ , we only have to simplify the first term.

For all sufficiently small $t>0$ and all $s \in \mathbb{R}$ with $|s|\leq 4$ , the function $g^{(4)}$ is differentiable on the open interval with endpoints 0 and st, so there is a point $\xi_{st}$ in this interval such that

\[ g^{(4)}(st) = g^{(4)}(0) + g^{(5)}(\xi_{st})st\]

by the mean value theorem. Inserting this into the first term of (7.10) and using the fact that $g^{(5)}$ and $\tilde H_1^u$ are bounded yields

(7.11) \begin{equation} \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) = t^5 g^{(4)}(0) \int_{-4}^4 \tilde H_1^u(s) {\text{d}} s+o(t^5)\end{equation}

as $t\downarrow 0$ .

Noting that $P_4$ integrates to 0 on each interval of unit length, another application of Fubini’s theorem, the refined partial integration formula [Reference Kiêu7, Lemma 4.1], and Lemma 3.1 gives

(7.12) \begin{align} \nonumber \int_{-4}^4 \tilde H_1^u(s) {\text{d}} s & =\int_{-4}^4 H_1^u(s) {\text{d}} s \\ & = {\mathbb{E}}_u \bigg(\bigg[K_2({-}4)+ \sum_{x\in X^u\cap [{-}4,4]} J_{{K_2}}(x) -K_2(4)\bigg]\big(K_1(0)-{\mathbb{E}}_u K_1(0)\big)\bigg).\end{align}

The arguments that lead to (7.4) in combination with (7.3) where $r=-4$ and $k=8$ imply ${\mathbb{E}}_u [K_2({-}4)K_1(0)]={\mathbb{E}}_u [K_2(4)K_1(0)]$ , and the two marginal terms in the last expression of (7.12) cancel. Hence, by the refined Campbell theorem [Reference Schneider and Weil10, Theorem 3.5.3] and the translation covariance of $ J_{{K_2}}$ ,

(7.13) \begin{equation}\begin{aligned} \int_{-4}^4 \tilde H_1^u(s) {\text{d}} s & = {\mathbb{E}}_u \sum_{x\in {X^u}\cap [{-}4,4]} J_{{K_2}}(x) \big(K_1(0)-{\mathbb{E}}_u K_1(0)\big) \\ & = {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{-4}^4 \big(K_1(x)-{\mathbb{E}}_u K_1(0)\big) {\text{d}} x = \sum_{j=-3}^{2} \theta_j + Q,\end{aligned}\end{equation}

where

\[ \theta_j= {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{s_j}^{s_{j+1}} \bigl(K_1(x)-{\mathbb{E}}_u K_1(0)\bigr) {\text{d}} x\]

and

\begin{equation*} Q = {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{-4}^{s_{-3}} \big(K_1(x)-{\mathbb{E}}_u K_1(0)\big) {\text{d}} x + {\mathbb{E}}^0 J_{{K_2}}(0) \int_{s_3}^{4} \big(K_1(x)-{\mathbb{E}}_u K_1(0)\big) {\text{d}} x .\end{equation*}

Here we have used the fact that $s_3 \leq 4$ , $s_5\geq 4$ , $s_{-3} \geq -4$ , and $s_{-5} \leq -4$ under $\mathbb{P}^0_u$ (the Palm distribution of $X^u$ ). Using Lemma 4.1, it is seen that $ J_{{K_2}}(0) = (1/12) h_1^3$ and consequently that $\theta_j$ evaluates to $ {\mathbb{E}}^0_u (1/144) h_1^3(h_{j+1}^3-h_{j+1}{\mathbb{E}}^0_u h_1^3)$ . As $h_{j+1}$ only depends on the perturbations $E_{j+1}$ and $E_j$ , we conclude by independence that $\theta_j =0$ for all $\vert j \vert >1$ . Moreover, a coupling argument shows that $Q = \theta_3 + \theta_4 = 0$ . The Palm expectation of

\[ J_{{K_2}}(0) \int_{-4}^{s_{-3}} (K_1(x)-{\mathbb{E}}_u K_1(0)) {\text{d}} x\]

is unchanged when we put $E_{-3}=E_5$ , $E_{-4}=E_4$ , and $E_{-5}=E_3$ . Under this coupling assumption, $s_{-3}=s_5-8$ , $s_{-4}=s_{4}-8$ , $s_{-5}=s_{3}-8$ , $h_{-3}=h_5$ , and $h_{-4}=h_4$ , and hence

\begin{align*} Q & = {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{-4}^{s_{5}-8} (K_1(x+8)-{\mathbb{E}}_u K_1(0)) {\text{d}} x + {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{s_3}^{4} (K_1(x)-{\mathbb{E}}_u K_1(0)) {\text{d}} x \\ & = {\mathbb{E}}^0_u J_{{K_2}}(0) \int_{s_3}^{s_5} (K_1(x)-{\mathbb{E}}_u K_1(0)) {\text{d}} x = \theta_3 + \theta_4.\end{align*}

Summarizing, we obtain from (7.13) that

\begin{align*} \int_{-4}^4 \tilde H_1^u(s) {\text{d}} s & = \tfrac1{144}\sum_{j=-1}^{1}{\mathbb{E}}^0_u h_{1}^3\big(h_{j+1}^3-h_{j+1}{\mathbb{E}}^0_u h_1^3\big)\\ & = \tfrac{1}{8}\big(2\mu_4+\mu_2\mu_4-\mu_2^3-\mu_3^2\big) ,\end{align*}

where the last equality follows from lengthy and tedious—but elementary—calculations. Inserting this into (7.11) yields the assertion (2.10).

The expression (2.11) of the remainder is found by different arguments which are outlined in the first author’s thesis [Reference Stehr11].

7.2. Cumulative sampling

Before turning to the proof of Theorem 2.5, we state in Lemma 7.3 below that the covariance function of the Peano kernel decreases exponentially, from which admissibility follows.

The unit-intensity scaled cumulative process $X^u$ is a stationary point process with i.i.d. holding times $\{\omega_i\}_{i\in\mathbb{Z}}$ . We assume that $\omega_1$ has cumulative distribution function F with density with respect to Lebesgue measure such that $F(0)=0$ . Moreover, since $X^u$ has intensity 1, the holding times satisfy ${\mathbb{E}} \omega_i=1$ . To explicitly construct the point process $X^u$ , the first point $X_0$ of $X^u\cap(0,\infty)$ is chosen with cumulative distribution function G,

\[ G(x)=\int_0^x \bar F(y) {\text{d}} y , \qquad x\geq 0,\]

where $\bar F(y)=1-F(y)$ ; see e.g. [Reference Asmussen1, Chapter V: Corollary 3.6]. Note that the distribution G has density $\bar F$ . Given $X_0$ , the last point $X_{-1}$ of $X \cap ({-}\infty,0)$ (i.e. the largest point) is chosen according to $X_{-1}= X_0 - \omega^*$ , where $\omega^*$ is the conditional distribution of $\omega_0$ given $\omega_0 > X_0$ . This ensures that $X_{-1} < 0$ , and corrects [Reference Ziegel, Baddeley, Dorph-Petersen and Jensen14], where $\omega_0$ was used instead of $\omega^*$ . Having chosen increments $\{\omega_i\}_{i\neq 0}$ independent of $X_{-1},X_0$ , and setting $x_0=X_0$ , $x_i = X_0 + \sum_{\ell=1}^i \omega_\ell$ , and $x_{-i} = X_{-1} + \sum_{\ell=1}^{i-1} \omega_{-\ell}$ for all $i\in \mathbb{N}$ , we obtain a realization $X^u=\{x_i\}_{i\in \mathbb{Z}}$ of the cumulative point process. This construction implies that the point interval containing the origin has the length-weighted distribution, as expected.

The following lemma is stated in terms of the scaled unit-intensity cumulative process $X^u$ , but it is easily seen that it might as well have been formulated in terms of the process X with intensity $1/t$ .

Lemma 7.3. Let $n\in \mathbb{N}$ and $m\in\mathbb{N}_0$ be given, and let the unit-intensity process $X^u$ be from the cumulative model such that ${\mathbb{E}} {\text{e}}^{\eta \omega_1} < \infty $ for some $\eta > 0$ , and such that Assumption 2.1 is satisfied. Then

(7.14) \begin{equation} H_m^u(s)=O({\text{e}}^{-\epsilon s}), \quad\quad s\to\infty,\end{equation}

for some $\epsilon > 0$ . In particular, $X^u$ and X are strongly n-admissible.

Proof. The admissibility claim obviously follows from (7.14).

In the case of the trapezoidal estimator we can state the theorem without the integrability assumption (2.2). This is because finite moments of the Peano kernel only require (2.2) to be true for $j\in\mathbb{N}$ . As we assume that the increments have exponential moments, in particular they have finite moments of any positive order, and hence all integrability results for the Peano kernels apply.

The proof relies on exponential decays in renewal theory, and we refer to [Reference Asmussen1, Chapter V] for an introduction. Moreover, for fixed $n\in\mathbb{N}$ and all $m\in\mathbb{N}_0$ , we will explicitly use the fact that $K_m(s)$ depends on n points of the underlying point process to each side of s.

Let $n\in\mathbb{N}$ and $m\in\mathbb{N}_0$ be given. Let $N=(N(s))_{s\geq 0}$ be a pure renewal process with increments $\{\omega_i\}_{i\in\mathbb{N}}$ , and let U be the corresponding renewal measure. Also, let $y_0=0$ and ${y_i=\sum_{\ell=1}^i \omega_\ell}$ for $i\in \mathbb{N}$ ; that is, $X^u\cap (0,\infty) = \{X_0 + y_i\}_{i\in\mathbb{N}_0}$ . Then $y_i \sim F^{*i}$ for all $i\in \mathbb{N}$ . Define $\psi:[0,\infty) \to \mathbb{R}$ by

\[ \psi(s) = {\mathbb{E}}[K_m(s;\ N) \boldsymbol{1}_{{y_{n-1}\leq s}}] = {\mathbb{E}}^0_u[K_m(s) \boldsymbol{1}_{{s_{n-1}\leq s}}],\]

and let $\tilde \psi:[0,\infty) \to \mathbb{R}$ be given by

\begin{equation*} \tilde \psi(s) = {\mathbb{E}}[K_m(s;\ N) \boldsymbol{1}_{{y_{n-1}\leq s}}\boldsymbol{1}_{{y_n > s}}] = {\mathbb{E}}^0_u[K_m(s) \boldsymbol{1}_{{s_{n-1}\leq s}}\boldsymbol{1}_{{s_n > s}}] .\end{equation*}

Then $\psi(s) = U*\tilde \psi(s)$ . This can be seen by a renewal argument obtaining the renewal equation $\psi = \tilde \psi + \psi*F$ , which has the desired solution. Another rather intuitive approach is to condition on the nth-to-last point of N prior to s happening at time $x\in[0,s)$ , which has probability $U({\text{d}} x) \bar F* F^{*(n-1)}(s-x)$ . Now we initialize a new independent pure renewal process at time x and we obtain, integrating over [0, s],

\begin{align*} \psi(s) & = \int_0^s {\mathbb{E}} [K_m(s-x;\ N)\mid y_{n-1}\leq s-x, \, y_n>s-x] \bar F* F^{*(n-1)}(s-x) U({\text{d}} x) \\[3pt] & = \int_0^s \tilde \psi(s-x) U({\text{d}} x) \\[3pt] & = U*\tilde\psi(s) .\end{align*}

The exponential moment assumption implies $\bar F(s) = O({\text{e}}^{-\eta s})$ , which in turn implies that also $\bar G(s) = O({\text{e}}^{-\eta s})$ , as $s\to\infty$ . Moreover, $\overline{G* F^{*i}}(s)=O({\text{e}}^{-\eta s})$ , $s\to\infty$ , for all $i\in\mathbb{N}$ . We consider

\[ {\mathbb{E}}_u K_m(0)K_m(s) = {\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0+y_{2n-2} \leq s}} + {\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0+y_{2n-2}>s}}.\]

An application of the Cauchy–Schwarz inequality yields

\begin{align*} {\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0 + y_{2n-2}>s}} & \leq [{\mathbb{E}}_u K_m^2(s)K_m^2(0)]^{1/2} \mathbb{P}(X_0 + y_{2n-2} > s)^{1/2} \\[4pt] & \leq C\, \big( \overline{G*F^{*(2n-2)}}(s) \big)^{1/2}\end{align*}

for some finite C. Hence ${\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0+y_{2n-2}>s}} = O({\text{e}}^{-\eta s/2})$ as $s\to\infty$ , and (7.14) therefore follows once we show that

(7.15) \begin{equation} {\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0+y_{2n-2}\leq s}}= ({\mathbb{E}}_u K_m(0))^2 + O({\text{e}}^{-\epsilon s})\end{equation}

for some $\epsilon>0$ , as $s\to \infty$ . We apply a renewal argument conditioning on the nth arrival in $X^u\cap (0,\infty)$ , that is, conditioning on the value of $X_0+y_{n-1}\sim G*F^{*(n-1)}$ , and then initializing a new independent pure renewal process,

\begin{align*} & {\mathbb{E}}_u K_m(s)K_m(0)\boldsymbol{1}_{{X_0+y_{2n-2}\leq s}} \\[3pt] &\quad = \int_0^s {\mathbb{E}}_u \big[K_m(0)K_m(s)\boldsymbol{1}_{{X_0+y_{2n-2}\leq s}} \mid X_0+y_{n-1}=v \big] \big(G*F^{*(n-1)}\big)({\text{d}} v) \\[3pt] &\quad = \int_0^s {\mathbb{E}} \big[K_m(s-v;\ N)\boldsymbol{1}_{{y_{n-1}\leq s}}\big] {\mathbb{E}}_u \big[K_m(0)\mid X_0+y_{n-1}=v\big] \big(G*F^{*(n-1)}\big)({\text{d}} v) \\[3pt] &\quad = \int_0^s \psi(s-v){\mathbb{E}}_u \big[K_m(0)\mid X_0+y_{n-1}=v\big] \big(G*F^{*(n-1)}\big)({\text{d}} v) .\end{align*}

Since $ {\mathbb{E}}_u K_m(0)={\mathbb{E}}^0_u J_{K_{m+1}}(0)$ by (4.3), an application of Fubini’s theorem yields

\begin{equation*}\begin{aligned} \int_0^\infty \tilde \psi(s) {\text{d}} s & = {\mathbb{E}} \int_{y_{n-1}}^{y_n} K_m(s;\ N) {\text{d}} s \\[2pt] & = {\mathbb{E}} K_{m+1}(y_{n-1}^+;\ N) - {\mathbb{E}} K_{m+1}(y_n^-;\ N) = {\mathbb{E}}_u K_m(0),\end{aligned}\end{equation*}

and consequently, by [Reference Asmussen1, Chapter VII: Theorem 2.10(iii)],

(7.16) \begin{equation} \psi(s) = U*\tilde \psi(s) = {\mathbb{E}}_u K_m(0) + O({\text{e}}^{-\epsilon' s})\end{equation}

for some $0<\epsilon' < \eta$ , as $s\to\infty$ . Furthermore, by another application of the Cauchy–Schwarz inequality, we conclude that

(7.17) \begin{equation}\begin{aligned} \int_0^s {\mathbb{E}}_u \big[K_m(0)\mid X_0+y_{n-1}=v\big] \big(G*F^{*(n-1)}\big)({\text{d}} v) \\ & = {\mathbb{E}}_u K_m(0) - {\mathbb{E}}_u K_m(0)\boldsymbol{1}_{{X_0+y_{n-1}>s}} \\& = {\mathbb{E}}_u K_m(0) + O({\text{e}}^{-\eta s/2})\end{aligned}\end{equation}

as $s\to\infty$ . Combining (7.16) and (7.17) yields (7.15).

As for the perturbed model, Theorem 2.5 is stated without Assumption 2.1. This is because the strong admissibility and the variance decomposition are satisfied for the trapezoidal rule, when (2.2) is assumed only for $j\in\mathbb{N}$ . This relaxed assumption is ensured by the finite exponential moments of the increments.

Proof of Theorem 2.5. Let $m\in\{0,1\}$ . As X is strongly admissible, Proposition 6.1 in combination with the decrease rate (7.14) yields the variance decomposition

\[ \mathbb{V}\text{ar} (\hat{V}_{1}({\kern1.3pt}f)) = ({-}1)^{m+1} t^{2m+2} 2 g^{(2m+1)}(0^+) H_m^u(0) + o(t^{2m+2}) .\]

From (7.8) and (7.9) and the fact that ${\mathbb{E}}^0_u h_1^j=\nu_j$ , we conclude (2.12) and (2.13).

8. An application in stereology

To illustrate the general theory, we describe a geometric application that was also the original motivation for this work. In stereology, the volume of a compact object $Y\subset \mathbb{R}^3$ can be approximated from sections with equidistant and parallel planes with joint normal direction $\nu$ in the unit sphere $S^2$ , if the area of each intersection profile is accessible; see [Reference Baddeley and Jensen3, Chapter 7].

Formally, if f(x) is the area of the intersection of Y with the plane $\{y\in \mathbb{R}^3\,:\, \nu^Ty=x\}$ positioned at a signed distance $x\in\mathbb{R}$ from the origin along $\nu$ , the integral $\int f {\text{d}} x$ coincides with the volume of Y by Fubini’s theorem. If f(x) is available at all points of the equidistant stationary point process $X=t(U+\mathbb{Z})$ , the volume-estimator (1.2) can be used and is called the (classical) Cavalieri estimator. When f(x) is known at the points of a stationary point process X with intensity $1/t$ , the so-called (generalized) Cavalieri estimator (1.1) can be used. However, as outlined above, the generalized Cavalieri estimator does not exploit the smoothness of f and thus has a suboptimal decrease rate as $t\downarrow 0$ . We therefore suggest employing Newton–Cotes estimators (of appropriate order) instead.

This section is devoted to Monte Carlo simulations illustrating the advantage of the new estimators when the sampling points are not equidistant. For illustration purposes we start by considering the Euclidean unit ball $Y=\{z\in\mathbb{R}^3\::\: \lVert{z}\rVert\leq 1\}$ . In this case, the measurement function is

\[ f(x)=\boldsymbol{1}_{{[{-}1,1]}}(x) \pi (1-x^2),\]

which is a $(1,\infty)$ -piecewise smooth function, as fʹ has jumps and is piecewise linear. Applying the classical Cavalieri estimator to such a function yields the extension term $\mathbb{V}\text{ar}_E (\hat{V}({\kern1.3pt}f))=\tfrac{\pi^2}{90}t^4$ from (2.5). Using sampling by the perturbed model or the model with cumulative errors we expect that the generalized Cavalieri estimator decreases at a rate of 3 or 1, respectively, whereas the trapezoidal estimator ( $n=1$ ) and Simpson’s estimator ( $n=2$ ) decrease at a rate of 4 in both point-models, an asymptotic behaviour visible in Figure 1 below. This figure shows the empirical variances of these three estimators based on 2000 Monte Carlo simulations as functions of the mean number of sections, that is $2/t$ , with the variance plot including the extension term of the classical Cavalieri estimator and the extension term of the trapezoidal estimator as given by the dominating terms in (2.9) and (2.13) for the perturbed and cumulative model, respectively. The variances in this and the following figures are shown in a double-logarithmic scale with $\alpha$ and $\hat{\alpha}$ being the theoretical and approximate rates of decrease ( $\hat{\alpha}$ has been found by the least squares method applied to the datapoints where $15 \leq 2/t \leq 40$ ).

Figure 1:

Empirical variance for the volume estimation of the unit ball in $\mathbb{R}^3$ based on perturbed sampling with $E_i \sim \text{Unif}(({-}s,s))$ and sampling with cumulative errors with $\omega_i \sim \text{Unif}((1-c,1+c))$ . We choose s and c such that the average relative deviation (the coefficient of error) of the point-increment from the ideal increment 1 is 5%. In both figures, the graph of the trapezoidal estimator (solid dark grey) is almost completely hidden by the graph of Simpson’s estimator (solid grey), and the trapezoidal extension term (dashed dark grey) is almost identical to the classical extension term (dashed black).

Figure 2:

Empirical variance for the volume estimation of a spindle-shaped body of revolution in $\mathbb{R}^3$ based on perturbed sampling with $E_i\sim \text{Unif}(({-}s,s))$ . We choose s such that the coefficients of error (CE) of the point-increments are 5% (left) and 20% (right).

The graphs of Figure 1 are characteristic for the behaviour of variances and extension terms for objects with (1,1)-piecewise smooth measurement functions. For instance, ellipsoids or, more generally, strictly convex bodies lead to the same variance behaviour, apart from the facts that intercepts of these curves may be shifted and the Zitterbewegung may differ.

For comparison, we therefore give another example, where the measurement function exhibits a higher order of smoothness. The measurement function

\[ f(x)=\boldsymbol{1}_{{[{-}1,1]}}(x) \frac{\pi}{2} (1+\cos \pi x)\]

is obtained from a spindle-shaped body of revolution, if all section planes are orthogonal to the rotation axis. The corresponding convex body is illustrated in [Reference Kiderlen and Dorph-Petersen6, Fig. 4]. The measurement function f is $(2,\infty)$ -piecewise smooth. Using this measurement function, the extension term of the classical estimator is $\mathbb{V}\text{ar}_E (\hat{V}({\kern1.3pt}f))=(\pi^6/60480)t^6$ . Figure 2 shows empirical variances based on perturbed sampling with the two extension terms included, where the extension term of the trapezoidal estimator is given as the sum of the dominating terms in (2.10) and (2.11). In Figure 2a we use small perturbations to illustrate the fact that the dominating term in (2.10) can be made arbitrarily small. Hence, a decrease rate of 6 for the variance of the trapezoidal estimator can be a good approximation with small perturbations, as the trapezoidal extension term is approximately given by $1.7 \cdot 10^{-4} t^5 + 3.0 \cdot 10^{-2}t^6$ here. Even when we consider $100 \leq 2/t \leq 200$ , we only obtain an approximate decrease rate of $\hat{\alpha} = 5.66$ . For comparison, Figure 2b gives a better illustration of the actual asymptotic rate of decrease which corresponds to the bound from Theorem 2.3, that is, $\alpha = 5$ . Here we use larger perturbations, which in turn gives an approximate trapezoidal extension term of $0.043 t^5 + 0.25 t^6$ . When we increase the number of intersecting planes to $2/t\leq 100$ , the actual rate becomes even more apparent, as we here obtain an approximate decrease rate of $\hat{\alpha} =5.19$ .

The last two simulations are meant to illustrate the findings in Theorem 2.3 for point process models where we do not have explicit formulae for the extension term. The first is the already discussed model with accumulated errors. To illustrate the wide range of point process models to which our results apply, we also simulated from the Matérn hard core process of type II; see [Reference Schneider and Weil10, Section 3.5, pp. 93-94], which satisfies the strong integrability assumption (2.2). The empirical variances for the aforementioned spindle-shaped body are depicted in Figure 3. It is worth noticing that the variance of the trapezoidal estimator under the Matérn model seem to satisfy the strong bound of Theorem 2.3, that is, (2.7). When the number of intersecting planes is increased to $2/t\leq 100$ the result is clearer, as we find approximate decrease rates of $\hat{\alpha} =4.94$ and $\hat{\alpha} =6.07$ for the trapezoidal estimator and Simpson’s estimator, respectively.

Figure 3:

Empirical variance for the volume estimation of a spindle-shaped body in $\mathbb{R}^3$ based on sampling with cumulative errors with $\omega_i\sim \text{Unif}((1-c,1+c))$ and sampling with a Matérn hard core process of type II with intensity 1 and a hard core distance of 0.4. The value of c is chosen so that the coefficient of error of the increment is 5%.