François Caron

Chris Holmes

Emmanuel Rio

On the sampling distribution of an 2 distance between Empirical Distribution Functions with applications to nonparametric testing

Keywords: two-sample test, nonparametric test, binary tree

We consider a situation where two sample sets of independent real valued observations are obtained from unknown distributions. Under a null hypothesis that the distributions are equal, it is well known that the sample variation of the infinity norm, maximum, distance between the two empirical distribution functions has as asymptotic density of standard form independent of the unknown distribution. This result underpins the popular two-sample Kolmogorov-Smirnov test. In this article we show that other distance metrics exist for which the asymptotic sampling distribution is also available in standard form. In particular we describe a weighted squared-distance metric derived from a binary recursion of the real line which is shown to follow a sum of chi-squared random variables. This motivates a nonparametric test based on the average divergence rather than the maximum, which we demonstrate exhibits greater sensitivity to changes in scale and tail characteristics when the distributions are unequal, while maintaining power for changes in central location.

Sur la distribution d'une distance 2 entre fonctions de distributions empiriques, avec application au test non paramétrique d'adéquation entre les distributions de deux échantillons.

Résumé : On considère la situation où un ensemble de deux échantillons {X 1 , . . . , X n (1) }, {Y 1 , . . . , Y n (2) }, d'observations réelles indépendantes est obtenu à partir de distributions inconnues {F X , F Y }, avec les fonctions de répartition empiriques { FX , FY }. Dans ce cas, il est connu que sous l'hypothèse nulle F X ≡ F Y la distance maximum ∞ || FX -FY || ∞ , a une distribution asymptotique de forme standard indépendantes de F . Ce résultat forme la base du test d'adéquation de Kolmogorov-Smirnov. Dans cet article, nous montrons que d'autres distances entre fonctions de répartitions existent pour lesquelles la distribution d'échantillonnage asymptotique est connue. En particulier, nous proposons une distance 2 pondérée dérivée d'une récursion binaire de R qui suit asymptotiquement la même distribution qu'une somme de variables aléatoires χ 2 . Ce résultat suggère un test non paramétrique basé sur la divergence moyenne plutôt que le maximum pour lequel nous montrons une plus grande sensibilité à des changements de variance et de queues de distribution, tout en conservant des résultats similaires pour détecter des changements du paramètre de location.

Mots-clés : test d'adéquation non paramétrique, récursion binaire 1 Introduction

The sampling distribution of the empirical distribution function defined as Fn (x) = k/n for real-valued x (k) < x < x (k+1) , where x (i) denotes the realised ith order statistic, is an important quantity in statistics in part as it allows for nonparametric hypothesis testing. For instance when two sets of samples have been obtained under different conditions {X 1 , . . . , X n (1) }, {Y 1 , . . . , Y n (2) } it is well known from the work of [START_REF] Kolmogorov | Sulla determinazione empirica di una legge di distribuzione[END_REF] that for a particular ∞ norm the distance

d ∞ ( FX , FY ) = || FX -FY || ∞ = max x∈R | FX (x) -FY (x)| has an asymptotic distribution of standard form independent of F X , F Y under the null hypothesis X i , Y i ∼ F X . This result
underpins the popular two-sample Kolmogorov-Smirnov test. In this paper we explore the sample distribution of other distance metrics d( FX , FY ). In particular we derive the distribution of an 2 distance calculated from a recursive binary partition of the real line R. We show that this statistic is independent of F X , F Y , and follows a weighted sum of χ 2 under the null hypothesis

X i , Y i ∼ F X .
This motivates a nonparametric test based on the weighted average divergence d w ( FX , FY ) rather than the maximum. We demonstrate numerically that this test appears to show greater sensitivity in detecting departures in the scale or tails of F X , F Y , while maintaining power to detect shifts in central location.

The rest of the paper is as follows. In Section 2 we outline our partitioning scheme of R and derive the asymptotic sampling distribution of the 2 distance under the null hypothesis. In Section 3 we provide some illustrations that the test can be more powerful than other conventional non-parametric test. In Section 4 we derive a similar test for k-sample test. Finally in Section 5

we provide a discussion of our work.

2 An 2 distance for two-sample test Suppose one has two sets of samples, {X, Y } of sample size {n (1) , n (2) } obtained say under different treatments. We write n = n (1) + n (2) . Let F0 (x) be the empirical distribution of the joint sample {X, Y }.

We shall first consider a dyadic (binary) tree that recursively partitions R into disjoint measurable sets such that at the mth level of the tree we find

R = ∪ 2 m -1 j=0 B (m) j where B (m) i ∩B (m) j = ∅
for all i = j. The jth junction in the tree at level i has associated set

B (i) j
and clearly

B (i) j = (B (i+1) 2j , B (i+1)
2j+1 ) for all i, j. It will be convenient in what follows to simply index the sets using base 2 subscript and drop the superscript so that, for example, B 000 indicates the first set in level 3, B 0011 the forth set in level 4 and so on. Such a recursive binary tree is represented in Figure 1. Let Π denote the partition structure defined by the collection of sets Π = (B 0 , B 1 , B 00 , . . .) and in particular consider a partition centered on the quantiles of F0 .

That is, at level m, where j * is the decimal representation of the binary number j and m is the level of the tree.

B j = [ F -1 0 (j * /2 m ), F -1 0 ([j * + 1]/2 m ))], (1)
We write resp. n

(1) j

and n

(2) j the number of observations in X and Y that lie in the subset B j of Π. We note n (12) j

= n

(1)

j + n (2) j . Hence n (1)
0 is the number of items from X in the first 50% percentiles of F0 , n

(1) 00 is the number of items from X in the first 25% percentiles, etc. We also have n

(k) j = n (k) j0 + n (k)
j1 , for k = 1, 2, 12, where 12 denotes the combined sample set.

