arXiv:1203.5010v2 [] 30 Nov 2016

Non-Parametric Approximations for
Anisotropy Estimation in Two-dimensional
Differentiable Gaussian Random Fields
Manolis P. Petrakis∗ and Dionissios T. Hristopulos†
Geostatistics Laboratory, School of Mineral Resources
Engineering, Technical University of Crete, Chania 73100,
Greece
Dated: March 28, 2018

Abstract
Spatially referenced data often have autocovariance functions with
elliptical isolevel contours, a property known as geometric anisotropy.
The anisotropy parameters include the tilt of the ellipse (orientation
angle) with respect to a reference axis and the aspect ratio of the
principal correlation lengths. Since these parameters are unknown a
priori, sample estimates are needed to define suitable spatial models
for the interpolation of incomplete data. The distribution of the
anisotropy statistics is determined by a non-Gaussian sampling joint
probability density. By means of analytical calculations, we derive
an explicit expression for the joint probability density function of
the anisotropy statistics for Gaussian, stationary and differentiable
random fields. Based on this expression, we obtain an approximate joint
density which we use to formulate a statistical test for isotropy. The
approximate joint density is independent of the autocovariance function
and provides conservative probability and confidence regions for the
∗
†

petrakis@mred.tuc.gr
Corresponding author; dionisi@mred.tuc.gr

1

anisotropy parameters. We validate the theoretical analysis by means
of simulations using synthetic data, and we illustrate the detection of
anisotropy changes with a case study involving background radiation
exposure data. The approximate joint density provides (i) a standalone approximate estimate of the anisotropy statistics distribution (ii)
informed initial values for maximum likelihood estimation, and (iii) a
useful prior for Bayesian anisotropy inference.

1

Introduction

Fast and accurate methods of anisotropy estimation are needed in various
fields to better model spatially extended processes and the properties of
heterogeneous materials (Guilleminot and Soize, 2012). The characterization and measurement of anisotropy in biological tissues, for example, is
important for diagnostic and medical reasons (Ranganathan et al., 2011;
Richard and Bierme, 2010). Significant changes in anisotropy over time may
suggest a crucial change in the underlying physical processes. For example,
an accidental release of radioactivity may significantly alter the anisotropy of
radioactivity patterns over the monitored area. Reliable and computationally fast detection of systematic changes in spatial distributions is crucial,
especially for automatic monitoring systems (Pebesma et al., 2011). Another
practical question is what constitutes a significant departure from isotropy to
necessitate the use of anisotropic autocovariance functions. Non-parametric
methods attempt to provide answers to such questions without requiring
knowledge of the autocovariance functions (henceforward, covariance function
for simplicity). Non-parametric isotropy tests are thoroughly reviewed in a
recent publication (Weller and Hoeting, 2015).
Two types of anisotropy are usually encountered in spatially extended
processes. Physical anisotropy implies tensor fields that represent directionally
dependent material properties such as transport coefficients in heterogeneous
media, e.g. (Adler, 1992). Statistical anisotropy characterizes scalar processes
(e.g., scalar permeability, pollutant concentrations), the correlation range of
which depends on the spatial direction. Geostatistical analysis employs two
types of statistical anisotropy: geometric and zonal (Zimmerman, 1993; Chilès
and Delfiner, 2012). Herein we focus on geometric anisotropy, which implies
SRFs with covariance functions that possess elliptical isolevel contours (see
Fig. 1). The estimation of anisotropy parameters is a topic of ongoing interest

2

in various engineering fields (Jiang, 2005; Okada et al., 2005; Feng et al., 2008;
Olhede, 2008; Le Bihan et al., 2001; Xu and Choi, 2009; Richard and Bierme,
2010; Wang and Leckie, 2012) and in data assimilation (Weaver and Mirouze,
2013). In geostatistics, the anisotropy is typically modeled by estimating the
empirical variogram in different directions and fitting anisotropic variogram
models (Chilès and Delfiner, 2012). For second-order stationary SRFs the
variogram is equivalent to the covariance function. However, the interpretation of such variogram analysis is not always straightforward (Weller and
Hoeting, 2015). Anisotropic modeling in the Bayesian framework has also
been investigated (Ecker and Gelfand, 1999, 2003). Recently, there is interest
in anisotropic models with locally varying parameters (Lillah and Boisvert,
2015). A study focusing on general characterizations of anisotropy beyond
the geometrical model appears in (Allard et al., 2015).
The mathematical framework for the study of anisotropy in spatial processes is based on spatial random fields (SRFs), also known as spatial random
functions (Adler, 1981; Christakos, 1992; Wackernagel, 1997; Lantuéjoul, 2002).
SRFs are used in several scientific and engineering disciplines that study spatially distributed processes (e.g., image processing, theory of transport in
heterogeneous media, wave propagation in random media, environmental
modeling). SRFs with Gaussian joint probability density function also provide the mathematical framework of Gaussian processes in machine learning.
Spatially referenced data are typically modeled as SRFs. The analysis of SRFs
based on data involves a number of distributional assumptions that need to
be validated. A common assumption is that of statistical stationarity which
states that the statistical properties are independent of the position. The less
strict second-order stationarity assumption is used in practice and requires the
expectation of the field to be constant and the covariance function to depend
only on the spatial lag. In the case of Gaussian random fields, second-order
stationarity is equivalent to strong stationarity. Isotropy is a stricter assumption that requires the covariance function to depend only on the magnitude
but not on the direction of the lag. For convenience, isotropic SRF models
are often used, even though many real data sets display anisotropic patterns.
In the case of two-dimensional SRFs that admit first-order derivatives
in the mean-square sense, a non-parametric and non-iterative method for
semi-analytic estimation of anisotropy parameters was proposed and studied
in (Hristopulos, 2002; Chorti and Hristopulos, 2008). This manuscript extends
the works above by investigating the joint dependence of the anisotropy
parameter estimates. We derive a non-parametric approximation of the
3

sampling joint probability density function (JPDF) of anisotropy statistics
for differentiable, stationary Gaussian SRFs. We prove this expression using
the Covariance Hessian Identity (CHI) (Swerling, 1962), the Central Limit
Theorem, Jacobi’s multivariate transformation theorem, and perturbation
analysis.
The non-parametric approximation yields a sampling JPDF which is
more dispersed in parameter space than the exact JPDF. This implies wider
probability regions for the anisotropy parameter statistics and confidence
regions for the anisotropy population parameters. Hence, if a sample is
classified as isotropic at confidence level p based on the approximate JPDF, it
is actually isotropic at p0 > p. The JPDF that we derive can also be used as
a prior in Bayesian model inference (Ecker and Gelfand, 1999, 2003; Schmidt
and O’Hagan, 2003; Zhang, 2012) or as a preliminary step in copula-based
spatial analysis (Kazianka, 2013).
This manuscript is structured as follows: In Section 2 we present essential
definitions and an overview of CHI. In Section 3 we derive a general expression
for the joint probability density fR̂,θ̂ (R, θ) for the anisotropy statistics (R̂, θ̂).
In addition, we obtain a relation for p-level probability regions of the anisotropy
parameters. In Section 4, we derive the non-parametric approximation of
fR̂,θ̂ (R, θ) and the corresponding probability region expression. In Section 5
we formulate a non-parametric test for isotropy. In Section 6, we validate the
theoretical results with numerical simulations and we illustrate the detection
of anisotropy changes with a case study involving radiation exposure data.
Finally, in Section 7 we review the main results obtained in this work, we
present our conclusions, and we outline directions for future research. Proofs
of theorems and lemmas are given in the Appendices.

2

Preliminaries

We use boldface symbols for vectors, matrices and tensors; the superscript
“t” denotes the vector or matrix transpose. D ⊂ R2 denotes the spatial
domain, |D| the enclosed area, s ∈ D the position vector in D, and ksk the
Euclidean norm of s. X(s, ω) represents a scalar SRF on the probability space
(Ω, F, P). The state index ω determines the field state and is suppressed in
the following for the sake of brevity. The events in F comprise the measured
SRF realization(s) or sample state(s). E [·] denotes the expectation over
the ensemble of states, and Cov (Z1 , Z2 ) = E [Z1 Z2 ] − E [Z1 ] E [Z2 ] is the
4

covariance of the random variables Z1 and Z2 . Realizations of an SRF X(s)
will be denoted by x(s).
We focus on wide-sense stationary Gaussian SRFs (GSRFs) with constant
mean m = E [X(s)] and covariance function c(r) = E [X(s) X(s + r)] − m2 .
We assume that the SRF is first-order differentiable in the mean square
sense for every s ∈ D, so that the partial derivatives ∂ 2 c(r)/∂ri2 in the
orthogonal directions i = 1, 2 exist at r = (0, 0). For Gaussian SRFs, mean
square differentiability essentially implies that the respective derivatives of
the sample states exist almost surely (Adler, 1981; Yaglom, 1987).
We assume
R
short-range correlations, i.e., with a finite correlation area dr |c(r)|. Such
correlation functions have a finite integral range.
The sample, xk = (x1 , . . . , xN )t comprises the values xk = x(sk ) of the
realization x(s), where sk , k = 1, . . . , N are sampling locations. We use the
following notation for the anisotropic parameters, illustrated in terms of the
anisotropic ratio: population parameters are marked by a star, i.e., R∗ . The
sampling function of R∗ is the random variable R̂. Specific numerical values
will be denoted by R. Sampling functions based on discrete approximations
of derivatives are denoted by Ř. The population anisotropy parameters are
illustrated in Fig. 1.
The Covariance Hessian Matrix H(r) (CHM) of a stationary, at least
first-order differentiable, SRF X(s) is defined as follows
∂ 2 c(r)
.
Hij (r) = −
,
∂ri ∂rj

