.1. In this experiment, the data under H 0 was generated by sampling from the mixture of both distributions. 100 trials were done. (NB: the conditions associated with TRank ref are specified in M. Depecker's PhD Thesis [Dep10].) . . . . . . .

ii First of all, I would like to thank my advisor, Dr. Frédéric Cazals, for trusting me and believing in my ideas, for giving me the necessary freedom to work on them and for supporting me with his experience, enthousiasm and positivity all along the process and, especially, in my darkest moments of doubt. I would like to express my gratitude to Prof. Dr. Peter Grünwald for fruitful discussions and important advice, in particular, with respect to the important optional stopping property of the tests proposed in this thesis, and for allowing me to expand my knowledge in the fascinating field of Statistical Learning. I would also like to thank again my Master's advisors Prof. Dr. Alfredo Viola and Dr. Gadiel Seroussi, for helping me in building my background in Algorithms, Coding and Information Theory, which was essential to this thesis. Additional thanks go to Prof. Dr. P. Grünwald and Dr. G. Seroussi, in their role of "Rapporteur," for their valuable comments and corrections, which greatly improved the quality of this manuscript.

Overall, I am grateful to Inria for providing me excellent conditions to work and, in particular, to cluster administrators for their efficient support when running my multiple experiments. I would also like to thank my labmates Andrea, Angeliki, Deepesh, Dorian, Noël, Romain, Simon and Tom for the discussions and shared fun moments. Special thanks to Tom for his software support.

I would like to thank my parents for supporting and encouraging me during all this time. Special thanks to my mother and Andrzej who made possible and easy our move to the wonderful city of Antibes. In this city, we met many new friends that I would like to thank for their support too. Last but not least, I will never finish to thank my wife Ana and my children Amélie and Yann, for their support and for enduring this long process that sometimes resembled an emotional rollercoaster.

iii

List of Figures

2.1 Random multiplexer generating pairs (label, position). . . . . . . . . . . . . . . .

2.2

Localizing the discrepancy: algorithm illustrated on a toy 1D example. (A, discrepancy estimation) We consider the height function defined by the discrepancy -δ(z), called the landscape for short. Practically, the estimate δ(z i ) is used at each sample z i (see text for details). (B, clusters as connected components) Focusing on discrepancies larger than a user defined threshold δ max , clusters are defined as connected components (c.c.) of the sublevel set of the landscape defined by δ max . On this example, four such c.c. are observed, two large ones and two small ones. Note that each cluster is charged to one local minimum of the function -δ(z). (C, smoothed clusters upon applying topological persistence) To get rid of small clusters, the landscape is smoothed using topological persistence (see text for details). In our case, this simplification yields two clusters. Note that the blue sample squeezed between C 1 and C 2 does not appear in the red cluster since its discrepancy is below δ max (and likewise for the red point with respect to C 2.14 Landscape simplification using topological persistence: simple strategy using Union-Find, versus refined strategy using the Morse-Smale-Witten (MSW) complex. (A,B) Two landscapes, with critical points of indices 0 (disks) and 1 (squares).

(C) The disconnectivity graph (DG) of both landscapes, namely the tree depicting the evolution of connected components of sublevel sets. Despite the differences between their MSW complexes, both landscapes share the same DG: upon passing the critical point e, the stable manifold of c merges with that born at a. (D,E) The MSW complexes of (A,B), respectively. In cancelling the pair of critical points (c, e), one does not know from the DG with which (a or b) the basin of c should be merged. But the required information is found in the MSW complex: on the landscape (A), c is merged with a; on the landscape (B), c is merged with b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.1

The random multiplexer used to generate pairs (label, position) . . . . . . . . .

Tree and partitions (A)

Partition of a 2D domain by a binary space partition tree, with two levels. (B) The corresponding full tree. To each node is associated a region; internal nodes also come with a splitting rule to define the regions of the two children.

Note that a tree T generates a set of partitions. Identifying a node with its associated region, these partitions are P T = {η 1 , {η 2 , n 3 }, {η 2 , η 4 , η 5 }}. In this thesis, we study problems related to learning and detecting multivariate statistical dissimilarity, which are of paramount importance for many statistical learning methods that are nowadays applied in an increasingly number of fields like biology, medical diagnosis, cheminformatics, fraud detection, computer vision, handwriting recognition, social networks, computational finance, recommendation systems, etc. 1 Statistical dissimilarity is, in essence, related to the problem of being able to distinguish one random variable from another (i.e., to distinguish their distribution). In particular, we consider the following problems:

• the two-sample problem: given two sets of observations, decide whether they come from the same distribution or not. For example, one may want to find out whether two classes of patients show different response to a drug, where the response can be described by a number of variables. We will also consider the case where these observations can be processed in a sequential manner which allows more flexibility in the sampling process and is of interest in data streaming applications.

• nonparametric multivariate effect size determination: in both previous cases, we consider unknown multivariate distributions without any assumption on their nature, i.e., the problems are nonparametric. Therefore, when a change is detected we are interested in determining the nature of the change. For example, if patients of the different classes show different response to the drug, one may want to identify where this difference is located in the space of variables considered to further infer if it is an intrinsic and important difference or eventually some random error due to the sampling process.

In Section 1.2, we introduce measures of statistical dissimilarity and the related problem of binary classification. In Section 1.3, we show how to learn these objects using data. In Section 1.4, we recall fundamental concepts of hypothesis testing, we present methods for nonparametric two-sample testing.

Statistical Dissimilarity

In this section, we define some important dissimilarity measures and the related problem of binary classification.

Dissimilarity Measures

Let us first present different ways of quantifying statistical dissimilarity. Let X and Y be two random variables defined on a continuous space Ω with densities f X and f Y and cumulative distributions F X and F Y respectively.

1 See, e.g., https://en.wikipedia.org/wiki/Machine_learning#Applications 1 Many useful measures of dissimilarity between distributions d(F X , F Y ) are not metrics but satisfy some of their properties, which are:

1. non-negativity : d(F X , F Y ) ≥ 0 2. identity of indiscernibles: d(F X , F Y ) = 0 ⇔ F X = F Y ⇔ f X = f Y a.e. 3. symmetry: d(F X , F Y ) = d(F X , F Y ) 4. triangle inequality : d(F X , F Z ) ≤ d(F Z , F Y ) + d(F Y , F Z ) .
Several dissimilarity measures have been introduced in the literature and some can be derived as special cases of others. We focus on three of them that are of special relevance for our work. We present their versions for continuous Ω since it is our focus, but notice that they can be defined for more general probability measures. f-divergences A common class of dissimilarity measures is the one of f-divergences. These divergences were introduced and studied in [START_REF] Csiszár | Eine informationstheoretische ungleichung und ihre anwendung auf den beweis der ergodizitat von markoffschen ketten[END_REF][START_REF] Morimoto | Markov processes and the h-theorem[END_REF][START_REF] Syed | A general class of coefficients of divergence of one distribution from another[END_REF] and are also know as Csiszár f-divergences, Csiszár-Morimoto divergences or Ali-Silvey distances. Intuitively, they can be seen as an average, weighted by a function f, of the density ratio f X (•) f Y (•) . For a strictly convex function f satisfying f(1) = 0, the f-divergence between f X and f Y is defined as

D f (f X f Y ) = Ω f f X (z) f Y (z) f Y (z)dz (1.1)
or, more generally, for probability measures P X and P Y such that P X is absolutely continuous with respect to P Y2 as

D f (P X P Y ) = Ω f dP X dP Y dP Y (1.2)
where dP X dP Y is the Radon-Nikodym derivative.3 An f-divergence is non-negative and satisfies the identity of indiscernibles, but they do not satisfy the two other properties of a metric.

The popular Kullback-Leibler (KL) divergence is an instance of this class, obtained with f(z) ≡ log(z) and denoted as D KL (f X f Y ). Other important dissimilarity measures are obtained with different functions f (see, e.g., [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF]).

As we will see in Section 1.4.2, several hypothesis tests are built upon f-divergences.

Jensen-Shannon Divergence

The KL divergence has two properties that may not be desired, namely, it is unbounded and it is not symmetric. Denoting the continuous entropy by H (f ) ≡ Ω f (z) log(f (z))dz and considering the average density f ≡ θ 0 f X +(1-θ 0 )f Y with θ 0 ∈ [0, 1], the Jensen-Shannon Divergence (JSD) [START_REF] Lin | Divergence measures based on the Shannon entropy[END_REF] is defined as follows

JS θ0 (f X f Y ) ≡ H f -θ 0 H (f X ) -(1 -θ 0 )H (f Y ) (1.3)
which can be equivalently expressed as

JS θ0 (f X f Y ) = Ω f (z) log f (z)dz -θ 0 Ω f X (z) log f X (z)dz -(1 -θ 0 ) Ω f Y (z) log f Y (z)dz (1.4) = θ 0 Ω f X (z) log f (z) f X (z) dz + (1 -θ 0 ) Ω f Y (z) log f (z) f Y (z) dz (1.5) ≡ θ 0 D KL f X f + (1 -θ 0 )D KL f Y f (1.6) since, by linearity, Ω f (z) log f (z)dz = θ 0 Ω f X (z) log f (z)dz + (1 -θ 0 ) Ω f Y (z) log f (z)dz.
Therefore, the JSD avoids infinite values by averaging the KL divergence of each density with respect to the average density. The JSD is non-negative and satisfies the identity of indiscernibles. When θ 0 = 1/2, the JSD is symmetric, bounded between 0 and 1 and its square root yields a metric [START_REF] Endres | A new metric for probability distributions[END_REF].

Interestingly, as will see in Section 1.2.2 and Chapter 2, the JSD is linked to binary classification.

Maximum Mean Discrepancy

Let F be a class of functions g : Ω → R, the maximum mean discrepancy

(MMD) is defined in [FM53, GBR + 12] as MMD[F, f X , f Y ] ≡ sup g∈F (E X [g(X)] -E Y [g(Y )]) . (1.7)
In words, the dissimilarity is characterized by a function g from a class F, whose expected value is large on f X and small on f Y . The class F must be "rich enough" to uniquely identify f X = f Y a.e.. Let us consider the space of continuous and bounded functions on Ω, i.e., F ≡ C 0 (Ω), then [GBR + 12, Lemma 1] implies the identity of indiscernibles. Nevertheless, this class of functions is too rich to be practical and the "richness" of the class affects the convergence rate of estimates [GBR + 12].

Reproducing kernels are a bridge from linearity to non-linearity for algorithms that can be expressed in terms of dot products, by mapping the original space to a higher (even infinite) dimensional one. Fortunately, this mapping does not need to be explicitly computed as only a reproducing kernel function is required to replace the dot product used by the algorithm (the so-called kernel trick ).

One way to obtain a class that is rich enough to identify any kind of dissimilarity but simple enough to have good convergence rates is to consider a reproducing kernel Hilbert space (RKHS) H , which is a Hilbert space (i.e., an abstract vector space with inner product •, • H ) of functions g : Ω → R equipped with the linear evaluation functional δ x : g → g(x) that is bounded on H. Then, the Riesz representation theorem ensures that ∀x ∈ Ω there exists a unique φ x ∈ H, called the feature mapping such that g(x) = g, φ x H . The reproducing kernel is the (positive definite) function

k : Ω × Ω → R defined as k(x, y) = φ x , φ y H .
Then, the unit ball (i.e., functions g such that g = g, g H ≤ 1) in a RKHS is a more restrictive class, but yet rich enough to uniquely identify whether f X = f Y [GBR + 12, Theorem 5] when a universal kernel is used. Furthermore, the MMD satisfies all the metric properties except non-negativity when H is a universal RKHS defined on a compact metric space, which is the case for the Gaussian

k(x, y) ≡ exp - x -y 2 2σ 2 (1.8)
and Laplace k(x, y) ≡ exp -xy σ (1.9)

cases. [GBR + 12, Lemma 6] provides an expression for squared MMD in terms of the reproducing kernel functions in a RKHS:

MMD 2 [F, f X , f Y ] = E X,X ′ [k(X, X ′ )] -2E X,Y [k(X, Y )] + E Y,Y ′ [k(Y, Y ′ )] (1.10)
where X ′ is an independent copy of X and Y ′ is an independent copy of Y . Another interesting special case of this measure is obtained for Ω = R, when considering in Eq. (1.7) the class of functions with bounded variation, i.e., functions g such that |∂ x g(x)|dx ≤ C for some constant C. Setting C = 1 yields the Kolmogorov-Smirnov metric [GBR + 12, Proposition 21] that yields the famous Kolmogorov-Smirnov two-sample test [START_REF] Kolmogorov | Sulla determinazione empirica delle leggi di probabilita[END_REF][START_REF] Smirnov | Table for estimating the goodness of fit of empirical distributions[END_REF].

As shown in [GBR + 12], the Monge-Wasserstein, or Mallows, distance can also be obtained as a special case of MMD using some specific function class in Eq. (1.7). Informally, if we see the two densities as different ways of piling up a certain amount of dirt over Ω, the Monge-Wasserstein distance is the minimum cost of turning one pile into the other, where the cost is the amount of dirt moved times the distance by which it is moved.

Binary Classification

As we will see throughout this work, statistical dissimilarity is intimately related to binary classification. For this task, it is assumed that samples are drawn from a mixture of the two distributions and they are labeled according to the distribution they were drawn from: the label 0 for f X and the label 1 for f Y .We denote by L the random variable corresponding to the label and θ 0 ≡ P (L = 0). The label alphabet is denoted as A ≡ {0, 1}. Therefore, a pair of random variables (Z, L), corresponding respectively to position (whose components are usually called features) and label, with joint density f Z,L , is considered. Since Z comes from a mixture of X and Y , it takes values on Ω as well. The density of each population can be expressed as

f X (z) = f (z|L = 0) and f Y (z) = f (z|L = 1). The density of Z is f Z ≡ θ 0 f X + (1 -θ 0 )f Y ≡ f .
Then, having seen a set of samples with their labels and the value of Z of a new sample, the goal is to predict its label L either in hard way by telling one value or in a soft way by estimating a probability distribution on its value while trying to minimize some loss function. Intuitively, the higher the dissimilarity between f X and f Y is, the easier the classification becomes.

Hard Classifiers and the 0/1 Loss A hard classifier is a function g : Ω → A. An error occurs if g(Z) = L. The 0/1 loss is defined as L 0/1 (g(z), l) ≡ 1 l =g(z) , where 1 denotes the indicator function. The Bayes error is given by the function g ⋆ that minimizes the probability of error, i.e.

L ⋆ 0/1 ≡ min g:Ω→A P (g(Z) = L) .
(1.11)

The mutual information I (Z; L) quantifies how much the knowledge of Z reduces the uncertainty about L. The JSD between f X and f Y turns out to be equal to I (Z; L) since, denoting by H (Z) the entropy of Z and by H (Z|L) the conditional entropy of Z given L, we have that

I (Z; L) = H (Z) -H (Z|L) (1.12) = H (f Z ) -θ 0 H (f X ) -(1 -θ 0 )H (f Y ) (1.13) ≡ JS θ0 (f X f Y ) .
(1.14) Furthermore, the JSD provides the following bounds on the Bayes error rate [START_REF] Lin | Divergence measures based on the Shannon entropy[END_REF]:

1 4 (h(θ 0 ) -JS θ0 (f X f Y )) 2 ≤ L ⋆ 0/1 ≤ 1 2 (h(θ 0 ) -JS θ0 (f X f Y )) (1.15)
where h(•) denotes the binary entropy (i.e., h(θ)

≡ -θ log 2 θ -(1 -θ) log 2 (1 -θ)).
As noted in [GBR + 12, Remark20][Fri04] and further discussed in [START_REF] Mark | Information, divergence and risk for binary experiments[END_REF], the optimal classifier g ⋆ (or some sequence of classifiers that approaches it), can be turned into a dissimilarity measure by considering

E X [g ⋆ (X)] -E Y [g ⋆ (Y )].
Soft Classifiers and the Logarithmic Loss A soft classifier (also known as probabilistic classifier ) is a function g : Ω → D A where D A is the set of all distributions over A (i.e., D A ≡ {p : l∈A p(l) = 1}). Of particular interest for our work is the logarithmic (log) loss, which is defined as L log (g(z), l) ≡log g(z)(l).

A vast literature in Information Theory exists for this loss (see, e.g., [START_REF] Cover | Elements of Information Theory[END_REF]). The log loss possesses the following important properties (see, e.g., [START_REF] Merhav | Universal prediction. Information Theory[END_REF] for more details):

• it is monotonically decreasing with the assigned probability to the observed label,

• technically, it is convenient to work with since it converts joint probabilities, or equivalently, products of conditional probabilities into cumulative sums, and

• it can be used for the purpose of testing the validity of statistical models (i.e., families of parametric distributions) by Dawid's prequential principle [START_REF] Philip | Present position and potential developments: Some personal views: Statistical theory: The prequential approach[END_REF] encompassed in Rissanen's Minimum Description Length principle [START_REF] Rissanen | Modeling by shortest data description[END_REF][START_REF] Rissanen | Fisher information and stochastic complexity[END_REF][START_REF] Rissanen | Stochastic complexity[END_REF][START_REF] Grünwald | The minimum description length principle[END_REF]. One of the reasons for this is that the loss is minimized, in the expected or in the almost sure sense, by the distribution that generates the samples, a property shared by very specific loss functions (see [START_REF] Merhav | Universal prediction. Information Theory[END_REF][START_REF] Miller | On loss functions which minimize to conditional expected values and posterior probabilities[END_REF]).

This last property is particularly important since it means that when learning a soft classifier one usually seeks to approximate the true conditional distribution g ⋆ (z)(•) ≡ P (L = •|Z = z). This makes the problem addressable by nonparametric regression with random design (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]), where given a random vector (Z, R), where Z ∈ R d and the response variable R ∈ R, the goal is to estimate the regression function defined as

m(z) ≡ E [R|Z = z] . (1.16) Therefore, setting R ≡ L yields m(z) = P (L = 1|Z = z) ≡ g ⋆ (z)(1).
(1.17)

Learning Dissimilarity Measures and Classifiers

In the previous section, we have introduced different measures of dissimilarity and the related problem of classification. In this section, we will address the problem of learning the involved objects, in particular, dissimilarity measures and classifiers (with emphasis on soft ones).

Settings

We first present some of the common settings for learning dissimilarity measures and classifiers.

Dissimilarity Measures

When learning dissimilarity measures, two i.i.d. samples are usually given, that is x 1 , . . . , x n0 and y 1 , . . . , y n1 whose corresponding random variables X i ∈ Ω and Y i ∈ Ω are i.i.d. with densities f X and f Y respectively.

Classification

For the classification problem, different settings can be considered. When the construction uses label information, it is called supervised learning. In this problem, samples (Z i , L i ) are usually assumed to be i.i.d. with joint density f Z,L .

Batch The batch (or off-line) learning setting is the classical case of supervised learning (see, e.g., [START_REF] Devroye | A probabilistic theory of pattern recognition[END_REF]). A training set of n labeled samples is given, that is (z, l) n ≡ (z 1 , l 1 ), . . . , (z n , l n ) from which a classifier ĝn : Ω → A is built. Then, this classifier is assessed on a different test set (z n+1 , l n+1 ), . . . , (z n+m , l n+m ). Note that we will also refer to the pair sequence (z, l) n as l n , z n .

On-line In the on-line setting, symbols are observed sequentially one by one. Based on previously observed pairs (z 1 , l 1 ), . . . , (z i-1 , l i-1 ), an on-line learning algorithm tries to build a classifier ĝi-1 : Ω → A predict the next label L i to come given its position z i . Therefore, an on-line learning algorithm can improve its performance as new samples are seen. This is also known as prediction with side information (see, e.g., [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF]). Notice that this is also a case of supervised learning.

Semi-supervised In the semi-supervised setting, the classifier has access to an additional set of unlabeled samples (z ′ 1 , . . . , z ′ u ). This can be of interest since it is common in practice that only few labeled samples are available for training while a larger quantity of unlabeled samples are available.

Roughly speaking, the semi-supervised learning is useful, when the unlabeled data helps elucidating the structure of the classifier. In other words, the knowledge on f Z acquired through the unlabeled data has to carry information that is useful in the inference of P (L = •|Z = z). This can be the case under certain assumptions (see, e.g., [CSZ + 06] for an in-depth treatment of the subject). For example, if we know that the points of each class tend to form clusters, then the unlabeled data can help in finding the boundary of each cluster more accurately. A simple method would be to identify clusters with some unsupervised method (without taking into account any label) and then assigning a label to each cluster according to the given labels.

Desired properties

Consistency When estimating dissimilarity measures, one usually seeks some sort of consistency, that is, that estimated measure of dissimilarity converges (in some sense) to the true dissimilarity measure as the number of samples used for the estimation tends to infinity.

For a soft classifier ĝn (z) to be consistent, it is usually required that it converges is some sense to P (L = •|Z = z). For example, using a nonparametric regressor m n (z) (i.e., an estimate of m(z) built using a set of n samples) the following strong sense known as strong pointwise consistency can be considered (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF])

m n (z) n→∞ ----→ m(z) a.s. (1.18)
for f Z -almost all z. Using the L 2 norm to quantify the difference or error between the estimate and the true regression function, two other notions of consistency are commonly used, which are the strong consistency

E Z |m n (Z) -m(Z)| 2 n→∞ ----→ 0 a.s. (1.19)
and the weak consistency

E L n ,Z n E Z |m n (Z) -m(Z)| 2 n→∞ ----→ 0.
(1.20)

These definitions translate straightforwardly to conditional probability estimation using the mapping of Eq. (1.17).

Bias-Variance Tradeoff

The bias-variance tradeoff is the dilemma that is faced when trying to simultaneously minimize two sources of error that prevent supervised learning algorithms from generalizing beyond their training set:

• the bias is the error generated by erroneous assumptions in the learning algorithm. High bias can cause underfitting, i.e., missing relevant relations between features Z and labels L.

• the variance is the error generated by the sensitivity of the learning algorithm to small fluctuations. High variance can cause overfitting, i.e., capturing noise or relations that do not generalize.

For example, in nonparametric regression, let m n be an arbitrary estimate, then for any Z ∈ R d the expected L 2 error of m n at Z can be written as [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]

E L n ,Z n E Z |m n (Z) -m(Z)| 2 = E Z E L n ,Z n |m n (Z) -m(Z)| 2 (1.21) = E Z [Var(m n (Z))] + E Z |Bias(m n (Z))| 2 (1.22)
where Var(m n (z)) is the variance of the random variable m n (z) and Bias(m n (z)) is the difference between the expectation of m n (z) and m(z), in both cases with respect to the random variables L n , Z n . As shown in [GK02, Sec. 2.3], the importance of these decompositions is that the expected variance and the expected squared bias depend in opposite ways on the "wiggliness" ("non-smoothness") of an estimate, yielding a bias-variance tradeoff. In other words, wiggly estimates allow to better fit the data and thus to reduce expected bias but, on the other hand, it is compensated by an increase in expected variance.