Consider the level j with given n (1) j ) whose probability function is given by Pr(n

(1) j0 = k|Π, n (1) j , n (12) j , n (12) j0 , H 0 ) = n (1) j k n (12) j -n (1) j n (12) j0 -k n (12) j n (12) j0 (2) if max(0, n (1) j -n (12) j + n (12) j0 ) ≤ k ≤ min(n (12) j0 , n
(1) j ) and 0 otherwise. n p is the usual binomial coefficient. The hypergeometric distribution is obtained when sampling without replacement n (12) j0 balls from an urn containing n

(1) j red balls and n (12) j -n

(1) j white balls. The conditional moments are then given by

µ j = E[n (1) j0 |n (1) j ] = n (1) j n (12) j0 n (12) j σ 2 j = E n (1) j0 -µ j 2 |n (1) j = n (1) j (n (12) j -n (1) j )n (12) j0 (n (12) j -n (12) j0 ) n (12) j 2 (n (12) j -1)
Definition 1 Let T be the weighted 2 distance defined by

T = j λ j z j ( 3
)
where

z j = n (1) j0 -µ j 2
and the sum goes over all elements j such that n

(12) j > 0. Let log 2 (a) =
log(a)/ log 2 and l(j) be the level of the tree associated to the binary number j. For all j with l(j) ≤ [log 2 (n)] -1, the weights λ j are defined by

λ j = n (12) (n (12) -1) 2 l(j) p j (1 -p j )n (12) j n (1) n (2) (4)
where

p j = n (12) j0 /n (12) j
, l(j) is the level of section j, and

λ j = 0 for l(j) ≥ [log 2 (n)].
The value z j measures the departure from the expected value at each level j, and there are 2 m such contributions at each level m of the binary tree. The realisation of T can then be used as a test statistic to quantify departures from the null hypothesis.

We now provide a result on the asymptotic distribution of the test statistic T under the null.

Proposition 1 Let W = j 2 -l(j) Y 2 j where Y j ∼ N (0, 1) for all j. Then inf{ T -V 2 : V Law = W } ≤ c 1 (n (1) n (2) ) -1 n 3/2 log 2 (n), (5)
where c 1 is some positive constant. It provides the rate of approximation O(n -1/2 log n) with respect to the minimal L 2 -distance as n tends to infinity, and consequently the rate of approximation

O(n -1/3 (log n) 2/3
) with respect to the Prokhorov distance. The proof is given in Appendix A.1.

Hence we can use T to test for the hypothesis F X ≡ F Y . Note we do not have the asymptotic distribution of the test statistic T under the alternative. Nonetheless, the next theorem provides a lower bound for T under H 1 .

Proposition 2 Assume that the ratio n 1 /n converges to α ∈]0, 1[ when n tends to ∞. Then under the alternative F X = F Y , we have

lim inf n→∞ n -1 T ≥ c 2 α 1 -α almost surely (6)
where c 2 > 0 is some strictly positive constant. The proof is given in Appendix A.2

Simulations

In Section 2 we showed how the test statistic T can be used to measure departures from the null hypothesis that the underlying distribution functions are the same. To examine the operating performance of the method we consider the following experiments designed to explore various canonical departures from the null.

a) Mean shift: Y (1) ∼ N (0, 1), Y (2) ∼ N (θ, 1), θ = 0, . . . , 3 b) Variance shift: Y (1) ∼ N (0, 1), Y (2) ∼ N (0, θ 2 ), θ = 1, . . . , 3 c) Mixture: Y (1) ∼ N (0, 1), Y (2) ∼ 1 2 N (θ, 1) + 1 2 N (-θ, 1), θ = 0, . . . , 3 d) Tails: Y (1) ∼ N (0, 1), Y (2) ∼ t(θ -1 ), θ = 10 -3 , . . . , 10
e) Skew: Y (1) ∼ N (0, 1), Y (2) ∼ SN (0, 1, θ), θ = 1, . . . , 10
f ) Lognormal mean shift: log Y (1) ∼ N (0, 1), log Y (2) ∼ N (θ, 1), θ = 0, . . . , 3
g) Lognormal variance shift: log Y (1) ∼ N (0, 1), log Y (2) ∼ N (0, θ 2 ), θ = 1, . . . , 3
where SN (0, 1, λ) is the skew normal distribution of skewness parameter λ. Comparisons are performed with n (1) = n (2) = 50 against the two-sample Kolmogorov-Smirnov and Wilcoxon rank test. To compare the models we explore the "power to detect the alternative". We evaluate numerically the value θ such that the statistical power is at least 80% for each of the three tests. The results, obtained from 10 6 independent simulations, are reported in Table 1. Globally, our approach gives similar results when detecting a shift in the median, but outperforms other methods for detecting shifts in variance or tails.

Extension to multiple sample test

Our partitioning approach can be readily extended to deal with multiple tests for multiple treatments or conditions. Consider now that we are given p samples y (1) ∼ i.i.d. p) , F (1) , . . . , F (p) unknown, and we wish to test the following hypothesis against the alternative that at least one of the distributions is different from the others. We consider the binary tree Π = (B 0 , B 1 , B 00 , . . .) constructed from the empirical distribution of y = (y (1) , . . . , y (k) ). We note n (k) j , j = {}, 0, 1, . . ., k = 1, . . . , p the number of items from y (k) in interval B j , and

