A test of Gaussianity based on the Euler characteristic of excursion sets

: In the present paper, we deal with a stationary isotropic random (cid:12)eld X : R d ! R and we assume it is partially observed through some level functionals. We aim at providing a methodology for a test of Gaussianity based on this information. More precisely, the level functionals are given by the Euler characteristic of the excursion sets above a (cid:12)nite number of levels. On the one hand, we study the properties of these level functionals under the hypothesis that the random (cid:12)eld X is Gaussian. In particular, we focus on the mapping that associates to any level u the expected Euler characteristic of the excursion set above level u . On the other hand, we study the same level functionals under alternative distributions of X , such as chi-square, harmonic oscillator and shot noise. In order to validate our methodology, a part of the work consists in numerical experimentations. We generate Monte-Carlo samples of Gaussian and non-Gaussian random (cid:12)elds and compare, from a statistical point of view, their level functionals. Goodness-of-(cid:12)t p (cid:0) values are displayed for both cases. Simulations are performed in one dimensional case ( d = 1) and in two dimensional case ( d = 2), using R .


Introduction
The question of the Gaussianity of a phenomenon is a historical and fundamental problem in statistical literature.This type of information can be crucial in many application problems: oceanography and waves behavior, hydrology and climatology, agronomy, neurology and spike behavior, insurance and finance, astrophysics.For instance, in this latter application, during the last decade a large number of researchers joined efforts to decide whether the Cosmic Microwave Background (CMB) temperature is Gaussian or not.The problem of determining whether an i.i.d random sample comes from a Gaussian distribution has been studied extensively.In the case where the mean and the variance of the random variable are known, one can use a classical goodness-of-fit test.However, if these parameters need to be estimated the test of Lillifors and a variant of the Cramer-Von Mises test, with estimated parameters, are well adapted.These tests are no more distribution-free and depend on the true value of the parameters.The p−values must be obtained by simulations.The situation becomes more complicated when the sample comes from a stationary stochastic process satisfying some mixing conditions.For this type of problem, some tests have been designed to determine whether the one dimensional marginals are Gaussian.We can cite, by way of illustration: the Eps test [12] based on the empirical characteristic function; the test built by Lobato and Velasco [22] that uses a test of symmetry and kurtosis; the test of Subba and Gabr [30] where the bi-spectral density is the basis for the test.A remarkable exception is the test built in Cuesta-Albertos et al. [11].There, the method for constructing the test uses a one dimensional random projection, after which the projected sample is subjected to a test that infers whether the one dimensional marginal is Gaussian.Actually, the random projection procedure allows one to test whether the whole distribution of the process is Gaussian or not, and hence is not limited to marginal distributions.
In the present article, we deal with a real valued stationary isotropic random field and we use the information given by level functionals of a single realisation of the process to build a test of Gaussianity.The level functional is the Euler characteristic (EC) of the excursion sets above some levels.Our first motivation comes from the article [1] Section 7 (Model identification).In the aforementioned paper, Adler suggests to use the EC of the excursion sets as a way to determine what is the actual distribution of the observed process.His words can describe better than ours the main goal "Suppose that we are given a real valued random field f on R 3 , which we observe on the unit cube [0, 1] 3 .Our problem is that we are not certain what the distribution of the field is.For simplicity, we shall assume that the choices are Gaussian and χ 2 , where both are assumed to be stationary and isotropic.Then one way to choose between the two is to calculate, from the data, an empirical Euler characteristic curve, ....." c.f. [1].If the data is Gaussian, then the theoretical curve is a very precise one (see Figures 1 and 9, left panels), depending on the second spectral moment of the process and some other invariant quantities.Otherwise, if it is a χ 2 -process, Kramer oscillator process or a shot noise process, a completely different curve appears (see Figures 6,8 left panel, 10 left panel, and 11).In what follows, we offer a methodology that tries to implement these ideas.
The idea of observing level functionals of a random field in order to infer some information on the distribution of the field is not new.In [20], Lindgren provides estimators for the second spectral moment of an univariate Gaussian process that are based on the number of upcrossings at various levels.In [8], Cabaña builds a test of isotropy, for a two dimensional random field, that is based on the perimeter and the area of excursion sets.In [31], the covariance function of a bivariate Gaussian field is inferred from the excursion sets Euler characteristic.The same idea has inspired many precursors working in materials science, see for instance [28], [24].In those papers, the modelling of images or slices of a two-phases materials (even more complicated) is achieved by using a two dimensional stationary Gaussian field that has been thresholded at a certain level.The observed data are the lengths of the respective two phases along any line extracted from the image.The aim of these studies is to identify the Gaussian covariance function.Let us also mention [2] where the authors start from the observation of a neurological space-time signal at some moderate levels and deduce some parameters that help in estimating the probability of exceeding very high values.Not far from this thema, one can find the question of exceedances or the study of extreme values, when considering high levels (see, e.g., the seminal work of Rice [27] and [35]).We will not go further in that direction and, at the opposite, stay with the observation of moderate levels.In all the mentioned papers, the field that is under study is assumed to be Gaussian.On the contrary, in the present paper, Gaussianity is not assumed but has to be inferred without knowing the spatial correlation.
To be more precise, we aim at proving that the function that associates u to the mean excursions EC at level u provides a kind of signature of the distribution of the random field under study.Although a so complex information as the knowledge of the whole distribution cannot be summarize in a single function, our guess is that the shape of its graph could be enough to discriminate between Gaussianity and non Gaussianity.Our main tool will be a Central Limit Theorem for the EC of an excursion set of a stationary isotropic Gaussian random field.This asymptotic normality takes place when the domain grows to R d .The result is proved in [13] with the help of a Breuer-Major Theorem [25].We will also need some generalisations, extensions and explicit computations of this result.In particular, in the present work, we extend the results in [13] by showing that the random variables ) for i = 1, . . ., m; k = 1, . . ., p; with φ(X, T i , u k ) the EC of a standard Gaussian field X associated to the level u k in the domain T i and E[φ(X, T i , u k )] the theoretical Gaussian mean of excursions EC at level u k , are asymptotically jointly distributed as Gaussian when the disjoint domains T i grow up and satisfy asymptotic independence (see Proposition 5 for further details).Furthermore, we provide a tractable expression of the associated asymptotic variances (see Propositions 2 and 8).
The aforementioned results are used to build a Gaussianity test for a standard random field X.Indeed, I. if the H0 hypothesis: "the random field X is Gaussian" holds true, then the test statistic based on the sum of the scaled Z i k 's follows an appropriate chi-squared distribution.This case is for instance, illustrated by the pretty QQ-plots alignments and high p−values of goodness-of-fit tests in Sections 4.1.3and 5.1.Conversely, II.if the underlying standard random field X is not Gaussian with a given second spectral moment (in particular if X is a χ 2 or a Kramer oscillator process), we deliberately center again, as in the I. case, the Z i k 's variables by using the (wrong) theoretical Gaussian mean of excursions EC at level u k (see Equation ( 6)) with the same second spectral moment.In this case, we obtain a very small goodness-of-fit p−values for the chi-squared distribution associated to the considered test statistics based on Z i k 's (see Sections 4.2.1 and 4.2.2).Then we are able to reject the H0 hypothesis under considered alternatives.The crucial point is that the theoretical mean of excursions EC has a very different shape in the Gaussian case with respect to the considered alternatives: χ 2 -process or a Kramer oscillator process.This evident difference allows us to easily discriminate between Gaussian and non Gaussian setting.
As mentioned above, we consider two alternative hypothesis versus the Gaussianity: χ 2 -process and Kramer oscillator process.These processes are chosen for two types of reasons.First, the mean excursions EC curve for the considered processes can be analytically known.Secondly, they are typical examples of models with specific properties: positive asymmetry or oscillation.This variety of behaviors can cover real-life situations in the applications.For the same types of reasons, we also consider a discontinuous shot noise process as an alternative.However, due to the non-smoothness property of this process, we separate its study in a specific section (see Section 6).A further discussion in this sense is postponed to the conclusion section.
Main contributions of the paper.Under the assumption of stationarity and isotropy for the Gaussian field X, we give a new explicit formula for the second moment of the excursions EC (see Proposition 1 and Proposition 2).We extend the results in [13] by showing a Central Limit Theorem for the joint excursions EC concerned with different levels and disjoint domains (see Proposition 5).We also show finiteness of the asymptotic variance of the excursions EC of a chi-square field (see Proposition 7).A numerically tractable formula for the asymptotic variance in the univariate Gaussian case is given in Proposition 8. We propose a statistical methodology to test the Gaussianity distribution and we implement that on simulated data-sets to evaluate the finite-sample performance of our strategy.
Outline of the paper.Section 1 contains the general setting and the definition of the observation tool, namely the Euler characteristic of excursion sets.In Section 2, we focus on the Gaussian hypothesis.Assuming that the field under study X is a stationary isotropic Gaussian field, we give explicit formulas for the first two moments of the excursions EC.We also recall the Central Limit Theorem satisfied by the excursions EC when the domain tends to R d (see Theorem 4) and give an extension for the joint excursions EC concerned with different levels and disjoint domains.It allows us to describe the situation in term of a statistical model with observations whose distribution is known under the null hypothesis H0.Section 3 is concerned with the study of two alternative distributions of the observed random field: χ 2 and Kramer oscillator.In both cases, we give an explicit formula for the mean EC of excursion sets.Section 4 is devoted to numerical illustrations for univariate processes.We generate trajectory samples of stationary processes, Gaussian and non Gaussian, and compare the theoretical mean function of the excursions EC and the empirical one.We also build some chi-square statistics and associated goodness-of-fit tests in order to quantify the deviation between Gaussian and non Gaussian case.At last, in Section 5, we go further in the numerical study by considering two dimensional random fields.We generate Gaussian and χ 2 samples and compare with the theoretical situation.All the random generations and numerical computations are performed with R. In the last Section 6, we consider a shot noise process as a potential alternative, giving a preliminar study that could lead to a more general test of Gaussianity for non smooth processes.In order to keep the methodological spirit of the present paper, we have reported the technical proofs to the appendix section.

