El Bouch

Michel J J Olivier
email: olivier.michel@grenoble-inp.fr

Comon Pierre
email: pierre.comon@grenoble-inp.fr

Multivariate Normality test for colored data

Keywords: statistics, multivariate kurtosis, detection of changes, dependence, Gaussianity

Performances of the Multivariate Kurtosis are investigated when applied to colored data, with or without Auto-Regressive prewhitening, and with or without projection onto a lower-dimensional random subspace. Computer experiments demonstrate the importance of taking into account the possible color of the process in calculating the power of the normality test, in all the scenarios.

I. INTRODUCTION

Detecting changes in the distribution of a stochastic process is a longstanding problem and a myriad of methods has been proposed; see e.g [START_REF] Basseville | Detecting changes in signals and systems-a survey[END_REF]. Our concern is the detection of non-Gaussian signals in a Gaussian background, i.e. Normality tests. The framework considered here is one in which timeseries are recorded on d sensors (typically d = 2 or 3 in many applications). Moreover, we are interested in the online problem of reacting to a change as quickly as possible after it occurs, also known as sequential detection problem [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF], [START_REF] Poor | Quickest detection[END_REF].

Normality tests belong to the the class of tests without alternative, unlike other approaches based on the Likelihood Ratio such as the CUSUM test [START_REF] Page | Continuous inspection schemes[END_REF], [START_REF] Hinkley | Inference about the change-point from cumulative sum tests[END_REF]. The case of i.i.d (scalar or d-variate) processes has received significant attention, see [START_REF] Henze | Invariant tests for multivariate normality: a critical review[END_REF], [START_REF] Moore | The effect of dependence on chi squared tests of fit[END_REF] for a detailed survey. On the other hand, few tests concern the case of dealing with time-series from multiple sensors that in practical applications cannot be considered to be i.i.d, a case we will refer to as n.i.d or colored. Moreover, the proposed approaches for dealing with multivariate colored processes, such as [START_REF] Hinich | Testing for Gaussianity and linearity of a stationary time series[END_REF] [START_REF] Moulines | Testing that a multivariate stationary time series is Gaussian[END_REF], can hardly be used in real-time applications. For this reason, we have derived a joint normality test for multivariate colored data [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF] that is implemented with low computational burden.

There has also been interest in testing the normality of unobserved regression residuals [START_REF] Jarque | A test for normality of observations and regression residuals[END_REF], also referred to as the process of innovations. More precisely, define the Multidimensional Auto-Regressive [START_REF] Lütkepohl | Introduction to multiple time series analysis[END_REF] (or Vector AR of order p denoted VAR(p)) model to describe the statistical behavior of the d × 1 vector of observation x(i) for i = 1, . . . , N :

x(i) = p k=1 A k x(i -k) + ϵ(i) (1)
where A k is a d × d matrix of unknown parameters and ϵ(i) is the ith unobservable residual assumed zero-mean and

This work has been supported by the MIAI chair "Environmental issues underground" of Institut MIAI@Grenoble Alpes (ANR-19-P3IA-0003). i.i.d. An additional assumption is that residuals are drawn from a normal distribution. The drawbacks of violating the latter assumption has been studied, for instance [START_REF] Hogg | Statistical robustness: One view of its use in applications today[END_REF] showed that the ordinary least squares method, which is usually used to estimate {A k } 1≤k≤p , is sub-optimal for heavy-tailed distributions. Thus, it is important to validate this assumption of normality in this linear model.

Moreover, since the residuals are estimated in practice, we expect that errors in the model specification and estimation will impact the whitening performance of the filter and the estimated innovation process ϵ could no longer be considered i.i.d. In this case, the normality tests designed for i.i.d processes become biased as shown in [START_REF] Moore | The effect of dependence on chi squared tests of fit[END_REF], highlighting the importance of deriving our joint normality test for variables that are not statistically independent.