i, j = 1, 2.

(1)

Let Xij (s) = ∂i X(s) ∂j X(s), i = 1, 2 be the gradient tensor, where
∂i X(s) = ∂i X(s)/∂si , i = 1, 2 are the partial derivatives of X(s). The mean
gradient tensor Q∗ , also known as the matrix of spectral moments (Adler,
1981), is defined as follows
.
Q∗ij = E [∂i X(s) ∂j X(s)] = E [Xij (s)] .
(2)
The matrix Q∗ is nonnegative definite as the covariance of the random gradient
∇X(s) = (∂1 X(s), ∂2 X(s))t . It satisfies the following theorem (Swerling,
1962):
Theorem 1 (Swerling’s CHI). Let X(s) be a statistically stationary SRF
with covariance function c(r) that admits partial derivatives ∂ 2 c(r)/∂ri2 at
r = (0, 0). Then
Q∗ = H(r)|r=0 .
(3)
5

y

A2

c(u) = φ(utVu)

ξ2

A1

u

ξ1
*

θ

x

Figure 1: Definition of elliptical anisotropy parameters. The ellipse denotes
an iso-level contour of an anisotropic covariance function c(·). The principal
correlation lengths are ξ1 , ξ2 ; R∗ = ξ2 /ξ1 is the anisotropy ratio, while the xaxis and A1 are the sides of the anisotropy angle θ∗ . An anisotropic covariance
function c(·) is obtained from a positive definite function ϕ(·) via a rescaling
V followed by a rotation by θ∗ .
To define the anisotropy parameters, consider a coordinate system aligned
with the principal axes of anisotropy, e.g., A1 and A2 (see Fig. 1). In the
principal system, c(u) = φ(ut Vu), where u = (u1 , u2 ) is the lag, V is a
diagonal 2 × 2 matrix, and ϕ(·) is a positive definite function. The principal
correlation lengths of X(s) are given by ξi−2 = −aσ −2 ∂ 2 c(u)/∂u2i |u=0 , for
i = 1, 2 where a is a positive constant (Chorti and Hristopulos, 2008). We
define the anisotropy ratio as R∗ = ξ2 /ξ1 , and the orientation (rotation) angle
θ∗ as the angle between the horizontal axis of the reference system and A1 .
The anisotropy parameters (R∗ , θ∗ ) satisfy the following theorem:
Theorem 2. Let X(s) be an SRF satisfying the conditions of Theorem 1,
and qd , qo represent the following ratios of gradient tensor elements Q∗ij :
1 + R∗ 2 tan2 θ∗
. Q∗22
,
= ∗ = ∗2
Q11
R + tan2 θ∗

(4a)

tan θ∗ (R∗ 2 − 1)
. Q∗
qo∗ = 12
=
.
Q∗11
R∗ 2 + tan2 θ∗

(4b)

qd∗

6

Then, the anisotropic ratio, R∗ and the orientation angle, θ∗ are given by


2qo∗
1
−1
∗
,
(5a)
θ = tan
2
1 − qd∗

−1/2
1 − qd∗
∗
R = 1+ ∗
.
(5b)
qd − (1 + qd∗ ) cos2 θ∗
The proof is based on Theorem 1 (Chorti and Hristopulos, 2008). Therein
the notation R = R2(1) = ξ1 /ξ2 was used, whereas above we defined R∗ = ξ2 /ξ1 .
The equations (4) and (5) follow from (Chorti and Hristopulos, 2008) by means
of the transformation R 7→ 1/R∗ .
Equations (4) are invariant under the pair of transformations tan θ∗ 7→
−(tan θ∗ )−1 , that is, θ∗ 7→ θ∗ ± π/2, and R∗ 7→ 1/R∗ . By restricting the
parameter space to R∗ ∈ [0, ∞) and θ∗ ∈ [−π/4, π/4), or equivalently to R∗ ∈
[1, ∞) and θ∗ ∈ [−π/2, π/2), ensures that the mapping (qd∗ , qo∗ ) 7→ (R∗ , θ∗ ) is
one-to-one, except for the point (1, 0) which maps to (1, θ∗ ) θ∗ being any angle
∈ [−π/2, π/2). Theorem 2 permits estimating the anisotropy parameters
without knowledge of the covariance function, if Q∗ can be estimated from
the data (Hristopulos, 2002; Chorti and Hristopulos, 2008).

3

Sampling Joint PDF of Anisotropy Statistics

Every realization x(s) yields a different estimate of Q∗ , leading to a probability
distribution for the statistics R̂ and θ̂. Below we derive the joint PDF fR̂,θ̂ (R, θ)
based on Jacobi’s theorems for the transformation of a multivariate probability
distribution under transformation of the respective variables (Papoulis and
Pillai, 2002) and the Central Limit Theorem (CLT).
We estimate Q∗ ij using the spatially averaged gradient tensor Q̂ij , where
i, j = 1, 2,
Q̂ij :=

N
N
1 X
1 X
Xij (sk ) =
∂i X(sk ) ∂j X(sk ).
N k=1
N k=1

(6)

The estimation of the field’s partial derivatives from the data is discussed in
Section 6. Replacing the expectation with the spatial average requires the
ergodic hypothesis. A necessary condition for ergodicity is that |D| → ∞ in
7

such a way that both ratios of domain length over the correlation length in
the respective direction tend to infinity. In practice, this means that for an
accurate estimate of Q∗ ij the domain length along each principal direction
should be considerably larger than the respective correlation length. In the
following, we assume that the asymptotic regime is defined by |D| → ∞ in
the sense defined above for ergodicity and N → ∞ (for application of the
CLT).

3.1

Joint PDF of Gradient Tensor Components

We define the following random vector
Q̂ = (Q̂11 , Q̂22 , Q̂12 )t =

N
N
N
1 X
1 X
1 X
X11 (sk ),
X22 (sk ),
X12 (sk )
N k=1
N k=1
N k=1

!t
,

(7)
that comprises the independent components of the fully symmetric gradient
tensor sampling function (Q̂12 = Q̂21 ). As we show below, Q̂ tends to follow
the joint Gaussian distribution in the asymptotic limit due to the Central
limit theorem.
PN
1
According to (6),
h Q̂
i ij = N k=1 Xij (sk ) and based on the definition (2)

it follows that E Q̂ij = Q∗ij , i.e., Q̂ is an unbiased estimator of Q∗ . By
definition, the covariance matrix CQ̂ is symmetric, namely Cij;kl = Ckl,ij ;
hence, it involves six independent entries.
Lemma 1 (Covariance matrix CQ̂ ). For a statistically stationary GSRF, the
six independent entries of CQ̂ are given by the following series
1 X
1 X
C
(r
)
=
[Hik (rnm )Hjl (rnm ) + Hil (rnm )Hjk (rnm )]
ij;kl
nm
N2 r
N2 r
nm
nm

1  ∗ ∗
=
Qik Qjl + Q∗il Q∗jk
N
1 X
+ 2
[Hik (rnm )Hjl (rnm ) + Hil (rnm )Hjk (rnm )] ,
(8)
N r 6=0

Cij;kl =

nm

for
(i, j, k, l) ∈ {(1, 1, 1, 1), (1, 1, 2, 2), (1, 2, 1, 2), (1, 1, 1, 2), (2, 2, 2, 2), (1, 2, 2, 2)}
8

where rnm = sn − sm is the lag vector between two locations sn and sm for
n, m = 1, . . . , N .
Proof. The proof is given in Appendix A. This is the only step in which we
employ the Gaussian assumption for the joint PDF in order to accomplish
the decomposition of higher than second-order moments based on the WickIsserlis theorem. However, the Gaussian assumption could be relaxed using a
variational Gaussian approximation.

The term N −1 Q∗ik Q∗jl + Q∗il Q∗jk in (8) is obtained from the summands
with rnm = 0 and leads to the non-parametric approximation of fR̂,θ̂ (R, θ)
as shown below. The sums over rnm 6= 0 include parametric corrections
that depend on the covariance function. In the approximate, non-parametric
expression we omit the parametric terms which are smaller. These terms have
an 1/N 2 prefactor, but they also involve N 2 summands. However, the products
of the covariance Hessians Hik (rnm ) Hjl (rnm ) that appear in the summands
decay very fast with krnm k. This is due to the fact that, according to (1), the
covariance Hessian decays in space proportionally to the second derivative of
the covariance function; assuming ergodic conditions, this decay is fast. Hence,
at most O(N ) of these terms, for which krnm k < min(ξ1 , ξ2 ), make a significant
contribution. Thus, the parametric correction is at most O(1/N ). On the
other hand, since it involves Hik (rnm ) at finite lag distances, the corrections
are smaller (in absolute value) than the non-parametric component. In the
isotropic case, H12 (rnm ) = 0 for every rnm .
Lemma 2 (Joint PDF of Q̂). Assume X(s) is a statistically stationary SRF
with short-ranged covariance c(r) whose spectral density satisfies C̃(k) ∼
O(kkk−3− ) for  > 0 as kkk → ∞. Then, the joint PDF of the vector
Q̂ which is defined by (7) tends asymptotically to the following trivariate
Gaussian
1
−1
∗ t
∗
e− 2 (Q−Q ) CQ̂ (Q−Q )
∗
fQ̂ (Q̂; Q , CQ̂ ) =
,
(9)
(2π)3/2 det(CQ̂ )1/2
h i
where E Q̂ = Q∗ and the covariance matrix CQ̂ is defined by (8).
Proof. The proof is given in the Appendix B. The condition C̃(k) ∼
O(kkk−3− ),  > 0 implies that for every kkk → ∞, there are  > 0 and
C∞ > 0, such that C̃(k) ≤ C∞ /kkk3+ . This is satisfied by most finite-range,
twice differentiable covariance functions, including the Gaussian, rational
9

