A note on the higher moments of the euler characteristic of the excursion sets of random fields

Ali Reza Taheriyoun^{1}

Shahid Beheshti University, IRAN

Khalil Shafie^{2}

University of Northern Colorado, USA

Mohammad Jafari Jozani^{3}

Université de Sherbrooke, CANADA

ABSTRACT

In this paper a useful and suitable form of the second moment of the Euler characteristic (EC) of excursion sets of a random field with respect to a special threshold is presented which can be used in numerical evaluations of the variance of EC. We also presented an extension to the higher moments of EC. Also a new method of simulating is presented and the results are compared with the three older well known methods of simulation based on the moment estimation and monte carlo idea.

Keywords: Euler characteristics; DT characteristic; Isotropic Gaussian fields; Monte carlo method.

1 Introduction

The maximum of a Gaussian random field was used to test activation at unknown points in Positron Emission Tomography (PET) of blood flow in the brain. PET is a functional, nuclear medical imaging technology allowing to display different metabolic states of the human body. It is common to estimate the number of the activation region using EC of the excursion sets. Activation can be expressed as crossing blood pressure from a special level like u (naturally the number of components of an active region is estimated by EC). Working with EC is not as simple as differential topology (DT) characteristic because of its geometric behavior. It is worth mentioning that EC is equal to the differential topology (DT) characteristic as long as the intersection of the boundary of the research area with the excursion set is a null set. So, it is common to use DT characteristic instead of EC in such cases.

Expected EC is a valuable parameter to evaluate some behavior of EC and it is a good representor of the number of active region’s components when the finite dimensional distribution of the random field is known and a sample of size two or more is available. Adler and Hasofer (1976) as well as Adler (1981) presented a good formulae to compute the expected DT. They also found a closed form of the expected DT characteristic of homogeneous Gaussian random fields. A useful and different method of computation of the expected EC is considered by Worsley (1994) with application to the expected EC of χ^{2}, F and t fields. Cao and Worsley (1999) also computed the expected EC of excursion sets for Hotteling’s T^{2} field using Worsley (1994) method.

1Department of Statistics, Faculty of Mathematical Sciences, Shahid Beheshti University, Tehran, Iran (e-mail:

a taheriyoun@sbu.ac.ir)

2Applied Statistic and Research Methods, (e-mail: khalil.shafie@unco.edu)

3Université de Sherbrooke, Département de mathématiques, Sherbrooke, QC, CANADA, J1K 2R1 (e-mail:

mohammad.jafari.jozani@usherbrooke.ca)

Many statistical inferences about EC, such as asymptotic properties of EC, are closely related to its second and higher moments. Higher moments of the EC can also be used to estimate the fractal dimension of a rough surface which is the subject of a work by the authors.

Increasing the roughness of realization of a random field increases the value of EC at special levels. Concerning this fact, it is important to determine the confidence boundary of EC for different levels of crossing by use of the variance of EC. Meanwhile, for the data of digital sky surveys, some cosmic waves are detected as 3-dimensional random fields and it is important to determine the EC of them in given levels. To this end, it is better to compute the EC and using the variance of it to construct a confidence boundary to make a criteria for classifying the complex of observations. A good entry point to recent statistical literature on this topic is discussed by Adler et al. (2006).

To our best knowledge there are few works related to characterizing the higher moments of EC. The factorial moment of the number of up-crossings for one dimensional random fields is considered by Cramer and Leadbetter (1967). Belyaev (1972a) gave an expression for the higher moments of EC without any proof. In Balyaev (1972b) details were given, but as mentioned by Adler (1981) there were some mistakes in that proof regarding to the appropriate use of the mean value theorem. Adler and Taylor (2007) present a very general theorem on the higher moments of a counting random field (Theorem 11.5.1 ,page 284) that may be used to obtain the factorial moments of the EC. In this work, based on a general method of Adler (1981) we give a direct proof for the exact and not asymptotic behavior of the higher moments of EC and provide some computing methods through a simulation study in deriving the second moment of EC for Gaussian fields.

To this end, in section 2 some required point-set definitions and concepts are represented. In section 3 a useful and suitable expression for the second moment of EC is presented and the result is extended to the higher moments of EC. Some numerical studies for evaluating the produced integrals are given in section 4.

2 Point-Set Representation

Our framework closely resembles the one introduced by Adler (1981), and several of the following definitions and notations are reproduced here for sake of completeness. Note that the excursion set inside a fixed set S is just the set of the points where the field exceeds a fixed and predefined threshold value. In other words if X(t) is a N-dimensional random field then for any fixed real number u the excursion set of X(t) above the level u in S is defined as

A_{u}(X(t), S) = {t ∈ S : X(t) ≥ u}.

In this paper we simply use A_{u} to denote the excursion set of X(t) above the level u in S.

Let X(t), t = (t_{1}, t_{2}, …, t_{N})′ ∈ ℜ^{N} be a homogeneous real-valued random field and S be a compact subset of ℜ^{N}. Throughout this paper we shall denote derivatives with respect to t_{i} by X_{i} and second order derivatives with respect to t_{i} and t_{j} by X_{ij}. Let D(t) be a (N −1)×(N −1) matrix of second order derivatives of X(t), with elements X_{ij}(t), i, j = 1, 2, …, N − 1. Index of D(t) is the number of negative eigenvalues of it. Under the suitable regularity conditions on X(t), DT characteristic of A_{u} is defined as

χ(A_{u}) = (−1)^{N−1} ∑_{j=0}^{N−1}(−1)^{j}χ_{j}(A_{u}),  (1)

where χ_{j}(A_{u}) is the jth type number, i.e., the number of points t ∈ S satisfying the following conditions:

a) X(t) = u,