Setting
Let us start by introducing some definitions and terminology.Let X be a real valued random field defined on R d .We will sometimes call it a multivariate process.We say that it is stationary if its distribution is invariant under translations in R d and isotropic if its distribution is invariant under rotations in R d .Hence, the distribution of a stationary isotropic random field is invariant under the action of the isometry group of R d .Assuming X is smooth enough, we write X ′ for its derivative, which is a random vector field from R d to R d .We also write X ′′ for its second derivative.For any t ∈ R d , we denote by ∇ 2 X(t) the 1  2 d(d + 1) random vector that contains the upper coefficients of the symmetric Hessian matrix X ′′ (t) and by X(t) the d + 1  2 d(d + 1) + 1 random vector (X ′ (t), ∇ 2 X(t), X(t)).Let us assume now that X is stationary, isotropic and centered.If almost every realisation of X is of class C 1 , then the covariance matrix of X ′ (0) has the following form for some λ > 0 usually named as second spectral moment.
All over the paper, we consider a real valued random field X defined on R d that satisfies the following assumption.

Assumption (A):
The random field X is stationary, isotropic, E(X(0)) = 0, Var(X(0)) = 1 and almost all realisations belong to C 3 (R d ).For any fixed t in R d , the covariance matrix of the random vector X(t) has full rank.
At last, the covariance function r is such that, where ψ(t) = max ) .
Remark: Assumption (A) could appear too strong for the task described in this paper.Actually, we are not interesting in working under the weakest assumption, our aim is rather to provide a methodology that applies for any dimension d.Since the theoretical tools that we use are mainly derived from [13], we chose to work with the same assumption setting as in this paper.

Notations.
• for any u ∈ R and any compact T ⊂ R d , we call "excursion set of X above the level u within the domain T " the following set {t ∈ T : X(t) ≥ u}, • p Z (.) denotes the probability density function of any random vector Z (assuming it exists), • | • | denotes without any ambiguity, either the absolute value, or the d-dimensional Lebesgue measure.