Universality Let us define, in the online case, the redundancy of a soft classifier ĝn (z) as the accumulated excess loss with respect to the true conditional distribution, i.e.,

R n (ĝ, g ⋆ ) ≡ n i=1 L log (ĝ i-1 (z i ), l i ) -L log (g ⋆ (z i ), l i ).
(1.23)

Then, a desired property is that the normalized redundancy 1 n R n (ĝ, g ⋆ ) vanishes to zero in some sense (e.g., in expected value or almost surely), as n tends to infinity [Grü07, Def. 2-3, p. 200]. Naturally, the convergence speed, which reflects the learning rate, is important too.

Note that in the literature, the term "universal" is also used in the non-stochastic case, where prediction is also done using distributions but the goal is to approach the performance of the best distribution of some given reference class (see, e.g., [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF]).

Effectiveness in High-Dimensions The traditional bane of nonparametric statistics is known as the curse of dimensionality [START_REF] Ernest | Adaptive control processes: a guided tour[END_REF], meaning that the prediction performance deteriorates dramatically as the dimension d of the data increases. In nonparametric regression, where m(x) is p times continuously differentiable, the best convergence rate of the expected error in norm L 2 for the worst case of m(x) (i.e., the minimax rate) is n -2p (2p+d) (see [START_REF] Charles | Optimal rates of convergence for nonparametric estimators[END_REF][START_REF] Charles | Optimal global rates of convergence for nonparametric regression[END_REF]). In other words, the number of samples required to attain a low expected error may be exponential in d, which can be difficult to satisfy even for moderate values of d. However, it is often the case that data which appears high dimensional has, in fact, some lower intrinsic dimensionality (in some sense), due to the fact that, roughly speaking, features may be correlated or redundant. This is the case in many real life data, e.g., in human body motion capture where 100 markers can be attached to the body and each marker measures position in three dimensions: dimension is 300 but motion is in fact constrained by joints and angles which yields data of much lower intrinsic dimension (see, e.g., [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF][START_REF] Verma | Which spatial partition trees are adaptive to intrinsic dimension?[END_REF]).

Therefore, methods whose error convergence rate depend only on the intrinsic dimension are of special interest.

Algorithmic Complexity Online methods ideally need to have constant processing time and memory in order to be effective in practice for potentially unbounded sample sizes. Nevertheless, this can be difficult to achieve when consistency and universality are sought.

Methods of Estimation

Now we present three classic and widely used classes of methods for building estimates for dissimilarity measures and for conditional probability, namely partitioning, nearest-neighbors and kernel based methods.

For conditional probability estimation via nonparametric regression, local averaging estimators of the following form can be used to estimate the regression function m(z) (Eq. (1.16)) (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]):

m n (z) ≡ n i=1 W n,i (z)R i (1.24)
where the weights W n,i (z) ≡ W n,i (z, Z 1 , . . . , Z n ) ∈ R depend on Z 1 , . . . , Z n .

In the case of the estimation of the f-divergence , a traditional approach is to first estimate the individual densities using some nonparametric method (e.g., kernel density estimation) and substitute in the definition formula. Nevertheless, the problem of estimating the individual densities is harder than estimating the divergence (see, e.g., [START_REF] Pérez-Cruz | Kullback-Leibler divergence estimation of continuous distributions[END_REF]), and thus, by Vapnik's advice [START_REF] Vapnik | The nature of statistical learning theory[END_REF], should be avoided.

Partitioning

This class of methods is based on considering a sequence of (finite or countably finite) partitions P n = {A n,1 , A n,2 , . . . } of R d , where A n,1 , A n,2 , . . . are Borel sets.

Classification In order to build a locally averaging estimator, consider the following functions

W n,i (z) ≡ 1 Zi∈An,j n l=1 1 Z l ∈An,j
for z ∈ A n,j .

(1.25)

The sequence of partitions is nested if the sequence of generated σ-algebras F(P n ) is increasing (i.e., it is a filtration).

For z ∈ R d , let us define A n (z) ≡ A n,j if z ∈ A n,j , .i.e., the cell of the n-th partition containing z. Denoting by λ the Lebesgue measure on R d and defining the diameter of a cell as

diam(A n (z)) ≡ sup u,v∈An(z) u -v (1.26)
the following result proves the convergence (see [GK02, Thm. 24.3 and Thm. 25.6])

Theorem. 1.1. Consider the estimator obtained from Eq. (1.24) and 1.25. If the sequence of partitions

P n is nested, diam(A n (z)) n→∞ ----→ 0 for each z ∈ R d and ∞ i=n λ(A i (z)) = ∞ (1.27)
for each z and n, then m n is weakly and strongly consistent. Moreover, if, for all z,

nλ(A n (z)) log n → ∞ (1.28)
and |R| < C for some C < ∞ then m n is also strongly pointwise consistent.

The difficult part of using partition based methods is to define how the space should be partitioned in order to satisfy conditions that makes the procedure consistent but also to be suitable for high-dimensional data with low-intrinsic dimension.

Regression trees proposed in [START_REF] Gordon | Consistent nonparametric regression from recursive partitioning schemes[END_REF][START_REF] Gordon | Almost surely consistent nonparametric regression from recursive partitioning schemes[END_REF] use a data dependent rule to split the space recursively yielding therefore a structure that can be represented as a tree , thus enjoying logarithmic time to find a cell when the tree is close to balanced.

Many extensions and variations to these idea have been proposed including the popular Breiman's Random Forests [START_REF] Breiman | Random forests[END_REF], online versions as in [START_REF] Denil | Consistency of online random forests[END_REF] and the more recent intrinsic-dimension adaptive Random Projection Trees [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF][START_REF] Kpotufe | A tree-based regressor that adapts to intrinsic dimension[END_REF]. In contrast to this intrinsic-dimension adaptive partitioning, axis-parallel splitting rules have been shown to be prone to the curse of dimensionality even when intrinsic dimension is lower [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF].

In [START_REF] Nobel | Histogram regression estimation using data-dependent partitions[END_REF], conditions that guarantee strong consistency of regressors based on data-dependent partitions are shown.

Divergence In [START_REF] Wang | Divergence estimation of continuous distributions based on data-dependent partitions[END_REF], the KL divergence is estimated by estimating the Radon-Nikodym derivative using the ratio of the empirical measures P X,n0 (A n,j ) and P Y,n0 (A n,j ) in each cell A n,j of each partition of the sequence, i.e.,

P X,n0 (A n,j ) ≡ 1 n 0 n0 i=1 1 Xi∈An,j (1.29)
P Y,n1 (A n,j ) ≡ 1 n 1 n1 i=1
1 Yi∈An,j .

(1.30)

The partitions are built using statistically equivalent blocks (Gessaman's partition [START_REF] Gessaman | A consistent nonparametric multivariate density estimator based on statistically equivalent blocks[END_REF][START_REF] Nobel | Histogram regression estimation using data-dependent partitions[END_REF]), with respect to the reference measure P Y , using the following procedure. First, the intended number of points per cell k n for the n-th partition must be chosen. Then, C n ≡ ⌊(n/k n ) 1/d ⌋ is the number of cells for the n-th partition. The iterative part starts with one cell containing all the points and then, for every coordinate i = 1..d, the following step is repeated:

• For each cell, project the samples along the i-th axis and split it into C n cells, using hyperplanes perpendicular to the i-th axis, such that all of them contain the same number of points (with the possible exception of the last cell)

The almost sure convergence of this method was proved under mild conditions. The authors propose some extensions that improve convergence rate by using non-uniform data-adaptive partitioning schemes and using a bias correction term. This method was further analyzed (for multivariate data) and extended in [START_REF] Silva | Universal consistency of data-driven partitions for divergence estimation[END_REF][START_REF] Silva | Information divergence estimation based on data-dependent partitions[END_REF] for more general partitioning schemes. More precisely, some asymptotic conditions on the maximum number of points per cell, the complexity of the partition rule (related to Vapnik-Chervonenkis complexity [START_REF] Vapnik | The nature of statistical learning theory[END_REF]) and the cells' diameters were shown to guarantee the almost sure consistency of the estimator.

Nearest Neighbors

This class of methods is based on considering a number of neighbor points to estimate the quantity of interest at each point z ∈ R d .

Classification In order to build a locally averaging estimator, consider the following weighting functions

W n,i (z) = 1/k n if Z i is among the k n nearest neighbors of z 0 otherwise (1.31)
where k n is a parameter that can depend on n and the notion of nearest depends on the chosen distance and some tie-breaking rule.

The following theorem gives conditions to guarantee strong pointwise consistency [GK02, Thm. 25.17]:

Theorem. 1.2. If |R| < C for some C < ∞, k n log n → ∞ and k n n → 0,
then the k n -NN estimate using Euclidean distance is strongly pointwise consistent.

In [START_REF] Kpotufe | k-nn regression adapts to local intrinsic dimension[END_REF], it was shown that for regression functions that are λ-Lipschitz, the rate of convergence of a k n -NN regressor depends only on the intrinsic dimension. Moreover, the authors propose a method to locally choose k, which guarantees nearly optimal minimax rate.

The drawback of the exact nearest neighbors approach is that the number of neighbors need to grow with n, thus making it inappropriate in an online setting.

Divergence In [WKV09, PC08], the KL divergence is estimated by using a k-nearest-neighbor approach similar to the one use for density estimation. For each sample x i , consider the volume V xi of the ball B centered at x i enclosing the k-nearest neighbors, then the density estimate is:

f0 (x i ) ≡ k n 0 V xi . (1.32)
For consistent density estimation a growing k is required. With respect to the Bias-Variance tradeoff for finite sample sizes, as discussed in [START_REF] Devroye | A course in density estimation[END_REF][START_REF] Bernard W Silverman | Density estimation for statistics and data analysis[END_REF], a smaller k leads to a lower bias and a higher variance.

In [START_REF] Pérez-Cruz | Kullback-Leibler divergence estimation of continuous distributions[END_REF], it is shown that a a fixed k guarantees almost sure consistency of the KL divergence estimation while empirically outperforming methods based on consistent density estimation. In [START_REF] Wang | Divergence estimation for multidimensional densities via k-nearest-neighbor distances[END_REF], it is shown that using an adaptive k (as a result of fixing a radius instead of fixing k) reduces the bias in the estimation.

Kernel

As a dot product, a kernel function can be seen as a similarity measure between points (see, e.g., [START_REF] Scholkopf | Learning with kernels[END_REF]). Generalizing this idea, a kernel function is a weighting function that can be used to define a weighted neighborhood of a given point.

Classification

The Nadaraya-Watson kernel estimate (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]) is a locally averaging estimator that uses the following weighting functions

W n,i (z) ≡ k n (z, Z i ) n j=1 k n (z, Z j ) (1.33)
where the kernel functions k n have the form k n (z, Z i ) ≡ K z-Zi hn where the bandwidth h n depends only on the sample size. Usually, K(x) is symmetric with respect to 0, nonnegative and is large when x is small. For example, the naive (or window ) kernel is defined as follows

K(x) ≡ 1 x ≤1 .
(1.34)

Stone's theorem [START_REF] Charles | Consistent nonparametric regression[END_REF] allows to prove weak convergence under general conditions on h n and K (see [GK02, Sec. 5.2]). The following theorem gives conditions that guarantee strong pointwise consistency when using the naive kernel [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]Thm. 25.11].

Theorem. 1.3. For the naive kernel, assume that h n n→∞ ----→ 0 and nh d n / log n n→∞ ----→ 0 and |R| < C for some C < ∞, then the kernel estimate is strongly pointwise consistent.

An important technique that is worth mentioning is the popular Support Vector Machine (SVM) which is meant for hard classification (see [START_REF] Scholkopf | Learning with kernels[END_REF]). Using the kernel trick, the method looks for the hyperplane that best separates the two populations in the higher-dimensional space which allows nonlinear separation in the original space. The distance of a new point to be classified to this hyperplane can be interpreted as a measure of confidence of the classification. After calibrating according to the margin (i.e., the separation of the two populations yielded by the hyperplane), these scores can be converted to values in [0, 1] using Platt Scaling [P + 99], which learns a logistic regression model on the scores. As discussed in [START_REF] Langford | Estimating class membership probabilities using classifier learners[END_REF], no guarantee on the quality of this estimate as a probability can be guaranteed, essentially because the margin is not a sufficient statistic.

Alternatively, in [START_REF] Langford | Estimating class membership probabilities using classifier learners[END_REF], a probing method to convert any hard classifier into a soft one with guarantees (i.e., a small error for hard classification implying an accurate estimation of the label probabilities) has been proposed.

In [START_REF] Kivinen | Online learning with kernels[END_REF], online kernels methods applied, in particular, to classification and regression have been proposed.

Divergence An important kernel related divergence measure that yields a powerful two-sample test is the MMD that we introduced in Section 1.2.1.

A biased estimate is obtained from Eq. (1.10) by substituting expected values with samples averages

[GBR + 12, Lemma 6] MMD b [F, X, Y ] ≡   1 n 2 0 n0 i,j=1 k(x i , x j ) - 2 n 0 n 1 n0,n1 i,j=1 k(x i , y j ) + 1 n 2 1 n1 i,j=1 k(y i , y j )   1 2 (1.35) which cost O (n 0 + n 1 ) 2 to compute. This biased estimate convergences in probability at rate O (n 0 + n 1 ) -1 2
to its population value. Alternative quadratic and linear unbiased estimators for MMD 2 are also presented in [GBR + 12]. Although experimentally of slower convergence, the linear alternative is of particular interest for streaming computations since it only requires O (1) memory, while the quadratic versions require O (n 0 + n 1 ) to compute all interactions.

[RW11] elucidates interesting connections between MMD and SVM, where they also show that the further apart the distributions are, the easier the discrimination problem becomes.

Universal Prediction Methods for Soft Classification

One way of building a universal soft classifier is by considering a sequence of nested parametric models, i.e., sets of parametric distributions on infinite sequences s ∞ ≡ s 1 , s 2 , . . . where ]). Let us define M ≡ ∪ ∞ i=1 M i If we denote by M ⋆ the model to which the true distribution p ⋆ ∈ M ⋆ belongs, there are two cases of special interest :

s i ∈ S ⊆ R d , such that M 1 ⊆ M 2 ⊆ . . . where M i = {p θ : θ ∈ Θ i ⊆ R d } (see, e.g., [ Grü07
• the parametric case: M ⋆ ⊆ M k for some k,

• the nonparametric but approximable case: if we define the information closure M as the set of distributions for infinite sequences on S that can be arbitrarily well approximated by elements of

M, i.e., M ≡ {p ⋆ | inf p∈M D (p ⋆ p) = 0}, we have that M ⋆ \ M = ∅ and M ⋆ ⊂ M .
In our case, one could set X ≡ Ω × A and use joint distributions. Nevertheless, as noted in [vEGdR12, p. 382], for the methods we consider next, it is equivalent to set X ≡ A and use conditional distributions instead.

Next, we present two methods to achieve universality.

Bayesian Distributions

In the Bayesian approach, a prior distribution represents one's "belief" on which of the distributions of M generates the data. Then, as data are seen, this prior is updated to obtain a posterior distribution. This posterior can then be used to predict future data (see, e.g., [START_REF] Grünwald | The minimum description length principle[END_REF]p. 73]).

In the case of a sequence of nested models, a prior w(•) on the models M i of the sequence and a prior π(•) on the parameter space Ω i of each model are considered in order to mix the models in the following way:

P Bayes(M) (s n ) ≡ i w(i)P Bayes(Ωi) (s n ) (1.36) with P Bayes(Θi) (s n ) ≡ θ∈Θi π i (θ)p θ (s n )dθ.
(1.37)

Using Bayes' formula, one obtains the following predictive formulas:

P Bayes(M) (s n+1 |s n ) = i w(i|s n )P Bayes(Θi) (s n+1 |s n ) (1.38) with w(i|s n ) = w(i)P Bayes(Θi) (s n ) i w(i)P Bayes(Θi) (s n ) (1.39)
and

P Bayes(Θi) (s n+1 |s n ) = θ∈Θi π i (θ|s n )p θ (s n+1 |s n )dθ (1.40) with π i (θ|s n ) = π i (θ)p θ (s n ) θ∈Θi π i (θ)p θ (s n )dθ
.

(1.41)

In order to show that P Bayes(M) is universal in the parametric case, we denote the true distribution as p θ ⋆ ∈ M k for some k and observe that

P Bayes(M) (s n ) ≡ i w(i)P Bayes(Θi) (s n ) ≥ w(k)P Bayes(Θ k ) (s n ) (1.42) which implies that -log P Bayes(M) (s n ) ≤ -log w(k) -log P Bayes(Θ k ) (s n ). (1.43)
Then, in [START_REF] Li | Iterated logarithmic expansions of the pathwise code lengths for exponential families[END_REF], it was shown that for p θ ⋆ belonging to an exponential family with k-parameters (here, we assume, for simplicity, that the index in M k correspond to the number of parameters), with p θ ⋆ -probability 1,

R n P Bayes(Θ k ) , p θ ⋆ ≡ -log P Bayes(Θ k ) (s n ) + log p θ ⋆ (s n ) (1.44) ≤ k 2 log n 2π -log π k (θ ⋆ ) + log det I(θ ⋆ ) + C n log log n + o(1) (1.45)
where I(•) denotes the Fisher information matrix and C n is a random variable such that with probability 1 is nonnegative for all n, is upper-bounded by k for all n, and that takes on a value ≥ 1 for infinitely many n. Combining Eqs. (1.43) and (1.44) yields that 1 n R n P Bayes(M) , p θ ⋆ is upper bounded by a quantity that vanishes to zero almost surely. It is important to note that the convergence rate depends on the number k of parameters of the true distribution.

In the non-parametric case, the universality of P Bayes(M) has been proved for some model classes under certain conditions (see [START_REF] Grünwald | The minimum description length principle[END_REF]Section 13.4.1]).

The Switch-Distribution

Universality of Bayesian distributions stems from the fact that their performance is asymptotically as good as any of the distributions under consideration. Nevertheless, they suffer from the catch-up phenomenon, i.e., their convergence rate is not optimal. In [START_REF] Van Erven | Catching up faster by switching sooner: a predictive approach to adaptive estimation with an application to the aic-bic dilemma[END_REF], it has been shown that a better convergence rate can be achieved by allowing models to change over time, i.e., by considering instead of a set of distributions M a (larger) set constituted by sequences of distributions of M. This can be explained intuitively by the fact that when a small number of samples has been seen, models with fewer parameters work better than more complex ones even when the true model is more complex. Therefore, it can be beneficial to consider distributions that switch from simpler distributions to more complex ones. The distribution has a still a Bayesian form but the mixture is done over sequences of models.

Detection: Hypothesis Testing

In this section, we recall fundamental concepts used in statistical hypothesis testing with a focus on nonparametric two-sample testing.

Hypothesis Testing

Statistical tests allow to make decisions using samples by means of analyzing some global property-a statistic-, which is a function of them, and estimating the probability that the observed property value has occurred by chance alone.

Hypotheses In classical hypothesis testing, the setting is asymmetric. The null hypothesis H 0 is meant to represent the belief that the statistician is trying to refute. The probability of observing the global property value by chance alone is therefore computed assuming the null hypothesis is true. The alternative hypothesis H 1 represents the class of situations a test would be sensitive to in order to refute the null hypothesis.

When the hypotheses have a parametric form (with a parameter θ), i.e.

H 0 : θ ∈ Θ 0 H 1 : θ ∈ Θ 1 , (1.46)
two types of hypotheses can be considered. Hypotheses that are completely specified, i.e., |Θ • | = 1, are simple, otherwise they are composite. Nonparametric hypotheses include unknown or infinite parameters and therefore are an extreme case of composite hypotheses.

Statistic A statistic T n is a random variable that is a function of the observed samples, aiming at revealing the evidence against the null hypothesis.

In order to provide evidence against the null hypothesis, the statistic must take extreme values (either large or small).

p-value A valid p-value p ∈ [0, 1] is a statistic such that P H0 (p ≤ α) ≤ α, ∀α ∈ [0, 1].
(1.47)

An exact p-value is uniform under the null hypothesis, i.e.,

P H0 (p ≤ α) = α, ∀α ∈ [0, 1]. (1.48)
An exact p-value can be interpreted as the probability of observing a statistic at least as extreme as the one actually observed, assuming that the null hypothesis is true, where "more extreme" depends on the way the hypothesis is tested.

Rejection under significance level α The null hypothesis is rejected with significance α if there exists a threshold T α such that P H0 (T n ≤ T α ) ≤ α (we assumed that a small value of the statistic represents stronger evidence for rejecting the null hypothesis).

When using a p-value, this is equivalent to rejecting when p is under the threshold α, since P H0 (p ≤ α) ≤ α. The range of values of T n for which the null is rejected is called the rejection region.

Error types Since these tests are based on a finite number of samples, it is possible that an incorrect answer is returned. A Type I error occurs if H 0 is rejected when it is actually true (i.e., a false positive). Conversely, a Type II error occurs if H 0 is not rejected when it is actually false (i.e., a false negative). In general, there is a trade-off between the two error probabilities.

The significance level α is therefore an upper bound on the type I error probability. Nevertheless, when the alternative hypothesis represents a large class of probability distributions, it can be difficult or impossible to provide guarantees for type II error. For this reason, when not enough evidence is found to reject the null hypothesis, we avoid saying that the null hypothesis is accepted.

Power and consistency

The power of a test is defined as the test's probability of correctly rejecting the null hypothesis (thus, 1 -P (Type II error)) which, in other words, is the test's ability of identifying the alternative hypothesis. A test is consistent if its power tends to one when the number of samples tends to infinity. When using a p-value, consistency can be expressed as

P H1 (p ≤ α) n→∞ ----→ 1.
(1.49)

The Neyman-Pearson Fundamental Lemma Ideally, a Uniformly Most Powerful (UMP) test is desired, i.e., a test that, for all the distributions considered by the alternative, has the highest power for a given level α.

The following fundamental theorem describes which tests are UMP when both the null and the alternative hypotheses are simple (see, e.g., [CB01, p.388]).

Theorem. 1.4 (Neyman-Pearson Lemma [START_REF] Neyman | On the Problem of the Most Efficient Tests of Statistical Hypotheses[END_REF]). Consider testing H 0 : θ = θ 0 versus H 1 : θ = θ 1 given the samples x n = x 1 , . . . , x n , where the densities are f θ0 (•) and f θ1 (•) respectively. Then, the likelihood ratio test that rejects H 0 when

T n ≡ λ n ≡ f θ0 (x n ) f θ1 (x n ) ≤ η (1.50) where P H0 (λ n ≤ η) = α (1.51)
is the most powerful test at level α for a threshold η ≥ 0.

When the alternative hypothesis is composite, it is often not possible to construct UMP tests (see, e.g., [START_REF] Casella | Statistical inference[END_REF]8.3.19]).

Sequential Testing

In sequential hypothesis testing (also known as sequential analysis), the sample size is not fixed in advance and samples are processed sequentially until conclusive results are observed in accordance with a pre-defined stopping rule. This allows to potentially reduce the number of samples that need to be collected in order to reach a conclusion, thus reducing costs of data collection (financial, human, etc).