b) X_{i}(t) = 0 i = 1, 2, …, N − 1,

c) X_{N}(t) > 0,

d) index of D(t) is equal to j.

It can be shown that the definition of χ(A_{u}) depends on particular coordinate system and changing this system or any rotation can change the value of χ(A_{u}). But when the excursion set does not touch the boundary of S (i.e., ∂S), that means ∂S ∩ A_{u} = ∅, then χ(A_{u}) is equal to EC which is invariant under the rotation of coordinates. This definition of χ(A_{u}) was arisen from Morse theorem and basically is the way of the counting the connected components of an excursion set in one dimensional random fields. Roughly speaking, the Euler characteristic counts the connected components of excursion set in one dimensional random fields. In two dimensions, it counts the number of connected components of A_{u} minus the number of ’holes’ that pass through the A_{u} and in three dimensions it is the number of connected components minus holes plus the number of ’hollows’ inside the A_{u}. Similar definition can be provided for the higher dimensions based on geometric viewpoint. It is natural to focus on DT characteristic instead of EC, because of its simpler form.

The expected DT characteristic is the alternative sum of expectations E[χ_{j}(A_{u})] which is derived by Adler and Hasofer (1976). A simple form of expected DT of homogenous zero-mean Gaussian fields was derived by Adler (1981).

3 Higher Moments of EC

In this section we find the second moment of DT characteristic of excursion sets of random fields and extend it to the higher moments through Theorem 2. To this end, lower and upper bounds for the expectation of the squared (and higher degrees of) DT characteristic are derived. This method is similar to the one used by Adler (1981) for deriving the expectation of DT. Note that using (1) we have

[χ(A)]^{2} = [(−1)^{N−1} ∑_{j=0}^{N−1}(−1)^{j}χ_{j}(A_{u})]^{2}

= ∑_{j=0}^{N−1} ∑_{k=0}^{N−1}(−1)^{k+j}χ_{k}(A)χ_{j}(A).  (2)

So, taking the expectation of (2) we have

E[χ(A)]^{2} = ∑_{j=0}^{N−1} ∑_{k=0}^{N−1}(−1)^{k+j}E[χ_{k}(A)χ_{j}(A)].

In the following lemma we derive an upper bound for E[χ_{k}(A)χ_{j}(A)]. Let ϕ be the joint distribution function of

[X(t), X(s), X_{1}(t), X_{1}(s), …, X_{N−1}(t), X_{N−1}(s), X_{N}(t), X_{N}(s), Z(t), Z(s)],

where Z(·) is the vector of the length of N(N−1)/2 obtained from D(·) by replacing the successive columns on and above the main diagonal of D(·) under one another.

Lemma 1 Let X be a homogeneous N-dimensional random field (N > 1) suitably regular with respect to a compact subset S ⊂ ℜ^{N} at level u and ∂S is a null set with respect to the Lebesgue measure. If the second-order partial derivatives of X have finite fourth moments, ϕ is continuous in each of its variables, and the conditional density of [X(t), X(s), X_{1}(t), X_{1}(s), …, X_{N−1}(t), X_{N−1}(s)] given X_{N}(t), X_{N}(s) and the second-order partial derivatives X_{ij}(t), X_{ij}(s) i, j = 1, 2, …, N − 1, are bounded above, then the expectation of χ_{k}(A_{u})χ_{j}(A_{u}) satisfies

E[χ_{k}(A_{u})χ_{j}(A_{u})] ≤ ∫_{S^{2}} ∫ ∫ ∫_{x_{N} >0} ∫_{y_{N} >0} x_{N}y_{N} (det D(t)) (det D(s))
× ϕ(u, u, 0, …, 0, x_{N}, y_{N}z, z′) dx_{N} dy_{N} dz dz′ dμ(t) dμ(s),  (3)

where the third and fourth integrals are taking over z ∈ ℜ^{N(N−1)/2} and z′ ∈ ℜ^{N(N−1)/2} for which the matrices D(t) and D(s) have k and j negative eigenvalues respectively.

Proof: For ϵ > 0, let σ(ϵ) is the ball of radius ϵ defined by σ(ϵ) = {x ∈ ℜ^{N} : ∥x∥ < ϵ}, and δ_{ϵ}(·) is a function on ℜ^{N} as δ_{ϵ}(·) = cI_{{x∈σ(ϵ)}} such that ∫_{σ(ϵ)} δ_{ϵ}(x)dx = 1. Define M(t) as follow,

M(t) = [ X_{1}(t) ; X_{2}(t) D(t) ; ⋮ ; X_{N}(t) X_{N1}(t) ⋯ X_{N N−1}(t) ] ,

with

D(t) = [ \begin{matrix}X_{11}(t) & X_{12}(t) & … & X_{1 N−1}(t) \\ X_{21}(t) & X_{22}(t) & … & X_{2 N−1}(t) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ X_{N−1 1}(t) & X_{N−1 2}(t) & … & X_{N−1 N−1}(t) \end{matrix} ] .

Now we can write

χ_{k}(A_{u})χ_{j}(A_{u}) = lim_{ϵ→0} lim_{ζ→0}[χ_{k}^{ϵ}(A)χ_{j}^{ζ}(A)],  (4)

such that

χ_{k}^{λ}(A_{u}) = ∫_{S} δ_{λ}(X(t)) I_{A′_{k}}(Y(t)) |det M(t)| dt, λ ∈ {ϵ, ζ},

where Y(t) = [X_{N}(t), X_{11}(t), …, X_{N−1 N−1}(t)], X(t) = (X(t) − u, X_{1}(t), …, X_{N−1}(t)) and A′_{k}, k = 0, 1, …, N − 1, is a subset of ℜ^{\frac{N(N−1)}{2} +1} such that for any Y(t) ∈ A′_{k}, X_{N}(t) > 0 and the index of D(t) is k. Then we have,