Euler characteristic.
The Euler characteristic of a compact domain K in R d can be heuristically defined in the case d = 1 as the number of connected components of K, or in the case d = 2 as the number of connected components minus the number of holes in K.In the case where K is an excursion set {t ∈ T : X(t) ≥ u}, with T a rectangle in R d and u a real number, there exists a rather tractable formula that uses the theory of Morse functions (see [3] Chapter 9, for instance).It states that the Euler characteristic of {t ∈ T : X(t) ≥ u} is equal to a sum of two terms.The first one only depends on the restriction of X to the interior of T , it is given by the quantity φ(X, T, u) defined in Equation (2) below.The second one exclusively depends on the behaviour of X on the l-dimensional faces of T , with 0 ≤ l < d.From now on, we focus on the term φ(X, T, u), named as "modified Euler characteristic" in [13], and we still call it Euler characteristic (EC).It is defined by the following, with T the interior of T and the "index" stands for the number of negative eigenvalues.
Special case 1: dimension one.When d = 1, writing [0, T ] instead of T for a while, Equation (2) becomes Morse's theorem says that this quantity is linked with the number of up-crossings, Taking expectation in both expressions and using stationarity yield the next formula that we will use in the forthcoming sections, namely Sections 3 and 4, Special case 2: dimension two.When d = 2, Equation ( 2) can be rewritten in the following way.
With the notations introduced within this equation, µ 0 (T, u) denotes the number of local maxima above u, µ 2 (T, u) denotes the number of local minima above u and µ 1 (T, u) the number of local saddle points above u.Hence, − #{local saddle points of X above u in T }.

Under Gaussian hypothesis
In this section, we assume that X is Gaussian and satisfies Assumption (A) described in Section 1.

First two moments of the Euler characteristic of an excursion set
Let T be a cube in R d .This section is devoted to explicit formulas for the first two moments of φ(X, T, u).They are based on the decomposition in (2) and on Rice formulas for the factorial moments of µ k (T, u) (see for instance [3] Chapter 11 or [5] Chapter 6).
In particular, using the stationarity of X, the expectation can be computed as follows where we recall that λ is the second spectral moment of X, see (1).Moreover, it is proved in [3] Lemma 11.7.1, through a regression and due to Wick's formula, that where H k stands for the Hermite polynomial of order k.Hence, the next formula holds In what follows, we will be particularly interested in the next function that yields E[φ(X, T, u)] = |T | C(u).Equation (7) shows that C(u) implicitly depends on X through its dimension parameter d and its second spectral moment λ.Whenever necessary in the next sections, we will emphasize this dependence by writing C(u) = C(u, λ).
For the second moment, a so nice formula as (6) seems to be out of reach.Nevertheless, in the next proposition, we express the second moment as an integral that can be numerically evaluated.Let us mention the paper [32] where a similar formula is established but concerning another modified Euler characteristic, namely the differential topology (DT) characteristic.In this reference, the authors also give the second moment of the excursions DT characteristic as an integral of conditional distribution functionals and they propose a Monte-Carlo simulation method to evaluate these functions.We will use the following functions, defined for u ∈ R and t ∈ R d , In the one dimensional case, explicit formulas for the expectation functions g(u) and G(u, t) are given in the numerical Section 4 (see Proposition 8).

Proposition 1. Assume that X is Gaussian and satisfies Assumption (A).
Then, for any u ∈ R, the map t → G(u, t) D(t) −1/2 is integrable on any compact set in R d and Proof.Integrability comes from [13] Proposition 1.1 since X ∈ C 3 .In order to compute the expectation of φ(X, T, u) 2 , let us start with (2).It yields The expectation of the first term is equal to where we have used Rice Formula to get the second line and stationarity as well as the independence between X ′ (0) and (X(0), X ′′ (0)) to get the third one.
For the second and third terms, we introduce, for k, l = 1, . . ., d and s, t ∈ R d , where D k denotes the set of symmetric matrices having index equal to d − k.Hence, Rice Formula for the second factorial moment allows us to obtain Given that D k ∩ D l = ∅ for k ̸ = l, one can adapt the proof of the second moment Rice formula to get Let us remark that It yields It remains to compute the probability density function of (X ′ (0), X ′ (t)).The covariance matrix of this vector is equal to . Hence the result.

Asymptotic variance
In the next proposition, we let the cube T grow to R d and we give a formula for the asymptotic variance of φ(X, T, u).Actually, we consider the image of a fixed cube T by the dilation t → N t and we let N grow to +∞.
Proposition 2. Assume that X is Gaussian and satisfies Assumption (A) and let T be a cube in R d .Then for any u in R, where C(u), g(u), G(u, t), D(t) have been defined in (7) and (8).
Proof.From Proposition 1, for a fixed cube T , we have where we have used the relation | dt to get the last line.Hence, the asymptotic formula can easily be derived using Lebesgue dominated convergence theorem conditionally to the fact that t → G(u, t) D(t) −1/2 − C(u) 2 belongs to L 1 (R d ).This point is the matter of the next lemma.

Lemma 3. The map t → G(u, t) D(t)
The proof of Lemma 3 is postponed in the appendix section.
Beyond the existence of a finite asymptotic variance as stated in the previous proposition, φ(X, T (N ) , u) satisfies a Central Limit Theorem.
Theorem 4 (Theorem 2.6 in [13]).Assume that X is Gaussian and satisfies Assumption (A) and let T be a cube in R d .Then for any u in R, the next convergence holds in distribution where N (0, V (u)) stands for the centered Gaussian distribution with variance V (u) (see Equation ( 9)).

Disjoint domains and various levels
We now consider two domains T 1 and T 2 that are disjoint and two levels u 1 and u 2 that can be equal or not.

Proposition 5. Assume that X is Gaussian and satisfies Assumption (A). (a)
Let T 1 and As N → +∞, ) converges in distribution to a centered Gaussian vector with diagonal covariance matrix ) where V (u i ) is prescribed by (9).

(b)
Let T be a cube in R d and let u 1 and u 2 belong to R. For any integer N > 0, we introduce As N → +∞, ) converges in distribution to a centered Gaussian vector with covariance matrix ) , where V (u i ) is prescribed by (9) and V (u 1 , u 2 ) is a finite real number.
Comment.In [13] Theorem 2.5, the covariance V (u 1 , u 2 ) is prescribed by a series expansion that, although convergent, is so awkward that it cannot be used to evaluate it in practise.

Proof of Proposition 5.
Concerning item (a), we first have to establish that the fields Z (N ) 1 and are asymptotically independent.Intuitively, it comes from the fact that the distance between the cubes T goes to infinity and the covariance function of X has a sufficient rate of decay at infinity due to Assumption (A).Precisely, the asymptotic decorrelation is given in the next lemma.