One classical example of such procedure is Wald's Sequential Probability Ratio Test (SPRT) [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF], which considers the following simple hypotheses:

H 0 : X has density f θ0 (•) H 1 : X has density f θ1 (•)
.

(1.52)

Then, given a sequence of i.i.d. observations x 1 , . . . , x i , . . . , the logarithm of the likelihood ratio is incrementally computed, i.e.,

S i ≡ S i-1 + log f θ1 (x i ) f θ0 (x i ) (1.53)
with S 0 ≡ 0. The stopping rule is based on two parameters -∞ < a < 0 < b < ∞ that control the desired Type I and Type II error probabilities and is defined as

     a < S i < b : continue S i ≥ b : accept H 1 S i ≤ a : accept H 0 . (1.54)
In practice, Wald proposes to use a ≡ log β 1-α and b ≡ log 1-β α , which guarantee an actual Type I error probability α ′ ≤ α and an actual Type II error β ′ ≤ β. Wald proved that this test finishes with probability one.

Wald's scheme relies on the perfect knowledge of the probability density function for each hypothesis. Recently, in [START_REF] Santiago-Mozos | On the uncertainty in sequential hypothesis testing[END_REF], this restriction has been somewhat relaxed for Bernoulli distributions, by allowing some uncertainty in the knowledge of the parameter p defining them.

For more material on statistical testing see, e.g., [START_REF] Casella | Statistical inference[END_REF][START_REF] Erich | Testing statistical hypotheses[END_REF].

The Two-sample Problem

Given samples from two distributions, a two-sample test aims at comparing two distributions by statistically asserting whether they are different in some aspect or not. For example, one could ask whether two distributions have the same variance or not, given a set of samples from each. It is also possible to consider more than two samples in which case it is referred as the k-sample problem.

The two-sample problem consists in determining whether the two distributions are equal or not. Parametric tests for the two-sample problem use some assumption on the nature of the distributions that generated samples (e.g., Gaussian), while nonparametric ones are meant to detect differences without assumptions. In the most general form, the null and the alternative hypotheses for the two-sample problem are

H 0 : f X = f Y a.e. H 1 : ¬H 0 .
(1.55)

Desired Properties and Negative Result

Type I Error Control In many cases, authors prove that the obtained tests are only asymptotically of level α. Nevertheless, for small sample sizes, an important feature of nonparametric two-sample tests is whether the distribution of the test statistic is distribution-free (i.e., independent of the sample's distribution) and whether an exact formula is known, in which case it is said that test is exactly of level α.

Another option is to obtain a bound for the threshold (which in terms of p-value would yield a valid p-value) which results in a conservative test with a Type I probability error actually lower than the nominal value α and consequently a lower power.

Type II Error and Consistency It is usually desired to have tests with proven consistency. In the nonparametric case, power and convergence speed depend on the nature of the dissimilarity, which makes no test superior than any other in every case.

In fact, even if a test for the two-sample problem is consistent, at a given fixed sample size, it is not possible to provide guarantees on the Type II error, without prior assumptions on the nature of the difference between the two distributions (see [GBR + 12] for a proving example).

Complexity Another important aspect, especially for high-dimensions and when the sample size is large, is the computational complexity involved in the calculation of the test statistic and its threshold.

Classical Nonparametric One-dimensional Two-sample Tests

The classical univariate Wald-Wolfowitz runs test is based on ordering the samples of the two populations and counting the number of runs of samples of the same population. A more powerful classical test is the one of Kolmogorov-Smirnov which is based on estimating the largest difference between the two unknown cumulative distributions. More precisely, the test statistic is D n0,n1 := max x | FX (x) -FY (x)|, where FX and FY are the empirical cumulative distributions.

The classical Mann-Whitney U test [START_REF] Henry | On a test of whether one of two random variables is stochastically larger than the other[END_REF] aims at detecting whether one random variable X tends to have larger values than another variable Y . More precisely, the test is consistent against the following alternative hypothesis: the probability of X exceeding an observation from Y (accounting for ties by adding 1 2 P (X = Y )) is greater than 1 2 [START_REF] Lehmann | Consistency and unbiasedness of certain nonparametric tests[END_REF], i.e., P (X > Y ) + 1 2 P (X = Y ) > 1 2 . This general alternative implies (see [START_REF] Fay | Wilcoxon-mann-whitney or t-test? on assumptions for hypothesis tests and multiple interpretations of decision rules[END_REF] for a discussion) more restrictive alternatives like X being stochastically larger than Y (i.e., the cumulative distribution functions satisfy F Y (a) > F X (a) for every a ) and locations shifts (i.e., F Y (a) = F X (a + δ), δ = 0).

The test can be one-sided, in which case the alternative is one of the previously discussed or two-sided, in which case the alternatives obtained by exchanging X and Y are also considered (i.e., Y tending to be larger than X).

The test is based on counting the number of times a sample from one population precedes a sample from the other and thus requires an ordering of the samples.

Next we review some multivariate tests for the two-sample problem against any alternative. It is important to remark that concepts like cumulative distribution functions, runs or ranks do not have a straightforward generalization for the multivariate case.

Multivariate Two-sample Tests Methods

Next, we discuss some approaches to design two-sample tests in R d and examples from the literature.

Extensions from One-dimensional Tests Although not straightforward, some extensions of 1D tests have been attempted in the literature.

For the Mann-Whitney test, the order can be obtained by looking at the sign of the differences between every pair of samples. Multivariate extensions of this test have been proposed in [START_REF] Möttönen | Multivariate spatial sign and rank methods[END_REF][START_REF] Oja | Multivariate nonparametric tests[END_REF], by extending the sign concept to the spatial sign function defined as

S(x) ≡ x -1 x, x = 0 0, x = 0 (1.56)
and then considering the pooled sample z n ≡ x 1 , . . . , x n0 , y 1 , . . . , y n1 and the spatial centered rank

R i ≡ 1 n n j=1 S(A x (z i -z j )) (1.57)
where A x is a data based transformation chosen to make the test affine invariant. Then,

Rk ≡ 1 n k n k j=1 R (j+ l<k n l ) (1.58)
is a measure of how much the k-th sample differs from the pooled one. The statistic is based on combining this measure for each population:

U 2 ≡ d c 2 x 1 k=0 n k Rk 2 (1.59) where c 2 x ≡ 1 n n i=1 R T i R i .
This test can be easily extended to the k-sample problem by extending the sum to k summands.

Notice that, contrarily to other tests based on interpoint distances that we will discuss below, this test depends only on the direction of the differences between points and, in this sense, it generalizes the Mann-Whitney test. Some other generalization of this type of tests have been proposed in [START_REF] Li | New nonparametric tests of multivariate locations and scales using data depth[END_REF] (along with some scale tests), using the data depth concept which provides an ordering of the sample data with respect to a generalized median point (see, e.g., [LPS99, ZS00, Ser06]).

Friedman and Rafsky propose the minimal spanning tree (MST) of the pooled sample set (considering some distance between the points) as a multivariate generalization of the univariate sorted list. Therefore, one of their MST based test generalizes the Wald-Wolfowitz runs test by defining a run as a subtree whose vertices belong to the same sample set. Therefore, the test statistic is the number of edges that connects points of X to points of Y . The computational cost of this method using Kruskal's algorithm for building the MST is O (n 0 + n 1 )

2 log (n 0 + n 1 ) .

An attempt to generalize Kolmogorov-Smirnov test is proposed in [START_REF] Peacock | Two-dimensional goodness-of-fit testing in astronomy[END_REF]. The basic idea is to consider the largest difference between the empirical cumulative distributions considering the 2 d ways of defining them. 4 For example, in 2D, given the usual coordinates x and y, we consider P (< x, < y),P (< x, > y), P (> x, < y) and P (> x, > y). Thus, it is necessary to look at every possible step of the empirical cumulative distributions, in every possible order. A step occurs when a point is encountered when cumulating along one of the coordinates. Thus, (m + n) d steps must be considered (corresponding to selecting one coordinate value from each possible point). This corresponds to analyzing each quadrant at each center that yields a different step: there are then 2 d (m + n) d possible quadrants. This procedure is of exponential complexity in d and therefore has been essentially proposed for 2 dimensions. Some variations and improvements of this test have been proposed (see [START_REF] Fasano | A multidimensional version of the kolmogorov-smirnov test[END_REF][START_REF] Lopes | The two-dimensional kolmogorov-smirnov test[END_REF]), but still, essentially, for 2 or 3 dimensions.

Divergence based Naturally, any statistic converging to a measure of dissimilarity as those discussed in Section 1.2.1 can be used for a two-sample test.

Two recent examples are the MMD test of [GBR + 06, GFHS09, GBR + 12] and the KL divergence based test of [START_REF] Pérez-Cruz | Kullback-Leibler divergence estimation of continuous distributions[END_REF].

Projection/Scoring/Ordering based A general method to build two-sample tests is to map the multivariate to a one-dimensional space using some kind of projection, scoring or ordering. Then, if the first step does not depend on the populations but just on the pooled sample, one can apply a classical one-dimensional test on the obtained data.

For example, in [START_REF] Friedman | Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests[END_REF], an ordering is obtained by some traversal of the minimum spanning tree of the pooled sample.

In [START_REF] Cuesta-Albertos | A sharp form of the cramer-wold theorem[END_REF], it is shown that, under suitable regularity conditions, just one randomly chosen onedimensional projection determines a distribution on R d . More precisely, if F X = F Y , the distributions of the projected variables obviously coincide but, if F X = F Y , the probability that they coincide is zero. Thus, the proposed method consists in projecting the sample points in a one-dimensional space spanned by some random vector h and then applying the univariate Kolmogorov-Smirnov test to the projected points. Because of the random nature of the procedure, the authors suggest to repeat it with several projections in order to increase the power, since given a sample, not all the projections are equally good at revealing the differences, because of the natural loss of information incurred when projecting. When combining the results, it is necessary to take into account the multiple testing problem (i.e., the increased probability of getting false positives) by, for example, correcting it via the Bonferroni [START_REF] Bonferroni | Teoria statistica delle classi e calcolo delle probabilita[END_REF] or the Benjamini-Hochberg-Yekutieli [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] procedures.

As remarked in [START_REF] Friedman | On multivariate goodness-of-fit and two-sample testing[END_REF], any soft classifier (or, more generally, any scoring function) trained on a separate data set can be used to build a two-sample test, again, using a one-dimensional test on the probabilities (or scores). If the same set is used for training and evaluating, then the procedure is no longer valid but one can resort to a permutation method as will see later.

In this spirit, there is the Area Under the Curve (AUC) based test [START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF] that splits the data into two sets. The first step consists in training a score function using the AUC maximization criterion on the first dataset. The second step consists in scoring the samples from the second dataset using the trained score function to finally apply a Mann-Whitney test.

Category based Another interesting approach suggested in [START_REF] Friedman | Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests[END_REF], consists in dividing the pooled samples into two mutually exclusive categories, based on any criterion not involving the X or Y labels. Then a test like Fisher's [START_REF] Fisher | On the interpretation of χ2 from contingency tables, and the calculation of p[END_REF] can be applied to detect a tendency of one of the two sample set to be in one of the categories. For instance, Friedman and Rafsky proposed (for comparison purposes) to use the following criterion for defining the categories: based on the MST of the pooled sample, they count the number of nodes of degree 1 for each label.

Of course, two categories are far from enough to detect any kind of difference. Nevertheless, the procedure can be extended to more categories using the G-test [START_REF] Kullback | Information Theory and Statistics[END_REF][START_REF] Kullback | Tests for contingency tables and markov chains[END_REF].

Interpoint Distances based Now, we review some interpoint distance based tests.

In [START_REF] Maa | Reducing multidimensional two-sample data to one-dimensional interpoint comparisons[END_REF], it is shown that, under some assumptions on the interpoint distance function h(•, •), a multidimensional two-sample test can be reduced to one-dimensional interpoint distance testing, i.e., testing whether the distributions of h(X, X ′ ), h(X, Y ) and h(Y, Y ′ ) are equal.

The MMD based test [GBR + 06, GFHS09, GBR + 12] also consider interpoint similarities by means of a kernel function.

Another related test is [START_REF] Baringhaus | On a new multivariate two-sample test[END_REF], which considers the following statistic:

T n0,n1 ≡ n 0 n 1 n 0 + n 1   1 n 0 n 1 n0 j=1 n1 k=1 X j -Y k - 1 2n 2 0 n0 j=1 n0 k=1 X j -X k - 1 2n 2 1 n1 j=1 n1 k=1 Y j -Y k   . (1.

60)

Local Coincidences based There is a long tradition of defining two sample tests based on certain local coincidences, i.e., with some definition of neighborhood, these tests combine the number of samples from each population to determine the global level of mixture of the two populations. For example, k-nearest-neighbor approaches have been proposed and analyzed in [FST73, Sch86, Hen88] . The statistic is usually defined as follows. Denoting by Z i the pooled samples, the following is defined

I i (r) ≡
1, if the rth nearest neighbor of Z i belongs to the same sample as Z i 0, otherwise .

(1.61)

Then, denoting n ≡ n 0 + n 1 , the statistic is defined as

T k,n ≡ 1 nk n i=1 k r=1 I i (r) (1.62)
which is simply the proportion of all k nearest-neighbor comparisons in which a point and its neighbor are member of the same sample. Under the null hypothesis, a low value is expected, since samples should be well mixed.

Another notable member of this category is Rosenbaum's Cross-Match test [START_REF] Rosenbaum | An exact distribution-free test comparing two multivariate distributions based on adjacency[END_REF] whose construction allows an exact characterization of its null distribution, as we will discuss later. First, a minimum distance non-bipartite matching is computed over the pooled sample. Supposing that n = n 0 + n 1 is even5 , a nonbipartite matching divides the n points into I = n/2 pairs of two points in such a way as to minimize the total of the I distances within the n/2 pairs. Then, the test statistic is the number of pairs containing a sample from both X and Y . This procedure has a time complexity of O (n 0 + n 1 )

3 .

In [START_REF] Biau | On the asymptotic properties of a nonparametric l1-test statistic of homogeneity[END_REF], a partition of the space is performed and the statistic is the sum, over all the cells, of the absolute value of the difference of counts between the two populations (i.e., the L 1 distance).

Note that many of the KL divergence estimators presented in Section 1.3.3 fall in this category.

Null Distribution

Given some statistic that aims at revealing differences between distributions if any, the next step is to characterize its distribution under the null hypothesis which must be distribution free, since nothing is assumed about f X and f Y .

Exact Distributions In some cases, it is possible to give an exact formula for the null distribution. Therefore, in these cases, the empirical type I error, as the number of trials tend to infinity, will be close to the specified α level.

A common way to obtain such exact formulas is when the statistic is based on ranks (but not directly on the sample values) to consider all distinct label permutations n0+n1 n0 for the given sample sizes. A classic example is the Mann-Whitney test where all possible permutations of labels are implicitly considered using a recursive formula. When the sample size is large it gets difficult to compute and one usually resorts to an asymptotic normal approximation.

Another example is Rosenbaum's Cross Match test [START_REF] Rosenbaum | An exact distribution-free test comparing two multivariate distributions based on adjacency[END_REF], where an exact formula is provided. Nevertheless, as for the Mann-Whitney test, it is difficult to compute for large sample sizes, therefore the author proposes an asymptotic approximation.

Note that there can be a mismatch between the empirical Type I error and the design level when ties are possible (i.e., when different permutations yield the same statistic value).

Non-asymptotic Bounds When exact distributions are difficult to derive, it is sometimes possible to find thresholds for which the type I error probability is bounded by α. In [GBR + 12], several bounds are provided for the MMD, but, as the authors note, those are distribution-free bounds that hold even in the worst case scenario and therefore these tests are conservative by design.

Asymptotic Distributions Asymptotic distributions are usually easier to obtain than exact ones and many authors propose to use them as approximations. Note that this is different from approximating a known exact distribution. In [GBR + 12], a test based on an asymptotic null distribution for an unbiased estimate of MMD is proposed, which yields a greater power than the tests based on bounds.

For the multivariate runs test proposed in [START_REF] Friedman | Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests[END_REF] and further analyzed in [START_REF] Henze | On the multivariate runs test[END_REF], a normal distribution allows to characterize the asymptotic null distribution. Normal based asymptotic distributions are also proposed for the k-nearest-neighbors based tests [START_REF] Jerome H Friedman | A nonparametric procedure for comparing multivariate point sets[END_REF][START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF][START_REF] Henze | A multivariate two-sample test based on the number of nearest neighbor type coincidences[END_REF] .

Bootstrapping Another approach is to approximate the unknown null distribution using the empirical distribution. The bootstrap method draws, with replacement, n samples from the empirical distribution defined by the pooled sample. Authors that use bootstrap generally prove that their tests are asymptotically of level α (see, e.g., [START_REF] Fromont | Kernels based tests with non-asymptotic bootstrap approaches for two-sample problem[END_REF] for a discussion).

In [START_REF] Baringhaus | On a new multivariate two-sample test[END_REF], the asymptotic null distribution of the test statistic is derived and shown to be the limit of the bootstrap distribution. In [GBR + 06, GBR + 12], a consistent bootstrap method is also proposed to estimate the threshold.

Permutation Methods The permutation method (see, e.g., [START_REF] Phillip | Permutation, parametric and bootstrap tests of hypotheses[END_REF]), also called without-replacement bootstrap, is an alternative approach to obtain a null distribution. It relies on the assumption that, under H 0 , samples are exchangeable. Instead of considering analytically all the possible permutations, permutation methods actually permute the labels and recompute the statistic each time.

Even for moderate sample sizes it can be impractical to consider all possible permutations. Therefore, authors usually propose to use a Monte-Carlo approach, i.e., to consider only a subset of all possible permutations. 6Contrary to the standard bootstrap, permutation methods with correct handling of ties (e.g., by breaking them at random) and proper formulas yield valid p-values and not approximations (see, e.g., [E + 04, PS10]). More precisely, let m the number of permutations performed (these can be sampled with or without replacement from the space of permutations), then denote by b the number of times that the statistic yielded by a permutation is more extreme than the statistic yielded by the actual observed data, then a valid p-value is obtained by the following formula

p = b + 1 m + 1 . (1.63)
When permutations are drawn without replacement and there are no ties, this formula gives an exact p-value, otherwise it provides an upper bound. In [START_REF] Phipson | Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn[END_REF], an exact formula is proposed for the latter case. This perfect control of the type I error probability is appealing. Nevertheless, it is usually not straightforward to prove consistency for permutation methods.

With respect to the number of permutations m, it has been shown in [START_REF] Jockel | Finite sample properties and asymptotic efficiency of monte carlo tests[END_REF] in the context of parametric tests, that there is a loss of power when using permutations but this loss decreases as m increases.

The permutation method is a quite general one that can be applied to any statistic. Many authors, in addition to other methods, include it in their works, for example [START_REF] Hall | Permutation tests for equality of distributions in highdimensional settings[END_REF][START_REF] Friedman | On multivariate goodness-of-fit and two-sample testing[END_REF][START_REF] Aslan | New test for the multivariate two-sample problem based on the concept of minimum energy[END_REF][START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF]. In [START_REF] Gretton | A fast, consistent kernel two-sample test[END_REF], the implementation called MMD boot uses a permutation method too.

Effect Size

The effect size is a general concept in statistics meaning a quantitative measure of the strength of a phenomenon, e.g.: the correlation between two variables, the regression coefficient, the mean difference, etc. See, e.g., [START_REF] Gail | Using effect size-or why the p value is not enough[END_REF] for a recent discussion on the practical importance of the effect size.

In the case of a two-sample comparison, if a significant difference is found, it might be slight or on the opposite, very large. When using real data, if large sample sizes are used, it is reasonable to expect to find significant differences due, for example, to some sampling variability. To know if an observed difference is not only statistically significant but also important (or meaningful), one needs to quantify the effect size.

For example, when considering restricted alternatives, like a difference in mean or variance, an estimate of the magnitude of the difference in these parameters is very important in addition to the statistical significance of the result.

The classical Mann-Whitney test can be complemented by the Hodges-Lehmann estimator ∆ [START_REF] Hodges | Estimates of location based on rank tests[END_REF], which estimates the median of the difference between a sample X and a sample Y by computing the median of the difference in the cartesian product of the values in the two samples. This quantity has the following critical property: if the values of one population are shifted by ∆, then, the difference between the two populations vanishes-i.e. there is no evidence to reject the null.

The multivariate extensions to the Mann-Whitney test presented in Section 1.4.2 propose analogues to the Hodges-Lehmann estimator based on the spatial centered rank (see, e.g., [START_REF] Oja | Multivariate nonparametric methods with R: an approach based on spatial signs and ranks[END_REF]). Nevertheless, for more general multivariate differences, the question is much harder and has not been addressed in the literature except for some mentions like in [START_REF] Friedman | On multivariate goodness-of-fit and two-sample testing[END_REF].

Contributions

The contributions of this thesis are organized in three chapters.

Chapter 2 In this chapter, we address the question that naturally comes after a two-sample test rejects the null hypothesis: why the two data sets do not look like coming from the same distribution? In other terms, we address the problem of expressing an effect size for the multivariate nonparametric two-sample problem. This chapter presents a two-stage method providing feedback on this difference.

The first step consists in estimating on each sample point a local discrepancy measure that is a decomposition of the Jensen-Shannon divergence. As discussed, in Section 1.2.1, this measure has several desirable properties when using θ 0 = 1 2 , in particular, it is symmetric and bounded between [0, 1]. Moreover, it can be expressed as the expected value over z ∈ R d with respect to the average density f of the KL divergence between the conditional distribution P (L = •|Z = z) and the Bernoulli distribution with parameter θ 0 = 1 2 . This KL divergence is thus taken as a measure of local discrepancy. The estimation of the conditional probability (i.e., a soft classifier) is done by means of a strongly pointwise consistent regressor based on k-nearest neighbors as explained in Sections 1.2.2 and 1.3.3.

In the second step, we use Morse theory which is concerned with the study of a height function h defined on a manifold. In particular, it studies the variation of the topology of sublevel sets, i.e., sets of the form h -1 (] -∞, λ]) where λ is some level. The height function studied is defined by the aforementioned discrepancy. Practically, we study a discretized version of this function using a nearest-neighbor graph connecting the sample points-the elevation of a sample being the estimated discrepancy. Only critical points of index 0 and 1 (corresponding, respectively, to minima and saddles) are considered, as in the Watershed Transform (used in image processing for segmentation). Then, Topological Persistence allows to identify persistent local minima of this height function (i.e., basins that last long before merging into others and thus represent more stable features) and to simplify it as clusters defining regions of points with high discrepancy and in spatial proximity.

Experiments are reported both on synthetic and real data (satellite images and handwritten digit images), ranging in dimension from d = 2 to d = 784, illustrating the ability of our method to localize discrepancies.

Chapter 3 In Section 1.2.2, we discussed how soft classification is intimately related to dissimilarity. In Section 1.4.2, we mentioned how the AUC based test trains a scoring function on a portion of the dataset and then scores the remaining samples using this function in order to obtain the final statistic.