quadratic, Bessel-J, and Matérn with ν > 1 covariance models (Lantuéjoul,
2002).

3.2

PDF of Gradient Tensor Ratios

Based on the joint PDF of Q̂, we derive the JPDF of the gradient tensor
ratios fq̂ (q; Q∗ , CQ̂ ), where q = (qd , qo )t .
Lemma 3 (PDF of gradient tensor ratios). For an SRF X(s) that satisfies the
conditions of Lemma 2, the joint density fq̂ (q; Q∗ , CQ̂ ) tends asymptotically
to the following non-Gaussian density
λ1

λ 2 e− 2
fq̂ (q; Q , CQ̂ ) =
8z15



√
2π (z22 +4z12 ) exp

z22
8z12



z2
erfc √
−4z1 z2 ,
2 2z1
(10)
where erfc(·) is the complementary error function, and z1 , z2 , λ1 , λ2 in the
above expression are given by the following expressions, where q0 t = (1, qd , qo )
∗







t

z12 (q; CQ̂ ) = q0 C−1
q0 ,
Q̂
∗t

(11a)

z2 (q; Q∗ , CQ̂ ) = −2 Q C−1
q0 ,
Q̂

(11b)

λ1 (Q∗ , CQ̂ ) = Q∗ t C−1
Q∗ ,
Q̂

(11c)

λ2 (CQ̂ ) = (2π)−3/2 [det(CQ̂ )]−1/2 .

(11d)

Proof. The proof is based on the transformation of the JPDF under the
change of variables Q̂ 7→ q̂ and is given in Appendix C.
We simplify (10) by explicitly showing the dependence of fq̂ (q; Q∗ , CQ̂ )
on N . First, note that as shown by (8) and the associated dimensional
analysis, CQ̂ ∝ O(1/N ). Since z12 > 0 for all correlated SRFs, we can define
√
√
y N = z2 /(2 2z1 ) and 2 N λ̃1 = λ1 . In light of z1 > 0 according to (11a)
and z2 < 0 according to (11b), it follows that y < 0. The JPDF is expressed
as follows in terms of y
fq̂ (q; Q∗ , CQ̂ ) =

 √ 
√ i
λ2 e−N λ̃1 h√
√ 3
π (2y 2 N +1) exp(y 2 N ) erfc y N −2 y N .
2z1
(12)
10

√
For y < 0 and N → ∞ we define x = y N , we use the identity erfc(x) =
2 − erfc(−x) and the asymptotic expansion of the complementary error
function (Abramowitz and Stegun, 1970, Eq. 7.1.23 and 7.1.24) to show that

2 
erfc(x) = 2 + e−x π −1/2 x−1 + O(x−2 ) .
Hence, to leading-order in N , the JPDF (12) is approximated as follows
√
h
i
2π λ2
2
2
∗
fq̂ (q; Q , CQ̂ ) ≈
(2y N + 1) exp (y − λ̃1 ) N .
(13)
z13

3.3

Joint PDF of Anisotropy Statistics

Theorem 3 (Joint PDF of anisotropy statistics). For an SRF X(s) that
satisfies the conditions of Lemma 2, the JPDF of the statistics R̂ and θ̂ is
given asymptotically by
2R |R2 − 1| fq̂ (q; Q∗ , CQ̂ )
fR̂,θ̂ (R, θ; Q , CQ̂ ) ≈
3 ,
R2 cos2 θ + sin2 θ
∗

(14)

where fq̂ (q; Q∗ , CQ̂ ) is given by (12).
Proof. The proof is given in Appendix D. It is based on the transformation
of the multivariate probability density function under the change of variables
q 7→ (R, θ)t .
The function fR̂,θ̂ (R, θ) is clearly non-Gaussian and depends on Q∗ and
CQ̂ , whereas q is expressed in terms of (R, θ) using (4a) and (4b). If the
rotation angle is measured in degrees instead of radians, fR̂,θ̂ (R, θ) should be
multiplied by π/180.

3.4

Probability Regions for Anisotropy Parameters

The probability region at a probability level p ∈ [0, 1] is the “volume” of
space which contains a proportion p of the anisotropy statistics, given the
true values (R∗ , θ∗ ). The probability region of the anisotropy parameters is
defined by the following equivalent equations
R
∗

RE dQ fQ̂ (Q11 , Q22 , Q12 ; Q , CQ̂ ),
p=
dqd dqo fq̂ (qd , qo ; Q∗ , CQ̂ ),
C0

R

dR dθ fR̂,θ̂ (R, θ; Q∗ , CQ̂ ),
C 00
11

Figure 2: Schematic illustrating the transformation of probability regions at
level p due to variable transformations.
where E ⊂ R3 , C 0 ⊂ R2 , and C 00 ⊂ [0, ∞) × [−π/4, π/4) represent the probability regions in the respective spaces. These equations represent the evolution
of the probability region under the variable transformations Q 7→ q 7→ (R, θ)
as shown schematically in Fig. 2.
Lemma 4 (Parametric equation of probability regions). For a SRF X(s) that
satisfies the conditions of Lemma 2, the probability region of the anisotropy
statistics corresponding to level p ∈ [0, 1] in (R, θ)-space is given by the
parametric equation
h
i
y 2 (q; Q∗ , CQ̂ ) − λ̃1 (q; Q∗ , CQ̂ ) N = ln(1 − p),
(15)
where y(·), λ̃1 (·) are defined in (12) and q 7→ (R, θ) by means of (4a)
and (4b).
Proof. The JPDF fQ̂ is given by the trivariate Gaussian (9). Hence, the
probability region of Q̂ is an ellipsoid whose surface satisfies the equation
(Q − Q∗ )t C−1
(Q − Q∗ ) = `p ,
Q̂

(16)

where `p = F −1 (χ2 = p, ν = 2) is the inverse of the chi-square cumulative
distribution function with ν = 2 degrees of freedom (Siotani, 1964). Under
the transformation Q 7→ q, the ellipsoid is projected onto an ellipse which is
deformed by the transformation q 7→ (R, θ) into an asymmetric convex curve
12

(see Fig. 2). Based on (C-2), the equation of the corresponding ellipsoid in
(u, qd , qo )-space is given by
z12 (q; CQ̂ ) u2 + z2 (q; Q∗ , CQ̂ ) u + λ1 (Q∗ , CQ̂ ) − `p = 0,
where the coefficients z1 (·), z2 (·), λ1 (·) are given by (11). The above quadratic
equation has a unique real solution u = Q11 for any q if the discriminant
vanishes, i.e.,


z22 (q; Q∗ , CQ̂ ) − 4z12 (q; CQ̂ ) λ1 (Q∗ , CQ̂ ) − `p = 0.
(17)
The equation above defines the probability region at level p. We can verify
using (11) that (17) represents an ellipse in the space of q, i.e., it is equivalent
to
t