Lemma 6. Under assumptions of Proposition 5 and using the same notation, it holds that
Lemma 6 is proved in the appendix section.By using this result and Theorem 4, we know that the covariance matrix of the random vector (Z ) tends to ) . One can use the same arguments as those in [13] (which were inspired by the Breuer-Major Theorem of [25]) to establish that any linear combination xZ has a Gaussian limit in distribution.Item (b) is proved in [13], Theorem 2.5.

Statistical model
We are now able to build a test for the H0 hypothesis: "the random field X is Gaussian".Actually, we assume that X is observed through the family (Y i k ) 1≤i≤m;1≤k≤p where for disjoint domains T 1 , . . ., T m that have the same large volume and are at large distance one from each other, and various levels u 1 ≤ . . .≤ u p .Under these conditions, by using Proposition 5 and Equation ( 6), one can write the following approximation valid under H0, where The deterministic value C(u k , λ) is given by Equation (7).It depends on the level u k and on the standard field X only through its second spectral moment λ.
) is given by ( 9), and for k ̸ = l, since no explicit formula is available for V (u k , u l ) in practise, it has to be estimated.Using the statistical model in ( 11)-( 12), we will focus on two particular cases.
(a) Diagonal case.We take m disjoint domains T 1 , . . ., T m and m levels u 1 , . . ., u m , such that each level u k is associated to a single domain T k .In this setting we have We take m disjoint domains T 1 , . . ., T m and p levels u 1 , . . ., u p , such that different levels u k are associated to the same domain T i .In this setting we have mp observations (Y i k ) i,k and their covariance matrix is given by (12).

χ 2 hypothesis
In this section, we deal with χ 2 distributions instead of Gaussian ones.Let us fix s, a non negative integer, as the degrees of freedom.We start with {X i (.)} s i=1 , an independent sample of standard stationary Gaussian fields on R d that satisfy Assumption (A) of Section 1.We denote by r X their covariance function and recall that r X (0) = 1.Consider now the following stationary fields Note that for any t ∈ R d , χ 2 s (t) is a chi-square random variable with s degrees of freedom.One get readily that Therefore, Z (s) also satisfies Assumption (A).Moreover, its second spectral moment is equal to λ with λ = −((r X ) 2 ) ′′ (0) = −2r ′′ X (0).
We are now interested in the expectation of the Euler characteristic of for a fixed cube T ⊂ R d and a fixed level u in R. A formula for the mean EC of excursion sets of χ 2 fields is given in [34], Theorem 3.5 (see also Theorem 15.10.1 in [3]).It applies in our context, although we have to handle carefully with the second spectral moment of the X i 's, which is equal to λ/2.It yields where P d,s (•) is a polynomial of degree d − 1 with integer coefficients (see Section 3.3 in [34]).In particular, we stress that P 1,s (u) = 1 and Let us recall that, in dimension one, the mean Euler characteristic of the excursion above the level u is equal to the mean number of upcrossings at level u, see Equation (4).With this point of view, Formula (13) can also be found in [29] for instance.Next proposition is concerned with the second moment in dimension one, in the same spirit as Proposition 2.
Proposition 7. Let d = 1.Let us assume that the X ′ i s are one dimensional i.i.d.Gaussian processes that satisfy Assumption (A) for i = 1, . . ., s.Then, for any u in R and T in [0, +∞), φ(Z (s) , [0, T ], u) admits a finite second moment.Morever, there exists v s (u) ∈ [0, +∞) such that lim The proof of Proposition 7 can be found in the appendix section.

Kramer oscillator hypothesis
We work in dimension d = 1.Let us consider the following system of stochastic differential equations, well known as Kramer oscillator system, where V (q) = a 0 q 4 − a 1 q 2 for positive constants a 0 and a 1 , σ and c are also positive constants, and W is a Brownian motion.The asymptotic properties of such a system have been studied for instance in [36].Besides, it is well known that the Markov process (Q, P ) has an invariant measure µ that can be written (up to a numerical constant factor) as ) dp dq.
From now on, we assume that (Q(0), P (0))'s distribution is proportional to µ, so that (Q, P ) is stationary.Here, we are interested in the process Q.It is stationary, centered, but certainly not Gaussian since its distribution is proportional to exp ) dq.Nevertheless, the distribution of its derivative process, Q ′ = P , is actually Gaussian with zero mean and variance equal to σ 2 2c .An application of Rice formula gives the mean number of upcrossings of Q in [0, T ] at any level u and hence, using (4), we get Moreover, let us remark that a suitable choice of the parameters σ, c, a 0 , a 1 allows us to get Var(Q(0)) = 1 and Var(Q ′ (0)) = λ, so that Q satisfies the same moments constraints as the generic process X in Section 1. Actually, it is sufficient to prescribe σ 2 2c = λ and to solve the following non linear equation

Univariate numerical illustrations
In this section, we focus on the one dimensional case and hence only deal with univariate processes.Both for Gaussian and alternative distributions, we first compare the theoretical formulas for the moments of excursions EC with the empirical ones, obtained by Monte Carlo simulations.Secondly, we perform the test statistics and calculate associated goodness-of-fit p−values.A graphical illustration of our test is also provided by drawing QQ-plots.

First and second moments: E[φ(X, T, u)] and V (u)
We rewrite Equations ( 6) and ( 9) of Section 2 in the case d = 1 under H0 hypothesis, i.e. when X is a stationary isotropic standard Gaussian process with covariance function r and second spectral moment and Var[φ(X, where D(t) = (2π) 2 (λ 2 − r ′′ (t) 2 ) and g(u), G(u, t) are given by ( 8).In the one dimensional case, it yields Our first task is now to provide formulas for g(u) and G(u, t) that can be numerically evaluated.They will obviously include Gaussian integrals and, therefore, we need to describe the Gaussian distributions that are involved.
On the other hand, the conditional distribution is a 4-dimensional centered Gaussian distribution with covariance matrix given by where C 11 (t) is the covariance matrix of the vector (X(0), X(t), X ′′ (0), X ′′ (t)), C 22 (t) is the covariance matrix of the vector (X ′ (0), X ′ (t)) and C 12 (t) is the matrix of the covariances between those two vectors.
We still need some extra notations: • ϕ denotes the standard Gaussian density, Φ the standard Gaussian distribution, and In Proposition 8 below, we finally give explicit formulas for the functions g and G.They will be useful to give a numerical evaluation of V (u) in (18) for various values of u.
Proposition 8. Let X be a Gaussian process that satisfies Assumption (A).
1.For any u in R, ) .
The proof of this proposition is postponed in the appendix section.
In order to illustrate ( 17) and ( 18), we generate a 300 Monte-Carlo sample of a stationary standard Gaussian process with covariance function r(t) = e −t 2 .Note that it implies the second spectral moment λ = 2.By using the Cholesky function chol in R programm we are able to numerically evaluate L(t) matrix, for all t.These quantities are necessary to build g, G and numerically approximate the integral that gives the asymptotic variance V (u) (see Equation ( 18)).In order to evaluate φ(X, T, u) on each realisation of X, we use Equation ( 3) for various values of u.Comparison between theoretical formulas and empirical counterparts are shown in Figure 1. ) with E(X(0)) = 0, Var(X(0)) = 1 and covariance function r(t) = e −t 2 .In this case λ = 2.