Within this framework, we concentrate on the Multivariate Kurtosis (MK) defined in [START_REF] Poor | Quickest detection[END_REF]. In [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF], we calculated the power of this test variable in the colored case, and in [START_REF] Elbouch | Joint Normality Test Via Two-Dimensional Projection[END_REF] we studied its performances when the observed multivariate process was projected onto an arbitrary subspace of low dimension (typically 1-D or 2-D), in particular for time series generated by colored copulas. The main contributions of the present paper are now the following:

• We compare the performances of the MK test with and without linear prewhitening, i.e. using x(n) or ϵ(n). • We observe the impact of projecting 3-D observations and their regression residuals onto an arbitrary plane (2-D projection) or direction (1-D projection). • For sake of time and memory effectiveness, the test statistics and the regression function are estimated recursively, by using both exponential averaging and the Recursive Least Squares (RLS) algorithm.

The paper is organized as follows: our recent results [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF] are summarized in Section II. The recursive implementation is described in Section III, and the complete algorithm in Section IV. Computer results are eventually reported in Section V.

II. ASYMPTOTIC STATISTICS OF THE KURTOSIS

Let x(n) = [x 1 (n), x 2 (n), . . . , x d (n)] T ∈ R d be a real- valued d-variate random variable. Let S(τ ) = E{x(n)x(n - τ )
T } be the covariance function whose entries are S ab (τ ). Also denote S(0) = S.

Let X def = {x(1), . . . , x(N )} be a stochastic stationary process of N random variables x(n). Our aim is to test:

H 0 : X ∼ n.i.d N (0, S) versus H0 (2)
Where x(n) are identically but not independently distributed (n.i.d). Following Mardia's definition of the Multivariate Kurtosis (MK) [START_REF] Mardia | Measures of multivariate skewness and kurtosis with applications[END_REF], our test statistic reads:

Bd (N ) = 1 N N n=1 x(n) T S -1 x(n) 2 (3)
with

S = 1 N N k=1 x(k)x(k) T (4)
A. Known results: Mardia's MK test for a i.i.d. process

The following results make use of Landau's notations o( 1 N ) and O( 1N ), to precise that the absolute approximation error is dominated by 1 N or is of the order of 1 N , respectively. Theorem 2.1: [START_REF] Mardia | Measures of multivariate skewness and kurtosis with applications[END_REF] Let X be an i.i.d. process of dimension

d. If X ∼ i.i.d N (0, S) then Bd (N ) ----→ N →∞ N (µ, σ 2 ), with: µ = d(d + 2) + o( 1 N ) (5)
σ 2 = 8d(d + 2) N + o( 1 N ) (6)
In [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF] we devised a similar theorem for multivariate data that are no longer considered i.i.d. In the following, we sketch the main assumptions to identify the statistics of Mardia's MK test in the multivariate n.i.d. case.

B. Joint normality test for a colored process

Proofs of the results presented in this subsection can be found in [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF], as well as theoretical details.

1) Mathematical Preliminaries: In order to guarantee convergence while relaxing the i.i.d assumption, we make the mixing assumption:

∞ τ =0 |S ab (τ )| 2 converges to a finite limit ∀(a, b) ∈ {1, . . . , d} 2 .
Let S = S + ∆, where ∆ is small compared to S, then we have the following lemma:

Lemma 2.1: The entries of matrix ∆ are of order O(1/ √ N ). Let G = S -1 . This lemma allows the following expansion of the precision matrix

G def = S -1
up to second order in ∆:

G = G -G∆G + G∆G∆G + o(1/N ). (7)
Finally, by injecting [START_REF] Moore | The effect of dependence on chi squared tests of fit[END_REF] in Bd (N ) defined in ( 3), an approximation of the test statistic in o(N -1 ) may be obtained, as Bd can be shown to converge to a normal variable. Thus we only need to calculate its mean and variance, whose expressions are given below.

2) Statistics for a colored scalar process:

E{ B1 } = 3 - 6 N - 12 N 2 N -1 τ =1 (N -τ ) S(τ ) 2 S 2 + o( 1 N ) (8)
Var{ B1 } = 24 N 1 + 2 N N -1 τ =1 (N -τ ) S(τ ) 4 S 4 + o( 1 N ) (9)
The dependence between time samples is taken into account in the covariance function terms S(τ ). In the case of i.i.d samples, we recover the expressions given in Theorem 2.1 above, for d = 1.

3) Statistics for a bivariate colored process:

E{ B2 } = 8 - 16 N - 4 N 2 N -1 τ =1 (N -τ )Q 1 (τ ) (S 11 S 22 -S 2 12 ) 2 + o( 1 N ) (10)
Var{ B2 } = 64 N + 16 N 2 N -1 τ =1 (N -τ )Q 2 (τ ) (S 11 S 22 -S 2 12 ) 4 + o( 1 N ) (11)
where Q i (τ ) are linear combinations of S ab (τ ) thus containing information about the time-dependence between samples. Despite the difficulty to derive such expressions, the computational burden involved by their implementation is low.

III. TIMES SERIES MODEL: MULTIDIMENSIONAL AUTO-REGRESSIVE MODEL

A. Definition

The VAR(p) model defined in (1) can be written in matrix format as:

Y = Zw + ϵ, (12)
where, ∀i ∈ {p, . . . , N -1}:

w = A 1 , A 2 , . . . , A p T ∈ R dp×d Y = x(p + 1), . . . , x(N ) T ∈ R (N -p)×d Z = z(p), z(p + 1), . . . , z(N -1) T ∈ R (N -p)×dp z(i) = x(i); x(i -1); . . . ; x(i -p + 1) ∈ R dp×1 ϵ = ϵ(p + 1), ϵ(p + 2), . . . , ϵ(N ) T ∈ R (N -p)×d
A Least Squares estimation can be performed to find ŵ that minimizes ∥Zw -Y ∥ 2 assuming that Γ = Z T Z is nonsingular:

ŵ = Γ -1 Z T Y (13)

B. Recursive Estimation of VAR(p) parameters

To estimate the regression parameters, a classical Recursive Least Squares methods (RLS) is proposed. This adaptive strategy is well known to be suitable for data showing some sort of non-stationarity, and has been widely used (see e.g. [START_REF] Plackett | Some theorems in least squares[END_REF] [17], [START_REF] Hallin | Mixed autoregressive-moving average multivariate processes with time-dependent coefficients[END_REF] and [START_REF] Grenier | Time-dependent arma modeling of nonstationary signals[END_REF] to only cite few).

Below, the main RLS equations are given for the present problem. Let λ 1 be the decay factor weighting the observations, and ŵk be the k th column of w defined in [START_REF] Hogg | Statistical robustness: One view of its use in applications today[END_REF], and

y k (t -1) = x k (1), x k (2), . . . , x k (t -1)
T .

y k (t -1) = Z ŵk (t -1) + ϵ k (t -1) (14) ŵk (t -1) = Γ -1 (t -1)Z T y k (t -1) (15)
Suppose we want to update the model with new observations x(t); a new row z(t) T is appended to Z in [START_REF] Hogg | Statistical robustness: One view of its use in applications today[END_REF], and a new observation x k (t) is appended to y k . In the RLS algorithm, Γ -1 (t) and ŵk (t) are recursively expressed for t > N as [START_REF] Hogg | Statistical robustness: One view of its use in applications today[END_REF] augmented with the row z(t) T , and ỹk = y k (t -1); x k (t) . Note that we can also vary λ 1 and the order p of the model to allow for more flexibility and adapt to local non-stationarities of the data, see e.g [START_REF] Giurcȃneanu | AR order selection in the case when the model parameters are estimated by forgetting factor leastsquares algorithms[END_REF], [START_REF] Cho | Time-frequency analysis using ar models with variable forgetting factors[END_REF]. This is not detailed in this communication.