In this chapter, instead of using the train/test paradigm, we propose a sequential procedure for nonparametric two-sample testing using online soft classifiers. More precisely, we show that any online soft classifier can be turned into a sequential two-sample test for which a valid p-value can be computed, yielding controlled type I error. Thanks to its sequential nature, our test benefits from the property of optional stopping meaning that the user can decide to stop at any time, always obtaining a valid p-value. This is in contrast to Wald's procedure (discussed in Section 1.4) that requires the rule to be fixed in advance and to the linear version of MMD [GBR + 12] which, even if processing is done in a sequential way, requires to fix in advance the sample size at which the rejection decision will be evaluated. This makes our test the first truly sequential nonparametric two-sample test.

We also show how to build consistent tests and, in particular, how to build them from a k-nearest neighbors nonparametric regressor. We also show that Bayesian mixtures and switch distributions can be used to increase power, while keeping consistency.

Chapter 4 Following the lines of Chapter 3, we build a sequential two-sample test using a spatial partition based sequential soft classifier instead of a nearest neighbors based one. The soft classifier used is semi-supervised meaning that, in an initial step, it uses an unlabeled set (e.g., a mix from the two populations being compared) to build the spatial partitions.

More precisely, our approach is based on Bayesian and Switch mixtures on an ensemble of partition based soft classifiers created by means of random projection trees [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF] (see Section 1.3.3). The partitions discretize the density functions that are being compared, to obtain corresponding probability mass functions.

If the chosen partitions reveal the differences between densities (i.e., probability mass functions are different), then the corresponding two-sample test is proved to be consistent.

The algorithm can be used in a streaming context since it has, per sample, a negligible memory footprint and logarithmic time complexity w.r.t. the unlabeled set size constant time complexity, which remains constant.

Importantly, this construction automatically adapts to the scale at which differences occur and does not require parameters to be set.

Chapter 2

Beyond Two-sample-tests: Localizing Data Discrepancies in High-dimensional Spaces 2.1 Introduction

Comparing datasets in high dimensional spaces

Datasets represented as point clouds are ubiquitous in science and engineering, used for applications in 2D and 3D space (e.g. to represent laser scans in mock-up design) as well as in high dimensional spaces (e.g. to represent images and documents, physical or biological phenomena, etc). In manipulating such data, several classes of questions are faced, such as matching, topological inference, or comparison. This latter endeavor, which is the topic of this work, calls for a discussion in three directions, namely statistics (two-sample tests), information theory and learning (divergence estimation), and geometry-topology.

From a statistical standpoint, a broad class of comparison methods, requiring tame assumptions on the data are nonparametric two-sample tests (TST), see, e.g., [GBR + 12] and the references therein. In a nutshell, a TST is a statistical hypothesis test checking whether the two sets can be seen as i.i.d. samples of an identical unknown distribution (the null hypothesis). In accepting or rejecting the null hypothesis, under a level of statistical significance α, a TST summarizes the body of information encoded in the points' coordinates into a single boolean value [START_REF] Friedman | On multivariate goodness-of-fit and two-sample testing[END_REF]. However, this boolean information is in general of limited interest, for several reasons. First, it is unlikely that two real life datasets come from exactly the same distribution. Since consistent TST detect any kind of difference of any size if enough samples are given, the rejection of the null is expected. Second, the magnitude (and the nature) of the differences, known as effect size, usually conveys more information than the mere presence of a difference. Therefore, a reject decision should be just a signal to examine further the data in order to understand. Developing a notion of effect size for non-parametric two-sample tests in high dimensions has not been explored yet, and is the goal of this work.

From an information theoretical and probabilistic standpoint, the comparison can be phrased as the problem of estimating the global Kullback-Leibler divergence between unknown distributions using the samples in hand. For example, a difference between images has been proposed [START_REF] Filipovych | Semi-supervised cluster analysis of imaging data[END_REF], by coupling a univariate Kullback-Leibler estimate (per pixel) and a decomposition of the discrepancy map thus defined using a watershed transform in image space. In a more general setting, there exist techniques to estimate this quantity that avoid density estimation in high dimensions (see, e.g., [WKV05, WKV09, PC08, SN07]) and that could be amenable to decomposition, so as to determine a contribution of individual points or of groups of points. Nevertheless, this divergence lacks important properties (symmetry and boundedness), which makes it more difficult to further process it.

Finally, the comparison can also be tackled from the geometric and topological perspectives. In geometric terms, one may compute some distance or matching (one-to-one, one-to-many, many-to-many) between the data, see [START_REF] Ding | FPTAS for minimizing earth mover's distance under rigid transformations[END_REF] and the references therein. While this procedure is informative for the two datasets in hand, the main difficulty consists in accommodating a probabilistic setting. In a more topological perspective, persistence theory [START_REF] Edelsbrunner | Computational topology: an introduction[END_REF], which aims at assessing the stability of topological features-generators of persistent homology groups, was recently used to compare persistence landscapes [START_REF] Bubenik | Statistical topological data analysis using persistence landscapes[END_REF]. Such comparisons are clearly important since oblivious to geometric transformations, but our focus is clearly on geometry dependent features.

Contributions

This work proposes, to the best of our knowledge, the first attempt to model the differences (the effect size, see above), between two datasets for which one has rejected the null hypothesis stipulating that they share the same underlying distribution. In a nutshell, we aim at clustering samples, based on two criteria, namely samples within a cluster should (i) contribute significantly to the difference between the two clouds, and (ii) form a connected region. We match these goals by two major steps. In the first stage, we model pointwise differences, which we call discrepancies, using the Jensen-Shannon divergence (JSD), which is symmetric and can be decomposed in terms of the conditional probability of belonging to one of the populations given a space location. This conditional probability can be naturally estimated using known techniques of non-parametric regression like the one based on k n nearest neighbours, which possesses strong asymptotic guarantees of consistency. In the second stage, using a nearest neighbor graph defined over the samples, we study the height function defined by the estimated discrepancy, and design a clustering procedure based upon topological persistence. Finally, we summarize the information by a bar plot giving the profile of the discrepancy aiming at representing a multivariate nonparametric effect size.

This work is organized as follows. Sections 2.2 and 2.3 respectively present the two major steps. Section 2.4 presents the construction of the bar plot and other plots that are useful to visualize the different steps of the method. Finally, Section 2.5 present experiments.

Estimating the discrepancy between datasets

We aim at modeling the discrepancy between two datasets x (n0) ≡ {x 1 , . . . , x n0 } and y (n1) ≡ {y 1 , . . . , y n1 }, in some fixed dimension Euclidean space R d . We view these data as coming from two unknown densities f X and f Y with corresponding cumulative distributions functions F X and F Y .

Jensen-Shannon divergence decomposition using conditional distributions

Let D KL (f g) be the Kullback-Leibler divergence (KLD) between two densities f and g defined as

D KL (f g) ≡ R d f (x) log f (x) g(x) dx (2.1)
with the conventions 0 log 0 = 0 and 0 log 0 0 = 0. The Jensen-Shannon divergence (JSD) defined in [START_REF] Lin | Divergence measures based on the Shannon entropy[END_REF], allows to symmetrize and smooth the KLD by taking the average KLD of f X and f Y to the average density f ≡ (f X +f Y ) /2, that is:

JS (f X f Y ) ≡ 1 2 (D KL (f X f ) + D KL (f Y f )) (2.2)
In addition to being symmetric, the JSD is bounded between 0 and 1 and its square root yields a metric. Note also that by taking the average, two random variables are implicity defined: a position variable Z with density f Z ≡ f and a binary label L that indicates from which original density (i.e. f X or f Y ) an instance of Z is obtained. Formally, considering the alphabet A = {0, 1} and X ∼ F X , Y ∼ F Y , the following pair of random variables is defined:

(L, Z) = (0, X) with prob. 1 2 (1, Y ) with prob. 1 2 .
(2.3)

In the sequel, we will consider the conditional and unconditional mass functions P (l|z) = P (L = l|Z = z) and P (l) = P (L = l) =1 2 respectively, as well as the joint probability density f L,Z . We will also use the notation f l to denote f X (resp. f Y ) when l = 0 (resp. l = 1). Before establishing lemma 2.1, a key property for our comparison problem, we recall the definition of the Kullback-Leibler divergence between two discrete distributions P and Q over A :

D KL (P Q) ≡ l∈A P (l) log P (l) Q(l) .
(2.4)

Lemma. 2.1. One has: JS (f X f Y ) = R d f Z (z)D KL (P (•|z) P (•)) dz.
(2.5)

Proof of lemma 2.1. Recall that the JSD can be expressed as follows:

JS (f X f Y ) ≡ 1 2 (D KL (f X f Z ) + D KL (f Y f Z )) = 1 2 R d f X (z) log f X (z) f Z (z) dz + 1 2 R d f Y (z) log f Y (z) f Z (z) dz .
By linearity of integration and by Bayes' Theorem, and noting that the conditional densities f (z|L = 0) = f X (z) and f (z|L = 1) = f Y (z) we have then

JS (f X f Y ) = R d l P (l)f (z|L = l) log f (z|L = l) f Z (z) dz = R d l P (l) f L,Z (l, z) P (l) log f L,Z (l, z) P (l)f Z (z) dz = R d f Z (z) l P (l|z) log P (l|z) P (l) dz ≡ R d f Z (z)D KL (P (•|z) P (•)) dz.
The previous lemma shows that the JSD can be seen as the average, over z ∈ R d , of the KLD between the conditional and unconditional distribution of labels. More formally, we define: Definition. 2.1. The discrepancy at location z is defined as the KL divergence:

δ(z) ≡ D KL (P (•|z) P (•)) .
(2.6) Note that δ(z) ranges between 0 and 1 and is 0 iff f X (z) = f Y (z). Note also that since P (l) is known but P (l|z) is not, the problem we consider now is the one of estimating P (l|z) at each given location z.

Conditional probability estimation via non-parametric regression

In order to estimate the conditional distributions, we can use random design 1 non-parametric regression.

Generic Framework

First, we define the basic concepts (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF] for more details). Definition. 2.2. Given a random vector (Z, R), where Z ∈ R d and the response variable R ∈ R, the regression function is defined as

m(z) = E [R|Z = z] .
(2.7)

In the regression problem, the goal is to build an estimator m n (z) of m(z) using a set of n i.i.d. realizations of (Z, R).

With respect to the guarantees that are provided for regressors, usually, the L 2 risk or mean squared error is considered, i.e.,

R d |m n (z) -m(z)| 2 µ(dz) (2.8)
where µ denotes the distribution of Z. Nevertheless, since our goal is to estimate the discrepancy δ(z), we seek pointwise guarantees for regressors. In particular, we will consider a strong form of consistency which is defined as follows.

Definition. 2.3. Denoting µ the distribution of Z, a sequence of regression function estimates {m n } is strongly pointwise consistent (s.p.c.)

if for µ-almost all z ∈ R d m n (z) n→∞ ----→ m(z) a.s.. (2.9)
In [GK02, Sec. 25.6], some s.p.c. regression estimates are presented. For example, regression estimates based on partitioning, kernel and nearest neighbors are s.p.c under certain conditions for their parameters, when the absolute value |R| < M , for some M . Now we describe the s.p.c. k n -nearest neighbor regression function estimate (see [GK02, Ch.6&25] for further details). Given the training data {Z i , R i } i=1,...,n , let us denote as R (i,n) (z) the response value corresponding to i-th nearest neighbor (with some tie-breaking rule) of z in Z (n) . Then, the k n -nearest neighbor (k n -NN) regression function estimate is defined by

m n (z) = 1 k n kn i=1 R (i,n) (z).
(2.10)

Then we have the following theorem [GK02, Thm. 25.17]:

Theorem. 2.1 (Strong pointwise consistency of k-NN). If |R| < C for some C < ∞, k n log n → ∞ and k n n → 0,
then the k n -NN estimate using Euclidean distance is strongly pointwise consistent.

Application to Conditional Probability and Discrepancy Estimation

In order to apply this framework to our problem, note that the correspondence R = L yields m(z) = P (1|z).

(2.11)

Then, we can use the following estimator for P (l|z)

Pn (l|z) ≡ |1 -l -m n (z)|.
(2.12) Note that it is required that 0 ≤ m n (z) ≤ 1, since we aim at estimating a conditional probability, and that it is satisfied by the k n -nearest neighbor regressor. Using Eq. (2.12), we finally obtain an estimator for δ(z) : δn (z) ≡ D KL Pn (•|z) P (•) .

(2.13) Proof. Let us write δn (z) = h(m n (z)), with h(x) ≡ x log (x/(1/2)) + (1x) log ((1x)/(1/2)). It is easy to see that h(x) is a continuous function (composition, sum and product of continuous functions). Then, we apply the continuous mapping theorem [START_REF] Billingsley | Convergence of probability measures (2nd Edition)[END_REF], on every z ∈ R d where m n (z) n→∞ ----→ m(z) to complete the proof.

Remark 2.1. Theorem 2.1, leaves many possible choices for k n . For example, in [Gyo81, CD14], convergence rates of nearest neighbor based regressors and classifiers are studied. Roughly speaking, when m(z) = P (1|z) satisfies an α-Holder like condition, it is shown optimal convergence rates, assessed by upper and lower bounds on classification rates, are obtained when k n = 2α 2α+d .

Joint distribution compatible sampling

In the regression framework, the samples must be i.i.d. from a joint distribution f L,Z . In our original problem, we have two sets of samples drawn independently from f X and f Y . In order to ensure this condition, we will use the random multiplexer depicted in Fig. 2.1. On each input it receives i.i.d. samples from each of the populations. Then, it generates an instance l of L. According to the value of l (0 or 1), it consumes the corresponding input (X, Y resp.) and outputs it along with l.

The following lemma shows that the output has the desired joint density.

Lemma. 2.2. An output pair from the random multiplexer has joint density f L,Z .

Proof. The joint density of an output pair (l, z) is

f (l, z) = f (z|L = l)P (l).
Since, by construction, f (z|L = l) = f l (z) and, by Eq. (2.3),

f l (z) = f (z|L = l), then f (l, z) = f (z|L = l) P (l) = f L,Z (l, z).
Remark 2.2. In practice, there is a finite set of i.i.d. samples of X and Y available. Then, at some point the multiplexer can have no more data to consume on one of the inputs while there is still data available on the other one. Therefore, some samples of the original sets would not be used and, thus, some loss of information is to be expected. To favor the inspection of all samples, we resample B times as follows:

1. For j ∈ 1..B do:

(a) Generate a sequence {z j i , l j i } i=1,...,mj using the random multiplexer, until one of the datasets is exhausted. Let m j be the number of samples collected.

(b) Build δj mj (z) using some consistent estimator trained on {z j i , l j i } i=1,...,mj 2. Define δ(z) ≡ median j∈1..B ( δj mj (z))

Notice that δ(z) is also a consistent estimator.

Remark 2.3. Although a random multiplexer that uses a Bernoulli distribution with parameter θ 0 = n 0 /n instead of 1/2 would exempt us from using the resampling procedure, this would require considering the Jensen-Shannon divergence with arbitrary weights (see [Lin91, Sec. IV]). The resulting effect size would lack the symmetry and normalization properties, therefore making the comparison between different effect sizes harder.

Localizing the discrepancy 2.3.1 Defining clusters from sublevel sets

Overview. We wish to identify groups of samples, called clusters, intuitively characterized by two properties (Fig. 2.2):

• first, the discrepancy of such samples should be significant;

• second, samples within a cluster should be associated with regions where the discrepancy peaks.

To meet the first goal, we assume the existence of a value δ max stipulating that below δ max , the discrepancy is not significant. As we shall see in Experiments, while this value is not unique in general, few clusters typically stand out, imposing a lower bound on δ max .

To meet the second goal, we use a strategy reminiscent of mode clustering [CM02, CadOS11], a general clustering strategy consisting of defining a cluster from the catchment basin of a local maximum of a density estimate, as well as watershed transforms used in image processing [START_REF] Roerdink | The watershed transform: Definitions, algorithms and parallelization strategies[END_REF]. Our approach uses the hyper-surface defined by the discrepancy -δ(z) (Eq. (2.6)). In the sequel, we refer to this hypersurface as the landscape of the height function, instinctively. In studying this landscape, clusters shall be defined from the catchment basins of selected local minima of the landscape.

In the sequel, we present the required background in the continuous setting, and further explain how to deal with the discrete datasets at hand.

Background.

As stated above, we wish to focus on samples of the landscape whose discrepancy is less than -δ max (Fig. 2.2). The region of the landscape matching this specification is called the sublevel set associated with -δ max . We define a cluster as a connected component (c.c.) of this sublevel set. We also note that each such cluster may be associated with the lowest local minimum of this connected component. However, in case of noisy estimates, many small clusters may be reported. To ease the interpretation of results, we therefore smooth the landscape, and focus on c.c. of the smoothed landscape.

The smoothing is carried out using topological persistence applied to the height function -δ(z). To explain this process intuitively, we borrow concepts from Morse theory [START_REF] Fomenko | Topological Modeling for visualization[END_REF]. The central theorem from Morse theory tells us that the topology of sublevel sets of a smooth function (of a Morse function, to be exact) changes upon passing critical points of the height function [START_REF] Milnor | Morse Theory[END_REF]. For example, imagine sweeping the landscape with a horizontal hyper-plane, in a bottom-up fashion: when the plane hits a local minimum, a c.c. of the sublevel set defined by the height of the minimum is created; this c.c. initially reduces to the minimum, and further grows to its catchment basin.

Since our focus is on c.c. of sublevel sets, two types of critical points matter. The first type consists of local minima, as each local minimum creates a c.c. of a sublevel set; we define the birth date of that c.c. as the height of the local minimum. The second type consists of index one saddles, since such a critical point joins the c.c. of (generically) two local minima, or creates a loop (a generator in the so-called order one homology). For a saddle joining two c.c., the height of that saddle defines the death date of the c.c. which was born with the highest local minimum of the pair-recall that a saddle generically links two local minima. Using the previous construction, one assigns a birth date and a death to connected components of sublevel sets, except that associated with the global minimum-which never dies. Representing these points in 2D, with the x and y axis respectively coding the birth and death dates, yields the so-called persistence diagram (PD) of connected components of sublevel sets [START_REF] Edelsbrunner | Computational topology: an introduction[END_REF]. In particular, the persistence of a local minimum is the difference between the death date and the birth date of its catchment basin, or equivalently, the height difference to the lowest saddle connecting this local minimum to a deeper local minimum.

Equipped with these notions, smoothing the landscape consists of retaining local minima whose persistence is above a user defined threshold ρ.

Algorithm. The significance threshold δ max and the persistence threshold ρ just discussed yield a partition of the PD into five regions defined by three lines (Fig. 2.3). A local minimum m of the landscape (and its catchment basin) actually gets qualified with respect to three criteria:

• Selected/rejected: m is selected provided that its birth date occurs before -δ max .

• Persistent/canceled: m is persistent if its persistence is ≥ ρ, a user defined threshold.

• Filtered/un-filtered: the catchment basin of m is filtered if the death date of m is larger than the threshold -δ max . Note that this step removes the samples whose discrepancy is not significant-with respect to the δ max threshold.

The possible combinations, illustrated on Fig. 2.3, are:

• m ∈ R 1 : rejected. Such a local minimum is rejected, since its discrepancy is less than δ max . No point of the catchment basin of m is found in a cluster reported.

• m ∈ R 2 : selected / canceled / un-filtered. Such a local minimum is selected, yet canceled by persistence. Because m dies before -δ max , all samples found in its catchment basin shall be part of a cluster reported.

• m ∈ R 3 : selected / canceled / filtered. A local minimum which is selected, yet canceled by persistence. However, because m dies after -δ max , only the portion of its catchment basin belonging to the sublevel set D ≤-δmax shall be found in a cluster reported.

• m ∈ R 4 : selected / persistent / un-filtered. Such a local minimum is selected, and is not canceled by persistence. Because m dies before -δ max , all the samples found in its catchment basin are found in a cluster reported.

• m ∈ R 5 : selected / persistent / filtered. This combination is similar to the previous case, except that samples whose discrepancy is less than δ max are discarded. Note that the cluster associated with the global minimum belongs to this region even though it is not found on the PD since the global minimum never dies.

From the previous discussion, one gets the following:

Observation. 1. The number of clusters reported is equal one plus the number of points found in the region R 5 of the persistence diagram. These clusters are contributed by a number of persistent local minima equal to one plus the number of points found in the region R 4 ∪ R 5 of the persistence diagram.

This observation also has the following practical implication: for a fixed number of clusters, increasing (resp. decreasing) the threshold δ max results in smaller (resp. larger) clusters, involving samples with high (resp. low) discrepancy.

Figure 2.2 Localizing the discrepancy: algorithm illustrated on a toy 1D example. (A, discrepancy estimation) We consider the height function defined by the discrepancy -δ(z), called the landscape for short. Practically, the estimate δ(z i ) is used at each sample z i (see text for details). (B, clusters as connected components) Focusing on discrepancies larger than a user defined threshold δ max , clusters are defined as connected components (c.c.) of the sublevel set of the landscape defined by δ max . On this example, four such c.c. are observed, two large ones and two small ones. Note that each cluster is charged to one local minimum of the function -δ(z). (C, smoothed clusters upon applying topological persistence) To get rid of small clusters, the landscape is smoothed using topological persistence (see text for details). In our case, this simplification yields two clusters. Note that the blue sample squeezed between C 1 and C 2 does not appear in the red cluster since its discrepancy is below δ max (and likewise for the red point with respect to C 3 and C 4 ).

(B)

-δ max (C) δ R d R d -δ max - δ - δ (A) R d z i δ(z i ) C 1 C 2 C 3 C 4 C 1 ∪ C 2 C 3 ∪ C 4 Figure 2.
3 Partition of the persistence diagram exploited to define clusters. The axes x and y correspond respectively to the birth and death dates. The partition of the domain y ≥ x is induced by three lines: (i) y = x + ρ which specifies the persistence threshold (ii,iii) x = -δ max , y = -δ max , with δ max the threshold on the significance of the discrepancy. See text for the specification of regions R 1 to R 5 .

y : Death

x : Birth y = x 0 0 -δ max -δ max R 1 R 3 R 5 R 4 R 2 y = x + ρ -1 -1

Data structures and algorithms

Having presented the clustering algorithm for a smooth landscape, we now transcribe the method to handle the sample points z i , each equipped with an estimate δ(z i ). Note that the following points need to be addressed:

• data structures to represent the landscape,

• identification of critical points from the sampling,

• computation of the catchment basins of the local minima,

• implementation of the persistence based simplification.

An elementary strategy. To represent the discretization of the landscape provided by the sample points, we build a nearest neighbor graph (NNG) connecting the samples, that is, we connect each sample to its k nearest neighbors-using for instance the Euclidean distance. Second, we lift this graph into R d+1 by adding a new coordinate to each sample, namely δ(z).