χ_{k}^{ϵ}(A_{u})χ_{j}^{ϵ}(A_{u}) = ∫_{S} δ_{ϵ}(X(t)) I_{A′_{k}}(Y(t)) |det M(t)| dt  (5)
× ∫_{S} δ_{ϵ}(X(s)) I_{A′_{j}}(Y(s)) |det M(s)| ds.

Now for any ϵ, ζ > 0 by convergence of χ_{i}^{ϵ}(A_{u}), as ϵ → 0 we have

E[χ_{k}^{ϵ}(A_{u})χ_{j}^{ζ}(A_{u})] = ∫ [∫_{S} ∫_{S} δ_{ϵ}(X(t)) δ_{ζ}(X(s)) |det M(t)| |det M(s)| dt ds]
×ϕ(·) dx dy dx′ dy′ dz dz′,

where x′ and y′ denote the vectors (x_{1}, …, x_{N−1}) and (y_{1}, …, y_{N−1}). Third integration is taking over all positive x_{N} and y_{N} and over all x, x_{1}, …, x_{N−1}, x_{N}, all x_{Nj}s and y, y_{1}, …, y_{N−1}, y_{N}, and all y_{Nj}s, j = 1, 2, …, N − 1 and over all x_{ij}s and y_{ij}s, 1 ≤ i ≤ j ≤ N − 1 such that the index of D(t) and D(s) be k and j, respectively. Now using Fubini’s theorem

E[χ_{k}^{ϵ}(A_{u})χ_{k}^{ζ}(A_{u})] = ∫_{S} ∫_{S} ∫ ∫ dz dz′ dx_{N} dy_{N} ∫ [δ_{ϵ}(x − u, x_{1}, …, x_{N−1})δ_{ζ}(y − u, y_{1}, …, y_{N−1})]

× |det M(t)| |det M(s)|ϕ(.) dx dy dx_{1} dy_{1} … dx_{N−1} dy_{N−1} dμ(t) dμ(t). (6)

Expanding detM(t) and detM(s) the inner integrals converge to

(−1)^{j+k}x_{N} y_{N} (det D(t)) (det D(s)) ψ(u, u, 0, 0, …, 0, x_{N}, y_{N}, z, z′),

as ϵ, ζ → 0, where ψ is the conditional density of [X(t), X(s), X_{1}(t), X_{1}(s), …, X_{N−1}(t), X_{N−1}(s)] given X_{N}(t), X_{N}(s) and the second-order partial derivatives X_{ij}(t), X_{ij}(s) i, j = 1, 2, …, N − 1.

Operating the dominated convergence theorem, the inner integral in(6) is bounded above by an integrable function. Note that the inner integral is not greater than

[B_{N}(ϵ)B_{N}(ζ)]^{−1}ψ_{2}(x_{N}, y_{N}, z, z′) ∫_{x=u−ϵ}^{u+ϵ} ∫_{y=u−ζ}^{u+ζ} ∫_{x_{1}=−ϵ}^{ϵ} ∫_{y_{1}=−ζ}^{ζ} … ∫_{x_{N−1}(t)=−ϵ}^{ϵ}

× ∫_{x_{N−1}(s)=−ζ}^{ζ} ∑_{k=1}^{N} ∑_{j=1}^{N} |x_{k}||y_{j}| |a_{k}||a_{k}| ψ_{1}(x, y, …, x_{N−1}, y_{N−1}|x_{N}, y_{N}, z, z′)

× dx dy … dx_{N−1} dx_{N−1},

where ψ_{1} and ψ_{2} are respectively the joint distribution function of [X_{N}(t), X_{N}(s), Z(t), Z(s)] and the conditional joint distribution function of [X(t), X(s), …, X_{N−1}(t), X_{N−1}(s)] given [X_{N}(t), X_{N}(s), Z(t), Z(s)], and B_{N}(ϵ) is the volume of σ(ϵ). So, for finite constants K_{N} and G_{N}, we have [B_{N}(ϵ)B_{N}(ζ)] = K_{N}ϵ^{2N} G_{N}ζ^{2N}. We know that there exists a finite C such that ψ_{1} < C. Thus the last expression is not larger than

CK_{N} G_{N}ψ_{2}(x_{N}, y_{N}, z, z′)

× (∑_{k=1}^{N−1} |a_{k}| + |x_{N}| |det D(t)| ) (∑_{k=1}^{N−1} |a′_{k}| + |y_{N}| |det D(s)| ) ,

where each terms of this expression is integrable. Now the proof will be completed using dominated convergence theorem, Fatou’s lemma and integrating out x_{Nj} and y_{Nj}, j = 1, 2, …, N − 1. We can easily show that

E[χ_{k}(A)χ_{j}(A)] ≤ lim_{ϵ→0} lim_{ζ→0} E[χ_{k}^{ϵ}(A)χ_{k}^{ζ}(A)]
= ∫_{S} ∫_{S} ∫ ∫ ∫_{x_{N} >0} ∫_{y_{N} >0} x_{N} y_{N} |det D(t)| |det D(s)|
ψ(u, u, 0, 0, …, 0, 0, x_{N}, y_{N}, z, z′)
× dx_{N} dy_{N} dz dz′ dμ(t) dμ(s),

where the third and fourth integrals are taking over z ∈ ℜ^{N(N−1)/2} and z′ ∈ ℜ^{N(N−1)/2} for which the matrices D(t) and D(s) have k and j negative eigenvalues respectively.

In what follows we show that the upper boundary obtained in Lemma 1, can also be a lower bound for expectation of χ_{k}(A_{u})χ_{j}(A_{u}). To do this, we need the following geometric concepts and notations.