Γ -1 (t) = λ 1 Γ(t -1) + z(t)z(t) T -1 (16) = λ -1 1 Γ -1 (t -1) -buu T (17) ŵk (t) = Γ -1 (t) ZT ỹk (18) = ŵk (t -1) -b(z T ( ŵk (t -1) + x k (t)u))u + x k (t)u, u = λ -1 1 Γ -1 (t -1)z(t) and b = (1 + z(t) T u) -1 . Note that Z is Z in

IV. CHANGE DETECTION ALGORITHM

A. The test in practice

First we choose the nominal level of the test:

α = P(choose H0 |H 0 is true). • It has been shown in previous sections that z = ( B d - E{ B d })/ Var{ B d } is asymptotically normal. • We reject the null hypothesis H 0 at a significance level α if: 2(1 -Φ(z)) < α
where Φ denotes the cumulative distribution function (cdf) of N (0, 1). • In practice, the time-structure is unknown and the covariance (matrix) function entries S ab (τ ) are replaced by their sample counterparts at time t.

B. Recursive normality test on regression residuals

The proposed test statistic is computationally efficient i.e it can be easily computed over a sliding window or by using the exponential weighting technique to test the normality of the estimated regression residuals available at time t assumed to follow a Gaussian distribution N (0, V (t)). Let 0 < λ 2 < 1:

V (t) = λ 1 V (t -1) + (1 -λ 1 )ε(t)ε(t) T (19) Bd (t) = λ 2 Bd (t -1) + (1 -λ 2 )(ε(t) T V -1 (t)ε(t)) 2 . ( 20
)
The algorithm for sequentially detecting changes reads:

Require: p ≥ 1, 0 < λ 1 , λ 2 < 1, α, δ Initialization: Γ(p) ← δI dp , Bd (p) ← 0, V (p) ← I d for p + 1 ≤ t ≤ N do Update Γ -1 (t) ▷ using (16) Update { A k (t)} 1≤k≤p ▷ using (18) Compute ε(t) = x(t) - p k=1 A k x(t -k) Compute z = ( Bd (t) -E{ Bd (t)})/ Var( Bd (t)) if 2(1 -Φ(z)) < α then Change is detected else No change
end if end for Forgetting factors λ 1 and λ 2 (for 2nd and 4th order statistics respectively) are usually chosen by a rule of thumb, depending on the time-scale of the change. To give a better intuition of this factor, one can calculate the length N eff of a uniform sliding window that would yield the same estimator variance, when data are i.i.d. normal. This leads to λ ∝ 1/N eff . More details will be provided in a subsequent paper.

For the RLS algorithm, a common choice is λ 1 = 0.99 which corresponds to using a virtual window's size N eff = 2 1-λ1 . As for the recursive estimation of the Kurtosis, to prevent the problem of slow convergence, we choose N eff = 1000.

V. COMPUTER RESULTS

A set of of Monte Carlo simulations is presented to compare the power of the proposed normality test when applied directly on data or on regression residuals. Then, we study the performance of the test statistic on a low-dimensional (2-D or 1-D) projection of the initial multivariate data. As a final illustration of the effectiveness of our method, we apply the change detection Algorithm presented in subsection IV-B on synthetic colored data undergoing an abrupt change in its distribution.

A. Applying the test on colored data 1) Directly on data: M = 2000 simulations are considered, each being based on a sequence of N = 1000 samples. First, we simulate 1-D AR(p) processes (for p ∈ {4, 14, 20}). The AR coefficients are computed such that the equivalent filter is low-pass, with band pass equal to .25 (normalized freq).

For each simulation a 2-D Gaussian (or Uniform) AR(p) process is constructed by time embedding : x(t) = {x(2t), x(2t + 1)}. Then, both Mardia's test derived for i.i.d. samples (denoted B1,i.i.d ) and the test for colored samples, whose statistics are defined by equations [START_REF] Hinich | Testing for Gaussianity and linearity of a stationary time series[END_REF][START_REF] Moulines | Testing that a multivariate stationary time series is Gaussian[END_REF] for colored samples (denoted B1 ) are applied on the marginals of the 2-D process, and compared. Finally, the statistics derived for n.i.d. bi-variate data B2 , described by equations [START_REF] Elbouch | A Normality Test for Multivariate Dependent Samples[END_REF][START_REF] Jarque | A test for normality of observations and regression residuals[END_REF] is also applied to the 2-D process.