q0 M q0 = 0, where
M = (C−1
Q∗ ) (C−1
Q∗ )t − (Q∗ t C−1
Q∗ − `p ) C−1
.
Q̂
Q̂
Q̂
Q̂
By inserting in (17) the functions y and λ̃1 used in (12) we obtain the
parametric equation 2(λ̃1 − y 2 ) N = `p , where by definition F (`p , ν = 2) = p.
Since F (x, ν = 2) = 1 − exp(−x/2) (Abramowitz and Stegun, 1970, Eq.
26.4.1), it follows that `p = −2 ln(1 − p), finally leading to (15).

4

Non-parametric JPDF and Probability Region

The expressions for fR̂,θ̂ (R, θ; Q∗ , CQ̂ ) and the probability regions of (R̂, θ̂)
above depend on the matrix CQ̂ , given by (8). CQ̂ involves the series (8) that
does not, in general, admit a closed form. If Hij (r) decays fast for increasing
(0)
krk we can use the explicit approximation CQ̂ ≈ CQ̂ , where
Q∗ 2
Q∗12 2
2  ∗11 2
Q12
Q∗22 2
=
N
Q∗11 Q∗12 Q∗12 Q∗22


(0)

CQ̂


Q∗11 Q∗12
.
Q∗12 Q∗22
2
1
(Q∗12 + Q∗11 Q∗22 )
2

(18)

Figure 3 illustrates this fast decay of Cij;kl (r) for isotropic (Fig. 3a) and
(0)
anisotropic Gaussian covariance (Fig. 3b) functions. We expect that CQ̂
13

(a) Isotropic: ξ = 1, σ 2 = 1.

(b) Anisotropic: ξ1 = 1, ξ2 = 2, θ = 30◦ , σ 2 = 1.

Figure 3: Plots of Cij;kl (r) for (a) isotropic (ξ = 1) and (b) anisotropic (ξ1 = 1,
ξ2 = 2 and θ = 30◦ ) Gaussian correlation functions. Cij;kl (r) essentially
vanishes outside a square of side a ≈ 2ξ (isotropic case) and a ≈ 3 max(ξ1 , ξ2 )
(anisotropic case).
14

will lead to a joint PDF with higher uncertainty, and hence more spread
out than the true PDF, because it does not incorporate spatial correlations.
We validated this intuitive argument by means of numerical simulations (see
Section 6.1).
Theorem 4 (Non-parametric JPDF). For an SRF X(s) that satisfies the
(0)
conditions of Lemma 2, the non-parametric JPDF approximation fR̂,θ̂ (R, θ;
R∗ , θ∗ , N ) of (R̂, θ̂) is given in the asymptotic regime by
(0)

(0)

fR̂,θ̂ (R, θ; R∗ , θ∗ , N ) = |det(Jθ,R )| fq̂ (R, θ; R∗ , θ∗ , N ),

(19a)

where
√
(0)
fq̂ (R, θ; R∗ , θ∗ , N )

≈


2π λ2;0
2
2y02 N + 1 eN (y0 −1/2) .
3
z1;0

The coefficients z1;0 , y0 , λ2;0 are given by the following expressions,
δθ = θ − θ∗ ,
 2

1
y0 = √
(R − 1)(R∗ 2 − 1) cos(2δθ) − (R2 + 1)(R∗ 2 + 1) ,
2 z1;0
√
3
2 (R∗ )3  2
2
,
(R
+
1)
−
(R
−
1)
cos(2θ)
λ2;0 =
π 3/2
2
z1;0
= (R2 − 1)2 (R∗ 2 − 1)2 cos(4δθ) − 4(R4 − 1)(R∗ 4 − 1) cos(2δθ)

(19b)
where

(20a)
(20b)
(20c)

+ (R4 + 1)(3R∗ 4 + 2R∗ 2 + 3) + 2R2 (R∗ 2 − 1)2 .
(0)

Proof. In (11) we replace CQ̂ with CQ̂ , defined by (18). Thus, z1 , z2 , λ1 , λ2
√
are replaced, respectively, by z1;0 , z2;0 , λ1;0 , λ2;0 ; then, y0 = z2;0 /(2 N z1;0 ).
(0)
Performing the calculations with CQ̂ we obtain (20a)-(20c). The asymptotic
result (13) of Lemma 3 is used in (19a) to obtain the non-parametric approximation (19b). Note that in the non-parametric approximation, the coefficient
λ̃1;0 in the exponent on the right hand side of (19b) is reduced to 1/2.
Numerical comparisons show that the absolute relative error between
(0)
the non-parametric JPDF fR̂,θ̂ (R, θ; R∗ , θ∗ , N ) calculated with (i) the exact
(0)

fq̂ (R, θ; R∗ , θ∗ , N ), obtained from (12) by inserting the approximate covari(0)

ance matrix CQ̂ , and (ii) the asymptotic limit given by (13), is less than
≈ 10−9 for N = 50 and ≈ 10−6 for N = 30.
15

(a) R∗ = 1.2, θ∗ = 20◦ , N = 100.

(b) R∗ = 1.2, θ∗ = 20◦ , N = 500.

(c) R∗ = 3, θ∗ = 10◦ , N = 100.

(d) R∗ = 3, θ∗ = 10◦ , N = 500.
(0)

Figure 4: Non-parametric JPDF fR̂,θ̂ (R, θ) for various anisotropy parameters
R∗ , θ∗ and sample size N .
Figure 4 demonstrates representative plots of the non-parametric JPDF
based on (19). Note the bimodal structure of the JPDF for N = 100 in Fig. 4a,
with one mode at R = 1.2 and the other (smaller) at R ≈ 0.8. This is due to
the considerable spread of θ̂, which results from the relatively small number
of sampling points and the degeneracy of the anisotropy vector, i.e., the fact
that the combination (R, θ) is equivalent to (1/R, θ − π/2); the degenerate
peak at (0.83, −70◦ ) is folded into the primary domain. On the other hand,
the smaller dispersion of θ̂ for R = 3 leads to a single mode even for N = 100.
Knowledge of the anisotropy JPDF allows the construction of probability
regions for the anisotropy statistics and confidence regions for the anisotropy
parameters. For an SRF X(s) that satisfies the conditions of Lemma 2, the
probability region corresponding to level p of the anisotropy statistics (based
on the simplifications of the non-parametric approximation), is given by the
16

following slight modification of (15)
y02 (R, θ; R∗ , θ∗ ) −

1
ln(1 − p)
=
,
2
N

(21)

where y0 is a function of the values R, θ, and the parameters R∗ , θ∗ as defined
in (20a).

5

Statistical Test of Isotropy

Theorem 5 (Isotropic ratio). Let X(s) be a statistically isotropic GSRF
(R∗ = 1) which is sampled at N points. Assume that the covariance c(r) is
short-ranged and its spectral density satisfies C̃(k) ∼ O(kkk−3− ) for  > 0
as kkk → ∞ as defined in Lemma 2. In addition, assume that the asymptotic
regime conditions hold. The probability interval of the anisotropic ratio at
probability level p (for N > 2`p ) is given by
!
p
p
N − 2 `p (N − `p ) N + 2 `p (N − `p )
,
(22)
,
N − 2`p
N − 2`p
where `p = F −1 (χ2 = p, ν = 2) = −2 ln(1 − p) is the inverse of the chi square
cumulative distribution function with two degrees of freedom.
Proof. For R∗ = 1 the angle dependent terms in the equations (20) vanish,
showing explicitly that the probability region is independent of θ. Plugging (20a) in (21) the following quadratic in R2 equation is obtained
N (R2 − 1)2 − 2`p (R4 + 1) = 0.
In fact, the probability region is reduced to a one-dimensional probability
interval whose endpoints coincide with the roots of the above equation. The
constraint N > 2`p is in practice satisfied for N → ∞ and ensures that the
roots of the above equation are positive real numbers. Based on the definition
of `p the constraint is equivalent to N > −4 ln(1 − p). For example, p = 0.95
implies `p ≈ 6 and N > 12.
Equation (22) is independent of c(r) and thus provides a non-parametric
approximation of the probability interval for R∗ . The JPDF (19) is indepen(0)
dent of θ and θ∗ for R∗ = 1. The PDF, fR̂ (R), of R̂ for R∗ = 1 and N = 100
17

3

(0)

fR̂ (R)

2.5
2
1.5
1
0.5
0
0.6

0.8

1

1.2

1.4

1.6

R
(0)

Figure 5: Non-parametric sampling PDF of the anisotropic ratio, fR̂ (R), for
an isotropic random field sampled at N = 100 points. Shaded area represents
the corresponding 95% probability interval (R− , R+ ) = (0.77, 1.29).
is shown in Fig. 5, including the 95% probability interval predicted by (22).
Note that the PDF has a node instead of a peak at R = 1. This is not an
artifact of the non-parametric approximation, since the complete JPDF (14)
also vanishes at R = 1. The node is due to the root of the Jacobian (D-1) at
R = 1, which reflects that the isotropic point (1, 0) in (qd , qo )-space is mapped
onto the straight line R = 1 in the (R, θ)-space. The node is also evident in
numerical simulations that do not use the Jacobian (see Figure 7 below).

6

Application to Simulated and Real Data

To apply the formalism developed above to data sets that comprise discrete
sets of values, we replace the partial derivatives by respective discrete operators
ˇ i (sk ), i = 1, 2. The respective estimates of Q∗ are denoted by Q̌. The
∂X
discretization introduces a bias that increases with the sparsity of the sampling
pattern. A “good” sampling pattern is characterized by a typical distance
â between nearest neighbors which is approximately uniform (ideally, a
regular lattice pattern is best) and â  min(ξ1 , ξ2 ), where ξ1 , ξ2 are the
ˇ i (sk )
principal correlation lengths. Different approaches for estimating ∂X
18

are investigated in (Chorti and Hristopulos, 2008). Herein, the centered
differences scheme is used for gridded data.
We denote average values of a statistic over different samples (repetitions)
by a bar over the respective symbol, i.e., Q̌ij . For simulated data, the
ensemble properties Q∗ and CQ̂ which are unknown a priori, are replaced
by the respective averages Q∗ ≈ (Q̌11 , Q̌22 , Q̌12 )t and CQ̂ ≈ CQ̌ . In the
(0)

non-parametric approximation, CQ̂ is obtained from (18) by replacing Q∗
with Q̌.

6.1

Simulated Scattered Data

We generate SRF realizations with specified (R∗ , θ∗ ) to validate the probability
region of the anisotropy parameters (21). Figure 6 and Table 1 investigate the
anisotropic case R∗ = 1.5, θ∗ = −30◦ , whereas the isotropic case is considered
TM
in Table 2 and Figure 7. A desktop computer with an Intel R Core i5-2500
(4 cores, 3.30 GHz) CPU running Matlab R R2015b under 64-bit Windows R
7, was used for all the simulations.
We simulate scattered data using the following method: First, a realization of an GSRF is generated on a regular grid. The Fourier Filtering
Method (Pardo-Igúzquiza and Chica-Olmo, 1993; Lantuéjoul, 2002; Hristopulos, 2005) is used on L × L square grids with lattice constant a = 1. We use
Gaussian, c(r) = σ 2 exp(−krk2 /ξ 2 ), and Matérn, c(r) = σ 2 21−ν Γ(ν)−1 ξ −ν ×
krkν Kν (krk/ξ), covariance functions (expressions correspond to the isotropic
case), where Γ(·) is the Gamma and Kν (·) the modified Bessel function of
order ν. In the Gaussian case, the correlation range is controlled by ξ whereas
in the Matérn case by both ξ and ν. The smoothness parameter ν adjusts
the differentiability of the SRF: ν = 1/2 corresponds to the non-differentiable
exponential function and ν → ∞ to the infinitely differentiable Gaussian.
For given ξ, the field is smoother for higher ν. To compensate for this effect
and to compare SRFs of similar spatial variability, we use rescaled correlation lengths ξ˜ = Ad ξ, where Ad is the
√ integral scale factor (Hristopulos
and Žukovič, 2011): In d = 2, A√
d = 2 πν for Matérn correlations whereas
for Gaussian correlations Ad = π. For equal rescaled correlation
lengths,
√
ξ˜Gauss = ξ˜Matérn , with d = 2, ν = 2, it follows that ξGauss = 2 2 ξMatérn .
We randomly choose a fraction of the grid points to mimic scattered data.
For a square lattice of side L a sample of N = (ρL)2 points are randomly

19

chosen from Gaussian and Matérn lattice SRFs. An estimate
of the mean
√
distance between N uniformly distributed points is L/ N = 1/ρ, thus ρ is
the mean sampling frequency. The samples respect the condition that the
correlation lengths exceed the mean distance between the points, as specified
in the first paragraph of this Section.
We employ the natural neighbor interpolation method (Fisher et al., 2005)
in Matlab R on an M × M square grid with M = 200. Natural neighbor
interpolation provides smooth surfaces and does not assume isotropy of the
data; however, it is defined only inside the convex hull of the data sites. Due
to the occasionally poor sampling near the domain boundaries, interpolation
artifacts appear (Bobach et al., 2009) as elongated islands, oriented vertically
along the left and horizontally along the bottom sides of the domain. Hence,
they tend to bias the anisotropy estimates towards higher or lower anisotropy
√
ratios and angles near zero. Thus, boundary strips of thickness L/ N
are discarded from the interpolation surface to minimize bias. The partial
derivatives are estimated via centered differences on the interpolated surface.
Finally, we perform anisotropy estimation for each sample and compute
the non-parametric probability region at p = 0.95 using ensemble averages.
Also, we compute confidence regions for each anisotropy estimate at several
confidence levels.
6.1.1

Anisotropic Scattered Data

Figure 6a demonstrates a realization of a zero-mean, unit variance anisotropic
GSRF with Gaussian covariance with R∗ = 1.5, θ∗ = −30◦ , ξ = 28.3 on a
600 × 600 grid. A randomly extracted set of N = 1296 points (ρ = 0.06)
is shown in Fig. 6b. The depicted smooth field is generated from the 1296
points by interpolation and is used to estimate (Ř, θ̌).
In Figure 6c the non-parametric probability region (red contour) at p =
0.95, defined by (21) is compared with CHI anisotropy estimates (blue crosses)
from 1000 SRF samples. For each sample, we estimate Q∗ by means of the
spatial average Q̌ and then calculate (Ř, θ̌) by applying (5). We estimate
(R∗ , θ∗ ) based on the Ř, θ̌, which are obtained from the ensemble average Q̌
by means of Theorem 2. The ensemble-based anisotropy estimate (Ř, θ̌), is
denoted by a small circle inside the cloud of the (Ř, θ̌) points. Figures 6d–6f
demonstrate the simulated scattered data probability region estimation for a
zero-mean, unit-variance Matérn covariance with R∗ = 1.5, θ∗ = −30◦ , ξ = 10,
ν = 2. The normality of Q̂, supported by CLT considerations as shown in
20

(a) Gaussian SRF, R∗ = 1.5, (b) 1296 random nodes and (c) Anisotropy estimates and
θ∗ = −30◦
interpolated field.
probability region

(d) Matérn SRF, R∗ = 1.5, (e) 1296 random nodes and (f) Anisotropy estimates and
θ∗ = −30◦
interpolated field.
probability region

Figure 6: Non-parametric probability region estimation for scattered data.
(a) and (d): Realization of zero-mean, unit-variance anisotropic Gaussian SRF
with ξ = 28.3 and Mátern SRF with ν = 2 and ξ = 10 on a 600 × 600 square
grid. (b) and (e): Random sample of N = 1296 points and interpolated field√on
a 200 × 200 grid using natural neighbors. Boundary strips of thickness L/ N
were discarded from the interpolated field to avoid interpolation artifacts
at the domain boundary. (c) and (f): Anisotropy estimates (crosses) are
generated from 1000 random samples of N = 1296 points; the continuous
curve corresponds to 95% non-parametric probability region calculated with
the ensemble-based anisotropy estimates (Ř, θ̌), which are denoted by a small
circle inside the cloud.

21

Lemma 2, was confirmed by normal probability plots (not shown here).
The non-parametric probability region (Theorem 4) extends beyond the
region obtained from the true JPDF (this is supported by Figures 6c and 6f
as explained below). We conducted numerical experiments (not shown here)
˜ and N to confirm that non-parametric probability
for several values of ξ/a
regions based on (19a) are more extended in parameter space than the regions
˜ → 0, i.e., as the spatial extent of the
based on the true JPDF (12). If ξ/a
correlations is reduced, the scatter cloud of (Ř, θ̌) expands and tends to fill
˜ increases,
the non-parametric probability region. On the other hand, as ξ/a
i.e., for dense sampling of the SRF, the scatter cloud tends to be confined
inside the smaller parametric region. These observations agree with our
earlier statement that the non-parametric approximation contains the true
probability region.
In Table 1 we validate the non-parametric anisotropy confidence region for
simulated scattered Gaussian (ξ = 28.3) and Matérn (ξ = 10, ν = 2) covariance functions with R∗ = 1.5 and θ∗ = −30◦ . We generate 1000 realizations
for different domain sizes (L = 600, 800, 1000, 1200) and mean sampling
frequencies (ρ = 0.04, 0.06) and we enumerate the number of simulations for
which the ensemble means Ř, θ̌ (as estimates of the population means) are
outside the non-parametric confidence region. The latter is computed for
each anisotropy estimate at different confidence levels (p = 0.95, 0.75, 0.68,
0.5, 0.25) using (21). If the true JPDF and the confidence regions of the
anisotropy statistics are known at the p levels above, the average number of
simulations for which the true confidence region does not contain the ensemble
means Ř, θ̌ is 50, 250, 320, 500, and 750 respectively. However, the number
of simulations for which Ř, θ̌ lie outside the non-parametric region (21) is
always less than expected for the true confidence regions. This observation
agrees with the proposition that the non-parametric confidence region (21)
contains the true confidence region.
6.1.2

Isotropic Scattered Data

We numerically validate the isotropy testing procedure by enumerating the
number of anisotropy estimates that fall outside the probability region at
p = 0.95 for 1000 realizations of simulated scattered data in different domain
sizes (L = 600, 800, 1000, 1200), mean sampling frequency (ρ = 0.04, 0.06),
and isotropic covariances (Gaussian with ξ = 28.3, Matérn with ξ = 10, ν = 2).

22

23
1.22, −28.9◦

391
680
1.14, −28.1◦

Nout,p=0.50

Nout,p=0.25

Ř, θ̌

608

209

Nout,p=0.68
289

121

78

150

Nout,p=0.75

4

1296

8

576

600

1.21, −29.8◦

601

324

159

106

10

1024

534

191

68

41

3

2304

1.30, −29.7◦

800

1.22, −30.1◦

592

288

141

83

1

1600

465

155

48

26

0

3600

1.31, −30.1◦

1000

1.16, −29.5◦

667

346

172

130

10

1024

572

257

120

69

3

2304

1.23, −29.8◦

800

1.16, −30.0◦

632

302

138

95

5

1600

515

206

78

42

2

3600

1.24, −30.3◦

1000

(b) Matérn, R = 1.5, θ = −30◦ , ν = 2, ξ = 10,ρ = 0.04, 0.06

Nout,p=0.95

N

L

Ř, θ̌

280
1.28, −29.8◦

380

Nout,p=0.50

102
560

194

Nout,p=0.68

68

3

1296

1.19, −29.0◦

132

Nout,p=0.75

600

658

12

Nout,p=0.95

Nout,p=0.25

576

N

L

(a) Gaussian, R = 1.5, θ = −30◦ , ξ = 28.3, ρ = 0.04, 0.06

1.17, −30.4◦

630

320

170

114

11

2304

450

166

47

20

1

5184

1.25, −30.7◦

543

227

86

51

4

5184

1.32, −30.9◦

1200

1.22, −30.4◦

574

274

128

85

6

2304

1200

Table 1: Numerical validation of the non-parametric anisotropy confidence region for simulated scattered
data. Anisotropy estimates were computed from 1000 random samples of N points from (a) Gaussian and
(b) Matérn lattice SRFs with R = 1.5, θ = −30◦ . For each lattice size L and sampling frequency ρ, Nout,p is
the number of anisotropy estimates for which the non-parametric confidence region computed at p = 0.95,
0.75, 0.68, 0.5, and 0.25 using Theorem 4, does not contain the ensemble mean (Ř, θ̌).

If the true JPDF and the probability regions of the anisotropy statistics are
known at p = 0.95 probability level, on average 50 out of the 1000 simulations
should fall outside the true region. In Table 2, Nout,iso is the number of
estimates that fall outside the p = 0.95 isotropy probability interval (R− , R+ )
using Eq. (22) while Nout is the number of samples that fall outside the
probability region calculated using the ensemble-based anisotropy estimate
(Ř, θ̌) and Theorem 4. The mean time t̄ for anisotropy estimation per
processor core is also shown with an error estimate of one standard deviation.
Figure 7 demonstrates the non-parametric probability regions and the
scatter cloud of anisotropy estimates for scattered data sampled from isotropic
Gaussian and Matérn lattice SRFs of increasing side L. The absence of
estimates near the R = 1 line agrees with the existence of a JPDF node at
R = 1 as discussed in Section 5. For smaller domains the anisotropy estimates
deviate from isotropy.
The computational complexity of natural neighbors interpolation is O((M +
N ) log N ) (Park et al., 2006), where M is the number of the interpolation
points. The complexity of derivative estimation using centered differences is
O(M ). Hence, O(2M ) operations are needed for computing ∂ˇi X(s), i = 1, 2
and O(3M ) operations for Q̂ij . Thus the total computation time tM,N is of
O(5M + (M + N ) log N ), from which we obtain tM 0 ,N 0 /tM,N = (5M 0 + (M 0 +
N 0 ) log N 0 )/(5M + (M + N ) log N ). For M = 2002 , tN =5184 /tN =576 = 1.28.
The time ratio obtained from Table 2a is t̄N =5184 /t̄N =576 = 1.36±0.06 and from
Table 2b is t̄N =5184 /t̄N =576 = 1.31±0.06. For M = 1002 , tN =5184 /tN =576 = 1.53
while the simulation times (average times per anisotropy estimation per
processor core) obtained for 1000 realizations of isotropic Matérn covariance
(not shown here) are t̄N =576 = 36.4 ± 1.5 msec and t̄N =5184 = 58.3 ± 2.6 msec,
giving t̄N =5184 /t̄N =576 = 1.60 ± 0.10.

6.2

Case Study: Radiation Exposure

We study anisotropy in two data sets of daily averages of radioactivity gamma
dose rates over part of the Federal Republic of Germany, which was provided
by the German automatic radioactivity monitoring network for the Spatial
Interpolation Comparison (SIC 2004) exercise (Dubois and Galmarini, 2006).
Dose rates are measured in nanosieverts per hour (nSv/h). The background
data set corresponds to typical radioactivity measurements (≈ 100 nSv/h),
which follow the Gaussian distribution (graph not shown), and thus their
skewness and excess kurtosis coefficients are close to zero. The emergency
24

25

9

1

1

5184

7
1.01, 43.6

11

0.990, −44.4
0.902, 1.11
129 ± 3

Nout

Ř, θ̌

R− , R+

t̄ (msec)

1.01, 42.7

12

12

1024

131 ± 4

130 ± 4

0.934, 1.07 0.926, 1.08

8

15

Nout,iso

1296

576

600

N

L

146 ± 5

0.950, 1.05
141 ± 5

0.940, 1.06

5

8

2304

143 ± 5

0.950, 1.05

0.993, −42.3

800

1.01, 44.2

2

5

3600

160 ± 5

0.960, 1.04

139 ± 6

156 ± 5

0.940, 1.06 0.960, 1.04

1.00, 43.6

10

10

1600

1000

(b) Matérn, ν = 2, ξ = 10, ρ = 0.04, 0.06

132 ± 3

0.926, 1.08

145 ± 5

0.950, 1.05

0.997, −37.6

12

10

2304

169 ± 7

0.967, 1.03

1.00, 42.8

4

7

5184

175 ± 6

0.967, 1.03

1200

150 ± 9

0.950, 1.05

131 ± 3

0

11

2304

129 ± 4

5

0

3600

t̄ (msec)

2

6

1600

0.934, 1.07

13

2

2304

1200

0.902, 1.11

3

12

1024

1000

R− , R+

12

Nout

4

1296

800

0.990, −42.2 0.986, −43.2 0.995, −41.6 0.992, −42.2 0.997, −43.0 1.00, 41.8 0.998, −42.5 1.00, 44.9

14

Nout,iso

600

Ř, θ̌

576

N

L

(a) Gaussian, ξ = 28.3, ρ = 0.04, 0.06

Table 2: Numerical validation of the non-parametric isotropy test at p = 0.95 for simulated scattered data.
Anisotropy estimates generated from 1000 random samples of N points from isotropic (a) Gaussian and (b)
Matérn lattice SRFs. For each lattice size L and sampling frequency ρ, Nout,iso is the number of simulations
that fall outside the non-parametric probability region estimated from Q̌ using Theorem 4. Nout is the
number of simulations that fall outside the non-parametric isotropy interval (R− , R+ ) based on (22). t̄:
average time per anisotropy estimation (per processor core) with an error estimate of one standard deviation.

(a) Gaussian, L = 600,
N = 576

(b) Gaussian, L = 800,
N = 1024

(c) Gaussian, L = 1000,
N = 1600

(d) Matérn, L = 600,
N = 576

(e) Matérn, L = 800,
N = 1024

(f) Matérn, L = 1000,
N = 1600

Figure 7: Non-parametric probability region estimation for isotropic scattered
data. The initial lattice SRF is defined over a square lattice with side
L = 600, 800, 1000. The continuous curve corresponds to 95% non-parametric
probability region calculated with anisotropy parameters estimated from
Q̌. (a)–(c) Anisotropy estimates (crosses) generated from 1000 random
samples obtained from a zero-mean, unit-variance isotropic Gaussian SRF
with ξ = 28.3. (d)–(f) Anisotropy estimates (crosses) generated from 1000
random scattered samples obtained from a zero-mean, unit-variance Matérn
SRF with ξ = 10.

26

Table 3: Summary statistics of radioactivity dose rate exhaustive data sets
(units are in nanosieverts per hour) and CHI anisotropy estimates. Abbreviations: min: minimum sample value; med: median sample value, max:
maximum sample value; std: sample standard deviation; skew: sample skewness coefficient; kurt: sample excess kurtosis coefficient; Ř, θ̌: estimates of
anisotropy parameters.
N = 1008

min mean med max

std skew kurt

Background 57.0 97.7 98.6 180.0 19.6 0.4

Ř

θ̌

0.6 1.18 7.36◦

Emergency 57.0 106.1 98.9 1528.2 92.5 11.3 144.1 0.45 −0.75◦
data includes a simulated local release of radioactivity which results in five
dose rate “measurements” around 10 times above background (exceeding
1000 nSv/h). These measurements are aligned in the East-West direction.
Table 3 summarizes the statistics of both data sets. The two rightmost
columns show the CHI-based estimates of anisotropy parameters. Since the
95% confidence interval for isotropy is (R− , R+ ) = (0.92, 1.08), this dataset
can be considered as slightly anisotropic. The direction of anisotropy is
different in the two sets: In the background set the axis A1 is tilted with
respect to the x-axis (which is aligned with the E-W direction) by 7.36◦ ,
while the dominant anisotropy axis is A2 since ξ2 = 1.18ξ1 . A2 is closer to
the y-axis, implying a dominant North-South anisotropy. In the emergency
set the axis A1 is slightly tilted with respect to the x-axis (by −0.75◦ ), and
the dominant anisotropy axis is A1 since ξ2 = 0.45ξ1 . Since A1 is closer to
the x-axis, this implies that the radioactive plume reverses the dominant
anisotropy direction to East-West.
We calculate the non-parametric JPDF and the 95% confidence regions of
the anisotropy statistics based on the estimated anisotropy parameters (i.e.,
by CHI anisotropy estimation on gridded values obtained by natural neighbor
interpolation) for both sets. The results are shown in Fig. 8. There is no
overlap of the two joint density functions, and the contours corresponding
to the 95% confidence regions do not intersect. These patterns suggest
statistically significant anisotropy difference between the background and the
emergency data due to the elevated values of the dose rate in the East-West
direction of the spreading plume which changes the orientation of the major
anisotropy axis.
27

Figure 8: Non-parametric joint PDF shaded surfaces (red and green online)
and 95% confidence regions (solid contours) for the radioactivity dose rate
data sets: background measurements (right) and emergency data simulation
(left).

7

Discussion and Conclusions

This work focuses on the estimation of geometric anisotropy in scattered or
grid-based two-dimensional data. We derive explicit expressions for the joint
PDF of the anisotropy statistics, given by equations (13)-(14), and for the
corresponding anisotropy probability regions at any level, i.e., equation (15).
The main assumptions used are that (i) the data are drawn from a jointly
Gaussian, stationary and differentiable random field and (ii) the covariance
function is short-ranged.
We also derive a non-parametric approximation for the joint PDF of
the anisotropy statistics, which can be used if the covariance function is
unknown a priori, or if estimation of the covariance is not desired. The
non-parametric approximation of the anisotropy joint PDF is given by (19).
The corresponding equation for the non-parametric approximation of the
probability region is given by (21). We also derive probability intervals for
the anisotropy ratio under the hypothesis that the sample comes from an
isotropic random field. These probability intervals are used to formulate a
non-parametric test of the isotropic hypothesis. We illustrate the application
28

of the joint PDF and the probability regions with simulated and real data.
The Gaussian assumption is used in the decomposition of the covariance
(0)
matrices CQ̂ and CQ̂ , i.e., to derive equations (8) and (18) by means of the
Wick-Isserlis theorem. This decomposition can be justified in principle, albeit
approximately, even for non-Gaussian densities, based on optimal variational
(Gaussian) approximations. Practical application of the derived formulas
requires the estimation of anisotropy statistics using CHI. Accurate estimation
based on CHI requires in addition to (i) and (ii) above the following: (iii)
a large sample size, N → ∞ and (iv) a sample domain that is large with
respect to the correlation area. The latter is difficult to satisfy in case of large
anisotropy (R∗  1 or R∗  1). In such cases, the CHI estimate tends to
underestimate the actual anisotropy. CHI anisotropy estimates include biases
due to (v) the finite step size of the grid and (vi) interpolation (in the case of
scattered data).
In summary, our approach consists of the following steps: (i) If necessary,
preprocess the data in order to remove trends and use transforms to reduce
deviations from the Gaussian distribution (e.g., Box-Cox transform). (ii)
Choose an interpolator which provides smooth interpolation surfaces. Several
interpolators were tested and compared in (Chorti and Hristopulos, 2008).
The interpolation grid should be dense to allow accurate approximation of the
spatial derivatives. (iii) Compute the partial derivatives of the interpolated
surface. (iv) Compute Q̂ij and use Theorem 2 to obtain anisotropy parameter
estimates. (v) The non-parametric JPDF is obtained from Theorem 4. The
probability regions are provided by Lemma 4 using the approximate nonparametric covariance matrix, i.e. (18). (vi) The isotropy test of Theorem 5
can be used to test for the presence of anisotropy. (vii) The differentiability
assumption can be tested a posteriori by determining the optimal anisotropic
variogram model using standard geostatistical procedures. In this step, the
CHI anisotropy estimates can be used to fix the anisotropy parameters or to
provide informed initial guesses for likelihood optimization.
Our approach provides a computationally efficient, albeit approximate,
method of geometric anisotropy estimation in two dimensions, because the
analytical expressions derived above can be evaluated with minimal computational cost. The most computationally intensive part is the interpolation of
scattered data onto regular grids in order to calculate derivatives. We use
natural neighbor interpolation which is computationally fast (its complexity
is essentially determined by Voronoi tesselation). For small datasets, the

29

computation time scales linearly with the number of nodes M of the interpolation grid, while for large datasets the computational cost is dominated
by O(M log N ) where N is the number of data points. A formulation of the
natural neighbor interpolation algorithm which directly provides the partial
derivatives of the interpolated surface is also available (Sambridge et al., 1995,
Appendix A1).
Our approach could be useful in estimating anisotropy in big data sets.
In addition, the non-parametric JPDF can be used as an anisotropy prior in
Bayesian and copula analyses (Kazianka, 2013). The method also provides
initial estimates for maximum likelihood estimation of spatial anisotropic
models (Pebesma et al., 2011). Furthermore, it furnishes an easily computable
indicator of physical change in spatially extended systems based on the
comparison of anisotropy probability regions.
Straightforward extension of this work is possible for the joint lognormal
distribution along the lines of (Chorti and Hristopulos, 2008). The global
statistical measures of anisotropy can be efficiently calculated for large domains
and can thus provide a useful statistic for large data sets. Local variations
of anisotropy can also be investigated using windowing methods. Capturing
such local variability has applications in the analysis of medical images,
e.g. (Richard and Bierme, 2010). Currently, the solution of the non-linear
CHI equations for d > 2 is not available in closed form. Hence, an analytical
expression of the anisotropy joint PDF in higher than two dimensions is not yet
feasible. Another path for future research is the development of an anisotropy
detection method which will involve local integrals of the field values. Such
an approach, if analytically tractable, will apply to non-differentiable random
fields as well.

Acknowledgment
This work was funded by the European Commission, under the 6th FP, by
the Contract N. 033811 with the DG INFSO, action line IST-2005-2.5.12 ICT
for Environmental Risk Management. The views expressed herein are those
of the authors and not necessarily of the European Commission.
We would like to thank Prof. Athanasios Liavas (School of Electronic and
Computer Engineering, Technical University of Crete) for a careful reading of
the manuscript and for suggesting improvements. In addition, we thank two
anonymous reviewers for their valuable input.
30

Appendix A: Proof of Lemma 1
Proof. Using the definition (7) we obtain
Cij;kl

N
N
1 X
1 X
= Cov
Xij (sn ),
Xkl (sm )
N n=1
N m=1
1 X
Cov (Xij (sn ), Xkl (sm )) .
= 2
N n,m

!

(A-1)

Due to the stationarity of X(s), the double series in (A-1) is reduced to a
single series over all (N 2 ) lag vectors rnm = sn − sm (n, m = 1, . . . , N ), i.e.,
1 X
Cij;kl = 2
Cov (Xij (s0 ), Xkl (s0 + rnm ))
N r
nm
1
1 X
=
Cov (Xij (0), Xkl (0)) + 2
Cov (Xij (0), Xkl (rnm )) . (A-2)
N
N r 6=0
nm

Covariance of the gradient tensor: Let r denote any lag vector (including
r = 0) between two points. Based on the definition of the covariance function
it follows that
Cov (Xij (0), Xkl (r)) = E [Xij (0)Xkl (r)] − E [Xij (0)] E [Xkl (r)] .

(A-3)

Note that
E [Xij (0)Xkl (r)] = E [∂i X(0) ∂j X(0) ∂k X(r) ∂l X(r)] .
For a differentiable and stationary SRF X(s), the gradient component
∂i X(s) is a zero-mean Gaussian SRF with covariance function given by (Abrahamsen, 1997; Yaglom, 1987)
E [∂i X(s) ∂j X(s + r)] = −

∂ 2 c(r)
.
∂ri ∂rj

(A-4)

Hence, E [Xij (0) Xkl (r)] can be calculated using the moment factorization
property of multivariate normal distributions (Isserlis, 1918; Wick, 1950)
E [Xij (0)Xkl (r)] = E [∂i X(0)∂j X(0)] E [∂k X(r)∂l X(r)]
+ E [∂i X(0)∂k X(r)] E [∂j X(0)∂l X(r)]
+ E [∂i X(0)∂l X(r)] E [∂j X(0)∂k X(r)]
= Hij (0)Hkl (0) + Hik (r)Hjl (r) + Hil (r)Hjk (r).
31

(A-5)

The last equality follows from Eq. (A-4) and the definition (1) of CHM. The
second term on the right-hand side of (A-3) is
E [Xij (0)] E [Xkl (r)] = E [∂i X(0)∂j X(0)] E [∂k X(r)∂l X(r)] = Hij (0)Hkl (0).
(A-6)
Thus, in light of (A-5) and (A-6), equation (A-3) becomes
Cov (Xij (0), Xkl (r)) = Hik (r)Hjl (r) + Hil (r)Hjk (r).

(A-7)

Equation (8) follows from (A-2), (A-7), and Theorem 1 for the zero-lag
CHM.

Appendix B: Proof of Lemma 2
Proof. To show that the JPDF of Q̂ tends asymptotically to the normal
distribution, we use the multivariate CLT theorem. The classical CLT for
scalar random variables is discussed in (Gnedenko and Kolmogorov, 1954;
Levy, 1954; Feller, 1971). The CLT extension to vector random variables is
as follows (Anderson, 1984):
Assume N independent and identically distributed vector variables Zk ,
k = 1, . . . , N with mean m and covariance matrix CZZ . Then, for N → ∞
the joint distribution of the random vector Z̄ = (Z1 + · · · + ZN )/N tends
to the multivariate normal distribution with mean m and covariance matrix
CZZ /N .
The above CLT is generalized to SRF averages. Loosely stated, an average
of a stationary random field with finite-range correlations over N → ∞
points tends to follow the joint normal probability distribution (Bouchaud
and Georges, 1990). Thus, the multivariate CLT applied to the random vector
Zk = (X11 (sk ), X22 (sk ), X12 (sk ))t leads to (9).
Next, we establish the condition for the SRFs to have finite correlation
range. The Xij (sk ) are stationary SRFs by virtue of the stationarity of
X(s). Hence, φijkl (r) := Cov (Xij (s), Xkl (s + r)) = Cov (Xij (0), Xkl (r)). Using (A-7), φijkl (r) = Hik (r)Hjl (r) + Hil (r)Hjk (r). The correlation range of
Xij (sk ) is determined by the integral


Z
1
Vc = max
dr φijkl (r), .
i,j,k,l
φijkl (0)

32

Based on (A-7), φijkl (0) = Qij Qkl + Qil Qjk and thus Rφijkl (0) has a finite
value if X(s) has finite correlation lengths. We calculate D dr φijkl (r) in the
asymptotic regime where |D| → ∞, and we express the integral in terms
of the Fourier transform of c(r). Any permissible covariance function c(r),
where r ∈ R2 , admits the following pair of transformations, where C̃(k) is
the spectral density:
Z
1
c(r) =
dk ek·r C̃(k),
2
(2π)
Z
C̃(k) = dr e−k·r c(r).
R
Based on the above, it follows that Hij (r) = (2π)−2 dk ki kj ek·r C̃(k), and
thus
Z
Z
1
dr φijkl (r) =
dk ki kj kk kl [C̃(k)]2 .
2
(2π)
√
In
k2 r2 is
the inner vector product, and
R the above,
R∞  =
R ∞ −1, k · rR = k1 r1R+
R 2π
∞
dk = −∞ dk1 −∞ dk2 or dk = 0 k dk 0 dφ in polar coordinates. The
existence of the above integral depends on the
R behavior of C̃(k) at kkk = 0
and kkk → ∞. Since c(r) is short-ranged, dr c(r) = C̃(0) is finite, and
thus the integrand is well-behaved at kkk = 0. At kkk → ∞, the integral
converges (using polar coordinates) if [C̃(k)]2 decays asymptotically faster
than kkk−6−2 , where  > 0. This ensures that φijkl (r) is short-ranged.

Appendix C: Proof of Lemma 3
Proof. The probability transformation Q 7→ q is performed as follows: Since
dim(q) = 2 < dim(Q) = 3, we append to q the dummy variable u = Q11 ≥ 0
and then integrate over u. Using definitions (4a) and (4b), the absolute value of
the Jacobian determinant for the transformation (Q11 , Q22 , Q12 ) 7→ (u, qd , qo )
is
∂(Q11 , Q22 , Q12 )
Jq =
⇒ |det(Jq )| = u2 .
∂(u, qd , qo )
The dummy variable u is integrated, leading to
Z ∞
∗
fq̂ (q; Q , CQ̂ ) =
du fQ̂ (u, qo u, qd u; Q∗ , CQ̂ ) u2 .
0

33

(C-1)

In terms of qd and qo , the exponent of the PDF fQ̂ (·), given by (9), becomes
(Q − Q∗ )t C−1
(Q − Q∗ ) = z12 (q; CQ̂ ) u2 + z2 (q; Q∗ , CQ̂ ) u + λ1 (Q∗ , CQ̂ ).
Q̂
(C-2)
By virtue of the above, (C-1) is expressed as follows
Z ∞
1
2
∗
fq̂ (q; Q , CQ̂ ) = λ2
du u2 e− 2 [(u z1 ) +u z2 +λ1 ] .
0

According to (11a), z12 > 0 because CQ̂ is a covariance matrix; hence CQ̂ as
well as C−1
are positive definite. Thus, the Gaussian integral above exists
Q̂
and its value is given by (10).

Appendix D: Proof of Theorem 3
Proof. Equation (14) follows from the transformation (qd , qo ) 7→ (R, θ) with
Jacobian matrix JR,θ . The transformed PDF is given by fR̂,θ̂ (R, θ; Q∗ , CQ̂ ) =
fq̂ (q; Q∗ , CQ̂ ) |det(JR,θ )|, where det(JR,θ ) is given by
det(JR,θ ) =

∂qd
∂R
∂qo
∂R

∂qd
∂θ
∂qo
∂θ

=

2R (R2 − 1)
R2 cos2 θ + sin2 θ

3 .

(D-1)

Restricting the parameter space to R ∈ [0, ∞) and θ ∈ [−π/4, π/4), or
equivalently R ∈ [1, ∞) and θ ∈ [−π/2, π/2), the transformation (qd , qo ) 7→
(R, θ) is one-to-one except at (1, 0) in (qd , qo )-space, which is mapped onto the
straight line R = 1 in the (R, θ)-space, in which the Jacobian (D-1) vanishes.
Finally, using Lemma 3, fR̂,θ̂ (R, θ) is given by (14).

References
Abrahamsen, P., 1997. A review of Gaussian random fields and correlation
functions. Technical Report 917. Norwegian Computing Center. Oslo,
Norway.
Abramowitz, M., Stegun, I.A., 1970. Handbook of Mathematical Functions.
Dover books on mathematics, Dover Publications. 1st edition.

34

Adler, P.M., 1992. Porous Media, Geometry and Transports. Butterworth
and Heinemann, Stoneham.
Adler, R.J., 1981. The Geometry of Random Fields. Wiley, New York. 1st
edition.
Allard, D., Senoussi, R., Porcu, E., 2015. Anisotropy models for spatial data.
Mathematical Geosciences , 1–24.
Anderson, T.W., 1984. An Introduction to Multivariate Statistical Analysis.
Wiley, New York. 3rd edition.
Bobach, T., Farin, G., Hansford, D., Umlauf, G., 2009. Natural neighbor
extrapolation using ghost points. Computer-Aided Design 41, 350–365.
Bouchaud, J.P., Georges, A., 1990. Anomalous diffusion in disordered media:
Statistical mechanisms, models and physical applications. Physics Reports
195, 127–293.
Chilès, J.P., Delfiner, P., 2012. Geostatistics: Modeling Spatial Uncertainty.
Wiley, New York. 2nd edition.
Chorti, A., Hristopulos, D.T., 2008. Nonparametric identification of anisotropic (elliptic) correlations in spatially distributed data sets. IEEE
Transactions on Signal Processing 56, 4738–4751.
Christakos, G., 1992. Random Field Models in Earth Sciences. Academic
Press.
Dubois, G., Galmarini, S., 2006. Spatial interpolation comparison (SIC)
2004: introduction to the exercise and overview of results, in: Dubois,
G. (Ed.), Automatic Mapping Algorithms for Routine and Emergency
Monitoring. Office for Official Publications of the European Communities, Luxembourg, European Communities. volume EUR-21595-EN, pp.
7–18. https://wiki.52north.org/pub/AI_GEOSTATS/EventsSIC2004/
EUR_SIC_2004_online.pdf. Accessed: 13 May 2016.
Ecker, M., Gelfand, A., 1999. Bayesian modeling and inference for geometrically anisotropic spatial data. Mathematical Geology 32, 67–82.

35

Ecker, M., Gelfand, A., 2003. Spatial modeling and prediction under stationary
non-geometric range anisotropy. Environmental and Ecological Statistics
10, 165–178.
Feller, W., 1971. An Introduction to Probability Theory. Wiley, New York.
Feng, L., Hotz, I., Hamann, B., Joy, K.I., 2008. Anisotropic noise samples.
IEEE Transactions on Visualization and Computer Graphics 14, 342–354.
Fisher, P., Ledoux, H., Gold, C., 2005. An efficient natural neighbour
interpolation algorithm for geoscientific modelling, in: Developments in
Spatial Data Handling. Springer Berlin Heidelberg, pp. 97–108.
Gnedenko, B.V., Kolmogorov, A.N., 1954. Limit Distributions for Sums of
Independent Random Variables. Addison Wesley, Reading, MA.
Guilleminot, J., Soize, C., 2012. Stochastic modeling of anisotropy in multiscale analysis of heterogeneous materials: A comprehensive overview on
random matrix approaches. Mechanics of Materials 44, 35–46.
Hristopulos, D., 2002. New anisotropic covariance models and estimation of
anisotropic parameters based on the covariance tensor identity. Stochastic
Environmental Research and Risk Assessment 16, 43–62.
Hristopulos, D., 2005. Spartan Gaussian random fields for geostatistical
applications: Non-constrained simulations on square lattices and irregular
grids. Journal of Computational Methods in Science and Engineering 5,
149–164.
Hristopulos, D., Žukovič, M., 2011. Relationships between correlation lengths
and integral scales for covariance models with more than two parameters.
Stochastic Environmental Research and Risk Assessment 25, 11–19.
Isserlis, L., 1918. On a formula for the product-moment coefficient of any order
of a normal frequency distribution in any number of variables. Biometrika
12, 134–139.
Jiang, X., 2005. On orientation and anisotropy estimation for online fingerprint
authentication. IEEE Transactions on Signal Processing 53, 4038–4049.

36

Kazianka, H., 2013. spatialcopula: A matlab toolbox for copula-based spatial
analysis. Stochastic Environmental Research and Risk Assessment 27,
121–135.
Lantuéjoul, C., 2002. Geostatistical Simulation: Models and Algorithms.
Springer, New York.
Le Bihan, D., Mangin, J., Poupon, C., Clark, C., Pappata, S., Molko, N.,
Chabriat, H., 2001. Diffusion tensor imaging: concepts and applications.
Journal of Magnetic Resonance Imaging 13, 534–546.
Levy, P., 1954. Théorie de l’ Addition des Variables Aléatoires. Gauthier
Villars, Paris.
Lillah, M., Boisvert, J.B., 2015. Inference of locally varying anisotropy fields
from diverse data sources. Computers & Geosciences 82, 170–182.
Okada, K., Comaniciu, D., Krishnan, A., 2005. Robust anisotropic Gaussian
fitting for volumetric characterization of pulmonary nodules in multislice
CT. IEEE Transactions on Medical Electronics 24, 409–423.
Olhede, S.C., 2008. Localisation of geometric anisotropy. IEEE Transactions
on Signal Processing 56, 2133–2137.
Papoulis, A., Pillai, S.U., 2002. Probability, Random Variables and Stochastic
Processes. McGraw-Hill, New York.
Pardo-Igúzquiza, E., Chica-Olmo, M., 1993. The Fourier integral method: an
efficient spectral method for simulation of random fields. Mathematical
Geology 25, 177–217.
Park, S.W., Linsen, L., Kreylos, O., Owens, J.D., Hamann, B., 2006. Discrete
Sibson Interpolation. IEEE Transactions on Visualization and Computer
Graphics 12, 243–253.
Pebesma, E., Cornford, D., Dubois, G., Heuvelink, G., Hristopoulos, D., Pilz,
J., Stoehlker, U., Morin, G., Skoien, J., 2011. Intamap: the design and
implementation of an interoperable automated interpolation web service.
Computers and Geosciences 37, 343–352.

37

Ranganathan, S.I., Ostoja-Starzewski, M., Ferrari, M., 2011. Quantifying
the anisotropy in biological materials. Journal of Applied Mechanics 78,
064501.
Richard, F., Bierme, H., 2010. Statistical tests of anisotropy for fractional
Brownian textures. application to full-field digital mammography. Journal
of Mathematical Imaging and Vision 36, 227–240.
Sambridge, M., Braun, J., McQueen, H., 1995. Geophysical parametrization
and interpolation of irregular data using natural neighbours. Geophysical
Journal International 122, 837–857.
Schmidt, A.M., O’Hagan, A., 2003. Bayesian inference for non-stationary
spatial covariance structure via spatial deformations. Journal of the Royal
Statistical Society: Series B (Statistical Methodology) 65, 743–758.
Siotani, M., 1964. Tolerance regions for a multivariate normal population.
Annals of the Institute of Statistical Mathematics 16, 135–153.
Swerling, P., 1962. Statistical properties of the contours of random surfaces.
IRE Transactions on Information Theory IT-8, 315–321.
Wackernagel, H., 1997. Multivariate Geostatistics. Springer Verlag, Berlin.
Wang, L., Leckie, C.A., 2012. Improved Gaussian process classification via
feature space rotation. Neurocomputing 83, 89–97.
Weaver, A.T., Mirouze, I., 2013. On the diffusion equation and its application
to isotropic and anisotropic correlation modelling in variational assimilation.
Quarterly Journal of the Royal Meteorological Society 139, 242–260.
Weller, Z.D., Hoeting, J.A., 2015. A review of nonparametric hypothesis tests
of isotropy properties in spatial data. arXiv preprint arXiv:1508.05973 .
Wick, G.C., 1950. The evaluation of the collision matrix. Physical Review
80, 268–272.
Xu, Y., Choi, J., 2009. Mobile sensor networks for learning anisotropic
Gaussian processes, in: American Control Conference, 2009. ACC ’09, pp.
5049–5054.

38

Yaglom, A.M., 1987. Correlation Theory of Stationary and Related Random
Functions I. Springer Verlag, New York.
Zhang, Z., 2012. Iterative posterior inference for bayesian kriging. Stochastic
Environmental Research and Risk Assessment 26, 913–923.
Zimmerman, D.L., 1993. Another look at anisotropy in geostatistics. Mathematical Geology 25, 453–470.

39

