Image Anal Stereol 2017;36:43-49 doi:10.5566/ias.1681
Original Research Paper
ESTIMATION OF SAMPLE SPACING IN STOCHASTIC PROCESSES
ANDERS RØNN-NIELSENB,1, JON SPORRING2 AND EVA B. VEDEL JENSEN3
1Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 København Ø,
Denmark; 2Department of Computer Science, University of Copenhagen, Universitetsparken 1, 2100 København
Ø, Denmark; 3Department of Mathematics, Aarhus University, Ny Munkegade 118, 8000 Aarhus C, Denmark
e-mail: arnielsen@math.ku.dk, sporring@di.ku.dk, eva@math.au.dk
(Received December 15, 2016; revised March 7, 2017; accepted March 7, 2017)
ABSTRACT
Motivated by applications in electron microscopy, we study the situation where a stationary and isotropic
random field is observed on two parallel planes with unknown distance. We propose an estimator for this
distance. Under the tractable, yet flexible class of Lévy-based random field models, we derive an approximate
variance of the estimator. The estimator and the approximate variance perform well in two simulation studies.
Keywords: Lévy-based modelling, Monte Carlo methods, random fields, stereology, section distance,
systematic sampling, variance, variogram.
INTRODUCTION
Systematic sampling is a wide-spread technique in
microscopy, stereology and in the spatial sciences in
general. In 3D, systematic sampling is often performed
in two steps. In a first step, a stack of equally spaced,
parallel plane sections through the spatial structure of
interest is generated. In a second step, observations are
usually made at a regular lattice of grid points within
each section, see, e.g., Baddeley et al. (2005, Chapter
7).
Real sampling procedures may involve technically
very challenging physical cutting. An example is
the generation of ultra thin sections for electron
microscopy where the intended distance h between
neighbour sections is in the order of nanometers. In
such cases, the actual section distance may deviate
from the intended distance.
In Ziegel et al. (2010; 2011), the effect of such
errors on the variance of a Cavalieri type of estimator
has been studied. (For a short account of the Cavalieri
estimator, see Baddeley et al. (2005, Section 7.1).)
Here, the aim is to estimate an integral Θ =
∫
R f (x)dx,
using the estimator Θ̂ = h∑k f (xk), where xk = u +
kh ∈ R, k ∈ Z, is a systematic set of points. A simple
geometric example concerns the estimation of the
volume of a bounded object in R3 in which case f (x)
may be the area of the intersection between the object
and a horizontal plane at height x ∈ R. Without errors,
xk− xk−1 ≡ h. In Ziegel et al. (2010; 2011), the effect
of errors in section positions is studied under different
spatial point process models for {xk}.
If the aim is to reconstruct the object rather
than estimating quantitative properties of the object,
it becomes important to know the actual realized
distances between neighbour sections. A concrete
example of this situation may be found in Sporring
et al. (2014) where ultra thin electron microscopy
sections are analysed.
In the present paper, we take up this problem. We
study the situation where a stationary and isotropic
random field {Xv : v ∈ R3} is observed on two parallel
planes L1 and L2 with unknown distance h. More
specifically, we observe the random field at a regular
lattice of grid points {v1i} on L1 and, similarly, at {v2i}
on L2 where ‖v1i−v2i‖= h for all i, see the illustration
in Fig. 1. The idea is now to obtain information about
the correlation structure of the field, using the observed
values within sections. This information is then used
in the estimation of the unknown distance h between
sections.
Fig. 1. Illustration of the sampling of the random field
on two parallel planes a distance h apart.
This procedure has been applied with success in
Sporring et al. (2014). In the present paper, we study
43
RØNN-NIELSEN A ET AL: Estimation of sample spacing in stochastic processes
the statistical properties of the resulting estimator of
h. Under the tractable, yet flexible class of Lévy-
based random field models, we derive an approximate
variance of the estimator. We also examine the utility
of the estimator and the variance approximation in two
simulation studies.
For an introduction to Lévy-based modelling, see,
e.g., Jónsdóttir et al. (2013b). In relation to variance
estimation, Lévy-based models have earlier been used
in circular systematic sampling (Jónsdóttir et al.,
2013a). Lévy-based models have also been shown to
be a useful modelling tool for Cox point processes and
growth (Hellmund et al., 2008; Jónsdóttir et al., 2008).
The paper is organised as follows. First, we
introduce the estimator and explore some of its
properties. Under the Lévy-based random field model,
we then derive an approximate variance of the
estimator and study its performance in two simulation
studies. The derivation of the variance formula and
some integral calculations relating to one of the
simulation studies are deferred to two appendices.
THE ESTIMATOR OF SECTION
DISTANCE
Consider a stationary and isotropic random field
{Xv : v ∈R3}. Since the random field is stationary and
isotropic, the variogram E(Xv1 − Xv2)2 only depends
on ‖v1−v2‖ and is thereby determined by γ : [0,∞)→
[0,∞), where
E(Xv1−Xv2)
2 = γ(‖v1− v2‖) .
We will assume that γ is strictly increasing. Note
that a strictly increasing γ is equivalent to a
strictly decreasing covariance function. A wide range
of covariance models has this property, including
the spherical, Gaussian, exponential and 3rd order
autoregressive covariance functions (Jónsdóttir et al.,
2013b).
Let L1 and L2 be two parallel planes in R3 with an
unknown distance h. Assume without loss of generality
that
L1 = {(x,y,0) | x,y ∈ R} ,
L2 = {(x,y,h) | x,y ∈ R} .
We furthermore assume that pairs of data points
(Xv1i ,Xv2i), i = 1, . . . ,n, are observed where, for all
i = 1, . . . ,n, v1i ∈ L1, v2i ∈ L2 and the points are on
the form
v1i = (xi,yi,0) , v2i = (xi,yi,h) ,
see also Fig. 1. The goal is to estimate h and to
find an approximate expression for the variance of the
estimator.
Since
S =
1
n
n
∑
i=1
(Xv1i−Xv2i)
2 (1)
is an unbiased estimator of γ(h), we propose to
estimate h by
ĥ = γ−1(S) . (2)
From a first–order Taylor expansion of γ−1 around
γ(h) we find
ĥ≈ h+(γ−1)′(γ(h))
(
S− γ(h)
)
= h+
1
γ ′(h)
(
S− γ(h)
)
.
From this it is seen that
E(ĥ−h)2 ≈ 1
γ ′(h)2
E
(
S− γ(h)
)2
,
so
Var(ĥ)≈ 1
γ ′(h)2
Var(S) .
Using a higher order Taylor expansion above would
lead to substantially more involved expressions for
the approximate variance of ĥ. A first–order Taylor
expansion is adequate if γ−1 is approximately linear
in the set, where S typically varies. In a concrete
application, this can be investigated, using the estimate
of Var(S) derived in the next section.
In most applications, the variogram γ is unknown.
In order to apply the estimator Eq. 2, it is therefore
needed to replace γ in Eq. 2 by an estimate,
based on the available observations within sections.
First, we recommend to check the stationarity and
isotropy assumptions by, e.g., estimating the variogram
separately in different areas and directions within
each plane. Here we use that the distances between
sample points within sections are known. The overall
variogram γ can be estimated by fitting a parametric
curve, induced by one of the known covariance
models, to the pooled empirical variogram. For a
recent application of this procedure see Jónsdóttir et
al. (2013b, Section 4).
Note that if h is large, then the observations on the
two planes may be almost uncorrelated, and it becomes
difficult or impracticable to apply the estimator. More
specifically, let us suppose that
Cov(Xv1 ,Xv2)→ 0 ,
as ‖v1− v2‖→ ∞. Then,
γ(h)→ 2Var(Xv) ,
44
Image Anal Stereol 2017;36:43-49
as h→ ∞. If h is so large that γ(h) ≈ 2Var(Xv), then
it may happen that S > 2Var(Xv) and in such cases ĥ
does not exist.
VARIANCE ESTIMATION UNDER
A LÉVY–BASED MODEL
THE LÉVY–BASED RANDOM FIELD
MODEL
It is possible to find an explicit expression for γ
and Var(S) under a Lévy-based random field model.
Assume that Z is an independently scattered infinitely
divisible random measure on R3. Such a measure is
called a Lévy basis (Barndorff-Nielsen et al., 2004)
and references therein for details. The measure Z
thereby satisfies that, for B⊆R3, Z(B) has an infinitely
divisible distribution and, for B1,B2 ⊆ R3 disjoint,
Z(B1) and Z(B2) are independent.
We assume that Z is stationary and isotropic and
that the field {Xv : v ∈ R3} is given by
Xv =
∫
R3
f (u+ v)Z(du) ,
where f is a kernel function, satisfying that f (v) =
f (‖v‖) and
∫
R3 f (v)dv < ∞. These assumptions make
the field {Xv : v ∈ R3} stationary and isotropic and
furthermore, it is possible to compute cumulants
and covariances for the field. The Lévy–based
models provide a very flexible and broad class of
models, that contains the important cases of Gaussian,
gamma, inverse Gaussian and normal inverse Gaussian
distributions. Under mild regularity conditions, we
have
Cov(Xv1 ,Xv2) = Var(Z
′)
∫
R3
f (u+ v1) f (u+ v2)du
= Var(Z′)
∫
R3
f (u) f (u+ v1− v2)du
= Var(Z′)K(v1− v2) ,
say, where Z′ is the associated spot variable, having the
distribution of Z(B), when B ⊆ R3 has volume 1. Due
to stationarity and isotropy, K(v) = K(‖v‖), say, and
the variogram is of the form
γ(h) = 2Var(Z′)(K(0)−K(h)) , h≥ 0 . (3)
Furthermore, the n’th cumulant of Xv can be expressed
as
κn(Xv) = κn(Z′)
∫
R3
f (u+ v)n du .
In Table 1 below, the distributions and the first four
cumulants of the spot variable are listed for the gamma,
the inverse Gaussian and the normal inverse Gaussian
bases. For more details see Jónsdóttir et al. (2013b,
Section 5).
THE VARIANCE OF S
In this section, we derive an expression for the
variance of S as defined in Eq. 1.
Proposition 1 Under the Lévy–based random field
model, we have
Var(S) =
1
n2
n
∑
i, j=1
[
κ4(Z′)ρ4(v1i,v1 j)
+2Var(Z′)2ρ2(v1i,v1 j)2
]
, (4)
where
ρ2(v1i,v1 j) =
∫
R3
(
f (u+ v1i)− f (u+ v2i)
)
·
(
f (u+ v1 j)− f (u+ v2 j)
)
du , (5)
ρ4(v1i,v1 j) =
∫
R3
(
f (u+ v1i)− f (u+ v2i)
)2
·
(
f (u+ v1 j)− f (u+ v2 j)
)2 du . (6)
Below, we indicate how this result is derived.
First, notice that
Xv1i−Xv2i =
∫
R3
[
f (u+ v1i)− f (u+ v2i)
]
Z(du) ,
i = 1, . . . ,n. For ease of notation, let us write
Yi = Xv1i−Xv2i , gi(u) = f (u+ v1i)− f (u+ v2i) .
Table 1. The spot variable and its cumulants for various Lévy bases.
Basis Gamma Inverse Gaussian Normal inverse Gaussian
Z(du) Γ(α du,λ ) IG(δ du,γ) NIG(α,β ,µ du,δ du)
Z′ Γ(α,λ ) IG(δ ,γ) NIG(α,β ,µ,δ )
κ1(Z′) α/λ δ/γ µ +δβ/(α2−β 2)1/2
κ2(Z′) α/λ 2 δ/γ3 δα2/(α2−β 2)3/2
κ3(Z′) 2α/λ 3 3δ/γ5 3δβα2/(α2−β 2)5/2
κ4(Z′) 6α/λ 4 15δ/γ7 3δ (α2 +4β 2)α2/(α2−β 2)7/2
45
RØNN-NIELSEN A ET AL: Estimation of sample spacing in stochastic processes
We find
Var(S) = E(S2)− (ES)2
=
1
n2
n
∑
i, j=1
E(Y 2i Y 2j )− (E(Y 21 ))2 . (7)
These moments can be calculated, using a well-known
result for the logarithm of the Laplace transform of an
integral with respect to a Lévy basis. The result holds if
E(eεZ′)< ∞ for some ε > 0. This condition is fulfilled
for the three Lévy bases listed in Table 1. Let for
λ1,λ2 ∈ R sufficiently close to 0 and i, j ∈ {1, . . . ,n}
Ki, j(λ1,λ2) = logE
(
eλ1Yi+λ2Y j
)
. (8)
Then, see, e.g., Hellmund et al. (2008, Eq. 10),
Ki, j(λ1,λ2) =
∫
R3
logE
(
e(λ1gi(u)+λ2g j(u))Z
′)
du . (9)
By differentiating Eqs. 8 and 9 with respect to λ1 and
λ2 and equating the two expressions at λ1 = λ2 = 0
with each other, we obtain the variance formula Eq. 4.
The details can be found in Appendix A.
For the Gaussian kernel function
f (x) =
1
(2πσ2)3/2
e−
1
2σ2
‖x‖2
, x ∈ R3 , (10)
with σ2 > 0, it is possible to calculate ρ2 and ρ4
explicitly, see Appendix B. We find
ρ2(v1i,v1 j)
=
1√
2
1
(2πσ2)3/2
e−
1
4σ2
‖v1i−v1 j‖2
(
1− e−
1
4σ2
h2
)
(11)
and
ρ4(v1i,v1 j) =
1
4(2πσ2)9/2
e−
1
2σ2
‖v1i−v1 j‖2
·
(
1−4e−
3
8σ2
h2
+3e−
1
2σ2
h2
)
.
(12)
These expressions for ρ2 and ρ4 may be inserted into
Eq. 4 to get an explicit expression for the variance of S
in the case of a Gaussian kernel function.
As we shall see in the next subsection, further
simplifications may be obtained by approximating the
double sum in Eq. 4 by integrals.
FURTHER APPROXIMATIONS
Assume that the points v1i = (xi,yi,0) in L1
constitute a regular lattice of grid points. Let A be
the subset of R2 spanned by the points (xi,yi), and
assume that A has area a. If the area of the fundamental
region of the lattice is sufficiently small, the formula
for Var(S) can be approximated as follows
Var(S) =
1
n2
n
∑
i, j=1
[
κ4(Z′)ρ4(v1i,v1 j)
+2Var(Z′)2ρ2(v1i,v1 j)2
]
≈ 1
a2
∫
A
∫
A
[
κ4(Z′)ρ4(v1,v2)
+2Var(Z′)2ρ2(v1,v2)2
]
dv1dv2 .
If f (x) ≈ 0 for ‖x‖ > d, where d is substantially
smaller than the width of A, we have the further
approximation
Var(S)≈ 1
a2
∫
A
∫
R2
[
κ4(Z′)ρ4(v1,v2)
+2Var(Z′)2ρ2(v1,v2)2
]
dv1dv2 .
Since both ρ2(v1,v2) and ρ4(v1,v2) only depend on
‖v1 − v2‖, the inner integral does not depend on v2.
Hence, we have
Var(S)≈ 1
a
∫
R2
[
κ4(Z′)ρ4(v,v0)
+2Var(Z′)2ρ2(v,v0)2
]
dv ,
(13)
where v0 ∈ A is a fixed point. By the same type
of approximation argument, the variance can also be
approximated by
Var(S)≈ 1
a
∫
A
[
κ4(Z′)ρ4(v,v0)
+2Var(Z′)2ρ2(v,v0)2
]
dv .
(14)
MONTE CARLO METHODS
If the integrals in Eqs. 13 and 14 are too
complicated to calculate, a solution could be to
estimate them by Monte Carlo methods. Assume for
convenience that A contains v0 = (0,0), and define
furthermore h̄ = (0,0,h) and v̄ = (v1,v2,0) for v =
(v1,v2) ∈ A. Then,∫
A
ρ4(v,v0)dv =
∫
A
∫
R3
(
f (u+ v̄)− f (u+ v̄+ h̄)
)2
·
(
f (u)− f (u+ h̄)
)2 dudv
and∫
A
ρ2(v,v0)2dv =
∫
A
∫
R3
∫
R3
(
f (u+ v̄)− f (u+ v̄+ h̄)
)
·
(
f (u)− f (u+ h̄)
)
·
(
f (ũ+ v̄)− f (ũ+ v̄+ h̄)
)(
f (ũ)− f (ũ+ h̄)
)
dudũdv .
46
Image Anal Stereol 2017;36:43-49
It follows that
∫
A ρ4(v,v0)dv can be estimated from the
empirical mean of repeated independent simulations of(
f (U +V̄ )− f (U +V̄ + h̄)
)2 ( f (U)− f (U + h̄))2
g(U)h(V )
where U = (U1,U2,U3) and V = (V1,V2) are
independent, U is simulated according to some
probability density g on R3 and V is simulated
according to some probability density h on A. The
value of
∫
A ρ2(v,v0)2dv can be estimated similarly.
SIMULATION STUDIES
GAUSSIAN KERNEL FUNCTION
Assume that f is the Gaussian kernel function
Eq. 10, implying a Gaussian covariance. Then, using
Eqs. 11 and 12, the approximation Eq. 13 becomes
Var(S)
≈ 1
a
∫
R2
[
κ4(Z′)ρ4(v,v0)+2Var(Z′)2ρ2(v,v0)2
]
dv
=
[
κ4(Z′)
1
(2πσ2)9/2
(
1
4 +
3
4 e
− 1
2σ2
h2− e−
3
8σ2
h2
)
+2Var(Z′)2
1
2
1
(2πσ2)3
(
1− e−
1
4σ2
h2
)2]
· 1
a
∫
R2
e−
1
2σ2
(v21+v
2
2) dv1dv2
=
1
a
1
(2πσ2)2
·
[
κ4(Z′)
1
(2πσ2)3/2
(
1
4 +
3
4 e
− 1
2σ2
h2− e−
3
8σ2
h2
)
+Var(Z′)2
(
1− e−
1
4σ2
h2
)2]
.
Furthermore, cf. Eq. 3 and Appendix B, we get with
h̄ = (0,0,h)
γ(h) = 2Var(Z′)
(∫
R3
f (u)2 du−
∫
R3
f (u) f (u+ h̄)du
)
= 2Var(Z′)
1
(4πσ2)3/2
(
1− e−
1
4σ2
h2
)
,
from which we find
γ ′(h) = Var(Z′)
1
(4πσ2)3/2
h
σ2
e−
1
4σ2
h2
.
Combining these results, we find
Var(ĥ)≈ 1
γ ′(h)2
Var(S)
≈ 16πσ
6
ah2
[
κ4(Z′)
Var(Z′)2
1
(2πσ2)3/2
·
(
1
4 e
1
2σ2
h2
+ 34 − e
1
8σ2
h2
)
+
(
e
1
4σ2
h2−1
)2]
.
(15)
In this example, we combine the Gaussian kernel
function with a normal inverse Gaussian (NIG) Lévy
basis for which the spot variable Z′ is NIG(α,β ,µ,δ )-
distributed. For such a Lévy basis,
κ4(Z′)
Var(Z′)2
=
3(α2 +4β 2)
δα2
√
α2−β 2
,
see Table 1. In Fig. 2, the resulting approximate
variance Eq. 15 of ĥ is plotted as a function of h for
a concrete choice of parameter values. The chosen
values of the parameters give a covariance structure
similar to that obtained in a concrete analysis presented
in Jónsdóttir et al. (2013b). In Fig. 2, the set A has been
chosen such that the effective support of the kernel
function is much smaller than the width of A.
0 2 4 6 8 10
0
40
80
12
0
h
V
ar
ĥ
Fig. 2. A plot of the approximate expression for Var(ĥ)
as a function of h. The kernel function is Gaussian
with σ2 = 5. The parameters of the NIG distribution
associated with the Lévy basis are α = 0.6, β = 0.4,
µ = 2.4 and δ = 2. Furthermore, A = [0,100]2.
Note that the variance increases substantially as a
function of h.
We have simulated from the model described
above for four values of h: 0.2, 0.4, 1, 2. In each
case, we made 1000 simulated values of ĥ each based
on simulation of two parallel fields with distance h.
From the simulated values of ĥ we calculated the
empirical mean Êĥ and the empirical variance V̂ar(ĥ).
The results of this can be seen in Table 2 together with
the approximate variance Var(ĥ), using the formula
Eq. 15, and estimates Ṽar(ĥ) of the approximate
variance, based on Monte Carlo simulation. Fig. 3
47
RØNN-NIELSEN A ET AL: Estimation of sample spacing in stochastic processes
Table 2. For four values of h, Êĥ is the empirical mean of 1000 values of ĥ simulated with a Gaussian kernel
function and a NIG Lévy basis. The empirical variance is denoted V̂ar(ĥ). Furthermore, Var(ĥ) is the approximate
variance, and Ṽar(ĥ) is the Monte Carlo estimate of the approximate variance. The parameters of the model are
as specified in Fig. 2.
h 0.2 0.4 1 2
Êĥ 0.1999647 0.3997637 1.001492 1.994434
V̂ar(ĥ) 6.512 ·10−5 2.638 ·10−4 1.640 ·10−3 7.692 ·10−3
Var(ĥ) 6.298 ·10−5 2.534 ·10−4 1.652 ·10−3 7.703 ·10−3
Ṽar(ĥ) 6.313 ·10−5 2.587 ·10−4 1.675 ·10−3 7.951 ·10−3
Table 3. For four values of h, Êĥ is the empirical mean of 1000 values of ĥ simulated with an exponential kernel
function and a NIG Lévy basis. The empirical variance is denoted V̂ar(ĥ). Furthermore, Ṽar(ĥ) is the Monte Carlo
estimate of the approximate variance. The parameter of the kernel function is σ = 0.56, while the parameters of
the NIG distribution are as in Fig. 2.
h 0.2 0.4 1 2
Êĥ 0.200 0.400 0.998 1.996
V̂ar(ĥ) 4.732 ·10−5 1.813 ·10−4 1.260 ·10−3 7.352 ·10−3
Ṽar(ĥ) 4.736 ·10−5 1.925 ·10−4 1.358 ·10−3 7.434 ·10−3
shows four jointly simulated parallel fields from the
model such that the last three fields are at distances
0.4, 1 and 2 from the first.
3.8
4.0
4.2
4.4
4.6
20 40 60 80 100
20
40
60
80
100
3.8
4.0
4.2
4.4
4.6
20 40 60 80 100
20
40
60
80
100
3.8
4.0
4.2
4.4
4.6
20 40 60 80 100
20
40
60
80
100
3.8
4.0
4.2
4.4
4.6
20 40 60 80 100
20
40
60
80
100
Fig. 3. Four jointly simulated parallel fields with NIG
Lévy basis and a Gaussian kernel function. The three
last fields are at distance 0.4, 1 and 2 from the first.
The parameters of the model are as specified in Fig. 2.
It is seen from Table 2 that the bias of ĥ is
negligible. Furthermore, the variance approximation
Eq. 15 and its evaluation by Monte Carlo methods
perform well.
EXPONENTIAL KERNEL FUNCTION
Now assume that f is an exponential kernel
function
f (x) =
σ3
8π
e−σ‖x‖,
with σ > 0. As shown, e.g., in Jónsdóttir et al. (2013b,
Section 3), this kernel function induces a 3rd order
autoregressive covariance function for the Lévy-based
random field. For this kernel function, expressions for
ρ2 and ρ4 in the variance approximations Eq. 13 and
Eq. 14 are not available. Accordingly, the approximate
variances need to be evaluated, using the Monte Carlo
method. An example is shown in Table 3 for a NIG
Lévy basis. Again, we have simulated 1000 values of
ĥ for the four choices of h: 0.2, 0.4, 1 and 2.
ACKNOWLEDGEMENTS
This research was supported by Centre for
Stochastic Geometry and Advanced Bioimaging,
funded by a grant from the Villum Foundation.
REFERENCES
Baddeley A, Jensen E B V (2005). Stereology for
statisticians. Boca Raton: Chapman & Hall/CRC Press.
Barndorff-Nielsen O, Schmiegel J (2004). Lévy-based
tempo-spatial modelling with applications to
turbulence. Uspekhi Mat Nauk 159:63–90.
48
Image Anal Stereol 2017;36:43-49
Hellmund G, Prokešová M, Jensen E B V (2008). Lévy-
based Cox point processes. Adv Appl Prob 40:603–29.
Jónsdóttir K Ý, Jensen E B V (2013a). Lévy-based error
prediction in circular systematic sampling. Image Anal
Stereol 32:117-25.
Jónsdóttir K Ý, Rønn-Nielsen A, Mouridsen K, Jensen E
B V (2013b). Lévy-based modelling in brain imaging.
Scand. J. Stat. 40(3):511–29.
Jónsdóttir K Ý, Schmiegel J, Jensen E B V (2008). Lévy-
based growth models. Bernoulli 14(1):62–90.
Kendall M G, Stuart A (1977). The advanced theory of
statistics, Vol. 1. Distribution theory, 4th Ed. London:
Griffin.
Sporring J, Khanmohammadi M, Darkner S, Nava N,
Nyengaard J R, Jensen E B V (2014). Estimating the
thickness of ultra thin sections for electron microscopy
by image statistics. In: Proc 2014 IEEE Int Symp
Biomed Imaging, Beijing. 157–60.
Ziegel J, Baddeley A, Dorph-Petersen K-A, Jensen E B
V (2010). Systematic sampling with errors in sample
locations. Biometrika 97:1–13.
Ziegel J, Jensen E B V, Dorph-Petersen K-A (2011).
Variance estimation for generalized Cavalieri
estimators. Biometrika 98:187–98.
APPENDIX A: DERIVATION OF
VAR(S)
In this Appendix, we derive the variance formula
Eq. 4 by differentiating Eqs. 8 and 9, and equating
the two obtained expressions at λ1 = λ2 = 0 with each
other. Recall that
Yi = Xv1i−Xv2i , gi(u) = f (u+ v1i)− f (u+ v2i) .
It is easy to differentiate Eq. 9 with respect to λ1 and
λ2. Notice that
Ki, j(λ1,λ2) =
∫
R3
logE
(
e(λ1gi(u)+λ2g j(u))Z
′
)
du
=
∫
R3
K
(
λ1gi(u)+λ2g j(u)
)
du ,
where K is the cumulant generating function of Z′. It
follows that
∂ k+`
∂λ k1 ∂λ
`
2
Ki, j(λ1,λ2)
∣∣∣
λ1=λ2=0
= κk+`(Z′)
∫
R3
gi(u)kg j(u)` du ,
(16)
where κk(Z′) is the k’th cumulant of Z′.
In order to differentiate Eq. 8, we use the theory
of joint cumulants (Kendall et al., 1977, Section 3.14,
Section 3.29). Since E(Yi) = E(Yj) = 0, we get
∂ 2
∂λ 21
Ki j(λ1,λ2)
∣∣∣
λ1=λ2=0
= E(Y 2i )
∂ 2
∂λ1∂λ2
Ki, j(λ1,λ2)
∣∣∣
λ1=λ2=0
= E(YiYj)
∂ 4
∂λ 21 ∂λ
2
2
Ki, j(λ1,λ2)
∣∣∣
λ1=λ2=0
= E(Y 2i Y 2j )−E(Y 2i )E(Y 2j )−2(E(YiYj))2 .
Combining these equations with Eq. 16, we find
E(Y 2i ) = κ2(Z′)
∫
R3
gi(u)2 du ,
E(Y 2i Y 2j ) = κ4(Z′)
∫
R3
gi(u)2g j(u)2 du
+κ2(Z′)2
[∫
R3
gi(u)2 du
]2
+2κ2(Z′)2
[∫
R3
gi(u)g j(u)du
]2
.
Inserting this into Eq. 7, we obtain Eq. 4.
APPENDIX B: DERIVATIONS
FOR THE GAUSSIAN KERNEL
FUNCTION
In order to derive Eqs. 11 and 12, we use that for
v ∈ R3∫
R3
f (u) f (u+ v)du =
1
(4πσ2)3/2
e−
1
4σ2
‖v‖2
and∫
R3
f (u)2 f (u+ v)2 du =
1
8(2πσ2)9/2
e−
1
2σ2
‖v‖2
.
Furthermore, for v1,v2 ∈R3 such that v1 ⊥ v2, we have∫
R3
f (u)2 f (u+ v1) f (u+ v1 + v2)du
=
1
8(2πσ2)9/2
e−
1
2σ2
‖v1‖2− 38σ2 ‖v2‖
2
and ∫
R3
f (u) f (u+ v1) f (u+ v2) f (u+ v1 + v2)du
=
1
8(2πσ2)9/2
e−
1
2σ2
(
‖v1‖2+‖v2‖2
)
.
Using these expressions in Eqs. 5 and 6 yields Eqs. 11
and 12.
49