A Siffer
email: alban.siffer@irisa.fr

P A Fouque

A Termier

Are your data gathered? The Folding Test of Unimodality

Keywords: Mathematics of computing → Multivariate statistics, Nonparametric statistics, Exploratory data analysis, Unimodality test, Multivariate statistics

Understanding data distributions is one of the most fundamental research topic in data analysis. The literature provides a great deal of powerful statistical learning algorithms to gain knowledge on the underlying distribution given multivariate observations. We are likely to find out a dependence between features, the appearance of clusters or the presence of outliers. Before such deep investigations, we propose the folding test of unimodality. As a simple statistical description, it allows to detect whether data are gathered or not (unimodal or multimodal). To the best of our knowledge, this is the first multivariate and purely statistical unimodality test. It makes no distribution assumption and relies only on a straightforward p-value. Through real world data experiments, we show its relevance and how it could be useful for clustering.

INTRODUCTION

Given a multidimensional dataset of numerical attributes, an important question is to understand the "grouping behavior" of the data points. This question is traditionally answered by clustering algorithms. "Grouping behavior" being an ill-defined concept, clustering algorithms answer a more precise question : they determine groups (clusters) of data points that are similar, while being different from the data points of the other groups. Sometimes, one may not need such a detailed result : a more basic "grouping behavior" question can be to determine if all the data points make one single, coherent group or not.

In statistical terms, this is formulated as determining if the distribution of the data points is unimodal or not. Per se, unimodality test is a first information about the data, as it tells that the data points cannot be obviously separated into distinct groups. In a medical experiment, this would for example tell that the patients all reacted in a similar way, and that no group of patients deviated significantly from the others. Unimodality can then be seen as a building block for clustering algorithms: first, it can be used as a low-cost test to determine if running a clustering algorithm is necessary or not. Second, it can be used to improve clustering results, for example to help in parameter estimation.

Unimodality tests exist in the literature (for example the dip test [START_REF] Hartigan | The dip test of unimodality[END_REF] or the Silverman test [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF]), however they are restricted to one dimensional data. Up to now, there exists no purely statistical unimodality test able to tackle multi-dimensional data without distribution assumption. One may be tempted to use unidimensional unimodality tests on each dimension of multidimensional data (in a way similar to [START_REF] Maurus | Skinny-dip: Clustering in a Sea of Noise[END_REF]), however the simple example of Figure 1 shows the limits of such an approach. On the bottom left side, one can see that the bidimensional data is bimodal. However, its projection of the X (top) and Y (right) axis are perfectly unimodal, preventing to detect the multimodality of that dataset with unidimensional unimodality tests. Our main contribution is to propose the first statistical unimodality test designed for multidimensional data. Our approach relies on a novel descriptive statistics for (multidimensional) data, that provides a measure of its unimodality character (level of unimodality). The proposed approach keeps the same elegant properties as well-known unimodality tests like dip: it is independant from the distribution of the data and has only one parameter (a p-value threshold to tune the confidence in the results). We provide an algorithm for computing our test and show that it can be adapted to the streaming context through an incremental version.

The outline of the paper is as follows: in Section 2, we present related work on unimodality tests. We then describe our folding test of unimodality in Section 3, and its experimental validation in Section 4. Section 5 concludes the paper.

Unimodality tests. Many statistical tests have been proposed in the literature and Stoepker gave recently a nice survey of them [START_REF] Stoepker | Testing for multimodality[END_REF]. There are two kinds of tests: testing unimodality versus bimodality and testing unimodality versus multimodality. Once we know that a distribution is k-modal, many algorithms can be used to learn this distribution from a few samples [START_REF] Daskalakis | Learning k-Modal Distributions via Testing[END_REF][START_REF] Daskalakis | Testing k-Modal Distributions: Optimal Algorithms via Reductions[END_REF].

To test the "unimodality character" of a distribution, the main approaches aim at estimating the distribution. They include histograms, kernel density estimates and mixture models.

Silverman's Test is a parametric test that uses the bandwidth of the kernel density estimate to test multimodality. A sufficiently large bandwidth gives a smooth unimodal estimate. If we need a large amount of smoothing to find an unimodal estimate, this indicates multimodality. This test has several drawbacks [START_REF] Hall | On the calibration of Silverman's test for multimodality[END_REF][START_REF] Stoepker | Testing for multimodality[END_REF] even if we know that the input distribution is a mixture of gaussians. Some corrections have been proposed in [START_REF] Hall | On the calibration of Silverman's test for multimodality[END_REF] and a procedure to use non-gaussian kernels is given in [START_REF] Hall | Bump hunting with non-Gaussian kernels[END_REF].