Definition 1 If s and t are two distinct points in ℜ^{N}, then the line segment joining s and t, L(s, t), is define as follows

L(s, t) = {u : u = θs + (1 − θ)t, 0 ≤ θ ≤ 1}.

In the following lemma, let L_{n}, n ≥ 1, be the lattice of points in ℜ^{N} whose components are integer multiples of 2^{−n}; i.e.,

L_{n} = {t ∈ R^{N} : t_{j} = i.2^{−n}, j = 1, …, N , i = 0, ±1 ± 2, …}.

and define Δ_{n}^{ϵ}(t) be hypercube centered on an arbitrary point t as follows

Δ_{n}^{ϵ}(t) = {s ∈ R^{N} : |s_{j} − t_{j}| ≤ (1 − ϵ)2^{−(n+1)}, j = 1, …, N},

with Δ_{n}(t) = Δ_{n}^{0}(t). Let also ω_{i} and ω_{ij} be the moduli of continuity of X_{i} and X_{ij} respectively which are defined as

ω_{i}(h) = sup_{∥t−s∥<h}|X_{j}(t) − X_{j}(s)|, and ω_{ij}(h) = sup_{∥t−s∥<h}|X_{ij}(t) − X_{ij}(s)|.

Lemma 2 Suppose that the assumptions of Lemma 1 hold, and S is convex and the moduli of continuity ω_{i} and ω_{ij} of X_{i} and X_{ij} satisfy the following condition for any ϵ > 0

P {max_{ij}[ω_{i}(h), ω_{ij}(h)] > ϵ} = o(h^{N}).  (7)

Then (3) holds with the inequality sign reversed.

Proof: For any i ∈ L_{n} let I_{ni}^{∗} = I_{χ_{k}(A∩Δ_{n}(i))χ_{j}(A∩Δ_{n}(i))>0}. An approximation (χ_{j}(A)χ_{j}(A))_{n} for χ_{k}(A)χ_{j}(A) is defined as

(χ_{k}(A)χ_{j}(A))_{n} = ∑_{i∈L_{n}} I_{ni}^{∗}.

By suitable regularity conditions for any n and i we have

P {χ_{k}(⋃ ∂Δ_{n}(i))χ_{j}(⋃ ∂Δ_{n}(i)) = 0} = 1.

Since ∂Δ_{n}(i) is in a N − 1 dimensional space, the derivative with respect to N^{th} element is zero for any points in it. Thus

(χ_{k}(A)χ_{j}(A))_{n} ⟶^{a.s.} χ_{k}(A)χ_{j}(A).

It is clear that I_{ni}^{∗} = I_{ni}^{j}I_{ni}^{k} where I_{ni}^{j} = I_{χ_{j}(A∩Δ_{n}(i))>0}. Since ∑_{i∈L_{n}} I_{ni}^{∗} is positive non-decreasing in n and (χ_{k}(A)χ_{j}(A))_{n} is non-decreasing in n, using monotone convergence theorem

E[χ_{k}(A)χ_{j}(A)] = lim_{n→∞} E[(χ_{k}(A)χ_{j}(A))_{n}]
≥ lim_{n→} ∑ P {χ_{k}(A ∩ Δ_{n}(i))χ_{j}(A ∩ Δ_{n}(j)) > 0}.  (8)

Now we should investigate P {χ_{k}(A ∩ Δ_{n}(i))χ_{j}(A ∩ Δ_{n}(j)) > 0}. For a given realization of the random field X define the function ω_{X}^{∗}(n) as

ω_{X}^{∗}(n) = max [ max_{1≤j≤N} ω_{j}(2^{−n}), max_{1≤i,j≤N} ω_{ij}(2^{−n})] .

Let M_{X(t)} = max [max_{1≤j≤N} sup_{t∈S} |X_{j}(t)|, max_{1≤i,j≤N} sup_{t∈S} |X_{ij}(t)|], and δ and ϵ be two arbitrary small real numbers and K is a positive number. If η = δ^{2}ϵ/[2^{N}N N!(K + 1)^{N−1}], then

P {ω_{X}^{∗}(n) > η} = o(2^{−Nn}),  (9)

as n → 0. For suitable n and t, s ∈ L_{n} such that Δ_{n}(s), Δ_{n}(s) ⊂ S assume that ω_{X}^{∗}(n) < η and the event G_{δK} happens, where

G_{δK} = {|det D(t)| > δ, |det D(s)| > δ, X_{N}(t) > δ, X_{N}(s) > δ, M_{X} < K,

index of D(t) = k, index of D(s) = j}.

Define t^{∗} and s^{∗} as the solutions of the following equations, when the unique solutions in fact exist:

X(v) = (v − v^{∗})M(v), v ∈ {t, s}.  (10)

We claim that if ω_{X}^{∗}(n) < η and for t^{∗} ∈ Δ_{n}^{ϵ}(t) and s^{∗} ∈ Δ_{n}^{ϵ}(s) the event G_{δK} occurs and for large enough n we have χ_{k}(A ∩ Δ_{n}(t))χ_{j}(A ∩ Δ_{n}(s)) > 0 with both χ_{k}(A ∩ Δ_{n}(t)) > 0 and χ_{j}(A ∩ Δ_{n}(s)) > 0. So, for large enough n in (9), we have

P {χ_{k}(A ∩ Δ_{n}(t))χ_{j}(A ∩ Δ_{n}(s)) > 0} ≥
P {G_{δK} ∩ [t^{∗} ∈ Δ_{n}^{ϵ}(t)] ∩ [s^{∗} ∈ Δ_{n}^{ϵ}(s)]} − o(2^{−nN}).  (11)