The obtained empirical rejection rates of hypothesis H 0 are computed as #Rejections M , for a test level α = 5%. The results are summarized in Tables I andII. of their variance in [START_REF] Hinkley | Inference about the change-point from cumulative sum tests[END_REF][START_REF] Moulines | Testing that a multivariate stationary time series is Gaussian[END_REF], we can see that that the variance of B1,i.i.d is underestimated, consequently, the latter over-rejects the hypothesis of Gaussianity. This is noticeable even for Gaussian AR(p) process (Under H 0 ), with rejection rates that surpass the nominal level 5%. • Both scalar tests B1,i.i.d and B1 perform poorly compared to the joint normality test statistic B2 . • When the correlation tails last longer (Table II), the overall performance of the test statistics tends to decrease. 2) On Regression residuals: We now study the power of the normality tests on estimated regression residuals. We utilize the Least Squares method to obtain ε. The simulation procedure and the tests are the same as those described in the previous paragraph V-A1. We study the case where the order p of the generated AR(p) process is known (we choose p = 20), and the case where the order is misspecified (p = 9). The empirical rejection rates are summarized in Table III. • If the model is perfectly known, then all test statistics perform well on the well estimated regression residuals. However, if the model order is under estimated, some important time correlation remain, and scalar tests perform poorly compared to the joint normality test, as expected.

B. Random Projections

We simulate a 3-D process VAR(p), p ∈ {5, 20}, of length N = 1000 following equation [START_REF] Basseville | Detecting changes in signals and systems-a survey[END_REF], where the inputs ϵ(t) are i.i.d with distributions either multivariate standard normal N (0, 1) or multivariate U(-2, 2).

Then M = 2000 different projections on an arbitrary plane (2-D projection) going through the origin of the 3-D space are computed, corresponding to as many 2-D time series. For comparison, the same set of observations is also projected M times on an arbitrary direction (1-D projection). Eventually, we run the same set of experiments on estimated regression residuals estimated by ordinary least squares method (OLS). The results are reported below.

3-D

On residuals VAR [START_REF] Giurcȃneanu | AR order selection in the case when the model parameters are estimated by forgetting factor leastsquares algorithms[END_REF] p = 10 1-D proj., B1 0.250 0.41 2-D proj., B2 0.580 0.9

• The test B1 applied directly on 1-D projections of either the observations or its residuals computed by OLS, performs poorly. This is in accordance with our observations from the preceding experiments. • Even with a misspecified order, the joint test statistic performs best (empirical power of 90%) on the 2-D low representation of regression residuals, as it is able to account for both temporal and spatial (between coordinates) dependences.

C. Change detection on synthetic data

Consider the case where a process ϵ is constituted of i.i.d samples following a standard normal distribution. The process undergoes a change at n c = 5000 in its distribution: samples in the interval [5000, 10000] are now following a uniform distribution U(-√ 3, √ 3). The change ends at n = 10000 and the samples are again normally distributed. The process ϵ is then (low pass) filtered using an AR(5) model and is now denoted x.

One realization of this process is given in Fig. 1. It is clear from this figure that the change in the distribution is unbeknownst to the human eye.

The change detection algorithm presented in subsection IV-B is applied to this realization by setting p = 5, λ 1 = 0.99, λ 2 = 0.998, α = 5% and δ = 1. A 2-D process is obtained by taking x(t) = {x(2t), x(2t + 1)}.

For comparison, the CUSUM algorithm [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF] is applied on the regression residuals by using its recursive form. The instantaneous log-likelihood ratio is computed as: In fact, the likelihood ratio in ( 21) is derived under the hypothesis that the process ϵ is i.i.d. As the latter's values are estimated recursively, they are more likely to have residual correlations between them and the hypothesis of independence no longer holds, explaining why the cumulative sum of L(t) keeps increasing.