The Excess Mass Test of Müller et al. [START_REF] Werner | Excess mass estimates and tests for multimodality[END_REF] tests unimodality versus bimodality. The idea is to estimate the excess mass of the modes which can be seen as the area between the density f (x) and a level L. The test compares the excess mass under unimodal and bimodal assumptions (through a statistics noted ∆). Actually, only a simple histogram is needed to estimate the density. Some properties about the distribution of ∆ have been discovered by Cheng and Hall [START_REF] Cheng | Calibrating the excess mass and dip tests of modality[END_REF] leading to a bootstrap procedure to compute p-values.

The Dip Test by Hartigan and Hartigan [START_REF] Hartigan | The dip test of unimodality[END_REF] measures the distance between the empirical cumulative distribution function (ecdf) and the set of the unimodal distributions U. The cdf of the set U are well characterized: they are convex on an interval ]-∞, x l ], then constant on [x l , x u ] and finally concave on [x u , +∞[. In short, the dip of an empirical cdf F corresponds roughly to the distance (with infinity norm) between F and the set U, i.e:

dip(F ) = d ∞ (F, U) = min U ∈U ∥F -U ∥ ∞ .
In [START_REF] Hartigan | The dip test of unimodality[END_REF], the uniform distribution is considered as the least favorable unimodal distribution. For an empirical cdf F n , the dip test compares dip(F n ) with the quantiles of dip(U n ) where U n is the ecdf of n iid uniform random variables (these quantiles can be pre-computed with Monte-Carlo simulations). The dip test reports a p-value p which is the probability to have an ecdf U n with dip(U n ) > dip(F n ). The common threshold α = 0.05 is used to make the decision: if p < α the ecdf F n is probably multimodal. This statistics and the test which results from it are rigorously detailed in [START_REF] Hartigan | The dip test of unimodality[END_REF] and a good analogy is given in [START_REF] Maurus | Skinny-dip: Clustering in a Sea of Noise[END_REF]. Even if the excess mass and the dip statistics are equivalent [START_REF] Cheng | Calibrating the excess mass and dip tests of modality[END_REF] the latter is the most commonly used statistics.

The main drawback of these tests is that they are restricted to real valued random variable (dimension 1). Moreover they need to estimate the distribution (density or cdf) and may require several passes over the data. Some multivariate extensions have been proposed like the RUNT test [START_REF] Hartigan | The runt test for multimodality[END_REF] or the MAP test [START_REF] Paul | The MAP test for multimodality[END_REF]. Unfortunately, these methods rely on the construction of several minimal spanning tree and the use of expensive optimization steps. So they are far more complex than our purely statistical approach.

Usage and related fields. Testing unimodality is not an end in itself. On the contrary, it may be used as a tool for related purposes. Particularly, some clustering algorithms use such approaches to find relevant groups of data. In [START_REF] Maurus | Skinny-dip: Clustering in a Sea of Noise[END_REF], the dip statistics is used to find clusters in noisy data (skinny-dip algorithm). The authors idea is to find all the bumps in the ecdf because they represent "high" density regions. This information comes precisely from the dip computation. As the dip cannot be computed in the multivariate case, the authors apply it on each dimension, with the issue we raised before.

In [START_REF] Kalogeratos | Dip-means: an incremental clustering method for estimating the number of clusters[END_REF], Kalogeratos and Likas present a combination of the dip statistics with the k-means algorithm (dip-means algorithm). More precisely, the dip test is performed on the pairwise distances of the points within a cluster. Indeed, a cluster is acceptable if it is unimodal (unimodality assumption). If k-means provides a cluster which is not unimodal, dip-means re-run the algorithm assuming an additional cluster (with an initial state carefully chosen).

Other approaches to estimate the k of k-means exist but they are likely to run the clustering algorithm several times with different numbers of cluster so as to optimize a predefined criterion [START_REF] Anil | Data clustering: 50 years beyond K-means[END_REF]. Examples of criteria are the Bayes Information Criterion (BIC), the Akaike Information Criterion (AIC), the Minimum Message Length (MML) [START_REF] Mario | Unsupervised learning of finite mixture models[END_REF], the Minimum Description Length (MDL) [START_REF] Mark | Model selection and the principle of minimum description length[END_REF] or the Gap Statistics [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF].

We may think that these tasks are similar with testing unimodality but they are slightly different. First, clustering is stronger because it provides the groups of data while unimodality testing decides whether data are gathered within a single group or not. However clustering needs more information: the number of clusters and/or some parameters to define the similarity between the observations. Thus, unimodality testing may be viewed as a macroscopic description of the data while clustering inspects it deeper.

Definition of unimodality. As presented above, unimodality is well-defined in dimension 1. In a few words, a distribution is unimodal iff the cdf is first convex, then flat, and finally concave. Unfortunately this definition does not generalize in higher dimensions. Thus, several unimodality notions have been developed but in nonequivalent ways. They are rigorously detailed in [START_REF] Dharmadhikari | Unimodality, convexity, and applications[END_REF].

In our work, unimodality refers to star unimodality (def 2.1 of [START_REF] Dharmadhikari | Unimodality, convexity, and applications[END_REF]). Simply speaking, a density f is unimodal about the mode m if f decreases with the distance to m.

CONTRIBUTION

This section presents our theoretical contribution. First we give the general intuition of our test of unimodality ( §3.1). The next part ( §3.2) deals with the mathematical formulation. It allows us to develop theoretical results ( §3.3), leading to our folding test of unimodality ( §3.4 and §3.5). Finally we deal with the incremental computation of the required statistics for the test ( §3.6).

Notations. In the next parts we will use the following notations: ∥ . ∥ denotes the euclidean norm. In the general case, X is a random vector of R d (d ∈ N * ), we note E[X ] ∈ R d its expected value. We assume that X is 3-integrable: it means that for all its components X i , E[X 3 i ] exists. We write Σ ∈ R d ×d its covariance matrix (we may use Var X in dimension 1). We assume that Σ is non-degenerated

(so it is invertible). If A ∈ R a and B ∈ R b are two random vectors, Cov (A, B) ∈ R a×b = E[(A -E[A])(B -E[B]) T ] denotes their co- variance. When X is a real valued random variable, we use also the common centered moments M k (X ) = E[(X -E (X )) k ].
The proposed algorithms may use the theoretical distribution of a random variable X , so we may use E[X ], Var X etc. but obviously they are valid on a given sample of n observations (drawn independently from the distribution of X ). In this case, the common estimators have to be used:

E[Z ] ≃µ Z = 1 n n i=1 Z i Var(Z ) ≃ 1 n n i=1 (Z i -µ Z ) 2 (dim. 1) Cov(A, B) ≃ 1 n n i=1 (A i -µ A ) • (B i -µ B ) T (so Σ(Z ) = Cov(Z, Z )) .
Thus, if f (X ) is a function applied on the theoretical random variable X , f (X 1 . . . , X n ) denotes its sample version (computed through the estimators given above).

General intuition

In this section, we present the general intuition of our approach aimed at testing whether a distribution is unimodal or not.

Univariate case.

Here we restrict the description of our approach to univariate distributions. The generalization will be made in the next paragraph.

Let us have a look at a bimodal density (figure 2a). In this case, the variance is likely to be "high" because the two modes make the data far from the expected value. Our idea is the following: if we fold up a mode to the other (with respect to the right pivot s * ), the resulting density (figure 2b) will have a far lower variance. Intuitively, this phenomenon will not appear for unimodal distributions (actually not with the same amplitude).

Let us sum up the approach: (1) find the right pivot s * , (2) fold up the distribution along s * , (3) compute the variance of the folded distribution and (4) compare it with the initial variance.

More formally, if we assume we get the pivot s * , the folding step is performed with the transformation X → |X -s * |, finally we will compute the folding ratio:

φ(X ) = Var |X -s * | Var X . folded variance initial variance
Higher dimensions. How can this approach be generalized? Actually, the transformation made on the distribution (the folding) is very similar except that the absolute value are replaced by the euclidean norm. Then we will consider Var ∥X -s * ∥ where X is a random vector of R d and s * ∈ R d is the pivot.

The figure 3 gives an empirical example of the folding mechanism in two dimensions. Given a trimodal distribution (figure 3a), and the right s * , the folding X → ∥X -s * ∥ turns the multi-modal bivariate distribution into a univariate distribution (figure 3b) which is likely to have a "low" variance. But obviously we have to mention about which reference this variance is low. The variance, which is (in dimension 1) the squared deviation of a random variable from its mean, is replaced by E ∥X -E[X ]∥ 2 in higher dimensions (it corresponds to the trace of the covariance matrix). Indeed, in the unimodal case this expected value will be much lower than in the multimodal case. Thus, the folding step will potentially have more impact in the latter. Finally, the folding ratio is generalized through:

φ(X ) = Var ∥X -s * ∥ E ∥X -E[X ]∥ 2 .
The pivot s * . Until this part, we have assumed to get the right pivot s * , i.e allowing the folding mechanism to significantly reduce the variance. Here we give more details about this pivot. Our goal is to find whether the variance may be significantly reduced by folding. So, the natural way is to find the pivot which cuts down the variance the most, which is (when it exists):

s * = argmin s ∈R d Var ∥X -s ∥.
This pivot is well-defined in dimension 1 but not in the general case (the minimum is not necessarily reached in higher dimensions) but we will try to get around this problem. An unimodal example in R 2 is given below (figure 4). In the figure 4a, we draw a sample from a multivariate normal distribution (which is unimodal). Moreover, we choose three different pivots s 1 , s 2 and s 3 and we plot the histogram of the folded observations ∥X -s 1 ∥, ∥X -s 2 ∥ and ∥X -s 3 ∥. In this case, we may notice that s 1 seems to reduce the variance the best (so it may be the best pivot s * ). Actually, in the unimodal case, the best pivot s * is likely to be close to the mode although in the multimodal case, it is likely to stand "between" the modes.

X -s 1 X -s 2 X -s 3 (b) Folded distribution

Formal approach

The previous paragraph introduced our main idea. Here, we develop it more formally. We introduce the function v X defined by

v X : R d → R s → Var ∥X -s ∥.
By definition, v X is lower-bounded by 0, so we can define the folding ratio as below

φ(X ) = inf s ∈R d Var ∥X -s ∥ E ∥X -E[X ]∥ 2 .
(1)

Thus, the computation of φ(X ) requires a minimization step, which is expensive. The ideal would be to have an expression of s * = argmin s ∈R d v X (s). Unfortunately, its existence is not clear because v X can be minimum to the infinity (we have noticed this divergence on some concrete and non absurd configurations).

In our approach, we aim to find a pivot even if the real one is not defined. Obviously this approximate pivot should produce a folding ratio which has intuitively the same properties as the theoretical one φ(X ). To circumvent the existence problem we rather use the following function:

v (2) X : R d → R s → Var ∥X -s ∥ 2 .
If we find the s which minimizes v

X , we are likely to have a good candidate to minimize v X . Moreover, the properties of v

(2) X are richer than the properties of v X . Lemma 3.1. For all s ∈ R d we have:

v (2) X (s) = 4 s T Σs -2 s T Cov ∥X ∥ 2 , X + Var ∥X ∥ 2 . Lemma 3.1 tells us that v (2)
X is a quadratic form of R d . Thus, thanks to the properties of the covariance matrix Σ, v

X is strictly convex and then it has a unique minimum, noted s * 2 . Theorem 3.1 gives its analytical expression.

Theorem 3.1. The function v (2)
X has a unique minimum s * 2 given by:

s * 2 = 1 2 Σ -1 Cov X, ∥X ∥ 2 . (2)
In particular, if X is a real random variable (dimension 1) then

s * 2 = E[X ] + 1 2 M 3 (X ) M 2 (X ) . (3)
With the previous result, we have always a pivot allowing to compute an approximate of the folding ratio

φ(X ) = Var ∥X -s * 2 ∥ E ∥X -E[X ]∥ 2 . ( 4
)
The previous expression leads to a straightforward algorithm to compute φ(X ) (algorithm 1). For simplicity reason, we will use "folding ratio" in the next paragraphs to express its approximate version.

Algorithm 1 Batch computation of φ(X )

Input: X Output: φ(X ) D ← E ∥X -E[X ]∥ 2 ▷ or D ← Tr(Σ) s * 2 ← 1 2 Σ -1 × Cov X,
O d 3 + n • d 2 s * 2 + n • d Var ∥X -s * 2 ∥

Unidimensional case study

In this paragraph, we treat the unidimensional case with some common distributions. We show that the approximate folding ratio φ(X ) is relevant to rank the distributions according to their "unimodal character". The expression of s * 2 is very convenient and allows us to calculate analytically its value (and also φ(X )) for some well-known distributions. We have summarized some of these results in table 1.

Distribution of X s * 2 φ(X ) Exponential E(λ) 2 λ 1 -4e -2 -4e -4 ≃ 0.385 Normal N(µ, σ 2 ) µ 1 - 2 π ≃ 0.363 Uniform U[a, b] a + b 2 1 4
Table 1: Analytical approximate folding ratio for some common unimodal distributions

The results above show two things: first we can notice that the folding ratios do not depend on the distribution parameters. This property is very interesting because it characterizes a whole distribution class through a single value. Second, if we merely compare the values, the distributions seem to be well-ranked according to their "unimodal character". Obviously, we concede this is not a well-defined notion (more a visual aspect) but we can sketch what we mean: a distribution is more unimodal when the relative density at its mode is higher.

Table 1 presents only unimodal distributions. Actually, analytical expressions for more complex distributions are more difficult to get. However, we may present a bimodal example: let δ ≥ 0 and X δ ∼ 1 2 N (0, 1) + 1 2 N (δ, 1). It means that X δ follows a bimodal normal distribution with a gap of δ between the two modes. In this case, the value of s * 2 is merely δ /2 but the folding ratio is given by:

φ(X δ ) = 1 - 1 δ 2 + 4 δ erf δ 2 √ 2 + 2 √ 2 √ π e -δ 2 8 2 ≃ δ ≥4 1 1 + δ 2 2
where "erf" is the classical error function. The figure 5a represents the function δ → φ(X δ ). We can notice that this function decreases when the gap between the modes δ increases. It naturally means that the folding step is more efficient when the modes are far (the initial variance is higher but the folding step cuts it down). We add the value of the uniform distribution to the figure (dotted line) and we highlight the gap making this Gaussian mixture "less unimodal" than the uniform density. The transition point corresponds to δ = 3.04, and the figure 5b presents the shape of the density for this gap. At this moment, we can observe that the distribution starts to become bimodal. It can be retorted because the density looks clearly bimodal however let us imagine we have only a sample of this distribution: there is a good chance the histogram will not be as accurate and the distribution will be considered as unimodal.

In this short study, we have shown that the approximate folding ratio is a relevant statistics to evaluate the "unimodality character". Furthermore, the uniform distribution seems to be the right reference to make a decision. In the next part, we will use this analysis to build our unimodality test.

The folding test of unimodality

In the previous parts, we have developed the approximate folding ratio: a relevant statistics to rank the distributions according to their unimodality character. Henceforth, we want to make a decision: is the distribution of X unimodal or not? To answer, we need to compare the folding ratio φ(X ) to a reference.

As claimed by Hartigan's work [START_REF] Hartigan | The dip test of unimodality[END_REF] and also observed in the paragraph above, the uniform distribution is likely to be this reference. Indeed, if we consider a sample of size n (from the uniform distribution), we need a larger n to assess its unimodality. More formally, Hartigan and Hartigan have shown that the dip statistics of the uniform distribution is asymptotically larger than other distributions ones, meaning that an empirical sample drawn from the uniform distribution need more observations to be considered as unimodal (worst case of unimodality).

In our work, we generalize this approach for all dimensions d ∈ N * . We consider the uniform distribution within the d-ball as the limit case of unimodality. We do not need to precise the radius of this d-ball because the folding ratio does not depend on it (see proposition 3.1). Proposition 3.1. Let R > 0 and B d (R) = x ∈ R d ∥x ∥ ≤ R the d-dimensional ball of radius R. The approximate folding ratio of the uniform random vector U d ∼ U (B d (R)) is:

φ d = φ(U d ) = 1 (d + 1) 2 .
The power of our method relies on its independence to distribution parameters, allowing to catch whole classes through a single value. As of now, we are able to build an unimodality test based on the comparison of folding ratios. Definition 3.1 (Folding test of unimodality). Let X be a 3integrable random vector of R d . We define the folding statistics Φ(X ) by

Φ(X ) = φ(X ) φ d = (1 + d) 2 φ(X ), (5)
leading to the folding test of unimodality: If Φ(X ) ≥ 1, the distribution of X is unimodal, if Φ(X ) < 1 it is multimodal.

Statistical significance

Theory. Obviously, we do not have the same level of confidence in the test whether we have 100 or 1 billion observations. Here we precise the decision bounds when the sample size is n.

Let us take a sample X 1 . . . , X n ∈ R d . We have seen how we can compute Φ(X 1 . . . , X n ). According to this value we can decide whether the distribution is unimodal or multimodal. Now, let us imagine that the sample U 1 . . . , U n is drawn from the uniform distribution: the test becomes very uncertain because the computed folding statistics would be close to 1. So we may understand that the test is more significant if we are far from the uniform case.

This classical statistical hypothesis testing problem leads us to provide p-values. For instance, let q > 0 and let us assume we have computed Φ(X 1 . . . , X n ) = 1 -q < 0, so the distribution is considered as multimodal. The probability to have a uniform example U 1 . . . , U n with a lower folding statistics is given by:

P (Φ(U 1 . . . , U n ) < 1 -q) = P (1 -Φ(U 1 . . . , U n ) > q) .
Conversely, we can imagine that we have computed Φ(X 1 . . . , X n ) = 1 +q. In this case, we can estimate the probability to have a uniform sample with a higher folding statistics:

P (Φ(U 1 . . . , U n ) > 1 + q) = P (Φ(U 1 . . . , U n ) -1 > q) .
Obviously we want these two probabilities to be low. It means that the p-value p = P (|Φ(U 1 . . . , U n ) -1| > q) must be as low as possible. In a nutshell, if we want a significant test (usually p ≤ 0.05), it implies to choose q high enough, making the uncertainty area 1 ± q wider. Conversely, if the test outputs a value Φ(X 1 . . . , X n ), it leads to a decision whose significance can be estimated by p (the lower it is, the more significant it will be).

Numerical quantiles. The ideal would be to know the distribution of Y d ,n = |Φ(U 1 . . . , U n ) -1|. However Y d ,n does not follow a common distribution, leading us to compute some quantiles with Monte-Carlo simulations. Table 2 presents the quantiles q according to d and n at a significance level p = 0.05 (10000 simulations for each tuple with Marsaglia's sampling method [START_REF] Marsaglia | Choosing a point from the surface of a sphere[END_REF]). For example, let us take a dataset of 1000 observations in 2 dimensions X 1 . . . , X n (n = 1000, d = 2). The distribution is considered significantly unimodal (at level p = 0.05) if Φ(X 1 . . . , X n ) ≥ 1.09 (and significantly multimodal if it is lower than 0.91). Given a set of similar tables for different values of p, we can infer the significance of any test output Φ.

Obviously we notice that the higher is n, the tighter are the bounds but the dimension increases the uncertainty. Unfortunately, we cannot provide reliable quantiles for higher dimensions with these low numbers of points because of the curse of dimensionality in the sampling. They can easily be computed in specific cases but the number of simulations must be high enough.

Incremental computation

In the paragraph 3.2 we proposed a batch (or offline) version for the computation of φ(X 1 . . . , X n ). However, in the streaming context, we can accelerate it with an incremental computation of s * 2 . As seen in the Section 3.2, we have an analytical expression for the pivot:

s * 2 = 1 2 Σ -1 Cov(X, ∥X ∥ 2 )
. Furthermore this expression can easily be computed incrementally or even updated over a sliding window. This is a very important property if we work on streaming data. The algorithm 2 details how to perform an update.

As an input we have a new incoming observation X new and the square of its norm R new . The idea is to compute their respective basic moments (cumulative sums S X and S R ) and their co-moment (line 4) to finally compute an estimate of Cov(X, ∥X ∥ 2 ) (line 7). About the inverse of the covariance matrix, we use the Sherman-Morrison formula [START_REF] Sherman | Adjustment of an inverse matrix corresponding to a change in one element of a given matrix[END_REF] to compute the inverse of the co-moment of X (line 5) and then Σ -1 (line 6). We recall:

SM(A -1 , u, v) = A + u × v T -1 = A -1 - A -1 × u × v T × A -1 1 + v T × A -1 × u .
Here we present an update, so we have to deal with the initialization. Basically, the variables n, S X , S R , V X ,R and C are set to 0. Unfortunately, the initialization of V inv X and Σ inv is not as simple because it requires several observations X i (at the beginning V inv X and Σ inv are singular matrices). In practice, we start from a small batch of data X init . Therefore, we can compute V X ← X T init × X init and then V inv X ← V -1 X (so Σ inv ). Finally, this update can easily be extended to sliding window requiring to store the quantities X i and

Algorithm 2 Incremental computation of s * 2 Input: X new , R new = ∥X new ∥ 2 Output: s * 2 1: n ← n + 1 2: S X ← S X + X new 3: S R ← S R + R new 4: V X ,R ← V X ,R + X new × R new 5: V inv X ← SM(V inv X , X new , X new ) 6: Σ inv ← n × SM(V inv X , -S X , 1 n S X ) ▷ Σ -1 7: C ← 1 n V X ,R -1 n S X × 1 n S R ▷ Cov(X, ∥X ∥ 2 ) 8: s * 2 ← 1 2 × Σ inv × C R i
, but here, we prefer not to overload the reader with additional formulas. Is it possible to go further by computing φ(X ) incrementally too? Unfortunately, this is not as straightforward. There is no problem for the denominator E ∥X -E[X ]∥ 2 because this is the trace of Σ. The main issue comes from Var ∥X -s * 2 ∥: a variance can easily be computed incrementally (like in the algorithm 2) but s * 2 changes at each iteration, so an update of all the distances ∥X i -s * 2 ∥ is needed. In the paragraph 3.2, we discussed about the complexity of the batch computation. In table 3, we compare its complexity with the incremental computation in 3 cases:

• single computation: only φ(X 1 , . . . X n ) • cumulative: φ(X 1 ), φ(X 1 , X 2 ) . . . φ(X 1 , . . . X n ) • sliding window: φ(X 1 , . . . X 1+w ), . . . φ(X n , . . . X n+w ) Computation(s) Batch Incremental Single d 3 + n • d 2 (+ n • d) Cumulative n • d 3 + n 2 • d 2 d 3 + n • d 2 + n 2 • d Sliding window n • d 3 + n • w • d 2 d 3 + n • d 2 + n • w • d
Table 3: Complexity of each method (in big O notation)

Without delving into the details, we may notice two phenomenons. First, the incremental computation of s * 2 decreases the dependency on the dimension d (actually the update of the matrix costs O(d 2 ) instead of O(d 3 ) in the batch version). Second, as the computation of Var ∥X -s * 2 ∥ cannot be done efficiently, the computation of the folding statistics cannot be cut down any more.

EXPERIMENTS

In this section, we highlight the relevance of the folding test of unimodality. We experiment it on real world multidimensional data to emphasize its correctness and its practical uses. Particularly, we show that it provides a paramount statistical description necessary for data analysis. We make our python3 code available for reviewers 1 .

Should I try to cluster Pokémon?

In this section, we will show that our test is able to avoid a useless clustering step. Particularly, we analyze the Pokemon Stats Dataset from kaggle 2 . This dataset gathers statistics of the Pokémon until the 6 th generation: it includes 21 variables per each of the 721 Pokémons. In this short experiment, we keep only the 6 basic fight statistics: Attack, Defense, HealthPoints, SpecialAttack, SpecialDefense and Speed.

A priori, in such a 6-dimensional space, this is not easy to claim whether some clusters arise. Either all the features pairs could be plotted, or a clustering algorithm could be run. Through the first method (figure 6), the distributions of the features and their pairs lead us to think that the whole distribution is rather unimodal.

When we compute the folding statistics to this dataset, we get Φ(X ) = 5.04. As guessed above, this result claims that the distribution is unimodal. To show that our approach is relevant, we compare the result of our test with the output of some well-known clustering algorithms. Unfortunately, most of them needs to set precisely the number of clusters, so we only test the approaches having some procedures to find it.

First we try to model the data with gaussian mixtures. The number of clusters is found by minimizing either the Akaike's or the Bayesian Information Criterion (AIC, BIC). Then we use the gap statistics [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF] of Tibshirani et al., which is a common statistics to estimate the number of clusters. We test also the Mean Shift algorithm [START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF] with different quantiles q ∈ [0.3, 0.8]. This parameter tunes the "similarity" of the points (it sets the bandwidth of the underlying gaussian kernel). Finally, we run the Affinity Propagation algorithm [START_REF] Brendan | Clustering by passing messages between data points[END_REF] with the euclidean distance as similarity. For these two last approaches, we show only "stable" results, so we do not present neither absurd parameter choices (the mean shift algorithm with q → 1 would make all the points similar) nor absurd outputs of some algorithms (the affinity propagation algorithm with the gaussian similarity did not output less than 49 clusters). Finally we compute the average silhouette score [START_REF] Peter | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF] to estimate the clustering quality (the closer to 1, the better).

The results are given in the table 4. We can notice that there is no consensus among the different methods: they do not output the same number of clusters. At the very least, this indicates that the dataset does not contain obvious clusters. Now regarding the returned clusters, the average silhouette scores are mostly low, even close to zero, meaning that these clusters are not well separated. The highest silhouette score is reached by Mean Shift, however the algorithm returned 3 clusters of size 717-3-1, so in practice it found only one cluster. These results are in favor of the unimodality of the data, agreeing with the result of our folding test of unimodality.

Before running clustering algorithms, our simple test can thus be a crucial step to decide whether clustering is relevant or not.

Stock market behaviors

The aim of this section is twofold. Firstly, we analyze the behavior of the folding statistics over a sliding window. Secondly, we show how the folding test of unimodality can help the clustering parametrization step.

Catching the behavior jumps. In this experiment we analyze a part of the DJIA 30 Stock Time Series dataset hosted on kaggle 3 . In particular we study the historical stock data from NYSE (between 2006-01-03 to 2017-29-12) of four major companies: Cisco (CSCO), Intel (INTC), IBM and Microsoft (MSFT). We consider the highest price daily reached, so we get a dataset with n = 3018 observations in dimension 4 (number of companies).

The figure 7 shows the evolution of the stock prices for these four companies. We may notice two different behaviors: IBM stocks vary more than Cisco, Intel and Microsoft ones. Now, we compute the folding statistics Φ over a sliding window of size 650 (about two years and half, because there is no data on Saturday and Sunday). It means that at every iteration, we take a snapshot of 650 observations in dimension 4, and we compute Φ on it. The results are presented in the figure 8. We add significance bounds at level 0.05 (q ≃ 0.15). One can roughly see three "moments of unimodality" (around iterations 500, 1500 and 1800).

What does it mean? If the distribution at these moments is unimodal, it means that the companies stocks are individually quite the same: their behavior is stationary within the window. At the To highlight this phenomenon, we show two windows on the time series (figure 9). One corresponding to the minimum of Φ (maximum of multimodality) and the other one corresponding to its maximum (maximum of unimodality). Indeed, we may notice that the most multimodal window is when the IBM stocks jump a lot unlike the others although the most unimodal window occurs when the four stocks are the most stationary.

If we look deeper at the scatter plots within these two windows (figure 10), we can confirm these behaviors (to ease the plot we have overlaid the scatter-plots IBM vs MSFT, IBM vs CSCO and IBM vs INTC). On figure 10a, we can observe about three clusters whereas no clusters really arises on the figure 10b. However we clearly see three little groups of data on this latter figure. In fact, they are not statistically substantial (not real modes). Running time. We compare the running time of two methods: computing the statistics on the complete window each time (batch computation), or using the incremental approach presented in section 3.6. The complete computation (2400 iterations) took 5.5 seconds with batch computation, and 4.5 seconds with the incremental algorithm (around 18% acceleration). Note that the variance computation part is identical in both methods: if we only compare the computation of s * 2 , the times are 3.4 and 0.9 seconds, giving a 73% decrease on running time. We can conclude that even on this relatively simple data (low dimension, few iterations), the incremental computations has a significant impact on running time.

Helping clustering algorithms. We have mentioned in the Pokémon experiment (section 4.1) that unimodality testing provides a different piece of information from clustering. In that case, it was relevant to avoid over-clustering. Here, we use it as a parametrization support for clustering.

Let us consider the DBSCAN algorithm [START_REF] Ester | A densitybased algorithm for discovering clusters in large spatial databases with noise[END_REF]. It requires two parameters: the neighborhood radius r and the minimum requested numbers of neighbors to be considered as a core point m s . Our idea is to tune the radius parameter so as to make DBSCAN output "1 cluster" only when Φ > 1. Obviously, we cannot ensure the quality of the clustering in the other regions however with this criterion we can narrow down the parameter research domain.

An example with m s = 80 (12.5% of the observations in the window) is presented on the figure 11. We plot the folding statistics (top) and the number of clusters output by DBSCAN with different neighborhood radius (r = 6, 7, 8, 9 and 10). Moreover we add all the unimodal areas (where Φ > 1 + q) and the first multimodal area (where Φ < 1 -q).

Let us analyze the unimodal area around iteration 1500. All the instances of DBSCAN output nearly a single cluster. But around it. 1800 (unimodal area too), they all return 2 clusters. In the first unimodal zone (around it. 500) DBSCAN(r = 7) is the only one not returning one cluster. Eventually, if we look at the first multimodal area (before it. 330), only the instances DBSCAN(r = 7) and DBSCAN(r = 8) output several clusters. We may extend the analysis to all the regions, making the instance DBSCAN(r = 8) the more faithful to the folding test of unimodality.

Figure 1 :

 1 Figure 1: Unidimensional test fails in 2D

Figure 2 :

 2 Figure 2: Folding mechanism for univariate distribution

Figure 3 :

 3 Figure 3: Folding mechanism in dimension 2

Figure 4 :

 4 Figure 4: Impact of the pivot location

∥X ∥ 2 return

 2 Var ∥X -s * 2 ∥/D We give some details about complexity. Let us consider a dataset of n observations in R d . The covariance matrix Σ needs O(n • d 2 ) operations and additional O(d 3 ) operations are required for its inversion. Computing the norm of the observations and the covariance Cov X, ∥X ∥ 2 costs O(n • d) operations. Finally, anoter O(d 2 ) operations are required for s * 2 and O(n • d) for the variance Var ∥X -s * 2 ∥. In a nutshell, the complexity of the batch computation is linear for the number of observation n :

  Distribution of X δ when δ = 3.04

Figure 5 :

 5 Figure 5: Analysis of a bimodal normal distribution

Figure 6 :

 6 Figure 6: Pairs-plot of the Pokémon six fight characteristics

Figure 7 :

 7 Figure 7: Stock evolution of the four companies

Figure 8 :

 8 Figure 8: Folding statistics within 650 days time windows moments when the distribution is multimodal, some stock prices are jumping.

Figure 9 :

 9 Figure 9: Extreme time windows

Figure 10 :

 10 Figure 10: Scatter-plots IBM vs Others

Table 2 :

 2 Quantiles q according to d and n at significance level p = 0.05

	n/d	1	2	3	4	5
	100	0.22 0.28 0.33 0.35 0.38
	200	0.15 0.2 0.23 0.25 0.27
	500	0.1 0.13 0.15 0.16 0.17
	1000 0.07 0.09 0.1 0.11 0.12
	2000 0.05 0.06 0.07 0.08 0.09
	5000 0.03 0.04 0.05 0.05 0.05
	10000 0.02 0.03 0.03 0.04 0.04
	20000 0.02 0.02 0.02 0.03 0.03

Table 4 :

 4 Number of clusters and average silhouette score according to the clustering method

	2 https://www.kaggle.com/alopez247/pokemon

git clone https://scm.gforge.inria.fr/anonscm/git/sharedcode/sharedcode.git

https://www.kaggle.com/szrlee/stock-time-series-20050101-to-20171231

We insist on the non-completeness of our test. We cannot claim that DBSCAN gives the correct clusters in the multimodal regions. Moreover, we know that other phenomenons should be taken into account for clustering like the influence of the second parameter and/or data normalization. We merely show that our lightweight test is able to provide some paramount information about data. The latter may be useful for a potential clustering stage.

CONCLUSION

This paper introduces the folding test of unimodality. To the best of our knowledge, this is the first time a multivariate and purely statistical unimodality test has been proposed. The test is based on a new folding statistics, which provides a score related to the "level of unimodality" of a distribution. This statistics does not depend on the parameters of the input distribution. The only parameter of the unimodality test is a natural p-value giving the desired significance level of the result.

Among the perspectives opened by this work, we envision to exploit our folding statistics to discover k values for the k-means algorithm. This was explored before in the dip-means approach, by using unidimensional unimodality tests on all pairwise distances. Using our folding unimodality test allows to immediately validate each cluster of k-means in its multi-dimensional space (unimodality assumption). This significantly reduces the complexity compared to dip-means.