To proof this claim it is enough to show that M(t) and M(s) are invertible for these realizations. For large enough n it is easy to show that

|detM(v)| > \frac{δ^{2}}{2} , v ∈ \{t, s\}.  (12)

Now, using (10), t^{∗} ∈ Δ_{n}^{ϵ}(t) and s^{∗} ∈ Δ_{n}^{ϵ}(s) are equivalent by

v − X(v).M^{−1}(v) ∈ Δ_{n}^{ϵ}(v), v ∈ {t, s}.  (13)

If κ and τ are any other points in Δ_{n}(t) and Δ_{n}(s) respectively, then under the conditions of lemma we have

|det D(τ) − det D(t)| < \frac{ϵδ^{2}}{2} ,

|det D(κ) − det D(s)| < \frac{ϵδ^{2}}{2} .

Since detD(t) > δ and detD(s) > δ, then for any τ ∈ Δ_{n}^{ϵ}(t) and κ ∈ Δ_{n}^{ϵ}(s) it follows that detD(τ) ≠ 0 and detD(κ). Also the index of D(τ) in Δ_{n}(t) and index of D(κ) in Δ_{n}(s) will be equal to k and j respectively. Similarly, we must have X_{N}(τ) > 0 through Δ_{n}(t) and X_{N}(κ) > 0 through Δ_{n}(s). If there exist one or more points τ ∈ Δ_{n}(t) at which X(τ) = (X(τ)− u, X_{1}(τ), …, X_{N−1}(τ)) = 0 then we have χ_{k}(A ∩ Δ_{n}(t)) ≠ 0 and similarly χ_{j}(A ∩ Δ_{n}(s)) ≠ 0. Similar to Adler (1981), the choice of suitable G_{δK} and ω_{X}^{∗}, t^{∗} ∈ Δ_{n}(t) and s^{∗} ∈ Δ_{n}(s) implies the existence of at least one τ ∈ Δ_{n}(t) and one κ ∈ Δ_{n}(s) at which X(τ) = 0 and X(κ) = 0. It means that

{ω_{X}^{∗}(n) < η, G_{δK}, t^{∗} ∈ Δ_{n}^{ϵ}(t)} ⇒ χ_{k}(A ∩ Δ_{n}(t)) > 0  (14)

Now from (11), (12) and (14) and performing the transformation (10) we have

P {χ_{k}(A ∩ Δ_{n}(t))χ_{j}(A ∩ Δ_{n}(s)) > 0} ≥ ∫_{G_{δK}∩{t^{∗}∈Δ_{n}^{ϵ}(t)}∩{s^{∗}∈Δ_{n}^{ϵ}(s)}} |det M(t)||det M(s)|
× ψ_{u}((t − t^{∗}) M(t), (s − s^{∗}) M(s), x_{N}, y_{N}, x_{ij}, y_{ij})
× dt^{∗} ds^{∗} dx_{ij} dy_{ij} − o(2^{−nN}),

where 1 ≤ i ≤ N and 1 ≤ j ≤ N − 1, and ψ_{u} is the joint density function of the vector [X(t) − u, X(s) − u, X_{i}(t), X_{i}(s), X_{ij}(t), X_{ij}(s)], 1 ≤ i ≤ N, 1 ≤ j ≤ N − 1. If we let n → ∞ it follows from the mean value theorem for the integral that the last expression multiplied by the number of terms in the summation (8), converges to

∫_{S} ∫_{S}(1 − ϵ)^{2N} ∫_{G_{δK}} |det D(t)| |det D(s)| x_{N} y_{N}ψ_{u}(0, 0, …, 0, 0, x_{N}, y_{N}, z, z′)
× dx_{N} dy_{N} dz dz′ dμ(t) dμ(s).

Now letting ϵ, δ and K^{−1} → ∞ and applying the monotone convergence theorem for the above expression we have

E[χ_{k}(A_{u})χ_{j}(A_{u})] ≥ (−1)^{j+k} ∫_{S} ∫_{S} ∫ ∫ ∫_{x_{N} >0} ∫_{y_{N} >0} x_{N}y_{N} |det D(t)| |det D(s)|
× ϕ(u, u, 0, 0, …, 0, 0, x_{N}, y_{N}, z, z′) dx_{N} dy_{N} dz dz′ dμ(t) dμ(s).

and it completes the proof.

Now, the following theorem is an immediate consequence of Lemma 1 and 2.

Theorem 1 Under the conditions of Lemma 2 we have

E[χ_{k}(A_{u})χ_{j}(A_{u})] = (−1)^{j+k} ∫_{S} ∫_{S} ∫ ∫ ∫_{x_{N} >0} ∫_{y_{N} >0} x_{N}y_{N} |det D(t)| |det D(s)|
×ϕ(u, u, 0, 0, …, 0, 0, x_{N}, y_{N}, z, z′) dx_{N} dy_{N} dz dz′ dμ(t) dμ(s),

where the third and fourth integrations are taken over all z ∈ ℜ^{N(N−1)/2} and z′ ∈ ℜ^{N(N−1)/2} such that the indexes of D(t) and D(s) are k and j and using (2) we can write

E[χ(A_{u})]^{2} = ∫_{S} ∫_{S} ∫ ∫ ∫_{x_{N} >0} ∫_{y_{N} >0} x_{N}y_{N} |det D(t)| |det D(s)|
× ϕ(u, u, 0, 0, …, 0, 0, x_{N}, y_{N}, z, z′) dx_{N} dy_{N} dz dz′ dμ(t) dμ(s).

where third and fourth integrations are over all z ∈ ℜ^{N(N−1)/2} and z′ ∈ ℜ^{N(N−1)/2}.