Estimation of the second spectral moment
It has been already quoted that the second spectral moment of X, denoted by λ, plays an important role in Equation (17).Actually, all the influence of X in E[φ(X, T, u)] is summarised through this parameter.From our statistical point of view, since we aim at inferring the distribution of X from observations of the excursions of X, parameter λ is unknown a priori.It has to be estimated from the observed level functionals X.
We estimate γ = λ 1/2 by using an unbiased estimator introduced by Lingren [20] for stationary zeromean Gaussian processes.It is based on p different levels u 1 < u 2 < . . .< u p by the following prescription, Actually, in the paper of Lingren, the number of up-crossings of level u k in the interval [0, T ], namely U (X, [0, T ], u k ), is used instead of φ(X, [0, T ], u k ) (recall the analogy described in Section 1).As a general rule, considering seems to be an acceptable choice (see discussion in [20]).In this case ).An illustration of this estimation procedure is given in Figure 2.This estimation of λ 1/2 is also used in Figure 1 (left), dashed line.The same estimation strategy could be developed in the alternative hypothesis framework as soon as we have at our disposal a closed formula linking the mean EC with the second spectral moment (see for instance Equation (21) in the univariate chi-square case).

Chi-square statistics
In the following we will consider particular sub-cases of the two cases presented in Section 2.4, diagonal case (a) and crossed case (b).In particular, we focus on the four next models.We will now associate chi-square statistics with each of the four aforementioned models.For the diagonal case (a.1), we consider ) is given by ( 17) and V (u) is given by (18).Furthermore, E and V ar in F a1 are respectively the empirical mean and the empirical variance on considered Monte-Carlo sample generations.Under H0 hypothesis, a consequence of ( 11)-( 12) is that both random variables F a1 and F a1 are approximately χ 2 m distributed, i.e. central chi-square with m degrees of freedom.We evaluate F a1 and F a1 on 300 Monte Carlo simulations.We choose m = 3 and u = 1.2.The QQ-plot comparison between the obtained empirical quantiles with the theoretical quantiles of a centered χ 2 m distribution is gathered in Figure 3 (first and second panels).

Consider now the diagonal case (a.2). Let
where E and V ar are respectively the empirical mean and the empirical variance on considered Monte Carlo sample generations.By using again the statistical model in ( 11)-( 12), both F a2 and F a2 are approximately χ 2 m distributed.An illustration is presented in Figure 3 where Z is the p-dimensional Gaussian vector given by Z E is the empirical mean on the Monte-Carlo simulations and the matrices Λ and Λ are defined below.Let Γ = (V (u k , u l )) 1≤kl,≤p be the (theoretical) covariance matrix of Z and let Λ stand for any square root of Γ.Similarly, let Γ be the empirical covariance matrix of Z evaluated on the Monte-Carlo sample generations, and let Λ be any of its square root matrix.Hence, F b1 and F b1 are both approximately χ 2 m p distributed, with m = 1.An illustration of the behaviour of F b1 is presented in Figure 4 (first panel), by choosing 300 Monte-Carlo simulations, m = 1, |T | = 200 and p = 3 different levels u 1 = −1.5, u 2 = 0, u 3 = 1.5.

Consider now the crossed case (b.2).
In this setting we have m disjoint domains T 1 , . . ., T m and p different levels u 1 < . . .< u p .Let where for any i ∈ {1, . . ., m}, Z i is the p-dimensional Gaussian vector given by and E is the empirical mean on the Monte-Carlo simulations.Let Γ = (V (u k , u l )) 1≤kl,≤p be the (theoretical) covariance matrix of Z i .Let Λ stand for any square root of Γ.Similarly, let Γ (i) be the empirical covariance matrix of Z i evaluated on the Monte-Carlo sample generations, and let Λ (i) be any of its square root matrix.Hence, Moreover, since for 1 ≤ i ̸ = j ≤ m, the Gaussian vectors Z i and Z j are independent, F b2 and F b2 are still centered χ 2 distributed with now mp degrees of freedom.
An illustration of the behaviour of F b2 is presented in Figure 4 2.

Alternative processes
In this section, we simulate trajectories of an alternative standard process Z.We compare it with a standard Gaussian process X having the same second spectral moment.
Similarly to Section 4.1.3,we now consider a test statistic.For sake of brevity, we only focus on the diagonal case (a.1) with a unique level u and m disjoint domains T 1 , . . ., T m , and we introduce where E[φ(X, T i , u)] is the mean Euler characteristic of the excursions of X given by ( 6) with d = 1 and V ar(φ(Z, T i , u i )) is the empirical variance on considered Monte-Carlo sample generations.In order to use this empirical variance, one has to be sure that the variance of the Euler characteristic of the excursions of Z is finite.We recall that this is the case when Z is chi-square distributed (see Proposition 7).Under H0 hypothesis, the D statistic is asymptotically distributed as a chi-square with m degrees of freedom.