The NNG representation provides a simple way to study our estimated discrepancy [START_REF] Chazal | Persistence-based clustering in Riemannian manifolds[END_REF]. To do so, let the star of a sample point be the set of its neighbors in the NNG. A local minimum is a sample such that all samples in its star have a higher elevation. For a sample which is not a local minimum, say p, one may estimate the negative gradient of the height function studied using the sample q maximizing the slope from p to q. Note that starting from a sample which is not a local minimum, the process of iteratively following this pseudo-gradient ends up at a local minimum. This process is called the flow process, and all samples flowing to a local minimum define its catchment basin.

Finally, to identify index one saddles, we proceed as in [START_REF] Chazal | Persistence-based clustering in Riemannian manifolds[END_REF], processing the samples by increasing elevation δ(z). In doing so, a saddle is a sample flowing to a local minimum, but having a neighbor in its lower star flowing to a distinct local minimum. (If several samples meet this criterion for the same two local minima, that with the lowest elevation only is retained.) Note that this construction aims at capturing the merge event between smooth sublevel sets, as described in the previous section.

Upon detecting a saddle, the persistence of the highest local minimum involved is computed, and compared against the threshold ρ. In case this persistence is less than ρ, a merge between catchment basins is performed using Union-Find [START_REF] Tarjan | Data Structures and Network Algorithms[END_REF]. Finally, the inspection of samples halts upon reaching the first sample whose elevation is above -δ max -or equivalently discrepancy below δ max .

A refined strategy. A computational drawback of the previous method is to focus on samples, while the construction of clusters uses the two criteria δ max and ρ. The latter criterion being inherently coupled to the persistence diagram, i.e. to the critical points rather than the non critical samples (Fig. 2.3), we now present a more efficient method fully exploiting critical points and their catchment basins. This refined algorithm is especially interesting when the number of samples and the number of critical points differ by several orders of magnitude, which is usually the case. The algorithm is more involved as it requires notions from Morse homology [START_REF] Banyaga | Lectures on Morse Homology[END_REF]. In the sequel, we therefore only sketch its main steps, referring the reader to the supplemental section 2.7.1 for more details. The main steps are:

• Step1: Pre-processing. The local minima and index one saddles are identified as explained above.

The connections between local minima and saddle yielded by the aforementioned flow process are used to define a graph connecting critical points. This graph aims at reproducing the so-called Morse-Smale-Witten complex used in Morse homology [START_REF] Banyaga | Lectures on Morse Homology[END_REF]. The catchment basin of each local minima is also built and stored with this local minimum.

•

Step 2: Landscape simplification. Using the MSW complex, the landscape is simplified using topological persistence, using the recursive procedure explained in [START_REF] Cazals | Reconstructing 3D compact sets[END_REF]. This simplification manipulates the MSW complex rather than the samples-since the catchment basin is stored with the local minimum, whence its efficiency. The output consists of local minima with persistence above the threshold ρ, together with merged catchment basins.

•

Step 3: Extraction of clusters. Finally, the connected components of sublevel sets are extracted from the simplified landscape.

Remark 2.4. The previous description focuses on the identification of samples with significant discrepancy. In studying the height function δ(z) instead of -δ(z), one can instead identify clusters of samples of low discrepancy.

Combining discrepancy estimation and localization

Qualifying the clusters

We decompose the JSD by clusters of points that are defined by the method described in Section 2.3. Then, the contribution of a cluster C reads as:

JS C (f X f Y ) ≡ 1 n 0 + n 1 z∈(x (n 0 ) ∪y (n 1 ) )∩C δ(z). (2.14)
Combining the analysis of sections 2.2 and 2.3 yields the workflow of Fig. 2.4.

Plots

The previous analysis are best exploited using the following plots:

• Raw data embedding: For samples embedded in 2D or 3D space, a plot of the points with a color to indicate the label (blue: 0; red: 1). For samples embedded in a higher dimensional space, a 2D embedding of these samples obtained using multi-dimensional scaling (MDS). In any case, the goal of this plot is to intuitively visualize the distributions of the two populations.

• Discrepancy shaded data embedding -aka discrepancy plot: A plot similar to the raw data plot, except that each sample is color coded using a heat palette from fully transparent white to red across yellow, as a function of the value of the estimated discrepancy δ(z). That is, a point with δ(z) equal to zero (resp. one) is colored fully transparent white (resp. non transparent red).

Figure 2.4 Workflow of the whole method. In blue: the parameters.

L ∼ B( 1 /2) x (n0) y (n1) {z ji , l ji } i=1,...,mj Repeat from j = 1 to B k n ′ b Datasets:two point clouds Output: estimated discrepancy δ(z i ), z i ∈ x (n0) ∪ y (n1) JSD decomp. by clusters {C i } JS Ci (f X f Y )
Output: divergence by cluster Output:

{C i } s.t. ∪C i = x (n0) ∪ y (n1) ∩C i = ∅
Clustering with topological persistence:

k-Nearest neighbor graph -Persistence: ρ -Filtering: δ max estimate for δ(z)

• Persistence diagram: A plot showing for each minima a red cross with coordinates (x, y) corresponding to its birth and death dates respectively, while analyzing the landscape whose elevation is the negative estimated discrepancy.

• Clusters: A plot similar to the raw data plot, with one color per cluster. The points not belonging to any cluster are colored in gray.

• JSD decomposition plot: A 1D plot presenting a synthetic view without relying on the MDS embedding. The x coordinate represents the sample space and always ranges from 0 to 1. The total estimated JSD is represented by the area under the dashed line. And the maximum possible JSD which is always 1 is represented by the area under the continuous line. One bar is depicted for each cluster plus another one (the last one) for the points not belonging to any cluster. The area of each bar represents the contribution of the corresponding cluster to the total JSD and its color corresponds to the proportion of samples 0 in the cluster.

Implementation

Our method is implemented in the scope of the Structural Bioinformatics Library (SBL), available from http://sbl.inria.fr/. The SBL is a generic C++/Python library mixing low level algorithmic and data analysis classes, as well as applications targeting specific modeling problems in computational structural biology. Practically, our method is made available in http://sbl.inria.fr/applications/ > Data analysis > Density difference based clustering.

Practically, four python scripts are made available (see also the three steps in the workflow of Fig. 2.4):

• sbl-ddbc-step-1-discrepancy.py: the script estimating the discrepancy.

• sbl-ddbc-step-2-clustering.py: the script computing the clusters from the watershed transform.

• sbl-ddbc-step-3-cluster-plots.py: the script analyzing the clusters computed at the previous step.

• sbl-ddbc-disc-based-colored-embedding.py: an additional script computing an 2D embedding of the data processed, based on multi-dimensional scaling.

The user manual2 precisely describes the parameters of these scripts as well as the formats of input and output files, making the results presented in this paper reproducible.

Experiments

We present results on four datasets featuring various difficulties: medium embedding dimension (real images), low intrinsic dimension (crenels), varying intensity of discrepancy (mixture of Gaussians), and real data embedded in high dimension (handwritten digits). For the sake of convenience, we refer to the two datasets to be compared as the blue and the red datasets. The number of neighbors was set to k n = n 2/3 . Note that the Maximum Mean Discrepancy two-sample test of [GBR + 12] rejects the null hypothesis in all the cases for a significance level α lower than 1%.

Model: Gaussian mixture

Specification. The goal is to assess the ability of the method to spot regions of different intensity of discrepancy.

Two Gaussian mixture models where randomly generated using MixSim R package [START_REF] Melnykov | MixSim: An R package for simulating data to study performance of clustering algorithms[END_REF] (see also supplemental section 2.8.1). The distributions for X and Y consist in two mixtures of four different two-dimensional Gaussians with equal weight and with some degree of overlap. In this example, n 0 = n 1 = 2000.

Results.The results are presented on Fig. 2.8 and 2.9 corresponding respectively to thresholds δ max1 = 0.13 and δ max2 = 0.25. The discrepancy plot clearly shows regions of different intensities. The persistence diagram (built with k = 6) highlights four persistent minima plus the global one, calling for simplification. On the other hand, the critical values associated with the local minima of the elevation are quite scattered. Setting δ max1 = 0.13 selects two large clusters, while setting δ max1 = 0.25 selects four smaller clusters, of which those identified by 0 and 2 exhibit higher discrepancies. This illustrates how δ max trades-off the number of clusters, their discrepancy and their size.

In any case, note that the regions do not necessarily coincide with original components of the mixture, i.e., they are a result of the comparison only.

Model: crenels

Specification. The goal is to assess the ability of the method to spot local differences, and to cope with data of low intrinsic dimension (one) in a high dimensional space. We create two one dimensional datasets, which differ by two crenels.

More precisely, in this dataset, each point corresponds to a vector in R m×m encoding the pixels of a square grayscale image (0 = black, ≥ 1 = white) of size m × m. The samples are the result of rotating the grayscale image i (Fig. 2.5). Therefore, taking m = 11 yields a dataset of intrinsic dimension one embedded in dimension d = 121.

Each sample of the blue population is an instance of a RV X obtained by rotating i with a random uniformly distributed angle A X , that is:

X = rotate(i, A X ), A X ∼ U(s, t).
(2. [START_REF]fJWWW +[END_REF] where the function rotate applies a bilinear filter to smooth the result (details in [POSaH]).

For the red population, consider two Bernoulli random variables B 1 ∼ B(p 1 ) and B 2 ∼ B(p 2 ), and two uniform variables U 1 ∼ U (a, b) and U 1 ∼ U (c, d). Each sample of the red population is an instance of a RV Y defined as:

Y = rotate(i, A Y ), A Y = B 1 (B 2 U 1 + (1 -B 2 )U 2 ) + (1 -B 1 )A X .
(2.16)

Note that the rotation used to obtain Y comes from the uniform distribution U (s, t) with probability 1p 1 , and that the discrepancy between both distributions are high in the angle ranges [a, b] and [c, d] (i.e. the crenels) and low everywhere else on the support, since, loosely speaking, points that are added to the crenels are missing from the rest of the support. Practically, we used:

• n 0 = 2000, n 1 = 2000 • s = -15, t = 15 • a = -4, b = -2, c = 9, d = 10
• p 1 = 0.3, p 2 = 0.5 The other crenel corresponds to the global minimum which does not appear in the plot since its death date is infinite. The clusters yielded by the simplification and filtering steps correspond to the crenels. In the divergence decomposition plot, we observe these two crenels with high discrepancy produced by a high proportion of red points and also a non negligible total discrepancy given by the rest of the points, which has a higher proportion of blue points.

Model: mixture of handwritten digits

Specification. The goal is to assess the ability of the method to spot local differences, and to cope with real life high dimensional data. This dataset is based on the MNIST dataset [START_REF] Lecun | The MNIST database of handwritten digits[END_REF] that contains examples of handwritten digits (Fig. 2.6). Each digit is represented as a 28 × 28 gray-scale image, yielding a representation in dimension d = 784. For our experiment, we used the digits 3,6 and 8 and we built our populations by sampling with replacement from these three populations. The following table summarizes the number of samples taken from each digit set for each population:

Statistical image comparison

When processing real images, two-sample tests are of mild interest since the null hypothesis is likely to be rejected. However, the JSD decomposition is still of interest to quantify the differences on a statistical basis.

Consider a digital image whose pixels use C color channels. For example, C = 1 in the monochrome case and C = 3 in the RGB color case (where the components of each vector correspond to the red/blue/green intensities). A digital image is a r × c matrix of pixels, such that each pixel takes values in [0, 1] C . (For a pixel, a value of 0 (resp. 1) represents the minimum (resp. the maximum) intensity for the corresponding color channel.) We follow the construction of [START_REF] Song | Relative novelty detection[END_REF] to build our samples, that is, by taking b × b pixel blocks yielding (rb + 1)(cb + 1) samples, each being a vector of dimension Cb 2 . Then, a discrepancy estimate δ(z) is computed on each sample z and assigned to the pixel located in the upper left corner of the corresponding block. (NB: two bands of width b -1 on the right and bottom side of the image are not assessed.) Using the satellite color images of [START_REF] Song | Relative novelty detection[END_REF] shown in Figure 2.7 and using the same block size b = 2, we compute the estimated discrepancy δ(z i ) for each sample z i . The results are shown in Figure 2.12. The 2D MDS embedding is performed on a random subsample of size 4000 of the pooled dataset. The discrepancy plot hints at two clusters corresponding to colors found in the oval and the rectangular fields present in the second picture. These clusters stand out in the persistence diagram (obtained from a nearest neighbor graph with k = 10 nearest neighbors), the global minimum not being represented since its death date is infinite. We select them by setting δ max = 0.1.

In the JSD decomposition plot, we observe a small JSD (depicted by the dashed line) contributed by two relatively small clusters consisting of a strong majority of label 1. The local minima of the two clusters correspond to blocks in the second image as shown in Figure 2.13. This experiment shows how our method can be used as a statistical image comparison tool, allowing to spot differences between two images by comparing densities in the color space. In this example, we expect the two green fields to be responsible of the statistical discrepancy and this is what we find in our analysis.

Conclusion

This work proposes the first method to model the difference between datasets given as point clouds, for which there is evidence showing that they do not have the same underlying distribution. The method relies on a pointwise estimation of an integrand related to the Jensen-Shannon divergence (JSD), a symmetric version of the Kullback-Leibler divergence. An estimate of the JSD is obtained for each sample using a non-parametric regression method relying on k n nearest neighbors estimates. Topological persistence is then used to gather samples in groups associated with local maxima of the JSD. All in all, our method delivers groups of samples with significant contribution to the JSD, and associated with local maxima of the JSD.

On the theoretical side, several questions are of major interest. A first goal will be to characterize the clusters returned by our procedure, based upon assumptions on the distributions underlying the data. This problem is related to the robustness of a recent clustering method combining mode seeking and topological persistence [START_REF] Chazal | Persistence-based clustering in Riemannian manifolds[END_REF], since under suitable conditions, persistent modes of the density and those defining clusters have been shown to match. Coming up with a similar line of argumentation is more challenging in our case since two densities are involved, and the magnitude of the JSD and those of these densities are independent quantities. A second goal will consist of generalizing the method to data associated with a (non Euclidean) metric space.

On the practical side, we believe that our method goes well beyond statistical analysis based on two-sample tests, which essentially summarizes the information contained in all coordinates into a single boolean value (accept or reject the null hypothesis). It should therefore prove of interest wherever twosample tests are used.

Acknowledgments. The authors wish to thank Tom Dreyfus for implementing the landscape analysis method used in Step 2.

Supplemental: Algorithms

A refined strategy to compute clusters

In this section, we provide further details on the refined strategy used to compute clusters, see section 2.3. Recall that the elementary strategy merely maintains the tree of merge events between catchment basins, also known as the disconnectivity graph, using a Union-Find algorithm. The refined strategy uses the Morse-Smale-Witten complex instead, namely the graph connecting critical points (Fig. 2

.14).

Step 1: construction of the MSW complex. From the critical points identified in the elementary strategy, we build the bipartite graph whose vertex sets are the local minima (index 0 critical points) and the saddles (index 1 critical points) respectively (Fig. 2.14(D,E)). We also build the stable manifolds of local minima (i.e. their catchment basins), and assign them to the local minima. Using this MSW complex, we also compute the persistence diagram associated with local minima.

Step 2: Simplifying the landscape. The general procedure to recursively simplify a landscape using the MSW complex has already been presented in the context of non manifold shape reconstruction [START_REF] Cazals | Reconstructing 3D compact sets[END_REF]. This procedure, which handles cancellation of non persistent critical points whatever their indices, is used here only for local minima and index one saddles. In a nutshell, the cancellation of a pair of critical points (a, b) whose indices differ by one consists of rerouting the connections of a and b in the MSW complex, and of redistributing the stable manifold of a [START_REF] Cazals | Reconstructing 3D compact sets[END_REF]. Note in particular that each remaining local minimum is endowed with two types of samples: those from its own SM and those from SM inherited from canceled local minima.

Step 3: sublevel set extraction. The previous simplification yields a partition of the landscape into the SM of the persistent minima. We remove from these SM the samples whose discrepancy is less than δ max , a task carried out in two steps. First, the samples from its own SM are filtered out. Second, the samples of basins inherited from the simplification are also filtered out, provided that such a basin was born before -δ max . In particular, inherited basins born after -δ max are ruled out in constant time. The persistent local minima and their remaining samples, if any, form the clusters.

Remark 2.5. In the previous explanation, the NNG connecting samples has been taken for granted for the step one. The NNG construction must indeed be carried out first, as filtering out samples with low elevation may deplete the neighborhood of selected samples, and jeopardize the identification of critical points.

Figure 2.14 Landscape simplification using topological persistence: simple strategy using Union-Find, versus refined strategy using the Morse-Smale-Witten (MSW) complex. (A,B) Two landscapes, with critical points of indices 0 (disks) and 1 (squares). (C) The disconnectivity graph (DG) of both landscapes, namely the tree depicting the evolution of connected components of sublevel sets. Despite the differences between their MSW complexes, both landscapes share the same DG: upon passing the critical point e, the stable manifold of c merges with that born at a. (D,E) The MSW complexes of (A,B), respectively. In cancelling the pair of critical points (c, e), one does not know from the DG with which (a or b) the basin of c should be merged. But the required information is found in the MSW complex: on the landscape (A), c is merged with a; on the landscape (B), c is merged with b. • n 0 = n 1 = 2000

• d = 2
• Rationale: A simple and easy to visualize dataset with regions of different intensities of divergence

The first population is drawn according to a mixture of four spherical Gaussians with equal probability, i.e.,

X = 4 i=1 1 {i=U } N 0 [i]
where U is a uniform discrete RV which takes values in {1..4} and

N 0 [i] ∼ N (µ 0 [i], σ 0 [i])
The randomly generated parameters by MixSim R Package were:

µ 0 [1] = (0.9556715, 0.3617815) µ 0 [2] = (0.6207539, 0.8498296) µ 0 [3] = (0.3077166, 0.2886823) µ 0 [4] = (0.1496965, 0.9773699)
distributions is effective in increasing power, as our tests outperform state-of-the-art ones on selected challenging datasets.

We note in passing that our contribution bears two main differences with Wald's sequential test [START_REF] Wald | Sequential tests of statistical hypotheses[END_REF]. First of all, this classical procedure works only for simple alternative hypotheses, since the probability of type II error must be kept under control. Although some extensions have been proposed (e.g., [START_REF] Santiago-Mozos | On the uncertainty in sequential hypothesis testing[END_REF]), they are not applicable to the nonparametric case. Another important difference is that in Wald's procedure the stopping rule must be fixed in advance in order to obtain a valid p-value for the whole procedure, while in our procedure optional stop is allowed.

The previous contributions are articulated as follows. Lemma 3.1 connects classical two-sample testing to probabilistic classifiers. Theorem 3.1 and Lemma 3.2 turn any sequential predictor (and, thus, any non-anticipating probabilistic classifier) into a sequential two-sample test allowing optional stopping. Theorem 3.2 gives sufficient properties on probabilistic classifiers to guarantee λ-consistency (a weaker sense of consistency) of the resulting two-sample test. Theorem 3.3 shows how such conditions can be obtained from strongly pointwise regressors.

Since regressors are in general parameterized and may be optimized, in Section 3.4, we use Bayesian mixtures of regressors, promoting the regressors with better past performance. Additionally, we suggest to use switch distributions to attenuate the power loss that is expected due to sequential learning.

Two-sample test based on sequential prediction

Problem statement

We state the problem as follows:

Problem. 3.1. Given a set of samples x 1 . . . x n0 and y 1 . . . y n1 whose corresponding random variables X i and Y i are i.i.d. with densities f X and f Y respectively, select one of the following hypotheses

H 0 : f X = f Y a.e. H 1 : ¬H 0 . (3.1)
In order to assess the strength of the evidence against H 0 , a random variable p is used, which is called a valid p-value (see, e.g., [CB01, Def. 8.3.26]) if 0 ≤ p ≤ 1 and

P H0 (p ≤ α) ≤ α, ∀α ∈ [0, 1]. (3.2)
Then, the lower is p the stronger is the evidence against H 0 . It is also possible, to set a threshold called significance level α, so that H 0 is rejected when p ≤ α.

Then, two types of errors must be considered. One faces a type I error when H 0 is rejected while it is actually true. One faces a type II error when H 0 is not rejected while it is actually false. The probability of Type I error is controlled by design and is upper bounded by α. Then, one usually considers the power of the test for a significance level α, which is P H1 (p ≤ α). The test is termed consistent for a given level

α when P H1 (p ≤ α) n0+n1→∞ -------→ 1.

Random labels framework

A sequential probability predictor (or predictor for short) processes sequentially input symbols belonging to some alphabet A. Before observing the next symbol in the sequence, it predicts it by estimating the probability of observing each symbol of the alphabet. Then, it observes the symbol and some loss is incurred depending on the estimated probability of the current symbol. Subsequently, it refines its model in order to better predict future symbols. The predictor can also be allowed to observe side-information to make better predictions.

Intuitively, if we shuffle the samples from both populations, and yet manage to predict the population each sample belongs to, then it is natural to think that there is some difference in the features, so that H 0 should be rejected.

In order to do this shuffling, we use a random device receiving samples from each of the populations(Fig. 3.1). The output corresponds to the input X with probability θ 0 or to Y with probability 1θ 0 , where 0 < θ 0 < 1 is a parameter that must be set. Formally, considering the alphabet A = {0, 1}, we define the following pair of random variables: (L, Z) = (0, X) with probability θ 0 (1, Y ) with probability 1θ 0 .

In a classical two-sample test setting, the inputs of the random device uniformly draws from each of the given finite populations until the selected input has no more samples available, in which case the paired sequence generated (L, Z)

N ends. Notice, that it can happen that N < n 0 + n 1 . In order to minimize the expected number of unused samples N -(n 0 + n 1 ), one should set θ 0 = n 0 /(n 0 + n 1 ).

Remark 3.1. If we are, for example, in a streaming setting where the inputs of the random multiplexer are plugged on streams with rates r 0 and r 1 (in symbols per time unit) it is natural to set θ 0 = r 0 /(r 0 + r 1 ) in order to match them. In this case, the sequence never ends.

Notations and problem reformulation

Our framework is based on considering two random variables Z (positions) and L (labels) representing, respectively, the pooled original samples, and the two populations these samples belong to. The following notations are used to describe the probabilistic properties of these random variables. The unconditional and conditional label probabilities are denoted P θ0 (l) and P θ(z) (l|z). One has: P θ0 (l) ≡ P (L = l) , P θ(z) (l|z) ≡ P (L = l|Z = z) (3.3) with θ 0 ≡ P (L = 0) and θ(z) ≡ P (L = 0|Z = z). The joint density and the joint density assuming independence are respectively denoted f θ(z) (z, l) and f θ0 (z, l). The mixture density for position is denoted

f (z) = l f θ(z) (z, l).
The entropy of random variables is denoted H (•), while the entropy of L conditioned on Z is denoted H (L|Z); finally, the mutual information between Z and L is denoted I (Z; L).

Remark 3.2. With the previous notations, one has f θ(z) (z, l) ≡ P θ(z) (l|z) f (z), and f θ0 (z, l) ≡ P θ0 (l) f (z).

(3.4)

In the setting of random labels, let us consider the following two-sample problem:

Problem. 3.2. Given a sequence of samples (l, z) n whose corresponding random variables (L i , Z i ) are i.i.d. with joint density f θ(z) (•, •), select one of the following hypotheses

H 0 : f θ(z) (z, l) = f θ0 (z, l) a.e. ∀l ∈ {0, 1} H 1 : ¬H 0 . (3.5)
The following lemma is a simple consequence of Bayes' formula applied to the joint densities.

Lemma. 3.1. The null hypotheses of Problems 3.1 and 3.2 are equivalent.

Proof. Consider the conditional density for position

f l (z) ≡ f X (z) if l = 0 f Y (z) if l = 1 .
Using Bayes' formula for the joint densities yields f θ(z) (z, l) = f θ0 (z, l) a.e. ∀l ∈ {0, 1} ⇔ f l (z)P θ0 (l) = f (z)P θ0 (l) a.e. ∀l ∈ {0, 1}

⇔ f X = f Y a.e..

Robust sequential p-value

Using the statement of Problem 3.2, we phrase our two-sample problem as a model selection problem and use a likelihood ratio test to obtain a p-value. The models we consider are distributions on length-n sequences, which can be obtained from sequential probability predictors. This approach has the advantage of providing a hypothesis test in which the sample size need not be fixed in advance as classical Neyman-Pearson does (see, e.g., [START_REF] Van Der Pas | Almost the best of three worlds: Risk, consistency and optional stopping for the switch criterion in single parameter model selection[END_REF][START_REF] Shafer | Test martingales, Bayes factors and p-values[END_REF]). More formally, Theorem. 3.1. Given some arbitrary distribution Q on infinite sequences l 1 , l 2 , . . . (Q (l n ) denoting its marginal distribution on the first n outcomes),1 a test that rejects H 0 at any index n when the likelihood ratio

P θ0 (l n ) Q (l n ) ≤ α (3.6)
has a Type I error probability less or equal than α for problem 3.2, i.e., P θ0 ∃n :

P θ0 (L n ) Q (L n ) ≤ α ≤ α. (3.7)
The likelihood ratio

P θ 0 (l n ) Q(l n ) is called a robust p-value.
Proof. We consider the i.i.d. sequence L n , Z n and the class of models

f θ(z) (•, •) to which f θ0 (•, •) belongs. Let us define p 1 (l n , z n ) ≡ Q (l n ) f (z n ), which is a function from R d•n × A n to R + 0 . Note that p 1 (l n , z n ) f θ0 (l n , z n ) = Q (l n ) f (z n ) P θ0 (l n ) f (z n ) = Q (l n ) P θ0 (l n ) .
Then, one has

E θ0 p 1 (L n , Z n ) f θ0 (L n , Z n ) = E θ0 Q (L n ) P θ0 (L n ) = L n ∈A n P θ0 (L n ) Q (L n ) P θ0 (L n ) = 1
where the last inequality stems from Q being a distribution. Equation 3.7 follows from Theorem 3.5(Appendix 3.9).

The following lemma shows that any sequential predictor assigning a probability to l i using past label data l i-1 , and possibly all the available position data z ∞ , complies with Thm. 3.1: Lemma. 3.2. Given a sequential probability predictor Q l i |l i-1 (i.e. l∈A Q (l|•) = 1), one can build the following function of sequences

l n ∈ A n Q (l n ) ≡ n i=1 Q l i |l i-1 , (3.8)
which is a distribution.

Proof. Let us prove by induction that it is a distribution using the distributive law to obtain the following expression

l n ∈A n Q (l n ) ≡ l1∈A Q l 1 |l 0 l2∈A Q l 2 |l 1 • • • ln∈A Q l n |l n-1
where l 0 is the empty sequence. The base case corresponds to right-most sum and, since Q is a sequential probability predictor, we have

ln∈A Q l n |l n-1 = 1.
Observing that each sum is a convex combination of an expression that is equal to one proves the claim.

Consistency via λ-pointwise universal distributions (λ-PUD)

We now define the requirements imposed on the distributions to obtain consistent tests, in a weaker sense that we call λ-consistency. To this end, we consider distributions depending on the position sequence z n , whence the notation Q (l n |z n ).

Definition. 3.1. Given 0 < λ ≤ 1, a distribution Q is λ-pointwise universal (λ-PUD) if -lim 1 n log Q (L n |Z n ) ≤ H (L|Z) -log λ a.s.
for any f θ(z) (•, •) generating the samples.

Before introducing the notion of λ-consistency, we prove the following lemma.

Lemma. 3.3. Under the alternative hypothesis, I (Z; L) > 0.

Proof. By definition,

I (Z; L) ≡ - z∈R d l∈A f θ(z) (z, l) log f θ0 (z, l) f θ(z) (z, l) .
Recall that for a strictly convex function h and a random variable X, Jensen's lemma states that (see, e.g., [CT06, Thm. 2.6.2])

E [h(X)] ≥ h(E [X]) (3.9)
and the equality in 3.9 implies X = E [X] with probability 1. We apply it to X = (Z,L) and to the strictly convex function h(X) =log(X). It remains to show that the random variable X is non-constant with non-null probability.

f θ 0 (Z,L) f θ(z)
Under H 1 , there exists a set S of non-null probability such that ∀z ∈ S ∀l ∈ A :

f θ0 (z, l) = f θ(z) (z, l) ⇔ P θ0 (l) f (z) = P θ(z) (l|z) f (z) (3.10) ⇒ P θ0 (l) = P θ(z) (l|z) (3.11)
If the ratio were constant we should have

f θ0 (z, 0) f θ(z) (z, 0) = f θ0 (z, 1) f θ(z) (z, 1) ⇔ P θ0 (0) P θ(z) (0|z) = P θ0 (1) P θ(z) (1|z) = 1 -P θ0 (0) 1 -P θ(z) (0|z)
⇔ P θ0 (0) = P θ(z) (0|z) (3.12) which contradicts Ineq. (3.11).

The following theorem introduces the λ-consistency property -from which one recovers the usual notion of consistency for λ = 1, and shows that this property is obtained from λ-PUD.

Theorem. 3.2 (λ-consistency). By Lemma 3.3, under the alternative hypothesis, I (Z; L) = ε > 0. Consider a λ-PUD Q such that λ > 2 -ε . Then, the test described in theorem 3.1 using that λ-PUD is consistent.

Proof. The probability of rejecting H 0 is P Z,L ∃n :

P θ0 (L n ) Q (L n |Z n ) ≤ α (3.13) = P Z,L ∃n : - 1 n log P θ0 (L n ) + 1 n log Q (L n |Z n ) ≥ - log α n . (3.14)
Since Q is λ-PUD and, by the Asymptotic Equipartition Property (see, e.g., [START_REF] Cover | Elements of Information Theory[END_REF]), we have lim -

1 n log P θ0 (L n ) + 1 n log Q (L n |Z n ) ≥ H (L) -H (L|Z) + log λ a.s. = I (Z; L) + log λ > 0
where the last inequality stems from I (Z; L) >log λ, which is a direct consequence of the assumptions. Therefore, there exist δ > 0 and n ′ (δ) and such that ∀n ≥ n ′ (δ), the left-hand side of the inequality in Eq. (3.14) is greater than δ almost surely. For the right-hand side of the inequality, there exists n ′′ (δ) such that ∀n ≥ n ′′ (δ), -log α n < δ. Thus, for any n ≥ max(n ′ , n ′′ ), the inequality holds, and thus the probability of rejecting is 1.

λ-Pointwise universal distributions via strongly pointwise consistent regressors

In this section, we construct λ-PUD thus yielding λ-consistent tests. The construction uses sequential predictors based upon strongly pointwise consistent regressors. These sequential predictors define a distribution via lemma 3.2, and therefore a robust sequential p-value, and yield a λ-PUD.

Sequential probability estimation from nonparametric regression. We first build a sequential predictor using nonparametric regression (see, e.g., [START_REF] Györfi | A distribution-free theory of nonparametric regression[END_REF]). Given a random vector (Z, R), where Z ∈ R d and the response variable R ∈ R, the regression function is defined as

m(z) = E [R|Z = z] . (3.15)
To obtain a sequential predictor, we consider the mapping R = L. This yields m(z) = P (L = 1|Z = z) .

(3.16)

Let m n (z) be an estimate of m(z) based on n i.i.d. realizations of (Z, R). Given some sequence of regression function estimators {m n } such that 0 ≤ m n (z) ≤ 1 ∀n, z, let us define the following sequential predictor:

Pi

l i |l i-1 , z i ≡ m i-1 (z i )1 li=1 + (1 -m i-1 (z i ))1 li=0 . (3.17)
Notice that in this case, Pi predicts l i based on l i-1 and z i and ignores future z samples, and thus is compliant with lemma 3.2. We are interested in the following sense of consistency since it allows to build λ-PUD. We call an s.p.c. sequence of regression estimates an s.p.c regressor. In [GK02, Sec. 25.6], some s.p.c. regressors are presented. For example, regressors based on partitioning, kernel and nearest neighbors are s.p.c under certain conditions for their parameters, when the absolute value |R| < M , for some M . Note that is our case since R ∈ {0, 1}. In our experiments, we use nearest neighbors regressors, which are described in Appendix 3.8.

Type I error and robust sequential p-value. We follow the construction of [START_REF] Algoet | Universal schemes for prediction, gambling and portfolio selection[END_REF][Sec. 2] in order to build a λ-PUD. For any distribution Q (l|z) and 0 ≤ λ < 1, let us define

Q λ (l|z) ≡ (1 -λ) 1 |A| + λQ (l|z) . ( 3
i ((L, Z) ∞ -∞ ) ≡ -log P λ -i L 1 |Z 1 , (L, Z) 0 -i+1 (3.24)
where

P-i L 1 |Z 1 , (L, Z) 0 -i+1 ≡ m -i (Z 1 )1 li=1 + (1 -m -i (Z 1 ))1 li=0 (3.25)
where m -i (z) is an estimate of m(z) based on (L, Z)

0 -i+1 . Then the claim is equivalent to lim n 1 n n-1 i=0 ĝλ i (T i (L, Z) ∞ -∞ ) = E -log P λ θ(z) (
P P-i L 1 |Z 1 , (L, Z) 0 -i+1 i→∞ ---→ P θ(z) (L 1 |Z 1 ) = 1 (3.28)
Remark 3.3. Compare this to the more rigid approach of splitting the sample in two and optimizing on one part and using the optimized value to test on the other (as for TreeRank [START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF] or optimized kernel MMD [GSS + 12, ZGB13]). Note also that our process keeps learning beyond the index of the sample corresponding to the training set size.

Remark 3.4. Another interest of mixtures is to use non universal distributions that could be tailored to some more specific situations on which they may exhibit a better power than the universal ones. For example, one could use a linear kernel support vector machine2 which would be expected to be more sensitive on distributions that differ by a shift.

If a better finite size performance is sought, other predictors could be plugged in our framework, which could be better for different kind of differences between distributions.

Model switching. The mixture of lemma 3.4 also allows model switching to avoid the catch-up phenomenon (see [START_REF] Van Erven | Catching up faster by switching sooner: a predictive approach to adaptive estimation with an application to the aic-bic dilemma[END_REF]). That is, even if we are under H 1 , when few samples are available it can be better to predict using P θ0 and then switch to Pn when more samples are available.

Let us define the following distribution that switches before time s from P θ0 to any predictor Pn :

P Sw {P θ 0 , Pn} (l n |s, z n ) ≡ P θ0 l s-1 Pn l n s |l s-1 , z n (3.34)
where l 0 is the empty sequence and P θ0 l 0 ≡ 1. Given a prior π(s) on the switching time s, we define the following switch distribution P Sw(π)

{P θ 0 , Pn} (l n |z n ) ≡ s π(s) P Sw {P θ 0 , Pn} (l n |s, z n ) . (3.35)
The following lemma follows from lemma 3.4.

Lemma. 3.5. Given a λ-PUD Pn , the switch distribution defined by Eq. (3.35) with a prior π such that log π(0

) = o(n) is λ-PUD.
A first option for the switch time prior is the horizon-free prior defined in [START_REF] Van Erven | Catching up faster by switching sooner: a predictive approach to adaptive estimation with an application to the aic-bic dilemma[END_REF] π S (s) ≡ 1 s(s + 1) .

(3.36)

Another option, when the horizon n is known, is the uniform prior

π U (s) ≡ 1 n .
(3.37) Note that Eq. (3.34) can be interpreted as a training phase (until s) followed by a test phase that keeps learning. Eq. (3.35) replaces the choice of s by a mixture of all possible values, weighted by a prior.

Implementation and complexity

Here, let n be the number of samples effectively processed until the process is stopped.

KNN The KNN predictor was implemented using exact brute force FLANN's implementation [START_REF] Muja | Flann, fast library for approximate nearest neighbors[END_REF] whose time complexity is linear in the number of samples seen so far. When a new sample z i is seen, its nearest neighbors are searched to compute the labels' probabilities as described in Appendix 3.8, then the pair (l i , z i ) is added to the training set. At the beginning, when the training set is empty, it predicts using P θ0 . Therefore, the time complexity is quadratic in n and the space complexity is linear in n.

Since the algorithms we describe next are independent on whether the predictors depend on z n or not, we omit z n from notation.

BMA The sequential formula for BMA is, by Bayes rule,

BMA {Q k } l i |l i-1 = k Q k l i |l i-1 µ(k|l i-1 ) (3.38)
where the weights correspond to the predictors' posterior probability, i.e.,

µ(k|l

i-1 ) = µ k Q k l i-1 k µ k Q k (l i-1 )
.

(3.39) Therefore, it suffices to sequentially maintain the total probability of each predictor using a log representation for probabilities to avoid numerical problems. More precisely, the following identities allow to multiply and to add two given probabilities x and y in log representation: log(xy) = log(x) + log(y) log(x + y) = log(x) + log(1 + exp(log(y)log(x))).

Practically, for the sum, if the difference in magnitude between x and y is too large3 , the result is simply the larger number.

Therefore, the BMA algorithm is linear in time and uses constant space.

Switching In order to perform sequential computation of the switch distribution we use the algorithm described in [E + 10] which is a special case of the Forward algorithm of [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]. In a nutshell, it is like BMA except that weights are updated in a different manner. Algorithm 1 describes it for our case. It takes as argument the prior π in the conditional form µ(s) ≡ π(T = s|T ≥ s), T representing the switching time, which verifies

π(s) = π(T = s) = π(T = s|T ≥ s)π(T ≥ s) = π(T = s|T ≥ s) s-1 i=1 (1 -π(T = i|T ≥ i)).
The standard and uniform priors presented in Section 3.4 yield respectively

µ S (s) = 1 s+1 µ U (s) = 1 n-s+1
, s ≥ 1.

As BMA, the switching algorithm is linear in time and uses constant space.

Experiments

Instantiations

Our constructions allow defining two-sample tests from λ-PUD in general (section 3.2), and s.p.c. regressors in particular (section 3.3). More precisely, we define the following λ-consistent two-sample tests:

1. KNN p ≡ P λ n obtained via KNN with k n = ⌈n p ⌉ and λ = 0.9999 2. BMA≡ BMA {KNNp} , p ∈ {.3, .5, .7, .9} with uniform prior µ

SW

π U ≡ P Sw(π U ) {P θ 0 ,BMA} 4. SW π S ≡ P Sw(π S )
{P θ 0 ,BMA} For tests 3. and 4., we consider two versions: the first one computes the likelihood ratio on the full sequence generated by the random device discussed in section 3.2.2 (whence the letter F for Full); the second one exploits the ∃ quantifier of Eq. (3.7), i.e., stops at the first index such that the likelihood ratio is below α (whence the letters OS for Optional Stop). Note that for tests 1. and 2., we only report results for the full version.

Algorithm 1 Sequential Switch Distribution Computation {Input: a sequence l n , a conditional prior µ(s) ≡ π(T = s|T ≥ s) } {Output:

P Sw(π) {P θ 0 , Pn} (l n ) } p ← 1 {Initialize weights:} {no switch before time 1} w θ0 ← 1 -µ(1) {switch before time 1} w θ(z) ← µ(1) for i = 1 . . . n do {Predict:} p ← p • w θ0 P θ0 l i |l i-1 + w θ(z) Pi-1 l i |l i-1 {Update posterior weights:} {no switch before time i + 1} w θ0 ← w θ0 P θ0 l i |l i-1 (1 -µ(i + 1)) {switch at time i or already switched before } w θ(z) ← w θ0 P θ0 l i |l i-1 µ(i + 1) + w θ(z) Pi-1 l i |l i-1 {normalize weights} w θ0 ← w θ 0/(wθ 0 +w θ(z) )
w θ(z) ← w θ(z)/(w θ 0 +w θ(z) ) return p

Contenders

We compare their performance against the following methods:

• MMD b : the Maximum Mean Discrepancy two-sample test presented in [GBR + 12] using the bootstrap approach to compute the rejection region. We used the MATLAB code available at http: //www.gatsby.ucl.ac.uk/ ~gretton/mmd/mmd.htm. 500 shuffles were used for the bootstrap.

• TRank: the AUC optimization based two-sample test presented in [START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF]. We used the R implementation available from CRAN Archive at http://cran.r-project.org/src/contrib/Archive/ TreeRank with the default parameters. Since we test our methods on the same distributions and sample sizes as in [START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF] we also include the performance stated in that paper, which we denote as TRank ref .

In those cases, we ran TRank with same split value as TRank ref .

• Kernel optimized MMD ( OPT, MxRt, MxMMD, L2, xvalc and med): these are the methods proposed in [GSS + 12]. They are all based on splitting the data in train/test sets, performing some kernel optimization on the train set then testing on the test set. The kernels considered are Gaussian ones with width σ ∈ {2 i } i=-15...10 .

Results

Setup. We set α = 0.05. Experiments were run with n 0 = n 1 since the implementations of selected contenders (MMD b ) require this balance. Therefore, we set θ 0 = 1/2 in order to match these proportions (see discussion in section 3.2.2). Following [CDV09, GSS + 12], our experiments assess the type I and type II errors of our tests and their contenders, under various conditions: Gaussian data in low and medium dimension, mixture of Gaussian data aiming at assessing the incidence of scale.

Gaussian data in dimension d = 4. Table 3.1 presents the results on data generated using the specification of [CDV09, Fig. 1], corresponding to different 4-dimensional Gaussians: A remarkable feature of this test is its semi-supervised nature, meaning that, in an initial step, it uses an unlabeled set for the spatial partitioning. Then, its uses a sequence of labeled samples as in Chapter 3. The unlabeled set can be, for example, a subset of the labeled one. We call our construction of Q Random Projection Context Trees since it uses the Random Projection Trees of [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF] to define a spatial partitioning whose cells are used as context to condition the probability of the labels.

The construction is better described in three steps corresponding to the following sections. The first one consists in a discretization of the R d using a spatial partition where in each cell a distribution is estimated. The second step consists in considering a hierarchical set of random partitions (obtained from the unlabeled set) each providing a discretized conditional distribution from the first step, and then mixing them using a weighting scheme. Practically, the hierarchical partition is based on a binary tree. The success of this construction depends on the ability of the random partitions to reveal cells with a proportion of label 0 different from the unconditional probability θ 0 . In the third step, we boost the probability of finding such cells by considering and ensemble of trees.

Complexity-wise, processing a new labeled sample has logarithmic time complexity w.r.t. the unlabeled set size, and negligible memory footprint, which makes it suitable for streaming data. Moreover, this construction exhibits automatic adaption to local scale and does not require parameters to be set.

Finally, we present experiments stressing the power of our test and its ability to cope with different scales.

Discretizing the sample space

θ Ci ≡ P (L = 0|Z ∈ C i ) = z∈Ci θ(z)f Z (z)dz z∈Ci f Z (z)dz , (4.1)
from which one defines the following vector characterizing the partition:

θ P = [θ C1 , . . . , θ C |P | ]. (4.2)
Let C P (z) denote the cell of P in which a point z falls, then we denote

P θ P (l|z) ≡ P (L = l|Z ∈ C P (z)) , (4.3)
and the associated joint distribution P θ P (l n |z n ) . (4.4)

Consistency via pointwise universal distributions on a revealing partition

In order to be able to detect a change, the partition P must reveal differences between the unconditional and the conditional distributions. Therefore, we define the following.

Definition. 4.1. A revealing cell C is a cell such that θ C = θ 0 . A revealing partition P is such that it contains at least one revealing cell.

Lemma. 4.1. Let P be a revealing partition, then

H θ P (L|Z) ≡ - C∈P P (Z ∈ C) l∈{0,1} P θ P (l|z) log P θ P (l|z) < H (L) (4.5)
where H (L) denotes the entropy of L.

Proof. Since conditioning reduces entropy on average, we have that H θ P (L|Z) ≤ H (L) with equality if and only if variables L and C P (Z) are independent, i.e., ∀l, c, P (L = l|C P (Z) = c) = P (L = l) which is not the case since P is a revealing partition and thus there exists a cell c such that the equality does not hold.

On the way to defining a distribution on infinite sequences (thus compatible with theorem 3.1) and also yielding a consistent test, we define: Definition. 4.2. A distribution Q is pointwise universal with respect to a partition P (denoted as

P -PUD) if lim n→∞ - 1 n log Q (L n |Z n ) ≤ H θ P (L|Z) a.s.. (4.6)
The following theorem guarantees consistency when P -PUD are used:

Theorem. 4.1 (Consistency). The test described in Theorem 3.1 with Q being a distribution on infinite sequences and a P -PUD is consistent.

Proof. The probability of rejecting H 0 is P Z,L ∃n :

P θ0 (L n ) Q (L n |Z n ) ≤ α (4.7) = P Z,L ∃n : - 1 n log P θ0 (L n ) + 1 n log Q (L n |Z n ) ≥ - log α n . (4.8)
Since Q is a P -PUD and, by the Asymptotic Equipartition Property (see, e.g., [START_REF] Cover | Elements of Information Theory[END_REF]), we have lim -

1 n log P θ0 (L n ) + 1 n log Q (L n |Z n ) = H (L) -H θ P (L|Z) a.s..
And by lemma 4.1, the last expression is strictly positive. Therefore, there exist δ > 0 and n ′ (δ) and such that ∀n ≥ n ′ (δ), the left-hand size of the inequality in Eq. (4.8) is greater than δ almost surely. For the right-hand side of the inequality, there exists n ′′ (δ) such that ∀n ≥ n ′′ (δ), -log α n < δ. Thus, for any n ≥ max(n ′ , n ′′ ), the inequality holds, and thus the probability of rejecting is 1. The corresponding full tree. To each node is associated a region; internal nodes also come with a splitting rule to define the regions of the two children. Note that a tree T generates a set of partitions. Identifying a node with its associated region, these partitions are P T = {η 1 , {η 2 , n 3 }, {η 2 , η 4 , η 5 }}. (C) The partition associated with the dashed cut consists of the two grayed cells.

(A)

η 2 η 3 (C) Cut η 1 η 2 η 3 η 4 η 5 (B)
Corollary. 4.1. Let P be a revealing partition. Then the test of Theorem 3.1 with the distribution P Jeffreys(P ) n is consistent.

Proof. Follows from P being a revealing partition, Lemma 4.2 and Theorem 4.1.

Remark 4.1 (Sequential computation of Jeffreys' Mixture). Let n i l (C (z)) denote the number of times the symbol l occurred in the sequence l i restricted to samples that fell in the cell containing z and let n i (C (z)) denote the total number of symbols of the sequence l i restricted to samples that fell in the cell containing z. Then, Jeffreys' mixture can be computed sequentially using Kritchevsy-Trofimov (KT) estimator [KT81], i.e,

P Jeffreys(P ) n (l n |z n ) = n i=1 n i-1 li (C (z i )) + 1 2 n i-1 (C (z i )) + |A| 2 . (4.21)
Therefore, the sequential computation only requires maintaining the counters n i l (C) for each cell C.

4.3 Using a hierarchical ensemble of partitions associated with a tree

Hierarchical partitioning based on random projection trees

Practically, we consider a class of partitions associated to full binary space partition trees.1 In such a tree, each internal node has two children, be they internal nodes or leaves, obtained by recursively applying a splitting rule (Fig. 4.1(A,B)). To define such partitions, let a tree cut be a simplification yielding a full binary tree. A partition consists of the cells associated with the leaves of a tree cut (Fig. 4.1(C)).

If we consider a perfect depth-K binary tree,2 there exists a doubly-exponential number of N K ≈ 1.5 2 K partitions embedded within (see, e.g., [START_REF] Kozat | Universal portfolios via context trees[END_REF]). Random Projection (RP) Trees create a data adaptive hierarchical partition P of R d by using projections and splittings along random lines [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF]. We use the RPTree-Max version, which splits the data at the median plus a random jitter. This construction is suitable for high-dimensional data that has lower intrinsic dimension. More precisely, given a cell C ∈ P of doubling dimension d ′ [START_REF] Assouad | Plongements lipschitziens dans R n[END_REF], with probability at least 1 /2, only O (d ′ log d ′ ) levels are required to halve the diameter of the data contained in C.

In our implementation, the trees are based on a set of N unlabeled samples and we split until there is only one point per cell. Although the unlabeled set used for tree construction can be an arbitrary one, in the experiments, we use all the positions occurring in the sequence to be considered. Notice that since no label information is used, this does not compromise the conditions of Theorem 3.1 and thus p-value validity. In the case Ex3, we observe a good performance of CTS2 that is comparable to MMDb, especially, when the sample size increases (Ex3b).

With respect to the higher dimensional cases of Table 4.2, we observe that CT*2 is prone to the curse of dimensionality as our KNN based methods, which is expected since they are based on spatial partitions.

In the subtle cases B2 and B3 of Table 4.3, we observe the ability of our method to automatically adapt to the right scale without requiring any scale related parameter to be set as it is the case for our KNN based mixtures and the Kernel optimized MMD versions.

Conclusions

This work introduces a sequential nonparametric two-sample test, based on sequential prediction of labels defining the two populations, with the novel characteristic of being semi-supervised.

We introduce the notion of revealing partition and use it to establish the consistency of the two-sample test. Our test is shown to be consistent when prediction is carried out using Context Tree Weighting and Context Tree Switching over a set of partitions containing a revealing one. These constructions do not require any parameter to choose in order to define the scale at which differences occur. Implementationwise, the random spatial partitions are constructed using random projection trees. The time complexity of processing a new observation only depends on the number of unlabeled samples used to build the trees initially, which therefore remains constant. This property together with a negligible footprint per observation makes the test suitable for streaming data. For future work, the problem of analyzing the probability of finding a revealing partition as a function of the number of trees used in the Bayesian mixtures is open. Also, accommodating data in a metric (non Euclidean) space, using e.g. metric trees, would widen the application spectrum of the method.

Chapter 5 Conclusion

Conclusions

This thesis was motivated by the idea of giving new insights and proposing new methods for two-sample comparison.

Although the natural framework was the one of statistical hypothesis testing, we did not want to propose "yet another two-sample test" since there are many of them (none uniformly better than the others) and it is relatively easy to invent new ones.

An important concern was the practical relevance of nonparametric multivariate tests. If a consistent test is used, any small difference is detected with probability arbitrarily close to 1, as long as enough samples are given. If data comes from random sources that are not identical in nature, it is natural to believe that slightly small differences actually exist and the test turns out to answer the question: "do we have enough samples so as to show that the sets come from different distributions?" which might not be of so much interest.

This straightforwardly led us to the problem of measuring an effect size which was an almost unexplored field for nonparametric multivariate two-sample tests. Using a novel combination of statistical learning, Morse theory and topological persistence, we made a progress in this direction by proposing a simplified decomposition of the Jensen-Shannon divergence, allowing the user to a have synthetic view, in the form of a bar plot, of the structure of the difference, even in high dimensions.

As we explored the literature, we found out how related dissimilarity measures, classification and two-sample testing are. This motivated the idea of building a two-sample test based on sequential soft classification, giving birth to the first nonparametric multivariate two-sample test that processes samples in a sequential manner and enjoys optional stopping, in contrast to classical two-sample testing that uses batch data.

In order to build sequential soft classifiers, we took two different directions. In the first one, we established an interesting connection between strongly pointwise consistent regressors and consistent tests. In the second one, we combined two powerful constructions, namely random projection trees and universal distributions based on context trees, to build sequential classifiers that exhibit automatic adaption to scale and yield powerful two-sample tests.

Discussion and Future Work

In many practical experiments, we observed that differences can often be spotted and more easily interpreted by univariate tests applied to each variable individually (in the same lines as the idea of one random projection being enough to detect any difference-see Section 1.4.2 and [CAFR07]).

For example, we compared two protein complexes datasets based on 19 energy measures1 : one set consisting of complexes found in nature and the other one consisting of complexes produced by docking algorithms from their individual components. In this case, the univariate Mann-Whitney test and the Hodges-Lehmann estimator where able to show strong differences for some of the energies. Another example consisted in comparing two datasets measuring cytokines' concentration2 in sick (rheumatoid arthritis) and healthy individuals [DFS + 15]. Here also, univariate analysis was able to detect and measure important differences in some individual cytokines' concentration.

Nevertheless, situations like the one depicted in Figure 5.1 are not detected with variable-wise testing. Therefore, we believe that, on a general perspective, the ability to provide nonparametric multivariate feedback (i.e., effect size) downstream two-sample testing may prove of ubiquitous interest in exploratory statistics and data science.

In order to continue this line of research, it would be interesting to characterize and provide theoretical guarantees on the cluster decomposition found by our procedure. Furthermore, it would be interesting to generalize the tool for data in (non-Euclidean) metric spaces in order to increase its range of applicability.

With respect to our sequential nonparametric two-sample tests, we believe that their novel features, namely sequentiality, suitability for streaming data and semi-supervised nature, clearly extend the range of applicability of classical two-sample tests, for example, towards big data streaming scenarios.

From a theoretical perspective, it would be interesting to find soft classifiers for data in (non-Euclidean) metric spaces that would provide consistency guarantees as well.

Another interesting line of research would be on the practical applied side where consistency might be of less importance: since any online soft classifier can be used to obtain a valid p-value, it is possible to use this framework to compare data as diverse as, for example, two streams of text messages. Our change detection test requires building hierarchical spatial partitions.

Space partitioning: Geometry The Computational Geometry Algorithms Library (CGAL) [cga] provides a large set of geometric operations, which are in particular useful for projecting, e.g, CGAL::Vector_d::operator* (const Vector_d< Kernel > &w)

and for defining hyperplanes or other oriented hypersurfaces for partitioning, e.g., Space partitioning: Hierarchy A hierarchical partitioning can be represented by a tree structure and an array containing all the sample points. Points do not need to be replicated in each node to represent the content of each cell. Instead, the array can be recursively sorted in a quicksort fashion by using as sorting criterion the side on the separating hypersurface on which the points lie.

Then, the elements contained in each cell can be simply represented by the boundaries delimiting the portion of the array that contains them.

Nonparametric Methods for Learning and Detecting Multivariate Statistical Dissimilarity

In this thesis, we study problems related to learning and detecting multivariate statistical dissimilarity, which are of paramount importance for many statistical learning methods nowadays used in an increasingly number of fields. This thesis makes three contributions related to these problems. The first contribution introduces a notion of multivariate nonparametric effect size shedding light on the nature of the dissimilarity detected between two datasets. Our two step method first decomposes a dissimilarity measure (Jensen-Shannon divergence) aiming at localizing the dissimilarity in the data embedding space, and then proceeds by aggregating points of high discrepancy and in spatial proximity into clusters.

The second contribution presents the first sequential nonparametric two-sample test. That is, instead of being given two sets of observations of fixed size, observations can be treated one at a time and, when strongly enough evidence has been found, the test can be stopped, yielding a more flexible procedure while keeping guaranteed type I error control. Additionally, under certain conditions,using nearest neighbors regression, when the number of observations tends to infinity, the test has a vanishing probability of type II error.

The third contribution presents a semi-supervised sequential nonparametric two-sample test, this time based on random spatial partitioning. The test also exhibits a vanishing type II error, under certain conditions on the partitions. This test automatically detects multi-scale differences. Processing a new observation has logarithmic time complexity w.r.t. the unlabeled training dataset size, and negligible memory footprint, which makes it suitable for streaming data.

Méthodes non-paramétriques pour l'apprentissage et la détection de dissimilarité statistique multivariée

Cette thèse présente trois contributions en lien avec l'apprentissage et la détection de dissimilarité statistique multivariée, problématique d'importance primordiale pour de nombreuses méthodes d'apprentissage utilisées dans un nombre croissant de domaines. La première contribution introduit la notion de taille d'effet multivariée non-paramétrique, éclairant la nature de la dissimilarité détectée entre deux jeux de données, en deux étapes. La première consiste en une décomposition d'une mesure de dissimilarité (divergence de Jensen-Shannon) visant à la localiser dans l'espace ambiant, tandis que la seconde génère un résultat facilement interprétable en termes de grappes de points de forte discrépance et en proximité spatiale.

La seconde contribution présente le premier test non-paramétrique d'homogénéité séquentiel, traitant les données issues de deux jeux une à une-au lieu de considérer ceux-ci in extenso. Le test peut ainsi être arrêté dès qu'une évidence suffisamment forte est observée, offrant une flexibilité accrue tout en garantissant un contrôle de l'erreur de type I. En utilisant des régresseurs basés sur les plus proches voisins, sous certaines conditions, la convergence vers 0 de la probabilité d'erreur de type II est établie.

La troisième contribution présente un test non-paramétrique d'homogénéité séquentiel semi-supervisé, basé sur des partitions aléatoires de l'espace. Le test a aussi une probabilité d'erreur de type II tendant vers zéro sous certaines hypothèses. Il détecte automatiquement des différences multi-échelles. Le traitement d'une nouvelle observation requiert un coût logarithmique en la taille du jeu de données non étiquetées utilisé pour l'apprentissage, et une taille mémoire négligeable. Ceci le rend particulièrement adapté aux flux de données.

  3 A P -PUD based on Jeffreys mixture . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Using a hierarchical ensemble of partitions associated with a tree . . . . . . . . . . . . . . 4.3.1 Hierarchical partitioning based on random projection trees . . . . . . . . . . . . . 4.3.2 Consistency via context tree weighting . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Increasing power using switch distributions . . . . . . . . . . . . . . . . . . . . . . 4.4 Using an ensemble of trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2 Choice of parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Conclusion 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Discussion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Resources: Algorithms and Software A.1 Resources for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Resources for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Resources for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  3 and C 4 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Partition of the persistence diagram exploited to define clusters. The axes x and y correspond respectively to the birth and death dates. The partition of the domain y ≥ x is induced by three lines: (i) y = x + ρ which specifies the persistence threshold (ii,iii) x = -δ max , y = -δ max , with δ max the threshold on the significance of the discrepancy. See text for the specification of regions R 1 to R 5 . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Workflow of the whole method. In blue: the parameters. . . . . . . . . . . . . . . 2.5 Rotated images (a) Original image i (b,c,d,e) Example rotated images . . . . . . . . . . 2.6 Subsets of handwritten digits used. Images cropped from http://www.cs.nyu.edu/ ~roweis/ data.html. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Original satellite color images from [STS09]. . . . . . . . . . . . . . . . . . . . . . . 2.13 Original satellite image (Fig. 2.7): the red and blue squares correspond to the local minima of the two clusters identified from the analysis presented on Fig. 2.12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Model: Gaussian mixture. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition . . . . . . . . . . 2.9 Model: Gaussian mixture. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition . . . . . . . . . . 2.10 Model: Crenels. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition . . . . . . . . . . . . . . . . 2.11 Model: Handwritten digits. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition . . . . . . . . . . 2.12 Model: Satellite images. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition . . . . . . . . . . . . .

  v

  (C) The partition associated with the dashed cut consists of the two grayed cells. . . . . . . . . . . . . . . . . . . . . . 5.1 Projection onto x or y axis does not allow to distinguish the datasets. . . . . . vi List of Tables 3.1 Gaussian data in dimension d = 4. The symbol H 0 or H 1 at the beginning of each line indicates whether the null is true or false. Numbers indicate the percentage of trials for which the null hypothesis was rejected, given α = 0.05. A total of 150 trials were done. . . 3.2 Gaussian data in dimension d = 10, 30. For conventions, see the caption of

Figure 2 . 1

 21 Figure 2.1 Random multiplexer generating pairs (label, position).L ∼ B ( 1 /2) X Y (L, Z) B: Bernoulli distribution

Figure 2 . 5

 25 Figure 2.5 Rotated images (a) Original image i (b,c,d,e) Example rotated images

Figure 2 . 6

 26 Figure 2.6 Subsets of handwritten digits used. Images cropped from http://www.cs.nyu.edu/ ~roweis/ data.html.

Figure 2 . 7

 27 Figure 2.7 Original satellite color images from [STS09].

Figure 2 .

 2 Figure 2.13 Original satellite image (Fig. 2.7): the red and blue squares correspond to the local minima of the two clusters identified from the analysis presented on Fig. 2.12.

Figure 2 .

 2 Figure 2.8 Model: Gaussian mixture. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition

Figure 2 . 9

 29 Figure 2.9 Model: Gaussian mixture. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition

Figure 2 .

 2 Figure 2.10 Model: Crenels. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition

Figure 2 .

 2 Figure 2.11 Model: Handwritten digits. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition

Figure 2 .

 2 Figure 2.12 Model: Satellite images. Figs, from top to bottom: Raw data embedding, Discrepancy plot, Persistence diagram, Clusters plot, Divergence decomposition

Figure 3 . 1

 31 Figure 3.1 The random multiplexer used to generate pairs (label, position)L ∼ B (θ 0 ) X Y (L, Z) B: Bernoulli distribution

Definition. 3

 3 .2. A sequence of regression function estimates {m n } is strongly pointwise consistent (s.p.c.) if m n (z) n→∞ ----→ m(z) a.s. (3.18) for µ-almost all z ∈ R d , µ denoting the distribution of Z.

  L|Z) a.s.. (3.26) Since P λ n is based on an s.p.c. regressor, we have that P-i l|z, (L, Z) → P θ(z) (l|z) a.s. (3.27) for all l ∈ A and f -almost all z ∈ R d . Since the measure of z values failing the previous convergence is null, one has

4. 2 . 1

 21 Spatial partition and discretized distributionA spatial partition P = C 1 , C 2 , . . . , C |P | of R dis a finite collection of cells of mutually disjoint interior, and whose union equals R d . Averaging θ(z) ≡ P (L = 0|Z = z) yields

Figure 4 . 1

 41 Figure 4.1 Tree and partitions (A) Partition of a 2D domain by a binary space partition tree, with two levels. (B)The corresponding full tree. To each node is associated a region; internal nodes also come with a splitting rule to define the regions of the two children. Note that a tree T generates a set of partitions. Identifying a node with its associated region, these partitions are P T = {η 1 , {η 2 , n 3 }, {η 2 , η 4 , η 5 }}. (C) The partition associated with the dashed cut consists of the two grayed cells.

Figure 5 . 1

 51 Figure 5.1 Projection onto x or y axis does not allow to distinguish the datasets.

••

  TreeRank (AUC based test) [CDV09]: TreeRank::TwoSample 2 Multivariate Cramer [BF04]: cramer::cramer.test • Nearest-Neighbors [Sch86]: MTSKNN::mtsknn • Mann-Whitney [MW47] (univariate): stats::wilcox.test and coin::wilcox test (the latter provides exact p-values in the presence of ties). • Kolmogorov-Smirnov [Kol33, Smi48] (univariate): stats::ks.test A.3 Resources for Chapter 4

  CGAL::Hyperplane_d CGAL::Sphere_d CGAL::Oriented_side Random directions can be generated using the function CGAL::Random_points_on_sphere_d<Point_d>(int dim, double r, Random &rnd=default_random) or by generating independent gaussian random numbers for each component (see, e.g., [Mul59, M + 72]) using, e.g., the boost class 3 boost::normal_distribution Fast exact and approximate algorithms for diameter computation (required, for example, for RPTree-Mean version of Random Projection Trees[START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF]) were proposed in[START_REF] Malandain | Computing the diameter of a point set[END_REF].4

  Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Statistical Dissimilarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Dissimilarity Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Binary Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Learning Dissimilarity Measures and Classifiers . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Desired properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Methods of Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Universal Prediction Methods for Soft Classification . . . . . . . . . . . . . . . . . 1.4 Detection: Hypothesis Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Estimating the discrepancy between datasets . . . . . . . . . . . . . . . . . . . . . . . . . Localizing the discrepancy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Combining discrepancy estimation and localization . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Qualifying the clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Model: Gaussian mixture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Model: crenels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.3 Model: mixture of handwritten digits . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.4 Statistical image comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Supplemental: Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sequential Non-parametric Two-Sample Test via Nearest Neighbors 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Two-sample test based on sequential prediction . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Random labels framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Notations and problem reformulation . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Robust sequential p-value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Consistency via λ-pointwise universal distributions (λ-PUD) . . . . . . . . . . . . 3.3 λ-Pointwise universal distributions via strongly pointwise consistent regressors . . . . . . 3.4 Increasing power using mixtures and switch distributions . . . . . . . . . . . . . . . . . . 3.5 Implementation and complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Instantiations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Contenders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Appendix: Nonparametric regression based on k n -nearest neighbors . . . . . . . . . . . . 3.9 Appendix: Optional Stopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 A Sequential Non-parametric Two-Sample Test via Random Projection Context Trees 59 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	1 Introduction
	1.1

1.4.1 Hypothesis Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 The Two-sample Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5

2 Beyond Two-sample-tests: Localizing Data Discrepancies in High-dimensional Spaces 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Comparing datasets in high dimensional spaces . . . . . . . . . . . . . . . . . . . . 2.1.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 2.2.1 Jensen-Shannon divergence decomposition using conditional distributions . . . . . 2.2.2 Conditional probability estimation via non-parametric regression . . . . . . . . . . 2.2.3 Joint distribution compatible sampling . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 2.3.1 Defining clusters from sublevel sets . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Data structures and algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 2.7.1 A refined strategy to compute clusters . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Supplemental: Data Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Gaussian mixture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i 3 A 4.2 Discretizing the sample space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Spatial partition and discretized distribution . . . . . . . . . . . . . . . . . . . . . 4.2.2 Consistency via pointwise universal distributions on a revealing partition . . . . . 4.2.

  Table 3.1. In this experiment, the data under H 0 was generated by sampling from the mixture of both distributions. 100 trials were done. (NB: the conditions associated with TRank ref are specified in M. Depecker's PhD Thesis [Dep10].) . . . . . . . . . . . . . . . . . . . . . 3.3 Lattice of Gaussians in dimension d = 2. For conventions, see the caption of Table 3.1. The data under H 0 was generated by sampling from the mixture of both distributions.

	Chapter 1
	Introduction
	1.1 Motivation

A total of 200 trials were done. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Gaussian data in dimension d = 4. The symbol H 0 or H 1 at the beginning of each line indicates whether the null is true or false. Numbers indicate the percentage of trials for

  Type II error and λ-consistency. Note that the constant term in Eq. (3.19) guarantees that the logarithm of Q λ is bounded. This allows applying Breiman's Extended Ergodic Theorem[START_REF] Breiman | The individual ergodic theorem of information theory[END_REF] to obtain the following theorem. -1 , Z -1 ), (L 0 , Z 0 ), (L 1 , Z 1 ), . . . . Let us consider the operator T i that shifts any sequence {. . . , s -1 , s 0 , s 1 , . . . } by i positions to the left and let us denote ĝλ

						.19)
	Plugging Pn into the previous equation yields P λ n , which, by lemma 3.2, can be turned into a distribution
	on length-n sequences and, thus, we have the following corollary.
	Corollary. 3.1. The test obtained with P λ n yields a robust p-value.
	Theorem. 3.3. P λ n based on a s.p.c. regressor yields a λ-PUD, thus a λ-consistent test. That is,
		lim n	1 n	log P λ n (L n |Z n ) = E -log P λ θ(z) (L|Z) a.s.	(3.20)
						≤ H (L|Z) -log λ	(3.21)
	Proof. Let us define the doubly infinite i.i.d. process (L, Z)
	The claim	lim n	1 n	log P λ n (L n |Z n ) = E -log P λ θ(z) (L|Z) a.s.	(3.22)
	is equivalent, by Lemma 3.2, to
	lim n	1 n	n i=1	log P λ i L i |Z i , L i-1 = E -log P λ θ(z) (L|Z) a.s..	(3.23)

∞ -∞ ≡ . . . , (L

Table 3 .

 3 1 Gaussian data in dimension d = 4. The symbol H 0 or H 1 at the beginning of each line indicates whether the null is true or false. Numbers indicate the percentage of trials for which the null hypothesis was rejected, given α = 0.05. A total of 150 trials were done. Table 3.2 Gaussian data in dimension d = 10, 30. For conventions, see the caption of Table 3.1. In this experiment, the data under H 0 was generated by sampling from the mixture of both distributions. 100 trials were done. (NB: the conditions associated with TRank ref are specified in M. Depecker's PhD Thesis [Dep10].) Case KNN .3 KNN .5 KNN .7 KNN .9 BMA SW (F )

	Case	KNN .3 KNN .5 KNN .7 KNN .9 BMA SW (F ) πU	SW (F ) πS	SW (OS) πU	SW (OS) πS	TRank TRank ref MMD b
	H 0 , Ex1	0	0	0	0	0	0	0	2	6	4.7	1	3.3
	H 1 , Ex2	5.3	100	100	100	100	100	100	100	100	100	99	100
	H 1 , Ex3a 0	0	0	0	0	18.7	5.3	20	9.3	40.7	45	90
	H 1 , Ex3b 0	0	1.3	0	0.7	48.7	16	65.3	21.3	76	73	98.7
							πU	SW (F ) πS	SW (OS) πU	SW (OS) πS	TRank TRank ref MMD b
	H 0 , d = 10 0	0	0	0	0	1	0	1	5	5		6
	H 0 , d = 30 0	0	0	0	0	0	0	0	3	6		6
	H 1 , d = 10 0	0	0	0	0	69	20	68	35	36	90	99
	H 1 , d = 30 0	0	0	0	0	22	5	29.5	10	25	85	95

Table 3 .

 3 3 Lattice of Gaussians in dimension d = 2. For conventions, see the caption of Table3.1. The data under H 0 was generated by sampling from the mixture of both distributions. A total of 200 trials were done.In this chapter, we propose a sequential two-sample test using sequential conditional probability estimation based on spatial partitioning instead of nearest neighbors as in Chapter 3. That is, we base again our test on the result of theorem 3.1 but we propose a different construction for Q based on a discretization of Z.

	Case H0, B1	KNN.3 0	KNN.5 0	KNN.7 0	KNN.9 0	BMA 0	SW (F ) π U 0	SW (F ) π S 0	SW (OS) π U 0.5	SW (OS) π S 5.5
	H1, B2	11	0.5	0	0	11	21.5	14	27.5		16.5
	H1, B3	78.5	62.5	0	0	85	93	89.5	96.5		91.5
		Case	TRank	MMD b	OPT	MxRt	MxMMD	L2	xvalc	med
		H0, B1	6	5	4	2.5	2.5	4	7	7
		H1, B2	100	6	18.5	15.5	15	14.5	15.5	5
		H1, B3	100	6.5	21	18.5	15	15	21	7.5

Table 4 .

 4 1 Gaussian data in dimension d = 4. The symbol H 0 or H 1 at the beginning of each line indicates whether the null is true or false. Numbers indicate the percentage of trials for which the null hypothesis was rejected, given α = 0.05. A total of 150 trials were done. Table 4.2 Gaussian data in dimension d = 10, 30. For conventions, see the caption of Table 4.1. In this experiment, the data under H 0 was generated by sampling from the mixture of both distributions. 100 trials were done. (NB: the conditions associated with TRank ref are specified in M. Depecker's PhD Thesis [Dep10].)

	Case	SW (OS) π U	SW (OS) π S	CTW2 (OS)	CTS2 (OS)	TRank TRank ref MMD b
	H 0 , Ex1	2	6	2.7	3.3	4.7	1	3.3
	H 1 , Ex2	100	100	100	100	100	99	100
	H 1 , Ex3a 20	9.3	32	61.3	40.7	45	90
	H 1 , Ex3b 65.3	21.3	84.7	97.3	76	73	98.7
	Case	SW (OS) π U	SW (OS) π S	CTW2 (OS)	CTS2 (OS)	TRank TRank ref MMD b
	H 0 , d = 10 1	5	1	2	5		6
	H 0 , d = 30 0	3	1	2	6		6
	H 1 , d = 10 68	35	16	47	36	90	99
	H 1 , d = 30 29.5	10	2	13	25	85	95

Table 4 .

 4 3 Lattice of Gaussians in dimension d = 2. For conventions, see the caption of Table4.1. The data under H 0 was generated by from the mixture of both distributions. A total of 200 trials were done.

		Case	SW (OS) π U	SW (OS) π S	CTW2 (OS)	CTS2 (OS)
		H 0 , B1 0.5	5.5		2.5	2.5
		H 1 , B2 27.5	16.5		100	100
		H 1 , B3 96.5	91.5		100	100
	Case	TRank MMD b OPT MxRt MxMMD L2	xvalc med
	H 0 , B1 6	5	4	2.5	2.5	4	7	7
	H 1 , B2 100	6	18.5	15.5	15	14.5 15.5	5
	H 1 , B3 100	6.5	21	18.5	15	15	21	7.5

P X is absolutely continuous with respect to P Y if for every measurable set A, P Y (A) = 0 implies P X (A) = 0.

The Radon-Nikodym derivative dP X dP Y is the function f such that, for every measurable set A, P X (A) = A f dP Y . For continuous random variables, the density is the Radon-Nikodym derivative of the corresponding (induced) measure with respect to the Lebesgue measure.

Note that, in fact, there are 2 d -1 independent orderings since one depends on all the others, for example, in 2D, P (> x, > y) = 1 -[P (< x, < y) + P (< x, > y) + P (> x, < y).

the odd case is treated by using a pseudopoint n + 1 having distance 0 to every point, then optimally matching the n + 1 points and discarding the pair containing the pseudopoint

Apparently, this was first introduced in[START_REF] Dwass | Modified randomization tests for nonparametric hypotheses[END_REF]

In contrast to fixed design regression, where [GK02, Sec. 1.9]: "one observes values of some function at some fixed (given) points with additive random errors, and wants to recover the true value of the function at these points."

http://sbl.inria.fr/doc/Density_difference_based_clustering-user-manual.html

When we write that Q is a distribution over infinite sequences l 1 , l

, . . ., we mean that Q is a probability distribution relative to some sample space Ω and associated σ-algebra F such that L 1 , L 2 , . . . are all random variables on the probability triple (Ω, F , Q) -see [Wil91, Section 4.5-4.8]. This automatically insures that the compatibility conditionl ′ ∈A Q (l n , l ′ ) = Q (l n ) holds.

using, e.g., Platt's scaling to obtain probabilities.

In our experiments, we check whether |log(y) -log(x)| > 100.

A full binary tree is a tree in which every node in the tree has either 0 or

children. 2 A perfect depth-K binary tree is a binary tree in which all interior nodes have two children and all leaves have the same depth K.

provided by Rosetta software, see http://rosettadesign.med.unc.edu/documentation.php.

Cytokines are a broad category of small proteins that are involved in cell signaling.

Available at www.gatsby.ucl.ac.uk/ ~gretton/mmd/mmd.htm and www.gatsby.ucl.ac.uk/ ~gretton/adaptMMD/ adaptMMD.htm.

Available at http://cran.r-project.org/src/contrib/Archive/TreeRank

Available at www.boost.org.

A C implementation can be found at http://www-sop.inria.fr/members/Gregoire.Malandain/diameter/.

Acknowledgments

σ 0 [1] = 0.04771839I 2 , σ 0 [2] = 0.02111844I 2 , σ 0 [3] = 0.03342965I 2 , σ 0 [4] = 0.07551961I 2 where I 2 is the 2 × 2 identity matrix.

The second population is drawn according to a mixture of four non-spherical gaussians with equal probability defined by the following randomly generated parameters (using analogous notation): Given two sets of samples x 1 , . . . , x n0 and y 1 , . . . , y n1 whose corresponding random variables X i ∈ R d and Y i ∈ R d are i.i.d. with densities f X and f Y respectively, a nonparametric two-sample test ambitions to determine whether f X = f Y . To this end, a statistic, i.e., a function aiming at revealing discrepancies of the data is defined. This function typically quantifies in a global way the local homogeneity of the mixture of the populations, this local homogeneity being assessed by methods as diverse as nearest-neighbors (e.g. [START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF][START_REF] Henze | A multivariate two-sample test based on the number of nearest neighbor type coincidences[END_REF][START_REF] Rosenbaum | An exact distribution-free test comparing two multivariate distributions based on adjacency[END_REF][START_REF] Pérez-Cruz | Kullback-Leibler divergence estimation of continuous distributions[END_REF]), spatial partitions (e.g. [START_REF] Biau | On the asymptotic properties of a nonparametric l1-test statistic of homogeneity[END_REF]) or kernels (e.g. [GBR + 12]). The way the data are processed allows classifying two-sample tests into two tiers. Tests of the first tier compute a statistic on the whole dataset, to reveal the discrepancy if any (e.g. classical tests like [START_REF] Bickel | A distribution free version of the smirnov two sample test in the p-variate case[END_REF][START_REF] Friedman | Multivariate generalizations of the wald-wolfowitz and smirnov two-sample tests[END_REF][START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF][START_REF] Henze | A multivariate two-sample test based on the number of nearest neighbor type coincidences[END_REF] and also more recent ones like [HT02, Ros05, GBR + 12]). Tests from the second tier use a learning perspective by splitting the data into a training and a test sets (e.g. [Fri04, CDV09, GSS + 12, ZGB13]).

In the training phase, the function aiming at revealing the discrepancies is learned and optimized (in the sense of its discrimination power). In the second phase, this function is evaluated on the complement of the training set to obtain the statistic. An appealing extension to this second tier is to use a sequential framework in which the function is optimized at each sample and then immediately used to test the next sample. In this sequential framework, it is natural and desirable to be able to stop at any time, i.e., as soon as a difference has been perceived. Whereas classical Neyman-Pearson null hypothesis testing requires the sampling plan or equivalently the stopping rule to be defined in advance to ensure the validity of the procedure (see, e.g., [START_REF] Wagenmakers | A practical solution to the pervasive problems of p values[END_REF]), Bayes factor model comparison makes this optional stopping possible (see, e.g., [START_REF] Shafer | Test martingales, Bayes factors and p-values[END_REF]). Note that selected tests, such as MMD l [GBR + 12], process samples sequentially, yet require a predefined sampling plan.

Contributions

We make three contributions. First, we design a framework for sequential nonparametric two-sample tests. The framework is based on sequential prediction of labels defining the two populations (or equivalently distributions on length-n sequences), and enjoys optional stopping. Second, under suitable conditions qualifying the difference between the tested distributions, consistency of the two-sample test is guaranteed when the sequential predictor is built from strongly pointwise consistent regressors, obtained in our case from k n -nearest neighbors (KNN) regressors. Third, we show that combing mixtures and switch and thus

We also have that ĝλ i is bounded between 0 and log 2 1-λ , so that the claim follows from Breiman's extended ergodic theorem [START_REF] Breiman | The individual ergodic theorem of information theory[END_REF].

The inequality 3.21 stems from P λ θ(z) (L|Z) ≥ λP θ(z) (L|Z).

As shown in [Alg92, Thm. 2], the following mixture is a pointwise universal distribution

where 0 ≤ λ k ր 1 and ∞ k=0 µ k = 1. Indeed, for each λ used in the mixture, one gets a λ-PUD by Lemma 3.4. Consequently, the mixture obtained for λ arbitrarily close to 1, yields a λ-PUD, approaching therefore the ideal case of λ = 1. Yet, this mixture requires an infinite sum, which motivates our use of λ-PUD -which yield λ-consistent tests.

Increasing power using mixtures and switch distributions

Let us first consider Bayesian Model Averaging (BMA). Given a set of distributions {Q k } k and weights µ k such that k µ k = 1, BMA produces the following mixture:

(3.33)

Pointwise universal mixtures. We first prove a lemma that allows building λ-consistent two-sample tests by mixing one λ-PUD with arbitrary distributions:

Lemma. 3.4. A mixture of distributions containing at least one λ-PUD Q 0 with weight µ 0 (n) such that log µ 0 (n) = o(n) yields a λ-PUD.

Proof. Consider a λ-PUD Q 0 and the arbitrary distributions Q k (l n |z n ) , k > 0. Assuming that k µ k = 1, we define the mixture:

An interesting application of mixtures is the following one. Since it is not clear which neighborhood size function k n is best for a KNN regressor, we consider a set of functions k n = n p , where p takes values in some set -finite for practical purposes. All these predictors yield λ-PUD if p < 1 (see Appendix 3.8). But using a mixture allows one to favor the best size function by updating the a posteriori weights according to past performance.

• Ex1 : two populations drawn from the same distribution, so that Type I error is assessed. Sample sizes: n 0 = n 1 = 1000.

• Ex2 : two populations drawn from two shifted Gaussians, so that Type II error is assessed. Sample sizes: n 0 = n 1 = 1000.

• Ex3 : corresponds to two subtly shifted Gaussians. Sample sizes: Ex3a: n 0 = n 1 = 3000; Ex3b:

Regarding Ex1 (Table 3.1, line 1), our tests using the full sequence are conservative -no type I error, an observation also valid in the other experiments. For the remaining experiments (Table 3.1, lines 2-4), we observe that the switching distribution with a uniform prior is the best amongst our predictors. Also, its power is comparable to that of TRank, yet worse than that of MMD b (Ex3b).

Gaussian data in dimension d = 10, 30. Table 3.2 presents the results on data generated using the same specifications of [START_REF] Clémençon | AUC optimization and the two-sample problem[END_REF][Fig. 2], which correspond to 10 and 30-dimensional gaussians shifted by ∆µ = 0.2 with n 0 = n 1 = 2000. Here we estimate Type II error probability.

Although individual predictors KNN i yield the worst performance, SW π U shows good performances, even outperforming TRank in dimension 10 -in our replica, which use default parameters. We notice, though, that performances degrade upon increasing the dimension, a likely consequence of distance concentration phenomena perturbing the choice of neighbors by KNN.

Lattice data and incidence of the scale. Table 3.3 presents the results on data corresponding to the specification of [GSS + 12]: two 5 × 5 grids of two-dimensional Gaussians (a.k.a. blobs) that differ in covariance (parameters: stretch=10, rotation angle=π/4, blob distance=15). We consider three cases: -B1: two populations drawn from the mixture of both blobs, so that Type I error is assessed. Sample sizes:

-B2: two populations drawn from each of the blobs, so that Type II error is assessed. Sample sizes:

-B3: Same as B2 but with larger sample sizes:

As noted in [GSS + 12], it is a prototypical example where MMD b fails and kernel width optimization is important, since differences occur at a smaller scale. For our predictors, this scale is captured by slower k n . It is important to emphasize, that predictors with larger k n are also consistent and, therefore, they would also detect the differences with more samples. Remarkably, when n 0 = n 1 ≥ 1500, SW π U (F or OS) outperforms all the MMD contenders.

Conclusion

This work introduces the first sequential nonparametric two-sample test, based on sequential prediction of labels defining the two populations. Our test is shown to be consistent when prediction is carried out by strongly pointwise consistent regressors. We show that mixtures of distributions based on KNN regressors are effective in favoring the best neighborhood size function. This update being carried out along the sequential process is more flexible than classical approaches splitting the data into a training and test sets. We also show that model switching increases the power, a fact related to the ability of automatically selecting the best splitting point in a train/test paradigm. Experimentally, while no test is expected to be the most powerful for all kinds of data, our best constructs outperform state-of-the-art ones on selected challenging datasets.

We foresee two outstanding questions. Complexity-wise, the regressors used rely on exact nearest neighbor queries, exhibiting linear complexity in the worst-case. Inferring whether our tests remain consistent when information of lower quality is used (e.g. approximate nearest neighbors) would allow using more efficient data structures. In addition, obtaining consistency guarantees under constant size memory requirements would be of special interest in a streaming context.

Then we have the following theorem [GK02, Thm. 25.17]:

then the k n -NN estimate using Euclidean distance is strongly pointwise consistent.

Appendix: Optional Stopping

In order to prove Theorem 3.1, we use the following result Theorem. 3.5 ([vdPG14, Thm. 3.1], special case of [START_REF] Shafer | Test martingales, Bayes factors and p-values[END_REF]). Given a sample space X and a set of parameters M , suppose we have data X 1 , . . . , X n i.i.d. ∼ P µ , X i ∈ X , µ ∈ M , and a test rejects the null hypothesis

where α is a constant and p 1 is a function from X n to R + 0 . If p 1 (x n ) is such that, for all n, the compatibility condition

holds and E µ0 p1(X n )

Pµ 0 (X n ) ≤ 1, then: Nevertheless, ML does not produce a distribution over {0, 1} n as required by Thm. 3.1. It can be corrected using Normalized Maximum Likelihood (NML) [START_REF] Shtar'kov | Universal sequential coding of single messages[END_REF][START_REF] Rissanen | Fisher information and stochastic complexity[END_REF] but it is hard to compute. Instead, we shall use the following Bayesian mixture, which is asymptotically equivalent to NML (when taking the normalized logarithm) :

with w(•) being Jeffreys' prior for the Bernoulli distribution [START_REF] Jeffreys | An invariant form for the prior probability in estimation problems[END_REF] defined, in this case, as

Using the gamma function Γ, it can be expressed as follows

Therefore, for any cell C of the partition P

Then, since C∈P n(C) = n, using the log sum inequality, we have that

Then, by adding inequality 4.19 for each cell of the partition and taking, in each cell C i , θ = θ Ci from θ P , we obtain the inequality of the statement.

The limit follows from the Asymptotic Equipartition Property (see, e.g., [START_REF] Cover | Elements of Information Theory[END_REF]).

Consistency via context tree weighting

Intuitively, the finer the partition the larger the chances of finding revealing cells with the largest difference θ Cθ 0 possible, which should be easier to detect. On the other hand, Eq. (4.16) shows that there is a cost on the convergence rate that is proportional to the size of the partition. For this reason, we are interested in distributions that optimize the number of parameters at the same time they optimize the distribution itself, i.e., twice-universal distributions. CTW. In order to build a twice-universal distribution, we use a Bayesian mixture of the distributions P Jeffreys(P ) n associated to each partition P ∈ P T using as weight w nat (P ) ≡ 2 -2|P |+1 , obtaining therefore

This can be efficiently computed using the recursive algorithm called Context Tree Weighting (CTW) [START_REF] Frans | The context-tree weighting method: Basic properties[END_REF] (with the generalization proposed in [START_REF] Frans | Context weighting for general finite-context sources[END_REF]), where z gives the context and each node can be split in only one way (accordingly to T ). This combination of spatial partitions and CTW has already been used in [START_REF] Kozat | Universal portfolios via context trees[END_REF] to build universal portfolios.

For any cell C partitioned into a set of cells {C i }, we have the following recursive formula

and for a cell C which cannot be split (leaf node)

Let C 0 be the cell corresponding to the root node of a tree T , then it can be proved by induction [WST95, Lemma 2] that P CTW(P T )

The next lemma shows that P CTW(P T ) n is asymptotically as good as any distribution P θ P such that P ∈ P T , the convergence rate depending on the size of P .

Lemma. 4.3. P CTW(P T )

n is pointwise twice-universal for the class of distributions defined by P T since it satisfies for every l n , z n and for any partition P ∈ P T :

-

n1 can be used to avoid estimating θ 0 which is already known.

Remark 4.3 (Sequential computation of CTW). At each internal node except the root (cf. Remark 4.2), by Eq. (4.23), CTW performs a Bayesian mixture of two models (Jeffreys and CTW from child nodes). This can be sequentially computed as in Eq. (3.38) and (3.39). Therefore, sequential computation of CTW only requires maintaining these weights for each internal node plus the counters mentioned in Remark 4.1. More precisely when a new observation (l, z) is processed, the tree is traversed from the root to a leaf node, following the path corresponding to the sequence of nested cells containing z, whose length is O (log N ), with N the number of unlabeled samples.

Increasing power using switch distributions

Since CTW is based on Bayesian mixtures it is also prone to the catch-up phenomenon (see Section 3.4). In order to avoid it and hopefully improve power, we use the recursive switch distribution called Context Tree Switching (CTS) proposed in [START_REF] Veness | Context tree switching[END_REF]. The intuitive idea is to allow using simpler models (i.e., smaller partitions) when few observations are available, and then switch to more complex models (i.e., finer partitions) when more data is available.

Like in CTW (see Eq. (4.23)), the recursive formula considers at each internal node two possible models, i.e., the one given by Jeffreys' estimation for the corresponding cell and the one combining the models corresponding to each child node. These two models are indexed by the symbols a, b respectively.

Let w C (•) be a prior over model index sequences i n(C) ≡ i 1 , . . . , i n(C) ∈ {a, b} n(C) at cell C, recursively defined by

where α n ≡ 1 n is the switch rate. The previous equation calls for two comments:

• the first line is needed, as cells void of points appear in the formula below, and

• the last line corresponds to the cases where one does not switch models, and switches models, respectively.

To define the formula replacing Eq. (4.23) in the context of CTS, we shall need:

• C (l n ): the subsequence of l n corresponding to the symbols that fell in C

• t C (k) : the smallest integer implicitly defined by the following equation:

(4.28)

The following equation replaces Eq. (4.23) and defines the CTS mixture at each internal node C (assuming two child nodes C 1 and C 2 ):

Using an ensemble of trees

When considering finite length performance, we can be unlucky and obtain a bad set of hierarchical partitions. Therefore, in order to reduce variance and, hopefully, to improve the probability of rejecting under the alternative hypothesis, we use an ensemble of J random trees. The corresponding partitions are denoted {P T1 , . . . , P T J }. Using them, we define the following mixture with the uniform prior w unif (J) ≡

Experiments

Implementation

The random projection context trees are represented by a binary tree structure in which each node representing a cell C stores the counters n l (C), l ∈ {0, 1}. Internal nodes also store a representation of the hyperplane used by the RP tree to split the cell into its two children. Therefore, when a new sample is processed the counters and the weights to be updated are those located in the path from the root node to the leaf node reached by going to the proper side at each internal node according to the side of the hyperplane the sample is located at.

Observation. 2. Processing a new observation has time complexity O(J log N ).

Observation. 3. Potentially, the counters could overflow, requiring an arbitrary-precision counter or, otherwise, approximate counting (see, e.g., [START_REF] Flajolet | Approximate counting: a detailed analysis[END_REF]).

Choice of parameters

We used J = 1000 trees and the whole sequence (since in this experiments they are known in advance) as unlabeled training set.

Results

We tested on the same datasets considered in Section 3.6 and we compared against the results presented there. The results are presented in Tables 4.1, 4.2 and 4.3.

Overall, we observe an empirical Type I error under its design value of α = 0.05, which is expected since the p-value is valid but not exact.

In all the cases, CTS2 exhibited a larger power than CTW2, therefore confirming the interest of switching to avoid the catch-up phenomenon. The two steps of our nonparametric effect size method are implemented using the following resources.

Exact and approximate nearest neighbors The first step consists in estimating a decomposition of the Jensen-Shannon divergence using k-nearest-neighbor based soft classification. Exact nearest neighbors soft classification is implemented in the knn3 function in the R package caret [fJWWW + 15]. Since this package does not allow incremental update of the training set, it is best suited for offline applications.

k-Nearest neighbor graph, watershed transform and topological persistence The second stage of our nonparametric effect size method combines a Morse theoretical analysis carried out on a k-nearestneighbors-graph (connecting the data points, each equipped with its estimated discrepancy), with a simplification resorting to topological persistence (aiming at identifying stable features). These tools are implemented in the Structural Bioinformatics Library ([CD15] and http://sbl.inria.fr), and will be released in the fall of 2015.

A.2 Resources for Chapter 3

Our sequential two-sample test relies on k n -nearest-neighbor online soft classification and it was compared against a number of tests.

Exact and approximate nearest neighbors For approximate and exact nearest neighbors in C++, one option is the Fast Library for Approximate Nearest Neighbors (FLANN) [START_REF] Muja | Flann, fast library for approximate nearest neighbors[END_REF]. This library allows one to define a custom point type which can be equipped with a label. It is then easy to compute probability estimates by counting each type of label among the nearest neighbor set returned. Note that this implementation allows incremental updates of the training set.