F (1) , y (2) ∼ i.i.d. F (2) , . . . , y (p) ∼ i.i.d. F (
H 0 : F (1) = F (2) = . . . = F (p) RR n°7931
n j = p k=1 n (k) j . Assuming n (k) j
> 0 for all k, then the vector

z j0 = (n (1) j0 , n (2) j0 , . . . , n (p-1) j0
) follows the multivariate hypergeometric distribution MvHypGeo((n

(1) j , n (2) j , . . . , n (p) j ), n j0 ) whose pdf is given by Π p k=1 n (k) j n (k) j0 nj nj0
This distribution is obtained when sampling without replacement from an urn with n (k) j balls of color k, k = 1, . . . , p. The two first moments µ j and Σ j are defined by

µ j = E [z j0 |z j ] = n (12) j0 n (12) j n (1) j , n (2)
j , . . . , n

(p-1) j Σ j (k, l) = cov(n (k) j0 , n (l)
j0 |z j ) = - n j0 n (k) j n (l) j n 2 j n j -n j0 n j -1 Σ j (k, k) = var(n (k) j0 |z j ) = n (k) j n j 1 - n (k) j n j n j n j -n j0 n j -1
and we can show that z j0 is asymptotically multivariate normal [START_REF] Fraser | A vector form of the Wald-Wolfowitz-Hoeffding theorem[END_REF] and so is asymptotically X 2 (p -1) distributed. We extend the test statistic T introduced in the above section by considering the similar test statistic T 2 defined by

S j = (z j0 -µ j ) T Σ -1 j (z j0 -µ j ) = n j -1 n j -n j0 p k=1 n (k) j0 -n (k) j nj0 nj 2 n (k) j nj0 nj RR n°7931
T 2 = j|#{k|n (k) j >0}≥2 2 -l(j) S j (7)
The same experiments as in Section 3 are performed. We consider a 4-sample test where we define F 3 = F 4 = F 1 . Our test is compared to the Kolmogorov-Smirnov test with pairwise comparisons, where the threshold value is set so that to have overall size α = 0.05. We also compare the test to the Kruskall-Wallis test. Results are given in Table 2. Similar to the twosample test, our test performs similarly to detect a shift in the median but outperforms other methods to detect a shift in the variance or tails.

Discussion

We have derived the sampling distribution for a weighted 2 distance between empirical sampling distributions drawn from the same probability law. This allowed us to derive a simple nonparametric test statistic for departures from the null hypothesis of no treatment effect. This method offers an alternative nonparametric procedure to the usual Kolmogorov-Smirnov two sample test. The Kolmogorov-Smirnov test is known to have relatively low power for against alternative which differ in scale [START_REF] Capon | On the asymptotic efficiency of the Kolmogorov-Smirnov test[END_REF][START_REF] Klotz | Asymptotic efficiency of the two sample Kolmogorov-Smirnov[END_REF]. Simulations conducted in Section 3 and 4 show that the new binary tree test has similar power against alternatives which differ in location while largely outperforming Kolmogorov-Smirnov test in scale/tail alternatives.

A Proofs

A.1 Upper bounds for the minimal L 2 -distance of T to the limiting law under H 0

In this subsection we prove Proposition 1. The first step is to approximate the test statistic T by a weighted sum of independent χ 2 (1)-distributed random variables. This will be done using the Tusnády's type lemma proved in Castelle and Laurent-Bonvalot (1998), which we now recall.

Lemma 1 (Lemma 2.5 in Castelle and Laurent-Bonvalot ( 1998)) -Let Y be a standard normal random variable and Φ be the distribution function of Y . Let G n,n1,n2 be the distribution function of the hypergeometric distribution H(n, n 1 , n 2 ). We denote by m the mean of H(n, n 1 , n 2 ) (recall

m = n 1 n 2 /n). Let p = n 1 /n, p = n 2 /n and δ = 2p -1, δ = 2p -1. Set σ = npp (1 -p)(1 -p ). ( 8
)
Then, for each positive η, there exists positive constants c and d such that, if |δδ | ≤ 1 -η, then

|G -1 n,n1,n2 • Φ(Y ) -m -σY | ≤ c + dY 2 . ( 9
)
Starting from Lemma 1, we now approximate the statistic T in L 2 by a weighted sum of squares of standard normal random variables. To achieve this approximation, we construct a random variable with the same distribution as T from a family of independent standard normal random variables indexed by the binary tree, in the same way as in [START_REF] Komlós | An approximation of partial sums of independent RV's and the sample DF i[END_REF]. So let B denote the binary tree and let (Y j ) j∈B be a collection of independent standard normal random variables. The random variables (n

(12) j
) j∈B are given (and the collection of Gaussian r.v.'s is independent of this collection). Suppose that the random variables n

(1) j , with the adequate multinomial distribution at each level (and consequently n

(2) j )), have already been defined up to level k from the Gaussian random variables (Y j ) l(j)<k . We want to define the random variables at the level k + 1. By definition (1),

n (12) j0 = [n (12) j /2]. ( 10
) At level k + 1, we set n (1) j0 = G -1 n 12 j ,n (12) j0 ,n (1) j • Φ(Y j )
Let m j denote the mean of n .

Set

s j = n (12) j p j (1 -p j )p j (1 -p j ).

RR n°7931

By Eq. ( 10),

δ j = 2p j -1 = 0 if n (12) j
is even, and

|δ j | = |2p j -1| = 1/n (12) j otherwise. It follows that |δ j | ≤ 1/3, provided that n (12) j
≥ 2. In that case, Lemma 1 applies with η = 1/3, and

|n (1) j0 -m j -s j Y j | ≤ c + dY 2 j . (11)
Let us now denote by n = n (1) + n (2) the global size of the sample, and write n in basis 2, that

is n = a 1 a 2 • • • a l with a 1 = 1.
We may without loss of generality assume that n (1) ≤ n (2) . At the level zero,

n = n (12) 1 = a 1 a 2 • • • a l and, at level l -2, n (12) j
≥ 2.

Hence we will stop the construction at level l -2 (note that [log 2 (n)] = l -1), so that all the numbers n (12) j satisfy n

(12) j ≥ 2 . The test statistic T is defined by

T = l-2 k=0 T k , where T k = {j:l(j)=k} λ j |n (1) j0 -m j | 2
for positive weights λ j to be specified later.

Let