First alternative: χ 2 process
In conformity with Section 3.1, we consider where the X i 's are independent copies of a stationary standard Gaussian process with covariance function r(t) = e −t 2 /2 .Hence, we get λ = −2r ′′ (0) = 2 and then the obtained process Z (s) has the same variance and the same second spectral moment λ as the previous Gaussian one presented in Section 4.1.
In the first and third panels of Figure 5, we display the boxplot for the ratio between the empirical mean of 300 Monte Carlo values of φ(Z (s) , T, u) and the theoretical mean given by (13) in the case d = 1, i.e., This expectation has to be compared with the expectation given by ( 17) when X is a stationary standard Gaussian process with second spectral moment λ = 2. Therefore, in the second and fourth panels of Figure 5, we display the boxplot for the ratio between the empirical 300 Monte Carlo mean value of φ(Z (s) , T, u) and the Gaussian theoretical expectation given in (17).
Remark.When the degrees of freedom s tend to infinity, the CLT implies that Z (s) tends in distribution to a stationary centered Gaussian process with covariance function equal to t → r(t) 2 , which implies a variance equal to 1 and a second spectral moment equal to 2 (exactly as the Gaussian process X considered in Section 4.1).On the other hand, using Stirling approximation for the Γ function, one can prove that the right-hand side of ( 21) tends to the right-hand side of (17), i.e., See also Formula (3.4) in [29] for the same remark.The comparison between the second and fourth panels in Figure 5, as well as Figure 6, illustrate this convergence.
In Figure 7 below, we show the QQ-plot of the test statistic D (see (20)) versus the quantiles of χ 2 m .Since the considered chi-square process Z (s) is not Gaussian, a deviation can be observed.In order and the theoretical mean given by Equation (21).Second and fourth panels: Boxplot for the ratio between the same empirical values and the theoretical mean in the Gaussian case given by Equation (17).In both cases λ = 2.The degrees of freedom s is chosen equal to 2 and 10 respectively.to quantify how we are really able to reject the null hypothesis H0 under the alternative, in Table 3 we display the generated p−values for the goodness-of-fit test of D. One can easily see that this test allows us to reject the fact that the considered underlying process Z (s) is Gaussian.This is particularly evident if one compares these p−values with those in Tables 1 and 2.

Boxplot using Chi−square and Gaussian Expectations, s=2
Figure 7 shows a huge deviation from the bisector of the first orthant in the case s = 2 whereas this deviation is less significant in the case s = 10.It emphasizes the convergence in distribution of Z (s) towards a Gaussian as s tends to infinity.7 versus a χ 2 m distribution with m = 6.The considered process Z (s) (•) is a chi-square univariate process (d = 1) with E(Z (s) (0)) = 0, Var(Z (s) (0)) = 1 and λ = 2.

Second alternative: Kramer oscillator process
In this section we generate a 300 Monte-Carlo sample of a Kramer oscillator process as defined in Section 3.2.In order to obtain a process Q with zero mean, unit variance and second spectral moment equal to 2, we solve Equation ( 16) and choose The generation procedure is the following.We define a discretized schema to simulate the solution (Q(t), P (t)) of the stochastic differential system in (14).Actually, we use the Metropolization of the Euler-Verlet schema with a sufficient small discretization step.The interested reader is referred to Algorithm 2.11 (Generalized Hybrid Monte-Carlo) in [19].
The comparison between the (theoretical) expectation given by ( 15) and the empirical one is shown in left panel of Figure 8 below.Furthermore, as we did in previous Section 4.2.1, we also consider the D test statistic.The QQ-plot of D versus χ 2 m distribution is shown in right panel of Figure 8.Since the considered process Q is not Gaussian, a strong deviation is observed.It allows us to reject the H0 hypothesis.To quantify this deviation we present the generated p−values for the goodness-of-fit test associated to D in Table 4.   8 (right) versus a χ 2 m distribution with m = 6.The considered process Q(•) is a Kramer oscillator process as defined in Section 3.2 with E(Q(0)) = 0, Var(Q(0)) = 1 and second spectral moment λ = 2.

Bivariate numerical illustration
In this section, we focus on two dimensional random fields.For both Gaussian and chi-square distributions, we apply the same methodology as in the previous section.In particular, we use the test statistic D introduced in (20), its definition being not specific to dimension one.

Under H0 hypothesis
We consider a stationary centered Gaussian random field X = {X(t) : t ∈ R 2 }.Its restriction to a finite regular grid included can be seen as a model for a grey level image.The modified Euler characteristic of an excursion set in the rectangle domain T above level u is given by Equation ( 5).On the other hand, Equation (6) in dimension d = 2 gives its expectation, In what follows, we generate a 300 Monte-Carlo sample of a bivariate stationary centered Gaussian random field X with covariance function r(t) = e −||t|| 2 , t ∈ R 2 .In that case, E(X(0)) = 0, Var(X(0)) = 1 and the second spectral moment λ is equal to 2. We use (5) in order to compute φ(X, T, u) for a fixed cube T and various values of u.The local extremum points of X are given by the R function extrema2dC in the EMD package.To identify the saddle points, we find all the stationary points of X in the considered domain T (i.e., all the points with an associated null gradient function) and we exclude the points previously identified as local extremum points.
In the left panel of Figure 9, the comparison between the theoretical mean given by ( 22) and the empirical mean based on the simulations is illustrated.
In the right panel, we consider the chi-square statistic D introduced in ( 20) for a unique level u and m disjoint domains T 1 , . . ., T m .The considered process is now the Gaussian bivariate process X.Furthermore, in Table 5, we display the goodness-of-fit test p−values for D. The obtained high p−values allow us to statistically accept the χ 2 m distribution for the test statistic D. In this case λ = 2.

Alternative: χ 2 bivariate process
As an alternative to hypothesis H0, let us consider where the X i 's are independent copies of a centered stationary Gaussian two dimensional field with covariance function r(t) = e −||t|| 2 /2 .As described in Section 3.1, for any fixed t ∈ R 2 , χ 2 s (t) has a central χ 2 distribution with s degrees of freedom.Furthermore, the field Z (s) is centered, stationary and its covariance function is given by r(t) = e −||t|| 2 .Hence its variance is equal to 1 and its second spectral moment is equal to λ = −2r ′′ (0) = 2.The expectation of φ(Z (s) , T, u) is given by ( 13) with We generate a sample of 300 realisations of the bivariate process Z (s) for s = 2 on a fixed cube T .We use Equation ( 5) to empirically compute φ(Z (2) , T, u) for various values of u.On the other hand, Equation ( 13) yields the following in the case In the left panel of Figure 10, we compare Equation ( 23) and its empirical counterpart based on the Monte Carlo simulations of Z (2) .In the right panel, we show the QQ-plot of the test statistic D versus the χ 2 m distribution, where D is computed with the realisations of Z (2) .Since the considered chi-square process Z (2) is not Gaussian, a large deviation can be observed.Associated generated p−value for the goodness-of-fit test of D is given in Table 6.This quantifies how we are really able to reject the null hypothesis H0 under the alternative.10 (right) versus a χ 2 m distribution with m = 3.The considered processZ (2) (•) is a chi-square bivariate process (d = 2) with covariance function r(t) = e −||t|| 2 .In this case E(Z (2) )(0) = 1, Var(Z (2) (0)) = 1 and λ = 2.