L(t) = -ln(2 √ 3) + 1 2 ln(2π) + 1 2 ε2 (t) (21)

VI. CONCLUDING REMARKS

A novel approach to testing the Gaussianity of colored data is proposed. In a first stage, data are whitened recursively.

Then the MK test is applied in a second stage on available regression residuals in an online manner. Computer experiments have been performed to assess the impact of the preprocessing stage on the power of the test. They have evidenced fact that the bivariate MK statistics we recently obtained (assuming both spatial and temporal dependence) perform always better than their scalar counterparts. The latter are indeed severely affected by serial correlation in residuals due to imperfect whitening in the first stage. The same conclusions can be drawn if data are randomly projected onto a low-dimensional subspace.

Future work will concern the validation of the change detection algorithm on real-world phenomena recorded by a network of d-variate sensors within a multiple hypothesis testing framework.

Fig. 1 .

 1 Fig. 1. One realization of a Gaussian AR(5) process that undergoes an abrupt change in the distribution of its excitation ϵ (from N (0, 1) to U (-√ 3, √3)). Affected samples are between red dashed lines.

Fig. 2 .

 2 Fig. 2. Evolution of the normalized test statistics B1 (in blue) and B2 (in orange). In red, the evolution of the cumulative sum s(n) = n t=1 L(t). Black horizontal dashed lines are the critical values ±1.96 corresponding to a test power of 95%. Red vertical dashed lines are the beginning and end of an abrupt change in the excitation statistics of an AR(5) process.

TABLE I EMPIRICAL

 I REJECTION RATES FOR 2000 SIMULATIONS FOR α = 5% SIGNIFICANCE LEVEL WITH THE B1,i.i.d , B1 , B2 TEST APPLIED DIRECTLY ON AR(p) DATA WITH p = 4, 14

			AR(4)		AR(14)
	Test Statistic	Gaussian Uniform Gaussian Uniform
	B1,iid	0.067	1.	0.123	0.512
	B1	0.052	0.99	0.045	0.456
	B2	0.06	1.	0.065	0.88

• It seems that B1,i.i.d has a better detection power than B1 . As a matter of fact, by comparing the formulas

TABLE II EMPIRICAL

 II REJECTION RATES FOR 2000 SIMULATIONS FOR α = 5% SIGNIFICANCE LEVEL WITH THE B1,i.i.d , B1 , B2 TEST APPLIED DIRECTLY ON AR(20) PROCESS

		AR(20)
	Test Statistic Gaussian	Uniform
	B1,iid	0.228	0.430
	B1	0.047	0.399
	B2	0.06	0.688

TABLE III EMPIRICAL

 III REJECTION RATES FOR 2000 SIMULATIONS FOR α = 5% SIGNIFICANCE LEVEL WITH THE TEST STATISTICS B1,i.i.d , B1 , B2 APPLIED ON ESTIMATED REGRESSION RESIDUALS USING OLS METHOD,

		WITH KNOWN AND MISSPECIFIED ORDER
		On Residuals of AR(20) On Residuals of AR(20)
			p = 20		p = 9
	Test Statistic	Gaussian Uniform	Gaussian	Uniform
	B1,iid	0.06	1.	0.064	0.582
	B1	0.05	1.	0.051	0.429
	B2	0.055	1.	0.06	0.850


Wolfram Kloppmann

Lise Leroux

Philippe Bromblet

Sophie Jugie

Pierre-Yves Le Pogam

Aleksandra Lipinksa

Alabaster, a millenary stone of European cultural heritage

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Powered by TCPDF (www.tcpdf.org)

Fig. 1

 1 Fig. 1 Principal alabaster deposits used for sculpture in Late Medieval and Early Modern times. Transdisciplinary network of the Materi-A-Net project.

Franco-German Materi-A-Net project co-funded by ANR and DFG.