T = l-2 k=0 T k , where T k = {j:l(j)=k} λ j σ 2 j Y 2 j .
Clearly

T -T 2 ≤ l-2 k=0 T k -T k 2 . Now T k -T k = {j:l(j)=k} λ j (|n (1) j0 -m j | 2 -σ 2 j Y 2 j ).
Recall that we have to bound up the L 2 -norm of this random variable. Now

E(|n

(1) j0 -m j | 2 | F k ) = σ 2 j .
Consequently

M k = E(T k -T k | F k ) = 0.
Next we bound up the conditional variance of (T k -T k ) given F k . From our construction the random variables (|n

(1) j0 -m j | 2 -σ 2 j Y 2 j )
j are independent at the scale l(j) = k conditionally to the σ-field F k generated by the random variables (n (1) j ) with l(j) = k. Hence the conditional variance is the sum over j of individual conditional variances. By the inequality (11), we have (|n

(1) j0 -m j | 2 -σ 2 j Y 2 j ) 2 ≤ (c + dY 2 j + (σ j -s j )|Y j |) 2 (c + dY 2 j + (σ j + s j )|Y j |) 2 . RR n°7931 Now σ j -s j = s j n (12) j n (12) j -1 -1 ≤ s j 2(n (12) j -1) ≤ s j n (12) j ≤ 1 4 .
Hence the conditional variance of T k -T k given F k , which we denote by V k , satisfies

V k ≤ c 2 0 {j:l(j)=k} λ 2 j (1 + σ 2 j ),
for some positive universal constant c 0 . Note that σ 2 j ≤ 2 -3 n (12) j

. Also, from the definition of

n (12) j , n2 -l(j)-1 ≤ n (12) j ≤ n2 1-l(j) ,
which ensures that σ 2 j ≤ n2 -1-l(j) . Since n2 -1-l(j) ≥ 1 for l(j) ≤ l -2, we get the deterministic upper bound

V k ≤ c 2 0 n {j:l(j)=k} λ 2 j 2 -k .
From the above bounds and the fact that T k -T k 2 2 = E(V k ), we now have:

T -T 2 ≤ c 0 √ n l-2 k=0 2 -k {j:l(j)=k} λ 2 j 1/2 . ( 12
)
Now, let

W = {j:l(j)<l-1} λ j E(σ 2 j )Y 2 j .
Clearly W is a weighted sum of χ 2 (1)-distributed independent random variables. We first give an explicit formula for W .

By definition, n

(1) j has the hypergeometric distribution H(n, n 1 , n

). Hence

E(σ 2 j ) = n (12) j0 n (12) j1 n (12) j (n (12) j -1) E(n (1) j n
(2)

j ) = n 1 n 2 n(n -1) p j (1 -p j )n (12) j
, which leads to the definition formula