Non smooth processes: a case study
In this section, we sketch a preliminary study concerning non smooth processes.For this purpose, the shot noise process is a convenient toy model.We will only consider the one dimensional case pointing out two recent references that deal with dimension two, [7] and [18].Let d = 1 and let us introduce the following shot noise process S, where a > 0 and Ψ is a homogeneous Poisson point process on R with intensity ν > 0. The process S is clearly stationary with zero mean and variance equal to νa.Moreover, its values almost surely belong to the discrete set {k − νa, k ∈ N}.
See associated goodness-of-fit test p−value in Table 6.
In what follows, we still aim at testing the Gaussianity hypothesis.If a whole trajectory of S could be observed, its discrete shape would clearly indicate the right choice between Gaussian and shot noise.We now assume that only excursion sets of S above a few levels are observed and we consider the number of "upcrossings".Note that the notion of upcrossings has to be properly defined in this context, since S is not continuous.Actually, it is proved in [6] (see also [16] for a different approach) that for any level u ∈ R \ {k − νa, k ∈ N} and any interval T , This formula, which obviously differs from (6), is an indicator that could enable to discriminate between a Gaussian or a shot noise model.
In the following, we generate 300 trajectories of such a process on a fixed interval T with a = ν = 1.Since aν = 1, then Var(S(0)) = 1.A comparison between the theoretical expectation and the Monte-Carlo empirical mean is presented in Figure 11.One can observe the evident different shape of the mean EC curve compared to the Gaussian case (see Figure 1, left).

Conclusion
In the present paper, we have proposed a methodology to test whether a random field defined on R d is Gaussian.The test statistics are computed from the observation of a single realisation, precisely the Euler characteristic of excursion sets at moderate levels is concerned.Our approach requires two theoretical ingredients that have to be valid under the Gaussian hypothesis: a Central Limit Theorem satisfied by the Euler characteristic of excursion sets when the domain tends to R d , and a closed formula that gives the mean Euler characteristic of the excursion sets.We have established both of them.We consider this methodological work as a first step, there are many open questions and possible future researches extending and going further this work.
There are two major domains where the Gaussianity hypothesis is really relevant and where our results could certainly be applied: neurology and cosmology.The neurological signals that are collected when studying brain activity are commonly observed through their level sets (see [35] and [21]).A natural question concerns the type of the signal distribution: Gaussian or Poisson?Gaussian or Oscillator?
The test that we have built, together with the alternatives that we have considered, could be an appropriate tool for studying this problem.On the other hand, in the analysis of the Cosmic Microwave Background (CMB) radiation, the question of Gaussianity has been tackled in a huge amount of publications (see [26] for a recent overview).In that context, the random field under study is defined on the two dimensional celestial sphere.Hence, in order to be applied, our methodology should be extended beyond the Euclidean case.Let us mention that recent studies, like [10], [23], [9] for instance, contain theoretical results on high energy behavior of spherical random fields that could provide the required background for a test of Gaussianity in the same spirit as the one presented here.
Furthermore, it is not difficult to imagine that the same question of Gaussianity is actually asked in other real-life situations, and that the potential alternatives take various shapes, like a chi-square when the process under study is always positive, or an oscillator when the process under study presents an almost periodic structure (sea waves) or a shot noise process when sudden peaks are observed.Our methodology could also serve as a goodness-of-fit test in the opposite way.For instance, in the area of geostatistics, the authors of [14] observe the Euler characteristic of excursion sets to detect abrupt changes on soil data sets.They use as an a priori model for the data a chi-square random field.A natural extension to this study could be a test of chi-square distribution versus not chi-square.
In order to go further in this direction, it would be necessary to establish Central Limit Theorems for the Euler characteristic of excursion sets for new classes of processes or random fields.Beside the present work, we implemented simulations that lead us to postulate that such a Central Limit Theorem could exist, in particular for the two alternative processes that we have considered.
At last, another appealing question concerns the hypothesis that the process under study is expected to satisfy, such as smoothness of the realisations and fast decay of the covariance function.As already mentioned at the beginning of the paper, Assumption (A) seems to strong at least in the one dimensional case.Indeed, let us remark that in that case, two important results of Section 2 remain valid under weaker assumptions.Actually, on the one hand, the Rice formula gives the mean number of upcrossings under the only assumption of a finite second spectral moment.On the other hand, it is well known that a sufficient condition for the number of upcrossings to have a finite second moment is the so-called Geman condition, which only requires that r ′′ , the second derivative of the covariance function of X, exists and satisfies some integrability condition near 0. Such considerations allows us to believe that the theoretical results that we have presented deserve to be studied under weaker assumptions, even in dimension greater than one.Another extension could concern non continuous random fields, like the shot noise processes.