An extensive generalization of Theorem 1 can be obtained using the same method; i.e., we can find the similar lower and higher boundaries for k multiplied type numbers of DT characteristic (with arbitrary indices) and use (2) for final expression. In the following theorem we consider this situation.

Theorem 2 Let X be a homogeneous N-dimensional random field (N > 1) suitably regular with respect to a compact and convex subset S ⊂ I_{0} ⊂ ℜ^{N} at level u where ∂S has zero Lebesgue measure. If the second-order partial derivatives of X have finite 2k^{th} moment, the joint density of vector

(X(t^{(1)}), …, X(t^{(k)}), X_{1}(t^{(1)}), …, X_{1}(t^{(k)}), …, X_{N}(t^{(1)}), …, X_{N}(t^{(k)}), Z(t^{(1)}), …, Z(t^{(k)})),

i.e., ϕ, is continuous in each of its variables, and the conditional density of [X_{1}(t^{(1)}), …, X_{1}(t^{(k)}), …, X_{N−1}(t^{(1)}), …, X_{N−1}(t^{(k)})] given X_{N}(t^{(1)}), …, X_{N}(t^{(k)}) and given the second-order partial derivatives X_{ij}(t^{(1)}), …, X_{ij}(t^{(k)}), i, j = 1, 2, …, N − 1, is bounded above and (7) holds, then for A_{u} = A_{u}(X, S) we have

E[χ_{j_{1}}(A_{u}) … χ_{j_{k}}(A_{u})] = (−1)^{∑_{i} j_{i}} ∫_{S^{k}} ∫ ∫_{x_{N} (t^{(1)})>0} … ∫_{x_{N}^{(k)}>0} x_{N}^{(1)} … x_{N}^{(k)}
×(det D(t^{(1)})) … (det D(t^{(k)}))
×ϕ(u, …, u, 0, …, 0, x_{N}^{(1)}, …, x_{N}^{(k)}z^{(1)}, …, z^{(k)})
×dx_{N}^{(1)} … dx_{N}^{(k)} dz^{(1)} … dz^{(k)} dμ(t^{(1)}) … dμ(t^{(k)}),

where the pre-last integration is taken over all z(t^{(i)}) ∈ ℜ^{N(N−1)/2}’s such that the index of D(t^{(i)}) is j_{i}, i = 1, 2, … k. So,

E[χ(A_{u})]^{k} = (−1)^{k(N−1)} ∫_{S^{k}} ∫ ∫_{x_{N} (t^{(1)})>0} … ∫_{x_{N}^{(k)}>0} x_{N}^{(1)} … x_{N}^{(k)} (det D(t^{(1)})) … (det D(t^{(k)}))
×ϕ(u, …, u, 0, …, 0, x_{N}^{(1)}, …, x_{N}^{(k)}z^{(1)}, …, z^{(k)})
×dx_{N}^{(1)} … dx_{N}^{(k)} dz^{(1)} … dz^{(k)} dμ(t^{(1)}) … dμ(t^{(k)}),

where the pre-last integration is taken over all z(t^{(i)}) ∈ ℜ^{N(N−1)/2}’s, without anymore restrictions.

Now, it is enough to calculate the integrals Theorem 2. However, there is no hope to find a closed form for the integrals, but the numerical methods can be used to evaluate them. The convergence of numerical methods for evaluation of multiple integrals are closely dependent to the geometry of the field. EC of excursion set of a random field on a given threshold is in the set ℜ^{+} ∪ {∞}. For example the mean of the EC of a homogenous gaussian random field is infinite as ||t|| ↓ 0 with the covariance function

r(t) = 1 − c||t||^{β} + o(||t||^{β}).

So, the methods of numerical calculation of suggested integrals which implies the boundary conditions on the exact value of integrals can not be used for the case. In the next section we provide numerical evaluation of the second moment of EC using some simulation methods.

4 Simulation Results

In this section we provide numerical evaluation of E[χ(A_{u})]^{k} for k = 2 through some simulation methods. To this end a new method of simulating E[χ(A_{u})]^{2} is presented and the results are compared with the two other well known methods of simulating based on monte carlo idea. For k = 2 and S = I_{0} which it is a unit cube in ℜ^{2}, X(t) is a two dimensional Gaussian random field and following Theorem 1 we have

E[χ(A_{u})]^{2} = ∫_{I_{0}} ∫_{I_{0}} ∫ ∫ ∫_{x_{2}>0} ∫_{y_{2}>0} x_{2} y_{2} |x_{11} y_{11}|ϕ(u, u, 0, 0, x_{2}, y_{2}, x_{11}, y_{11})
×dx_{2} dy_{2} dx_{11} dy_{11} dt ds.  (15)

• First Method (M_{1})

Let w_{i} = x_{2i}y_{2i}|x_{11i}y_{11i}|ϕ(u, u, 0, 0, x_{2i}, y_{2i}, x_{11i}, y_{11i}), where (t_{i}, s_{i}, x_{2i}, y_{2i}, x_{11i}, y_{11i}), i = 1, 2, …, n are samples taken from related integration areas given in (15). It is easy to see that a solution of (15) is simply given by the mean of w_{i}s and will be devoted by M_{1}.

• Second Method (M_{2})

This method is based on generating a sample of size n from the joint distribution of t, s, X_{2}(t), X_{2}(s), X_{11}(t) and X_{11}(s), which given t, s, ϕ(u, u, 0, 0, X_{2i}, X_{2i}, X_{11i}, X_{11i}) is a multivariate Gaussian joint density with means zero. Note that the covariance function of X(t), which has been used in numerical methods is assumed to be the Gaussian covariance function as follows

r(t, s) = exp(−(t_{1} − s_{1})^{2} − (t_{2} − s_{2})^{2}).