W = n 1 n 2 n(n -1) {j:l(j)<l-1} λ j p j (1 -p j )n (12) j Y 2 j (13) Now we bound up T -W 2 . Clearly T -W 2 ≤ l-2 k=0 {j:l(j)=k} λ j (σ 2 j -E(σ 2 j ))Y 2 j 2 . Now σ 2 j -E(σ 2 j ) = p j (1 -p j ) n (1) j n (2) j -E(n (1) j n (2) j ) n (12) j Recall n (2) j = n (12) j -n (1) j . Therefrom T -W 2 ≤ l-2 k=0 ( A k 2 + B k 2 ), with A k = {j:l(j)=k} λ j p j (1 -p j )n (12) j n (12) j -1 (n (1) j -E(n (1) j )Y 2 j and B k = {j:l(j)=k} λ j p j (1 -p j )n (12) j n (12) j -1 (n (1) j p j -E(n (1) j p j )Y 2 j .

Now (n

(1) j ) j:l(j)=k depends on the random variables (Y k ) l(k)<l(j) , which ensures that this random vector is independent of the Gaussian random vector (Y j ) j:l(j)=k . Furthermore, it can be proven that (n (1) j ) l(j)=k is a negatively associated random vector. It follows that, for distincts j and j with l(j) = k, cov (n (1) j , n

(1) j ) ≤ 0 and cov (n

(1) j p j , n (1) j p j ) ≤ 0. Consequently A k 2 2 ≤ 3 {j:l(j)=k} λ 2 j p j (1 -p j )n (12) j n (12) j -1 2 Var n (1) j ≤ 3 4 {j:l(j)=k} λ 2 j Var n (
j p j ) ≤ 4Var n (1) j . It follows that T -W 2 ≤ 3 2 √ 3 l-2 k=0 {j:l(j)=k} λ 2 j Var n (1) j 1/2 . Since n (1) j has the hypergeometric distribution H(n, n 1 , n (1)
, we have:

Var n

(1)

j = n 1 n 2 n (12) j (n -n (12) j n 2 (n -1) ≤ n (12) j n 1 n ≤ 2 1-k n 1 ≤ 2 -k n, since n 1 ≤ (n/2). Hence T -W 2 ≤ 3 √ n l-2 k=0 2 -k {j:l(j)=k} λ 2 j 1/2
.

RR n°7931

Both the above bounds and ( 12) then imply that

T -W 2 ≤ (c 0 + 3) √ n l-2 k=0 2 -k j:l(j)=k λ 2 j 1/2 . ( 14
)
We now choose the coefficients λ j in such a way that it will be feasible to approximate the density of W . To this aim, we define λ j by the equation

λ j p j (1 -p j ) = n(n -1) 2 l(j) n (12) j n 1 n 2
For this choice of λ j , W has the same distribution as

W = l-2 k=0 2 -k Z k ,
where (Z k ) k≥0 is a sequence of independent random variables with respective laws χ 2 (2 k ). Now p j belongs to [1/3, 2/3], from the definition of n (12) j0 . Hence p j (1 -p j ) ≥ (2/9), and consequently

λ j ≤ 9n/(n 1 n 2 ).
Setting c 1 = 9(c 0 + 3), applying ( 14) and noticing that l

-1 ≤ log 2 (n) = log(n)/ log 2, we get that T -W 2 ≤ c 1 (n 1 n 2 ) -1 n 3/2 log 2 (n), (15)
which provides the rate of approximation O(n -1/2 log n) with respect to the minimal L 2 -distance as n tends to infinity, and consequently the rate of approximation O(n -1/3 (log n) 2/3 ) with respect to the Prokhorov distance.

A.2 Lower bounds for T under H 1

In this subsection, we give some lower bounds for the orders of magnitude of T under H 1 .

Throughout we assume that n tends to ∞ and (n 1 /n) converges to α in ]0, 1[. Let F = αF X +

(1 -α)F Y . We set

G 1 = F X • F -1 and G 2 = F Y • F -1 .
By the strong law of large numbers, for any fixed j at a fixed level m,

lim n→∞ n -1 (n 1 j0 -µ j ) = α(G 1 (2 -m-1 (2j + 1)) -(G 1 (2 -m j) + G 1 (2 -m (j + 1))/2).
Now lim n→∞ p j = 1/2. Moreover, it can be proven that lim n→∞ n -1 2 l(j) n (12) j = 1.

RR n°7931

The above facts ensure that the weights λ j satisfy lim n→∞ nλ j = 4/(α -α 2 ).

Consequently, a.s. lim n→∞ n -1 λ j z j = 4α 1 -α (G 1 (2 -m-1 (2j + 1)) -(G 1 (2 -m j) + G 1 (2 -m (j + 1))/2) 2

In order to write the statistic T in a different way, we now assume that F (1) and F (2) have a density w.r.t. the Lebesgue measure. Then G 1 has a locally integrable density g 1 on [0, 1].

We now recall the definition of the Haar system in L 2 ([0, 1]). Let e 0 = 1, and, at, for m ≥ 0, at scale m, for j = 0, 1, . . . 2 m -1, the functions ẽm,j be defined by ẽm,j (x) = 1 for x ∈]j2 -m , (j + 1/2)2 -m ], ẽm,j (x) = -1 for x ∈](j + 1/2)2 -m , (j + 1)2 -m ],

and ẽm,j (x) = 0 otherwise.

With this definition, lim n→∞ n -1 λ j z j = α 1 -α 1 0 (g 1 (u) -1)ẽ m,j (u)du (here M is a fixed integer such that A M > 0 and consequently the term on right hand is positive and does not depend on n).

RR n°7931

  j ) 2 is a 2-Lipschitz function of p j . Consequently Var ((p j ) 2 ) ≤ 4Var (p j ) and therefrom Var (n

  ẽ m,j (u)du. Since the Haar system is a total system in L 2 ([0, 1]), for any function g with null integral over [(u) -1) 2 du. Now, under H 1 , G 1 (x) = x on some subset of [0, 1] with positive Lebesgue measure,

Table 1 :

 1 Value θ needed to achieve 80% power for the binary tree, the Kolmogorov-Smirnov and the Wilcoxon two-sample tests. Lower values indicate better performances.

		Tree Test K-S Wilcoxon
	Gaussian: mean	0.67	0.67	0.58
	Gaussian: variance	2.19	2.73	-
	Gaussian: mixture	1.42	1.64	-
	Gaussian: skewness	1.13	1.14	0.94
	Gaussian:tail	2.02	2.97	-
	Lognormal: mean	0.69	0.68	0.58
	Lognormal: variance	2.20	2.78	-

Table 2 :

 2 Value θ needed to achieve 80% power for the binary tree, the Kolmogorov-Smirnov and the Kruskal-Wallis four-sample tests. Lower values indicate better performances.

		Tree Test K-S Kruskal-Wallis
	Gaussian: mean	0.67	0.63	0.57
	Gaussian: variance	2.25	3	-
	Gaussian: mixture	1.44	1.70	-
	Gaussian: skewness	1.01	1.05	0.83
	Gaussian:tail	2.02	3.16	-
	Lognormal: mean	0.67	0.63	0.56
	Lognormal: variance	2.27	3.12	-

RR n°7931

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Zoran Tiganj

Mamadou Mboup

Sylvain Chevallier

Emmanuel Kalunga

Online frequency band estimation and change-point detection

Detecting change-points in time series is at the heart of numerous applications, as abrupt changes in signal properties are quite common in natural and industrial processes.

In this paper we propose an algorithm for a particular changepoint detection problem where the frequency band of the signal changes at some points in the time axis. Apart from detecting the change-points, the proposed algorithm is also able to estimate the frequency bands. The main idea of the algorithm is to consider a simple local bandlimited model to represent the input signal in each sliding time window. The local model consists of a sum of two exponentials that, in the frequency domain, reads as the transfer function of a second order bandpass filter. Relying on the operational calculus, we obtain an explicit estimation of the parameters, that indicate the cut-off frequencies associated to each time window. The implementation is done in discrete time domain where the cut-off frequencies are computed as a combination of outputs of several FIR filters, providing a low computational cost online estimation. Experimental results show the efficiency and the stability of this algorithm even in presence of a moderate amount of noise.

I. INTRODUCTION

Change-point is a term commonly used to refer to the points in time where the properties of otherwise stationary process abruptly change. Processes with abrupt changes are common in natural and industrial systems. To illustrate the change-point detection problem we mention electrocardiogram (ECG) as a good example, since it is characterized with specific deflections which differ significantly from the rest of the signal. Another example is electroencephalogram (EEG) during epileptic seizures, which contains so-called spike and wave waveforms.

Detection of such change-points stands as an important research problem. Depending on the particular application, sometimes the detection is done after the recoding is finished (offline detection), but often the detection is required during the recording (online detection). A large amount of literature is devoted to the problem, in its general form. We refer to [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF] for an introductive presentation containing many application examples and also [START_REF] Mboup | A Volterra filter for neuronal spike detection[END_REF]- [START_REF] Tiganj | An algebraic method for eye blink artifacts detection in single channel eeg recordings[END_REF] for studies more related to the work presented here.

In this paper we address a very particular change-point detection problem. We consider a piecewise noisy input signal, such that each piece is characterized with several frequency components, belonging to some frequency band. The band changes from one piece to its neighbor so that two consecutive pieces always occupy different frequency bands. Here we propose an algorithm for online estimation of the cut-off frequencies of the bands. This also gives us an estimation of the time instants when the bands change.

Potential applications of such algorithm, same as for the change-point detection algorithms in general, are very diverse. They span form industrial to biomedical applications. For example, the epileptic spikes and wave complex mentioned above are characterized by an abrupt change in frequency. This change is assumed to be caused by sudden activation of a very large number of neurons whose synchronous cumulative activity brings new frequency components into the EEG spectrum [START_REF] Litt | Prediction of epileptic seizures[END_REF]. Fast detection of epileptic attacks is important to reduce the possible consequences for the patient, thus the online performance of the algorithm is often required. Another example based on the EEG signal is related with braincomputer interface (BCI) applications, where so called event related desynchronization/synchronization (ERD/ERS) caused by imagination of a movement is used to give a command to the computer. The ERD/ERS phenomenon is also caused by simultaneous activation of large neural circuits in the motor cortex during imagination of a movement of some body part [START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF]. The manifestation of the phenomenon is again the change of several frequency components in the signal. Since BCI applications are always in real time, the online detection of ERD/ERS is required. Industrial applications include faults detection and diagnosis in electro-mechanic controlled plants [START_REF] Kiltz | Fault diagnosis on a magnetically supported plate[END_REF], [START_REF] Ali | Fault diagnosis without a priori model[END_REF], in chemical processes [START_REF] Himmelblau | Fault detection and diagnosis in chemical and petrochemical processes[END_REF], detection of structural damages [START_REF] Salawu | Detection of structural damage through changes in frequency: A review[END_REF], quality control [START_REF] Himmelblau | Process analysis by statistical methods[END_REF], etc.

The main idea of the proposed algorithm is to consider a simple local bandlimited model to represent the input signal in each sliding time window of a given length. The local model consists of a sum of two exponentials which, in the frequency domain, reads as the transfer function of a second order bandpass filter. Two exponentials are enough to estimate the cut-off frequencies. Of course, we could use more than two, the price would be payed by increasing complexity of the algorithm. Relying on the operational calculus, we obtain an explicit estimation of the parameters that indicate the cutoff frequencies. To estimate them we take advantage of the algebraic framework of [START_REF] Fliess | An algebraic framework for linear identification[END_REF] (see also [START_REF] Neves | Algebraic Parameter Estimation of Damped Exponentials[END_REF], [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF], [START_REF] Ushirobira | Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data[END_REF]). Finally we slide the time window along the signal and compute the cut-off frequencies for each time instance. This is in time domain realized using a convolution integral. The solution is implemented in discrete time domain as a combination of FIR filters, providing an efficient online algorithm.

The proposed algorithm is detailed in section II, first for a lowpass model. for which the discrete time implementation is shown, and then transferred into bandpass model. In section III, results are shown on synthetic signals, ranging from low number of sine components to bandpass signal. We also provide a parameter analysis. The section IV concludes this paper.

II. ESTIMATION OF THE CUT-OFF FREQUENCIES

Let x(t) denotes the unobserved noise-free signal: x(t) = y(t) -(t). Let T > 0 be given. To any τ 0, we associate the interval I T τ = [τ, τ + T ). Thus the restriction of the signal x(t) in I T τ is defined as:

x τ (t) = H(t)x(t + τ ), t ∈ [0, T ) (1)
We propose to model the input signal with a local model that consists of a sum of two exponentials which, in the frequency domain, reads as the transfer function of a second order bandpass filter. We consider a second order lowpass filter form, which will be transformed later into a bandpass form. Signal x τ (t) reads in the operational domain as:

x τ (s) = b 0 (τ ) + b 1 (τ )s 1 -a 1 (τ )s -a 2 (τ )s 2 (2)

A. Finding coefficients of the lowpass model

To find the bandpass cut-off frequencies, only a 1 (τ ) and a 2 (τ ) must be characterized, as we will see later b 0 (τ ) and b 1 (τ ) are irrelevant for this problem. Thus we annihilate them by taking the second order operational derivative d 2 /ds 2 of the equation [START_REF] Mboup | A Volterra filter for neuronal spike detection[END_REF]. The equation then becomes: (2) . (3) Since there are now two unknowns (a 1 (τ ) and a 2 (τ )) and only one equation, we take another d/ds to obtain a system of two equations. The second equation is:

x τ (s) (2) = 2a 1 (τ )x τ (s) (1) + a 1 (τ )sx τ (s) (2) + 2a 2 (τ )x τ (s) + 4a 2 (τ )sx τ (s) (1) + a 2 (τ )s 2 x τ (s)
x τ (s) (3) = 3a 1 (τ )x τ (s) (2) +a 1 (τ )sx τ (s) (3) +6a 2 (τ )x τ (s) (1) + 6a 2 (τ )sx τ (s) (2) + a 2 (τ )s 2 x τ (s) (3) . (4)
By the very classical rules of operational calculus, multiplication by s corresponds to time derivation:

sx τ ¡ d dt x τ (t).
As numerical differentiation is difficult and ill-conditioned, equations ( 3) and ( 4) are divided by a s ν , for an integer ν > 2.

Only negative powers of s, and hence integral operators, will intervene in the resulting equation.

The unobserved signal x τ can now be replaced by its noisy observation counterpart y τ . Instead of a 1 (τ ) and a 2 (τ ), only their estimates ã1 (τ ) and ã2 (τ ) are now available.

The equations ( 3) and ( 4) can be reformulated as

P Θ = Q (5)
where Θ = ã1 (τ ) ã2 (τ ) and:

P = 2xτ (s) (1) s ν + xτ (s) (2) s ν-1 2xτ (s) s ν + 4xτ (s) (1) s ν-1 + xτ (s) (2)
s ν-2 3xτ (s) (2) s ν

+ xτ (s) (3) s ν-1 6xτ (s) (1) s ν + 6xτ (s) (2) s ν-1 + xτ (s) (3) s ν-2 , Q = xτ (s) (2)
s ν xτ (s) (3) s ν .

B. Expressing the system in the time domain

The numerical estimates are now obtained by expressing the system (5) in the time domain. For this, let us first recall that differentiation with respect to s corresponds to multiplication by -t: d ds xτ ¡ -tx τ (t), and division by s m , m > 0 maps to the m th -iterated time integration:

s -m x ¡ 1 (m-1)! t 0 (t -α) m-1 x(α)dα.
With these rules in mind and using integration by parts, the time domain analog of the system (5) is expressed as:

P (t, τ ) Θ = Q(t, τ ) (6)
where:

Q(t, τ ) = d 1 d 2 d 1 = - t 0 (t-τ ) ν-1 τ 2 (ν-1)! x τ (λ)dλ d 2 = t 0 (t-τ ) ν-1 τ 3 (ν-1)! x τ (λ)dλ P (t, τ ) = d 3 d 4 d 5 d 6 d 3 = t 0 2(t-λ) ν-1 λ (ν-1)! -(t-λ) ν-2 λ 2 (ν-2)!
x τ (λ)dλ

d 4 = t 0 -2(t-λ) ν-1 (ν-1)! + 4(t-λ) ν-2 λ (ν-2)! -(t-λ) ν-3 λ 2 (ν-3)!
x τ (λ)dλ

d 5 = t 0 -3(t-λ) ν-1 λ 2 (ν-1)! + (t-λ) ν-2 λ 3 (ν-2)!
x τ (λ)dλ

d 6 = t 0 6(t-λ) ν-1 λ (ν-1)! -6(t-λ) ν-2 λ 2 (ν-2)! + (t-λ) ν-3 λ 3 (ν-3)!
x τ (λ)dλ .

Redefining the intervals I T τ by their causal counterparts, I T τ = [τ -T, τ ], the system from equation ( 6) is considered for t = T and with the change of variable λ → T -λ. By defining

h 1 (λ) = -T 2 (1 -λ) 2 λ ν-1 h 2 (λ) = T 3 (1 -λ) 3 λ ν-1 h 3 (λ) = 2T (1 -λ)λ ν-1 -(1 -λ) 2 λ ν-2 (ν -1) h 4 (λ) = -2λ ν-1 + 4(1 -λ)λ ν-2 (ν -1) -(1 -λ) 2 λ ν-3 (ν -1)(ν -2) h 5 (λ) = T 2 -3(1 -λ) 2 λ ν-1 + (1 -λ) 3 λ ν-2 (ν -1) h 6 (λ) = 6T (1 -λ)λ ν-1 -6T (1 -λ) 2 λ ν-2 (ν -1) +T (1 -λ) 3 λ ν-3 (ν -1)(ν -2) .
the elements of the coefficients P and Q could be rewritten as convolution integrals:

d i = ∞ 0 h i (λ)x(t -λ)dλ where i = 1, . . . , 6.

C. Discrete-time implementation

We turn now to the discrete-time setting where only time samples of the observation, y m = y(t m ), are available. Assuming a uniform sampling, let t m = mT s , m = 0, 1, . . ., with a sampling period T s . Let T = M T s and ν > 2 be given. Set τ = nT s and compute the finite sequence for i = 1, 2, ..., 6 {h i,m = h i (t m )} M m=0 . Then, by a numerical integration method with abscissas t m and weights W m , m = 0, • • • , M , we obtain the approximation of d 1 , d 2 , d 3 , d 4 , d 5 and d 6 .

d i (τ ) ≈ d i,n = M m=0 g i,m y n-m , where g i,m = W m h i,m .
In all the sequel, we use the trapezoidal method for which W 0 = W M = 0.5, and

W m = 1, m = 2, • • • , M -1. As τ = τ n = nT s varies (sliding windows I M n ), d 1 , d 2 , d 3 , d 4 , d 5
and d 6 read as the outputs of the FIR filters, with impulse response sequences {g i,m } M m=0 for i = 1, 2, ..., 6. Accordingly, solving the system from equation ( 6) we obtain discrete approximations of a 1,n ≈ a 1 (nT s ) and a 2,n ≈ a 2 (nT s ).

D. Transferring lowpass into bandpass model

The general bandpass form is obtained by transferring the low-pass general form:

h(s) = b(τ ) a(τ ) + s (7)
into band-pass:

h(h1(s)) = b(τ ) a(τ ) + s 2 +ω1(τ )ω2(τ ) (ω2(τ )-ω1(τ ))s . (8)
Equation ( 8) can be written in similar form as equation (2):

h(h1(s)) = (ω2(τ )-ω1(τ ))b(τ )s ω1(τ )ω2(τ ) 1 + (ω2(τ )-ω1(τ )) ω1(τ )ω2(τ ) a(τ )s + 1 ω1(τ )ω2(τ ) s 2 . (9)
Now we chose a(τ ) = 1 and from here ω1,n and ω2,n , the discrete approximations of ω1 (τ ) and ω2 (τ ), can be expressed as:

ã2,n ω2 2,n -ã1,n ω2,n + 1 = 0, (10)
ω1,n = -1 ã2,n ω2,n . (11)
Thus for each n (as I M n slides), the cut-off frequencies of the model ω1,n and ω2,n can be estimated. Since ω2,n is a solution of a square equation and ω1,n is computed using ω2,n , we will have two solutions for each estimate of the cut-off frequencies.

III. RESULTS

This section shows the results obtained for different signal complexity, ranging from few sine functions to bandpass signal, and proposes a parameter analysis of the algorithm. The evaluations are always performed on signals composed of three pieces, each one having different spectral properties. The length of the signal was 90s with sampling frequency of 100Hz giving 9000 samples. We assess the ability of the algorithm to accurately estimate the frequency band occupied by each piece of the signal.

In the simulations to follow, the estimated values of both solutions for ω1,n and ω2,n are first processed by a median filter (the length of the median filter was set to 35 samples) and then smoothed by a moving average filter (the length of the moving average filter was set to 45 samples). Such processing was necessary to reduce the perturbations of the estimations caused mostly by the numerical instability of the matrix inverse P -1 . The lengths of the two filters were chosen experimentally, by observing the stability of the results. Parameter ν was set to 3 and parameter M was set to 10. The influence of these values will be discussed later.

The first test signal (top row of Fig. 1) is composed only from two sine functions in each of the three pieces. Each piece displays different frequencies for the sine functions, as shown in the middle row of Fig. 1 (here and on the remaining figures, each frequency plot refers to the corresponding part of the signal). Such signal is well suited for the second order model used in this algorithm. Both solutions for ω1,n and ω2,n are given on the plot in the bottom row of Fig. 1. Even though they are not fully stable, they indicate well the spectral structure of the corresponding parts of the signal.

Figure 2 shows the estimation obtained for an input signal consisting in a sum of 5000 sinusoids, each with random phase and amplitude. The signal is again composed of three pieces: for the first one the frequency components were randomly drawn from the interval 15 -18Hz, for the second from the interval 12 -15Hz and for the third from the interval 9 -12Hz. The oscillations of the cut-off frequency estimations are now larger due to the very high degree of undermodelling. Nonetheless, the general trend of cut-off frequency values still indicates the occupied band, and the change-points between the bands are clearly recovered as well.

Instead of having only a finite number of frequency components, each of the three pieces of the signal is now constituted as a full bandpass signal. The results, shown on Fig. 3, are similar to those obtained with a large but finite number of sinusoids. This suggests the robustness of the algorithm to the accuracy of the model, given that the model is only a second order.

In the previous experiments, the signal did not contain any frequency component except of those in the bands mentioned, resulting in a noise-free analysis. We now add different amounts of noise by generating white noise and limiting its spectrum to the 1 -25Hz band. A case for SNR = 6dB is shown on figure 4. Figure 5 shows the effect of different levels of the noise on the mean values and standard deviations of the estimated cut-off frequencies for the signal in figure 3. It is clear that the noise perturbations significantly affect the accuracy and stability of the results. However, for relatively high SNR the results could still be acceptable for different applications.

To get an idea on how the performance of the proposed parametric method compares to standard non-parametric approaches, we compared our algorithm with STFT. We applied STFT on the noisy input signal shown on the first plot of figure 4. We set the window length for STFT to 90 samples, the same as the total length of the filters used with the proposed method. The last plot on figure 4 displays the result of STFT. It is clear that, even though the price for good time resolution is payed by lower spectral resolution, the non-parametric approach gives useful results as well. Change points appear clearly (no false negatives), which is however not the case for the cut-of frequencies. In on-line applications, large errors in estimation of the cut-of frequencies could cause false alarms in detection of change points (false positives). The results suggest that this is as well the case for the proposed approach.

We explore now the influence of the free parameters M and ν on the estimation accuracy. Increasing M means increasing the impulse length of the filters h i,m . This results in higher sensitivity of the algorithm so that even the frequency components of the lower power can be detected. In general the estimated values of the frequency bands are wider than with lower M . This is demonstrated on figure 6, which is created using the same signal as figure 5, but with M set to 30 instead of 10. The increase of ν has a similar effect: it results in wider effective impulse length of the filters (the number of the samples significantly different from 0). The results for ν = 6 are shown on figure 7.

IV. CONCLUSION

Even though we modeled the signal by a simple second order model, the estimated cut-off frequencies were meaningful even for the bandpass signal. Moreover, the estimations were meaningful even in the presence of a moderate amount of noise. We also showed that the sensitivity of the algorithm can be experimentally tuned by adjusting the free parameters M and ν.

ACKNOWLEDGMENT

This work was supported by the foundation EADS and PHC COGITO project.

  of three consecutive segments each with 30s duration. The signal in each segment is a sum of 2 sine components. for each of the two estimations of the cut-off frequencies ω 1 and ω 2 after filtering with median and moving average filters.

Fig. 1 . 1 and ω 2 after 2 Fig. 2 .

 11222 Fig. 1. Estimation results for piecewise signal with each segment composed of 2 sine functions.

1 and ω 2 after

 12 of three consecutive segments each with 30s duration. Each segment consists of a bandpass signal. for each of the two estimations of the cut-off frequencies ω filtering with median and moving average filters.

Fig. 3 . 1 and ω 2 after

 312 Fig. 3. Estimation results for piecewise bandpass signal.

Fig. 4 .

 4 Fig. 4. Estimation results for noisy piecewise bandpass signal. SNR = 6dB.

  third segment of the signal

Fig. 5 .

 5 Fig. 5. Mean values of the estimated cut-off frequencies for each of the three segments for different SNR (in dB). Vertical bars denote the standard deviation of the estimated cut-off frequencies.

Fig. 6 .

 6 Fig.6. Analysis of influence of parameter M . Simulation setup is similar to Fig.5, but the value of M is set to 30 instead of 10.

  third segment of the signal

Fig. 7 .

 7 Fig.7. Analysis of influence of parameter ν. Simulation setup is similar to Fig.5, but the value of ν is set to 6 instead of 3.