A.1. Proof of Lemma 3.
The level u is fixed, so we do not refer to it in the following lines.For t ∈ R d \ {0}, we define h(t) = G(u, t) D(t) −1/2 − C(u) 2 .We have to establish that h is integrable on R d .We already know that h is integrable on any compact set of R d , so it remains to study the behavior of h(t) for large ∥t∥.
We write K = 1 2 d(d + 1) + 1 and we introduce the function where det(z) stands for the determinant of the d × d symmetric matrix whose upper coefficient are given by the 1 2 d(d + 1) dimensional vector z.Recall that we denote by X(.) the stationary and isotropic random vector field (X ′ (t), ∇ 2 X(t), X(t)), which has dimension Let Σ be the covariance matrix of X(0) and, for any t ∈ R d , let Γ X (t) be the covariance matrix of (X(0), X(t)).The latter can be written by blocks, each of one having dimension D × D, as follows ) .
Moreover we know that all the terms of matrix Γ X (t) are uniformly bounded in absolute value by ψ(t).Using Theorem 3.1 of [33] we get, for any where • the sum runs over all D × D matrices J = (J jk ) 1≤j,k≤D with non negative integer entries and for any D-dimensional index k and any z ∈ R D , with ϕ Σ the probability density function of any D-dimensional Gaussian vector N (0, Σ), We can bound |h(t)| in the following form where where the last inequality is a consequence of expanding det(∇ 2 X(0)) as a multivariate polynomial function of degree d evaluated at the coordinates of the matrix X ′′ (0).We now concentrate on H A J ((0, ∇ 2 X(0), X(0)), Σ).Following [33], we have where Z is any D-dimensional centered Gaussian random vector with covariance matrix Σ. Remember that Σ is the covariance matrix of X(0) = (X ′ (0), ∇ 2 X(0), X(0)), so it can be factorized as Σ = ( and Z can be expanded as Z = (Z 0 , Z 1 ) where Z 0 is a d-dimensional centered Gaussian vector with covariance matrix λI d and Z 1 is a K-dimensional centered Gaussian vector with covariance matrix Σ 1 , Z 0 and Z 1 being independent.Then Σ −1 ((0, x) ) , where the two blocks of coordinates are independent.So where the D-dimensional vector A J is equal to (A (0) J ) with blocks of respective size d and K.For the first term, we note that H A (0) where H k (.) is the usual multidimensional Hermite polynomial of multi-order k.Concerning the second term, we compute E thanks to a formula for multivariate Gaussian integrals that can be found in [17] for instance (see equation (35) in this reference),

J
where diag(B) denotes the vector containing the diagonal coefficients of any square matrix B. Finally, since H A (0) where we have used that |A (0) Coming back to inequality (27) yields The same inequality can be established for D(B J ).Hence, from (26) we get The series in the right hand is convergent for ∥t∥ large enough.Indeed we can choose A such that for ∥t∥ > A, we have 2D∥Σ −1 ∥ ψ(t) ≤ η < 1 and so the series is bounded by (1 − η) −1 .
The result follows since on the other hand, by assumption (A), ∫ ∥t∥>A ψ(t)dt < +∞.

A.2. Proof of Lemma 6.
Without loss of generality, we can assume that |T .Hence, by stationarity It is proved in [13] Proposition 1.3 that, for a cube T in R d and a level u in R, the following expansion holds both a.s. and in L 2 (Ω).
where D = d + 1 2 d(d + 1) + 1, H n denotes the D dimensional Hermite coefficient of multi-order n ∈ N D and Y (t) = Λ −1 X(t) with Λ any square root of the covariance matrix of X(0).Moreover the Hermite coefficients a(n, u) are such that for any q ≥ 1, where C is some positive constant that only depends on d, Λ and u.Hence, thanks to (28) and to the orthogonality of Hermite polynomials, we have for any where One can adapt Arcones inequality, for instance following the proof of [4] Lemma 1 step by step, to establish that there exists a constant κ > 0 depending on the covariance function of X such that Let now η > 0 be fixed.Thanks to Assumption (A), we can choose N large enough such that κψ(τ ) ≤ 1/2 as soon as ∥τ ∥ ≥ N δ, ∫ ∥t∥≥N ψ(t) dt < η and N δ > 2. Hence, using (29), we deduce that for ∥τ ∥ ≥ N δ, Note that C stands here for a positive constant that may change from a line to another.Coming back to the covariance of Z (N ) 1 and Thus, Lemma 6 is proved.
We now focus on the asymptotic variance as T goes to infinity.We will actually prove that lim Following (32), we now expand J(t, x) as J(t, x) := 4 u σ 2 ε (t) In first place, let us study the limit when t → ∞ of the (k = 0)-term in the integrand.We introduce the following notation, ) sin s−2 ψ dψ.
We have given that d 2 u0 (t) and using (31).Therefore, recalling the identity This last equality holds true if we show that the difference appearing into the integrand is bounded by an L 1 (R) function, outside of a compact interval.Actually, under (A), it is easy to prove that for t large enough, |I 0 (t) − K s (u) 2 | ≤ C| du(t) √ uσε(t) | ≤ Cψ(t).Hence, the limit (34) is established.
In the sequel we are going to study the asymptotic behaviour of the remaining terms.Let us introduce
Summing up, we have Integrating by parts the second integral, we get that it is equal to λϕ(u)(2Φ( λu The second term is equal to ) 1/2 u). Finally, ) 1/2 u).

Computation of G(u, t).
Since Γ(t) is the covariance matrix of (X(0), X(t), X"(0), X"(t)) conditioned to (X ′ (0) = X ′ (t) = 0), one can write From now on, we remove the dependence on t.We begin by writing Z = L Y, where L is a lower triangular matrix such that LL T = Γ(t).Denoting by l ij for i ≥ j the elements of L, we have and the expectation (36) can be written as By expanding the second sum, the integral can be written as the sum of the following terms (starting with index i = 4 term)

For
the crossed case (b.1), let us now consider different levels u 1 < . . .< u p and one single domain T , i.e. m = 1.Let us define

2 Fig 5 .
Fig 5.  First and third panels: Boxplot for the ratio between the empirical 300 Monte Carlo values of φ(Z (s) , T, u) and the theoretical mean given by Equation(21).Second and fourth panels: Boxplot for the ratio between the same empirical values and the theoretical mean in the Gaussian case given by Equation(17).In both cases λ = 2.The degrees of freedom s is chosen equal to 2 and 10 respectively.

Fig 6 .
Fig 6.Boxplots of the empirical 300 Monte Carlo values of φ(Z (s) , T, u) for the considered chi-square univariate process Z (s) and different values of u.Red points represent the theoretical means given by Equation (21) for the same values of u; blue ones are the Gaussian means given by Equation (17).In this case λ = 2, |T | = 200, s = 2 (left panel) and s = 10 (right panel).

Table 1
(third and fourth panels), by choosing 300 Monte Carlo simulations, m = 3, |T i | = 200 and different levels u 1 = −1.2,u 2 = 0, u 3 = Goodness-of-fit test p−values associated to test statistics F a1 , F a1 , F a2 and F a2 in Figure3..2.Furthermore, in Table1we display the goodness-of-fit test p−values for F a1 , F a1 , F a2 and F a2 in Figure3.Obtained p−values allow us to statistically accept the χ 2 m distribution for these test statistics.Moreover, we highlight how p−values increase when one evaluates the EC curve for several levels u i (see F a2 and F a2 cases). 1

Table 3
Goodness-of-fit test p−values associated to test statistics D in Figure

Table 4
Goodness-of-fit test p−value associated to test statistic D in Figure

Table 6
Goodness-of-fit test p−value associated to test statistic D in Figure