Suppose that t and s are independently and uniformly distributed on I_{0}. Given t and s the joint distribution of X(t), X(s), X_{1}(t), X_{1}(s) is known. So, it is easy to generate a sample from the joint distribution of X_{11}(t), X_{11}(s) given X(t), X(s), X_{1}(t), X_{1}(s) and t, s. Now, using the method introduced by Robert (1995) a suitable sample from the positive part of the joint distribution of X_{2}(t), X_{2}(s) given t, s, and X_{11}(t), X_{11}(s), X(t), X(s), X_{1}(t), X_{1}(s) is generated. Finally, the above method is repeated n times following the idea of monte carlo method to construct the values of the variable Z_{i} as follows

Z_{i} = X_{2i}X_{2i}|X_{11i}X_{11i}|ϕ(u, u, 0, 0, X_{2i}, X_{2i}, X_{11i}, X_{11i}).

It is easy to see that the mean of Z_{i}s is a solution of (15) which will be denoted by M_{2}.

• Third Method (M_{3})

In this method the main strategy is decomposing the right hand side of (15) to the three parts in which each of them can be easily computed using monte carlo method. To do this note that we can easily show

ϕ(u, u, 0, 0, x_{2}, y_{2}, x_{11}, y_{11}) = ϕ(u, u, 0, 0)ϕ_{1}(x_{11}, y_{11}|u, u, 0, 0)
×ϕ_{2}(x_{2}, y_{2}|x_{11}, y_{11}, u, u, 0, 0),

where ϕ_{1} is the joint conditional density function of (X_{11}(t), X_{11}(s)) given (X(t), X(s), X_{1}(t), X_{1}(s)) and ϕ_{2} is the joint conditional density function of (X_{2}(t), X_{2}(s)) given (X_{11}(t), X_{11}(s), X(t), X(s), X_{1}(t), X_{1}(s)). So, we have

E[χ(A_{u})]^{2} = ∫_{I_{0}} ∫_{I_{0}} ϕ(u, u, 0, 0) ∫ ∫ x_{11} y_{11} ϕ_{1}(x_{11}, y_{11}|u, u, 0, 0)
∫_{x_{2}>0} ∫_{y_{2}>0} x_{2} y_{2} ϕ_{2}(x_{2}, y_{2}|x_{11}, y_{11}, u, u, 0, 0)
dx_{2} dy_{2} dx_{11} dy_{11} dt ds.  (16)

To compute the first integral using monte carlo method it is required to generate a sample of size n_{2} × n_{11} × n from the positive part of the conditional related bivariate Gaussian distribution given (X_{11}(t), X_{11}(s), X(t) = u, X(s) = u, X_{1}(t) = 0, X_{1}(s) = 0), which implies constructing suitable samples of X_{11}(t)s and X_{11}(s)s. To do this n_{11} × n observations should be taken from the conditional related bivariate Gaussian distribution given (X(t) = u, X(s) = u, X_{1}(t) = 0, X_{1}(s) = 0) and hence it is necessary to generate a sample of size n of t and s from a uniform distribution on I_{0} × I_{0}, i.e., t_{1}, t_{2}, …, t_{n} and s_{1}, s_{2}, …, s_{n}. According to this sample the observations x_{11}(i, 1), x_{11}(i, 2), …, x_{11}(i, n_{11}) and y_{11}(i, 1), y_{11}(i, 2), …, y_{11}(i, n_{11}) are constructed corresponding to t_{i} and s_{i}, i = 1, 2, …, n. Note that the first double integral in (16) can be rewritten as

∫_{x_{2}>0} x_{2}ϕ_{2}(x_{2}|x_{11}, y_{11}, u, u, 0, 0) ∫_{y_{2}>0} y_{2}ϕ(y_{2}|x_{11}, y_{11}, u, u, 0, 0, x_{2})dy_{2} dx_{2}.  (17)

So, for each pair of (x_{11}(i, j), y_{11}(i, j)) the values x_{2}(i, j, k), i = 1, 2, …, n, j = 1, 2, …, n_{11} and k = 1, 2, …, n_{2} are generated from the positive part of the related Gaussian distribution. Then for each x_{2}(i, j, k) a sample of size n_{2} of y_{2}s is generated and first integral in (17) will be simply the mean of y_{2}s. Multiplying these means to the related x_{2}(i, j, k) and taking an average over k will approximate the first double integral (17), m(i, j) say. Hence, multiplying the m(i, j)s in related x_{11}(i, j)s and y_{11}(i, j)s and taking an average over all j indices, yields a monte carlo solution to the second double integral in (16), m(i) say. Finally the resulted approximation of (16), say M_{3}, is simply the average of all m(i) ϕ_{X(t_{i}),X(s_{i}),X_{1}(t_{i}),X_{1}(s)_{i}}(u, u, 0, 0), i = 1, 2, …, n.

Table 1 shows numerical values of the approximation of the E[χ(A_{u})]^{2} applying the above mentioned methods with (n_{11}, n_{2}) = (100, 100) for n ∈ {2, 4, 8, 16, 32}. For each combinations of n_{11}, n_{2} and n the simulation method is repeated r times r ∈ {100, 200, 300, 500} and approximation of E[χ(A_{u})]^{2} is given by the mean of the obtained values. For each method the standard error of the approximated values in r replication is given.

Using conditional expectations, it is easy to show that the generated estimators in the second and third methods are unbiased for E[χ(A_{u})]^{2} while the generated estimator in the first method is biased. Regarding the obtained results in table 1 it seems that the standard errors of M_{3} tend to zero as r is growing. So M_{3} converges to E[χ(A_{u})]^{2} faster than M_{1} and M_{2}. Table 1 shows the robustness of M_{3} against the gridding changes while the moment estimator is not robust and has a strange convergence to zero which is (E[χ(A_{u})])^{2}. For a given level of gridding the moment estimator is increasing in r which shows the exact value of E[χ(A_{u})]^{2} is greater than what has been shown by the moment estimator specially in higher grids.

Theoretically, consistency of M_{3} and the moment estimator of E[χ(A_{u})]^{2} is not obvious, but the property of their unbiasedness ensures their convergence to the parameter of interest for high level of iteration. Of course the numerical results of our simulation for the moment estimator and M_{3} are not noteworthy different but the exact relation between the values of r, n and n_{11} with the order of the convergence of these two methods is not discussed here.

To operate the procedure of computing of E[χ(A_{u})]^{2} on a certain image, one may use the maximum likelihood estimator of covariance function on relative grids and use the numerical gradient of estimated covariance function to generate the required vectors in M_{3}. In other hand, it is possible to estimate \frac{∂^{a+b}}{t_{i}^{a}t_{j}^{b}} r(t) using the maximum likelihood estimator of covariance function of numerical gradient of the image. Indeed there is no idea to show that, problem is robust under these to ways. Having no idea for generating a 2 × 2 Gaussian field the first four rows of Table 1 and 2 is empty for the moment estimator.

Table1. Estimations using n_{11} = n_{2} = 100
n r  M_{3}  M_{2}  M_{1}  MME
2 100 0.050887(0.093890) 1.77e-5(6.63e-6 ) 2.11e-8(1.55e-7)  ——(——)
2 200 0.047358(0.077103) 2.06e-5(9.89e-5) 1.05e-8(1.10e-7)  ——(——)
2 300 0.043963(0.065895) 1.96e-5(8.92e-5) 1.57e-8(1.63e-7)  ——(——)
2 500 0.043524(0.057481) 2.14e-5(9.91e-5) 2.27e-8(1.16e-7)  ——(——)
4 100 0.045924(0.078972) 8.68e-6(1.69e-5) 3.14e-8(3.19e-7) 0.031836(0.032083)
4 200 0.048488(0.085760) 6.09e-5(0.521e-3) 1.51e-8(2.22e-7) 0.031914(0.028660)
4 300 0.047634(0.077614) 4.35e-5(0.426e-3) 1.02e-8(1.84e-7) 0.032213(0.030366)
4 500 0.044597(0.064852) 4.19e-5(0.373e-3) 1.43e-8(1.55e-7) 0.036375(0.027624)
8 100 0.037594(0.022187) 1.30e-5(1.98e-5) 2.07e-7(2.07e-6) 0.014575(0.010019)
8 200 0.039537(0.026331) 3.81e-5(0.218e-3) 1.04e-7(1.46e-6) 0.014381(0.010578)
8 300 0.041579(0.029636) 2.97e-5(0.179e-3) 6.94e-8(1.19e-6) 0.012012(0.009034)
8 500 0.041390(0.027861) 4.41e-5(0.138e-3) 1.01e-7(9.19e-7) 0.013203(0.009454)
16 100 0.041696(0.041834) 2.36e-5(5.07e-5) 4.17e-9(4.05e-8) 0.007242(0.003849)
16 200 0.043560(0.034256) 1.86e-5(0.463e-3) 3.14e-9(9.14e-8) 0.007754(0.003482)
16 300 0.041970(0.027774) 2.11e-5(0.485e-3) 9.07e-9(9.14e-8) 0.008246(0.004387)
16 500 0.040107(0.024471) 2.74e-5(0.417e-3) 1.17e-8(8.63e-8) 0.008639(0.004301)
32 100 0.044741(0.033502) 2.82e-5(7.73e-5) 4.47e-9(3.74e-8) 0.005738(0.001761)
32 200 0.043184(0.026085) 4.01e-5(0.154e-3) 3.92e-9(3.06e-8) 0.005903(0.001875)
32 300 0.043199(0.023747) 3.61e-5(0.133e-3) 3.14e-9(2.55e-8) 0.006108(0.002045)
32 500 0.040239(0.020818) 4.59e-5(0.106e-3) 2.82e-9(2.19e-8) 0.006467(0.001882)

Acknowledgements

The authors are thankful to two anonymous referees whose comments and suggestions improved the presentation of our results.

References

1. Adler, R. J. (1981). The Geometry of Random Fields. John Wiley, New York.

2. Adler, R.J., and Hasofer, A. M. (1976). Level crossing for random fields. Annals of probability, 4, 1–12.

3. Adler, R. J., and Taylor, J. E. (2007). Random Fields and Geometry. Springer, New York.

4. Adler, R. J., Taylor, J. E. and Worsley, K. J. (2006). Random Fields and Geometry: Applications. In preparation.

5. Belyaev, Y. K. (1972a). Point processes and first passage problems. Proc. Sixth Berkeley Symp. Math. Statist. prob. , 2, 1–17. Univ. of California Press, Berkeley.

6. Belyaev, Y. K. (1972b). The general formulae for the mean number of crossings for random processes and fields (In Russian). In Bursts of Random Fields, Moscow Univ. Press Moscow 38–45.

7. Cramer, H., and Leadbetter, M. R. (1967). Stationary and Related Process. John Wiley, New York.

8. Robert,C. P. (1995). Simulation of truncated normal variables. Statistics and Computing, 5, 121–125.

9. Worsley, K. J. (1994). Local maxima and the expected Euler characteristic of excursion sets of χ^{2}, F and t fields. Advances in Applied Probability, 26, 13–42.

10. Worsley, K. J., Evans, A.C., Marett, S., and Neelen, P. (1992). A three dimensional statistical analysis for CBF activation studies in human brain. Journal of Cerebal Blood Flow and Metabolism, 12, 900–918.