Keywords: Classification, Mélanges gaussiens, Sélection de variables, Critère pénalisé, Heuristique de pente, Classification de courbes

 pour sélectionner le nombre de composantes du mélange et le sous-ensemble des variables significatives pour la classification. Ce critère dépend de constantes multiplicatives inconnues qui doivent être évaluées en pratique. Une méthode heuristique dite "de la pente" est proposée et expérimentée pour résoudre ce problème. Des exemples numériques sur données simulées, un exemple de classification de courbe et une application génomique mettent en évidence l'intérêt de cette procédure.

Introduction

Model-based clustering methods consist of modelling each cluster with a parametric distribution and considering the mixture of these distributions to describe the whole dataset. They provide a rigorous framework to assess the number of mixture components and to take into account the variable roles.

Currently, cluster analysis is more and more concerned with large datasets where observations are described by many variables. This large number of predictor variables could be beneficial to data clustering. Nevertheless, the useful information for clustering can be contained into a variable subset and some of the variables can be useless or even harmful to choose a reasonable clustering structure. Several authors suggest variable selection methods for Gaussian mixture clustering which is the most widely used mixture model for clustering multivariate continuous datasets. These methods are "wrapper" methods since they are included into the clustering process. [START_REF] Law | Simultaneous feature selection and clustering using mixture models[END_REF] have introduced the feature saliency concept. Regardless of cluster membership, relevant variables are assumed to be independent of the irrelevant variables which are supposed to have the same distribution. [START_REF] Raftery | Variable selection for model-based clustering[END_REF] recast variable selection for clustering into a global model selection problem. Irrelevant variables are explained by all the relevant clustering variables according to a linear regression. The comparison between two nested variable subsets is performed using Bayes factor. A variation of this method is proposed in [START_REF] Maugis | Variable selection for clustering with Gaussian mixture models[END_REF] where irrelevant variables can only depend on a relevant clustering variable subset and variables can have different sizes (block variables). Since all these methods are based on a variable selection procedure included into the clustering process, they do not impose specific constraints on Gaussian mixture forms. On the contrary, [START_REF] Bouveyron | High-dimensional data clustering[END_REF] consider a suitable Gaussian mixture family to take into account that data live in low-dimensional subspaces hidden in the original space. However, since this dimension reduction is based on principal components, it is difficult to deduce from this approach an interpretation of the variable roles.

In this paper, a new variable selection method for clustering is proposed. It recasts variable selection and clustering problems into a model selection problem in a density estimation framework. Suppose that we observe a sample from an unknown probability distribution with density s. A specific collection of models is defined: Each model S (K,v) corresponds to a particular clustering situation where the cluster number is K and v is the relevant clustering variable subset. A density t belonging to S (K,v) is decomposed into a Gaussian mixture density with K components on the relevant clustering variable subset v and a multidimensional Gaussian density on the other variables. Definitions of models S (K,v) are precised in Section 2. The problem can be formulated as the choice of a model among the collection since this choice automatically leads to a data clustering and a variable selection. Thus, a data-driven criterion is needed to select the "best" model among the model collection. This article is the companion of [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF] where a penalized likelihood criterion is proposed. The results obtained in this previous paper allow to specify the general shapes of the criterion penalties but additional work is necessary to use these penalties in prac-tice. The aim of this paper is to describe the practical use of these theoretical results. The so-called "slope heuristics" method is applied to calibrate penalties on the data.

First, our methodology is applied to a curve clustering problem. Curve clustering deals with the problem of identifying homogeneous groups in a set of functional data. This situation occurs in many areas of sciences, for instance in genetics, neuroscience, economics and engineering. Many methods for making curve clustering are based on different versions of the k-means algorithm. A widely used technique consists of finding a convenient projection of the functional data into a finite dimensional subspace, and next of applying a k-means procedure on the finite dimensional data obtained. In this context, B-spline bases are currently used, see for instance [START_REF] Abraham | Unsupervised curve clustering using B-splines[END_REF] and [START_REF] García-Escudero | A proposal for robust curve clustering[END_REF]. An other approach proposed by [START_REF] Tarpey | Clustering functional data[END_REF] is to adapt the kmeans algorithm for functional spaces. With a different point of view, [START_REF] Ma | A data-driven clustering method for time course gene expression data[END_REF] use mixture models on B-splines coefficients, as in the works of [START_REF] James | Clustering for sparsely sampled functional data[END_REF] for sparsely sampled functional data. In most of the cited works, each curve is described with a coefficient vector and in practice, the number of these coefficients can be of the same order as the curve number. This high dimensional context makes our method desirable to solve curve clustering problems. We illustrate the application of our method for curve clustering on the study of an oil production curve sample.

Next, we show that our method can be applied to the transcriptome data clustering in a particular context. During these last years, biologists are interested in determining biological functions of genes. In this aim, clustering methods such as hierarchical clustering or k-means algorithm are commonly applied to find clusters of co-expressed genes (see for instance [START_REF] Sharan | Cluster analysis and its applications to gene expression data[END_REF]Jiang et al., 2004, and references therein). Since the experiment number increases in available transcriptome datasets, variable selection is more and more considered in order to improve the clustering and its interpretation to help biologists. In the model-based clustering context, [START_REF] Maugis | Variable selection for clustering with Gaussian mixture models[END_REF] apply their variable selection method for a transcriptome dataset analysis for instance.

The paper is organized as follows: Section 2 presents the collections of Gaussian mixture model used in this paper. In Section 3, the general framework of model selection for density estimation based on Kullback-Leibler contrast is first recalled. Next, we present a penalized criterion to select a model into the considered collections of Gaussian mixture models, and we recall the results obtained on this criterion in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF]. Section 4 is devoted to the description of slope heuristics and its practical use in our specific context. Simulations and applications for a curve clustering problem and for a genomics study are presented in Section 5.

2 Gaussian mixture models

Multivariate Gaussian models and clustering

Centered observations y = (y 1 , . . . , y n ) with y i ∈ R Q are assumed to be a sample from a probability distribution with unknown density s. This target distribution is proposed to be INRIA estimated by a finite mixture model in a clustering purpose. Model-based clustering consists of assuming that the data come from a source with several subpopulations. Each subpopulation is modelled separately and the overall population is a mixture of these subpopulations. The resulting model is a finite mixture model. When the data are multivariate continuous observations, the component parameterized density is usually a multidimensional Gaussian density. The general form of a Gaussian mixture model with K components is

K k=1 p k Φ(•|η k , Λ k )
where the p k 's are the mixing proportions (∀k = 1, . . . , K, 0 < p k < 1 and K k=1 p k = 1) and Φ(•|η k , Λ k ) denotes the Q-dimensional Gaussian density with mean η k and variance matrix Λ k . The mixture model is an incomplete data structure model: The complete data are ((y 1 , z 1 ), . . . , (y n , z n )) where the missing data are z = (z 1 , . . . , z n ) with z i = (z i1 , . . . , z iK ) such that z ik = 1 iff y i arises from component k. The z's define an ideal clustering of the data y associated to the mixture model. After an estimation of the parameter vector thanks to the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], a data clustering is deduced using the maximum a posteriori principle (MAP rule):

ẑik = 1 if pk Φ(y i |η k , Λk ) > pl Φ(y i |η l , Λl ), ∀l = k 0 otherwise.

Definitions of the considered Gaussian mixture models

Currently, statistics deals with problems where individuals are explained by many variables.

In principle, the more information we have about each individual, the better a clustering method is expected to perform. Nevertheless, some variables can be useless or even harmful to obtain a good clustering of data. Thus, it is important to take into account the variable role in the clustering process. In this goal, we propose to consider Gaussian mixtures with a specific form, based on the following idea. On irrelevant variables, individuals have an homogenous behavior around the null mean (centered data) allowing not to distinguish different subpopulations. Hence data are assumed to have a spherical Gaussian joint law with null mean vector on these variables. On the contrary, the different component mean vectors are free on relevant clustering variables. Variance matrices restricted on relevant variables are either taken completely free or are chosen in a specified set of definite positive matrices depending on the considered situation. This idea is now formalized. Let V be the set of the nonempty subsets of {1, . . . , Q}. A Gaussian mixture family is characterized by its number of mixture components K ∈ N and its relevant variable index subset v ∈ V whose cardinal is denoted α. In the sequel, the set of index couples (K, v) is M = N × V. Consider the decomposition of a vector x ∈ R Q into its restriction on relevant variables x [v] = (x j1 , . . . , x jα ) and its restriction on irrelevant variables x [v c ] = (x l1 , . . . , x l Q-α ) where v = {j 1 , . . . , j α } and v c = {l 1 . . . , l Q-α } = {1, . . . , Q}\v. On relevant clustering variables, a Gaussian mixture f is chosen among the following family

L (K,α) = w ∈ R α → K k=1 p k Φ(w|µ k , Σ k ); ∀k, µ k ∈ [-a, a] α , (Σ 1 , . . . , Σ K ) ∈ D + (K,α) 0 < p k < 1, K k=1 p k = 1
where 0 < a and D + (K,α) denotes a family of K-uples of α × α symmetric definite positive matrices which is related to the Gaussian mixture shape specified hereafter. On irrelevant variables, a spherical Gaussian density belonging to the following family is considered

G (α) = u ∈ R Q-α → Φ(u|0, ω 2 I Q-α ); ω 2 ∈ [λ m , λ M ] .
Thus, the family of Gaussian mixture associated to (K, v) ∈ M is defined by

S (K,v) = x ∈ R Q → f (x [v] ) g(x [v c ] ); f ∈ L (K,α) , g ∈ G (α) .
(1)

The dimension of the model S (K,v) is denoted D(K, v) and corresponds to the free parameter number of Gaussian mixtures in this model. It only depends on the number of components K, the Gaussian mixture shape and the number of relevant variables α.

In this paper, four collections of Gaussian mixtures are considered but the same notation S (K,v) is used for the four model collections to make easier the reading of this article. The Gaussian mixture notation for those collections is taken from [START_REF] Biernacki | Model-based cluster and discriminant analysis with the MIXMOD software[END_REF].

• For the [L k B k ] collection, the variance matrices are assumed to be diagonal and free.

Thus, the variance matrices have the following form

Σ k = diag(σ 2 k1 , . . . , σ 2 kα )
where the eigenvalues σ 2 kj are assumed to be in the interval [λ m , λ M ]. The associated dimension of model S (K,v) is equal to D(K, v) = K(2α + 1).

• For the [L k C k ] collection, the variance matrices are assumed to be totally free. Thus, the variance matrices are α×α positive definite matrices whose eigenvalues are assumed to belong to the interval [λ m , λ M ]. The associated model dimension is

D(K, v) = K[1 + α + α(α+1)

2

].

• For the [LB k ] collection, the variance matrices are assumed to be diagonal and to have the same volume i.e.

∀k = k , |Σ k | 1 α = |Σ k | 1 α .
The variance matrices are decomposed into Σ k = βB k where the common volume β belongs to [β m , β M ] and B k is a diagonal matrix with a determinant 1 and with diagonal coefficients in the interval [λ m , λ M ]. Here, the model dimension is equal to D(K, v) = 2Kα + 1.

• For the [LC] collection, the variance matrices are all equal to a free positive definite matrix Σ whose eigenvalues are assumed to be in the interval

[λ m , λ M ]. The model dimension is D(K, v) = K(1 + α) + α(α+1)

2

.

INRIA

In this paper, for each of the four possible model collections, the variables can be assumed to be ordered or not ordered. If variables are ordered, the relevant variable subset is v = {1, . . . , α} and can be assimilated to its cardinal α. Moreover, note that a density of S (K,v) can be written as a global Gaussian mixture with mean vectors η k = (µ k , 0, . . . , 0) and block-diagonal variance matrices Λ k with diagonal-blocks Σ k and ω 2 I Q-α . Consequently, a data clustering can be deduced from the MAP rule given in Section 2.1.

These Gaussian mixture families allow to recast clustering and variable selection problems in a global model selection problem. A criterion is now required to select the best model according to the dataset.

A new penalized likelihood criterion 3.1 Model selection principle

Density estimation deals with the problem of estimating the unknown distribution of a sample y 1 , . . . , y n . In many cases, it is not obvious to choose a model of adequate dimension. For instance, a model with few parameters tends to be efficiently estimated whereas it could be far from the true distribution. On the opposite situation, a more complex model easily fits data but estimates have larger variances. The aim of model selection is to construct data-driven criterion to select a model of proper dimension among a given list. A general theory on this topic, with a non asymptotic approach is proposed in the works of Birgé and Massart (see Massart, 2007, for an overview). This model selection principle is now described in our setting.

Let S be the set of all densities with respect to the Lebesgue measure on R Q . The contrast γ(t, •) = -ln{t(•)} is considered, leading to the maximum likelihood criterion. The corresponding loss function is the Kullback-Leibler information. It is defined for two densities s and t in S by KL(s, t) = ln s(x) t(x) s(x) dx if sdx is absolutely continuous with respect to tdx and +∞ otherwise. Noticing that s is a minimizer of the Kullback-Leibler information over S, s is also a minimizer over S of the expectation of the empirical contrast, defined by

γ n (t) = - 1 n n i=1 ln {t(y i )} .
A countable collection of models {S (K,v) } (K,v)∈M is considered and let ŝ(K,v) be a minimizer of the empirical contrast γ n over the model S (K,v) . Substituting the empirical criterion γ n to its expectation and minimizing γ n on S (K,v) are expected to lead to a sensible estimator of s, at least if s belongs (or is close enough) to model S (K,v) . The model we want to select is the one presenting the smallest risk

(K , v ) = argmin (K,v)∈M E[KL(s, ŝ(K,v) )].
However, the function ŝ(K ,v ) , called oracle, is unknown since it depends on the true density s. Thus, the aim is to find a data-driven criterion to select an estimator such that its risk is as close as possible as the oracle risk. The model selection via penalization procedure consists of considering some penalized criterion

crit(K, v) = γ n (ŝ (K,v) ) + pen(K, v) (2)
where pen is a penalty function pen :

(K, v) ∈ M → pen(K, v) ∈ R + . Then the se- lected model ( K, v
) is a minimizer of the penalized criterion (2) and the associated selected estimator is ŝ( K,v) . The final purpose of a non asymptotic approach is to obtain a penalty function providing an oracle inequality. This oracle inequality allows to compare the risk of the penalized maximum likelihood estimator (MLE) ŝ( K,v) with the benchmark inf

(K,v)∈M E[KL(s, ŝ(K,v) )],
for a fixed number n of observations.

Theoretical results

From a theoretical point of view, the problem of defining a convenient penalized likelihood criterion for our specific Gaussian mixture model collections has been treated in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF]. This work has been made possible thanks to the use of a general model selection theorem for MLE, proposed by [START_REF] Massart | Concentration inequalities and model selection[END_REF]. The application of this theorem requires the control of bracketing entropy of Gaussian mixture families. It allows to obtain penalty function forms and associated oracle inequalities.

In the sequel, the norm √ f -√ g 2 between two nonnegative functions f and g of L 1

is denoted d H (f, g). We note that if f and g be two densities with respect to the Lebesgue measure on R Q , d 2 H (f, g) is the squared Hellinger distance between f and g

d 2 H (f, g) = R Q f (x) -g(x) 2 dx.
In the following, d H (f, g) is improperly called Hellinger distance even if f and g are not density functions. The following theorem summarizes the theoretical results established in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF]. These results are valid for the four Gaussian mixture models at hand. They are proved for the [L k B k ] shape and the [L k C k ] shape in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF] and additional proofs for the two other shapes [LB k ] and [LC] are available in the Appendices.

Theorem 1 [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF]). For the four Gaussian mixture collections, 1. If the variables are ordered, there exists two absolute constants κ and C such that, if

pen(K, v) ≥ κ D(K, v) n 2 A + ln 1 1 ∧ D(K,v) n A + 1 INRIA then the model ( K, v) minimizing crit(K, v) = γ n (ŝ (K,v) ) + pen(K, v) on M exists and E d 2 H (s, ŝ( K,v) ) ≤ C inf (K,v)∈M [KL(s, S (K,v) ) + pen(K, v)] + 1 n .
(3)

2. If the variables are not ordered, there exists two absolute constants κ and C such that, if

pen(K, v) ≥ κ D(K, v) n 2 A + ln 1 1 ∧ D(K,v) n A + 1 2 ln 8 exp(1)Q (D(K, v) -1) ∧ (2Q -1) then the model ( K, v) minimizing crit(K, v) = γ n (ŝ (K,v) ) + pen(K, v) on M exists and E d 2 H (s, ŝ( K,v) ) ≤ C inf (K,v)∈M [KL(s, S (K,v) ) + pen(K, v)] + 2 n . (4)
In the four cases, A is a function of parameters λ m , λ M , a, Q and also β m and β M for the

[LB k ] shape such that A = O( √ ln Q) as Q goes to infinity.
The penalty functions take into account the model complexity through D(K, v) and the richness of model family. Actually the number of models with the same dimension is larger in the non-ordered variable case, the associated penalty functions have an additional logarithm term, depending on the dimension.

The other logarithm term, common to both cases, is probably not necessary to define efficient penalties. As explain in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF], the reason for that is certainly that the general model selection theorem for MLE is stated in a local version whereas we only manage to apply the global version in our framework. Logarithm terms are not detected in practice as shown in Section 5.1 and thus only the preponderant term in D(K,v) n is retained in the penalty form.

Contrary to classical criteria for which Q is fixed and n tends to infinity, our result allows to study cases for which Q increases with n. For specific clustering problems where the number of variables Q is of the order of n or even larger than n, the oracle inequality is still significant.

The point we want to stress here is that Theorem 1 gives the general form of penalty functions but is not totally explicit since the results depend on absolute unknown constants and mixture parameters are not bounded in practice. Consequently a method is to be applied to calibrate the penalty function for a practical use of these results.

Slope heuristics

The aim of this paper is to show how the theoretical results of Section 3.2 can be applied in practice. Since the lower bounds on penalty functions in Theorem 1 are defined up to an unknown multiplicative constant, this theorem does not provide directly a usable model selection criterion. Last years, some efforts have been paid to overcome such a difficulty. [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] propose a practical method based on a mixture of theoretical and heuristic ideas for defining efficient penalty functions from the data. This heuristics is only proved in [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] in the framework of Gaussian regression with a homoscedastic fixed design and more recently generalized by [START_REF] Arlot | Slope heuristics for heteroscedastic regression on a random design[END_REF] in the heteroscedastic random-design case. Nevertheless applications of this method are developed in many other frameworks: For instance, in multiple change points detection by [START_REF] Lebarbier | Detecting multiple change-points in the mean of Gaussian process by model selection[END_REF], in genomics applications by [START_REF] Villers | Tests et sélection de modèles pour l'analyse de données protéomiques et transcriptomiques[END_REF] and in Gaussian Markov random fields by [START_REF] Verzelen | Adaptive estimation of regular Gaussian Markov random fields[END_REF]. This section first describes the main ideas of this heuristics, called the "slope heuristics", and next details its practical use in our framework.

Rationale for the slope heuristics

In many situations, the considered model collection contains several models with the same dimension. In order to penalize each model of dimension D in the same way, a new collection (S D ) D∈D is considered such that S D is the union of all the models S (K,v) having the same dimension D(K, v) = D. It is recalled that γ and γ n are the Kullback-Leibler contrast and its associated empirical contrast respectively (see Section 3.1). Moreover s D and ŝD denote a minimizer of KL(s, •) and γ n (•) on S D respectively.

As for criteria due to [START_REF] Mallows | Some comments on Cp[END_REF] and [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF][START_REF] Akaike | A new look at the statistical model identification[END_REF], Birgé and Massart criterion is based on an unbiased risk estimation. The ideal model to estimate s is the one minimizing the risk E [KL(s, ŝD )]. Nevertheless, it is impossible to choose this optimal model since s is unknown. A solution is to find a penalty function, called optimal penalty such that the empirical risk is as close as possible to the risk.

The following classical decomposition of the risk of each estimator ŝD is considered

E [KL(s, ŝD )] = b D + E(V D )
where the bias term b D := KL(s, s D ) and the variance term is defined by V D := ln(s D /ŝ D )s.

Note that the bias b D decreases whereas the variance term V D tends to increase when the dimension D increases. Among the model collection D, the selected model D is the one minimizing the criterion D → γ n (ŝ D ) + pen(D).

(5)

Defining bD := γ n (s D ) -γ n (s) and V D := γ n (s D ) -γ n (ŝ D ), the selected model dimension is also a minimizer of

γ n (ŝ D ) -γ n (s) + pen(D) = bD -V D + pen(D) = KL(s, ŝD ) + ( bD -b D ) -(V D + V D ) + pen(D). ( 6
)
Because of the law of large numbers, it is reasonable to assume that bD -b D ≈ 0. Furthermore, concentration arguments allow to suppose that KL(s, ŝD ) is close to its expectation which is the risk of ŝD . In order to make the quantity in (6) close to the risk E [KL(s, ŝD )],

INRIA the optimal penalty is defined by

pen opt (D) = V D + V D .
Next, the main hypothesis of this heuristics is to assume that V D ≈ V D . An argument to justify this hypothesis is that, in the expressions of V D and V D , the probability measure and the corresponding empirical measure play a similar role. If one permutes these measures inside the definitions of V D and V D , and also in the definitions of s D and ŝD , then V D is changed in V D and reciprocally. Finally, this hypothesis leads to pen(D) = 2 V D . Turning back on the expression of V D , it can be written

V D = bD + γ n (s) -γ n (ŝ D ).
For large dimensions, the bias term stabilizes itself since the approximation of the model cannot be appreciably improved. Thus, the behavior of V D according to the model dimension is known for large dimensions via -γ n (ŝ D ). In our framework, penalty functions could be regarded as proportional to the dimension (see remarks after Theorem 1) and if the slope Ĉ of the linear part of -γ n (ŝ D ) is known, the final penalty is pen(D) = 2 ĈD.

Using the slope heuristics

This section details how the slope heuristics is applied to select a Gaussian mixture model among a family S (K,v) (K,v)∈M with M := {(K, v); 2 ≤ K ≤ K max , v ∈ V} where the maximum number of mixture components K max and the mixture shape are fixed by the user. The heuristics makes use of three steps:

1. Estimation step: The MLE is computed for each model S (K,v) . According to (1), the mixture parameters and ω 2 can be independently estimated. The estimated mixture parameters p1 , . . . , pK , μ1 , . . . , μK , Σ 1 , . . . , Σ K are computed with the Expectation Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] using Mixmod software [START_REF] Biernacki | Model-based cluster and discriminant analysis with the MIXMOD software[END_REF]) and ω2 = 1

n n i=1 y i[v c ] 2 . Thus, the MLE is ŝ(K,v) (x) = K k=1 pk Φ x [v] |μ k , Σ k × Φ(x [v c ] |0, ω2 I Q-α ).
2. Penalty determination step: First, models are grouped according to their dimension in order to obtain the model collection (S D ) D∈D . The function D → -γ n (ŝ D ) is plotted and a threshold D 0 is chosen by the user such that the function has a linear behavior for D ≥ D 0 . Using a robust regression [START_REF] Huber | Robust Statistics[END_REF] of -γ n (ŝ D ) on D, the slope Ĉ of the linear part is estimated. Thus the optimal penalty function is defined by pen opt (D) = 2 ĈD.

The robust regression procedure allows to mitigate the influence of possible estimation errors of the first step.

Applications

This section is devoted to the application of our method on simulated and real datasets. We show that our method allows us to determinate the variable role to improve the clustering. Moreover, we check that the penalized estimator mimics the oracle. The slope heuristics is compared with the classical criteria used for Gaussian mixture model selection: AIC, BIC and ICL. They are respectively defined by

crit AIC (D) = γ n (ŝ D ) + D n , crit BIC (D) = γ n (ŝ D ) + D ln(n) 2n and crit ICL (D) = crit BIC + ENT n
with the entropy term ENT = -

n i=1 K k=1 z ik ln(t ik )
where z is given by the MAP rule and t ik = pk Φ(yi|μ k , Σk ) P K l=1 pl Φ(yi|μ l , Σl ) . The interested reader is respectively referred to [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF][START_REF] Akaike | A new look at the statistical model identification[END_REF], [START_REF] Schwarz | Estimating the dimension of a model[END_REF] and [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] for more details on these criteria.

The method is applied on simulated datasets in Section 5.1 and then on real datasets. In Section 5.2, our procedure is carried out on a curve clustering example for oil production profiles. In Section 5.3, a transcriptome dataset is studied with our method to obtain coexpressed gene clusters.

Simulated datasets

Comparison to the oracle and to other criteria. The aim of this first example is to compare the slope estimator with other penalized estimators and to study its behavior with respect to the oracle. The dataset consists of n = 2000 points described by Q = 32 variables. The data are simulated according to a mixture of four equiprobable Gaussian distributions N (µ k , Σ k ) where µ 1 = (3, 2, 1, 0.7, 0.3, 0.2, 0.1, 0.07, 0.05, 0.025), µ 2 = 0 10 , µ 3 = -µ 1 , µ 4 = (3, -2, 1, -0.7, 0.3, -0.2, 0.1, -0.07, -0.05, -0.025), and Σ 1 = Σ 3 = Σ 4 = I 10 and Σ 2 = diag(2, 1.9, 1.8, . . . , 1.1).

The vector 0 p denotes the null vector of length p. Twenty two independent variables, sampled from a N (0, 1), are appended. Consequently, the true density belongs to the model S (K0,v0)

INRIA

where K 0 = 4 and v 0 = {1, . . . , 10} (α 0 = 10) and the variables are ordered. Note that the discriminant power of the relevant variables decreases with respect to the variable index.

In other words the four subpopulations of the mixtures are progressively gathered together into a unique Gaussian distribution, as shown in Figure 1. This procedure is repeated 1000 times with new simulated dataset each time. The distribution used for simulations is the same as before except that the last ten variables have been removed in order to reduce the computation times. Consequently, the "true model" still corresponds to (K 0 , α 0 ) = (4, 10) but the total number of variables is now Q = 22. These simulations allow to compare the behavior of the slope estimator with the oracle and with the behaviors of three other estimators given by AIC, BIC and ICL criteria. Since the true density is known, a Monte Carlo procedure gives the following oracle model estimation K oracle = 4 and α oracle = 9. Note that even if the true density belongs to the density collection, the oracle model is not equal to the corresponding true model. The results are summarized in Table 5.1. The two criteria BIC and ICL select a model with 4 components and most of the times with 6 relevant variables. It is shown in [START_REF] Keribin | Consistent estimation of the order of mixture models[END_REF] that a model selection procedure using BIC is consistent to find the number of components of a Gaussian mixture when the component densities are bounded. But as far as we know, there no consistency result for such a variable selection and clustering problem. The results of Table 5.1 show that the model selected by BIC is not the true model. In this context, even if BIC tries to find the true model, this could be only done for unrealistic large samples. The behavior of the ICL method is not surprising since the aim of this criterion is to provide a mixture model leading to a sensible partitioning of the data. The AIC method selects too many components and too many relevant variables since it underpenalizes models in the mixture context. From a clustering point of view, BIC, ICL and the slope method have similar performances. The interest of this first simulated example is to illustrate the different behaviors of criteria. As expected, the slope method selects a model close to the oracle model. Waveform dataset The waveform dataset, available at the UCI repository [START_REF] Blake | Uci repository of machine learning databases[END_REF], is composed of three groups based on a random convex combination of two of three wave functions sampled at the integers from 1 to 21, with noise added. A detailed description is available in [START_REF] Breiman | Classification and regression trees[END_REF]. The dataset consists of 5000 observations described by 40 variables. The last nineteen are noisy variables, sampled from a N (0, 1) density. By construction, the first and the twenty-first variables have also the same distribution N (0, 1). Consequently they are both irrelevant for the clustering and thus there is 19 variables which are potentially relevant for clustering. First, the data have been centered. Contrary to the previous example, the variables are not ordered. Ideally, the model collection should be based on all the possible relevant variable subsets v. Nevertheless, the selection among this model family is impossible because of the large cardinal of this collection and the resulting computation times. To get round this problem, the variables are ordered by decreasing order of their variances. With this ordering, the last twenty-one variables are the variables sampled from the N (0, 1) density.

The model selection is proceeded with the two mixture shapes [L k B k ] and [L k C k ]. The plots of D → -γ n (ŝ D ) for the estimation of Ĉ are given on the top of Figure 4 for these two collections of models. To compare the two model collections, both corresponding fittings of the dimension model surfaces on maximum loglikelihood surfaces are presented on the bottom of Figure 4. As expected, we observe that the fitting is dramatically better for the model collection associated to the [L k C k ] collection. Indeed the relevant variables are dependent by construction. The [L k B k ] model collection is not rich enough for this problem. This [L k C k ] collection leads to select a model with K = 3 clusters and α = 19 relevant variables. Despite the simulated data do not follow a Gaussian mixture, the procedure provides a stable and sensible solution. As to other criteria, they all select α = 19 with respectively K = 2, 3 and 10 for ICL, BIC and AIC.

It has been said before that a clustering based on an ill-chosen collection of significant variables can lead to a large clustering error rate. As a matter of fact, Table 2 shows that the three true clusters are found with an error rate of 14.3%. Figure 3 plots the error rate in function of the number α of significant variables. Each curve corresponds to a fixed number of components in the mixture. Choosing a model with more significant variables deteriorates the clustering performance.

Curve clustering

An oil field production profile is the curve of oil production versus time. In the following, the term reserves (or ultimate reserves) denotes the amount of oil that is produced during the exploitation of an oil field. The reader is referred to [START_REF] Babusiaux | Oil and gas exploration and production, reserves, costs, contracts[END_REF] for more details about the exploration and the production of oil. It is well known by the oil industry that production profiles of large fields have a different shape than production profiles of little fields. Indeed, little fields tend to produce their reserves in a short time and early pass the production peak. On the contrary, large fields slowly produce their reserves and their production presents a plateau during several years at the top level. Figure 5 illustrates this behaviour with productions of three fields of different size. In order to compare the production profile shapes, we consider production profiles normalized by the amount of reserves contained in each field. The study's aim is to validate that a clustering of normalized production profiles is consistent with the values of the reserves variable.

The database is composed of several hundred of oil production profiles corresponding to hydrocarbon layers in the North Sea1 . The data used in the procedure are obtained from the original curves as follow. First, each production profile is normalized by the reserves of the corresponding field2 . Ideally, it is desirable to proceed to the clustering on complete production profile. This is impossible since most of the fields are still in production nowadays. Figure 5 suggests that the beginning of the production curve is sufficient to distinguish different shapes in the curve family. Thus, we only consider the subsample composed of 180 fields which have started their production more than 64 months ago. Next, a discrete wavelet transform (DWT) is proceeded on each of these normalized curves. This decomposition has the advantage of giving information on each curve at different resolution levels. This transformation has already been used in curve classification (see for instance [START_REF] Berlinet | Functional classification with wavelets[END_REF]. The reader is referred to [START_REF] Percival | Wavelet methods for time series analysis[END_REF] for details on the DWT. Let W i be the wavelet coefficient vector of the i th curve. Since the length of each curve is 64, the dimension of W i is also 64. The vector W i is defined by

W i = (V i6 , W i6 , . . . , W i1 )
where W ij is a vector of length 64/2 j which is composed of all the wavelet coefficients corresponding to the scale j. The coefficient V i6 is equal to the mean of the curve i divided by √ 64. The hierarchical structure of the DWT suggests a natural order of the wavelet coefficient variables according their resolution. Indeed, V i6 and W i6 give informations about the general shape of the curve i whereas W i1 and W i2 give informations about details on it. We do not use the coefficients in W i1 and W i2 since they correspond to the finer resolution. We will see that the remaining coefficients are sufficient to propose a sensible clustering. Moreover, all the wavelet coefficient variables are centered and scaled to unit variance to make easier the fitting of the multidimensional Gaussian distribution N (0, ω 2 I Q-α ) on the coefficient vectors which are not used for the clustering. These new coefficients are denoted V i6 and W ij where j ∈ {3, . . . , 6}. The procedure is performed on the sample y where y i = V i6 , W i6 , . . . , W i3 for an ordered model collection [LB k ]. This mixture collection allows to avoid estimation problems when the variances are too small.

Figure 6 clearly shows the expected linear behavior of D → -γ n (ŝ D ) in large dimensions. The selected model minimizing the penalized criterion deduced from the slope heuristics has K = 3 components and α = 20 clustering variables. Finally, the MAP rule gives a clustering of the curves. The clusters contain 31, 140 and 9 curves respectively. 7 displays the mean cluster of the normalized production profiles in each cluster. Boxplot of the logarithm of the reserves variable for each cluster are displayed in Figure 8. The second cluster mainly corresponds to the large fields whereas the first and the third clusters contain fields of medium size and small size respectively. The shape of normalized production profiles can be explained by the reserves variable. The reader is referred to [START_REF] Michel | Modélisation de la production d'hydrocarbures dans un bassin pétrolier[END_REF] for more details.

Analysis of a transcriptome dataset

Currently, biologists are interested in gene functional analysis. It is usually considered that co-expressed genes are often implicated in the same biological function and consequently they are potential candidate to be co-regulated genes. Thus biologists try to extract groups of co-expressed genes according to transcriptome datasets in order to characterize more precisely their biological functions. Moreover an experiment selection for the clustering is desirable to improve the clustering and its interpretation with a biological point of view.

Here we study a transcriptome dataset of Arabidopsis thaliana extracted from the database CATdb [START_REF] Gagnot | CATdb: a public access to Arabidopsis transcriptome data from the URGV-CATMA platform[END_REF]. To build this database, an identical statistical analysis for all transcriptome experiments has been performed to remove the technical biases (normalization) and to determine the gene significantly differentially expressed (differential analysis) between two conditions. In this differential analysis, we test if a gene is non-differentially expressed or not in the experiment j. For this test, a test statistic corresponding to the normalized differential expression and a p-value adjusted by the Bonferroni method are determined. Then a gene is declared to be differentially expressed when its Bonferroni p-value is lower than 0.05. The reader is referred to [START_REF] Lurin | Genome-wide analysis of arabidopsis pentatricopeptide repeat proteins reveals their essential role in organelle biogenesis[END_REF] for a description of such an analysis.

We focus on 305 genes of Arabidopsis thaliana studied on ten experiments which correspond to mutant conditions or different stress situations. These genes are declared differentially expressed in the two last experiments and non-differentially expressed in five experiments. Table 3 gives the number of differentially expressed genes per experiment. Each gene is described with a vector y i ∈ R 10 , the component y ij corresponding to the test statistic calculated in the experiment j for the differential analysis. Since there is not a natural way to order the variables, our procedure for non-ordered variables is performed with the [LC] mixture collection. The maximum number of components is fixed to K max = 40. After the estimation step, we notice that the function D → -γ n (ŝ D ) has a linear behavior for D ≥ 220 (see Figure 9), thus the slope method can be applied. The procedure selects a clustering with K = 8 clusters based on the seven variables v = {3, 5, 6, 7, 8, 9, 10}. The eight clusters have different size (see Table 4) and the clustering shows some interesting similar behaviors of expression profiles (see Figure 10). A similar clustering can be found if all the variables are considered (α = 10 fixed) but with the variable selection, the interpretation of the clustering is made clearer.

First, we note that the two benchmark experiments (9 and 10) where all genes are differentially expressed are selected. Moreover, the three variables which are not selected for the clustering are three variables where all genes are non-differentially expressed. The average behavior of genes per cluster is the same in the irrelevant experiments 1, 2 and 4 since it is concentrated around zero. On the contrary, genes of the cluster 2 have a particular behavior in experiments 7 and 8. Their expression difference decreases between the two experiments (7 and 8) whereas the genes of the other clusters have the same expression in these two experiments (see Figure 10). This remark may explain why the two experiments 7 and 8 where all genes are non-differentially expressed are selected for the clustering while experiments 1, 2 and 4 are not. This clustering can help biologists to find gene biological functions. For instance, 12 genes for which biologists know any biological function are clustered with other genes for which biologists have some information.

Discussion

In this paper, a methodology has been proposed to take into account the variable role in a clustering process in a model-based cluster analysis setting. The interest of our approach is to recast these two problems into a model selection problem where model collections are indexed by two quantities (K and v). The practical use of the penalized likelihood criterion, proposed in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF], is based on a slope heuristics method allowing to calibrate the multiplicative constant in the penalty term. The slope Ĉ is estimated on the restriction of D → -γ n (ŝ D ) to D ≥ D 0 , namely where this function has a linear behaviour. A robust regression is used allowing to mitigate the influence of possible estimation errors. The estimated slope Ĉ which depends on the user choice threshold D 0 can be confirmed by plotting the estimated slope in function of different values of the threshold. In a neighbourhood of Ĉ, the slope has to be stable. The behaviour of this slope heuristics method has been studied on simulated and real datasets. It has been compared with standard criteria used in Gaussian mixture clustering context. BIC, ICL, AIC and our criterion have different goals thus it is natural to observe different behaviours in the examples. In particular, our criterion tends to select the oracle model.

Our method can be efficiently applied when the variables could be ordered as illustrated in the curve clustering study in Section 5.2. When the variables cannot be ordered in a natural way, our method is difficult to use when the number of variables is too large. An exhaustive research of the best model becomes untractable. A possible way to circumvent this problem is to find a convenient strategy allowing to run the estimation step only for a model subset with a reasonable size.

In all cases, the user has to check that the linear behaviour of D → -γ n (ŝ D ) is observed and also that the dimension model surface fits the maximum likelihood surface on high dimensions (see Figure 4). If it is not the case, several explanations can be given. First, the model dimension can be too low, namely the maximum number of components K max has to be increased. The problem can be also related to the model family choice. Roughly speaking, the family model leads to a stabilization of the bias in large dimension if the family model efficiently approaches the true density. Otherwise, the collection of models has to be changed in order to obtain a better fitting between the two surfaces.

Our theoretical and practical results could be extended for most of the twenty-eight Gaussian mixture shapes proposed by [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF]. In particular, without variable selection (α = Q), our method allows us to select the number of clusters which is the fundamental problem in model-based clustering. It could be envisaged to adapt our works to select also the Gaussian mixture type.

with s α ∈ {-N α , . . . , N α } where N α = a(1+δ) α+z 2 - S J 4 √ λ 1-α M βM c1δ
and c 1 := 1 -2 -1 4 . For a given couple (z, J), this insures that for all vectors µ ∈ [-a, a] α , there exists a vector ν (z,J) of this lattice such that

(1 + δ) z β M λ M α-1 q=1 ν (z,J) -µ 2 q (1 + δ) j(q)+1 2 + ν (z,J) -µ 2 α (1 + δ) -S J 2 +α λ α M ≤ c 1 α δ 2 . (7)
For a couple (z, J) and a vector ν (z,J) defined as before, the two following associated functions are considered

l(x) = (1 + δ) -2α Φ x|ν (z,J) , (1 + δ) -1 4 β z+1 B l J u(x) = (1 + δ) 2α Φ x|ν (z,J) , (1 + δ) β z B u J . (8)
Second, we check that for all x ∈ R Q , l(x) ≤ Φ(x|µ, βB) ≤ u(x). According to Proposition 6 of Appendix C in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF] which allows to upper bound the ratio of two Gaussian densities, we get

Φ(x) u(x) ≤ (1 + δ) -2α |(1 + δ)β z B u J | |βB| exp 1 2 (ν (z,J) -µ) {(1 + δ)β z B u J -βB} -1 (ν (z,J) -µ) ≤ (1 + δ) -3α+1 4 exp 1 2δ (ν (z,J) -µ) (β z B u J ) -1 (ν (z,J) -µ) and l(x) Φ(x) ≤ (1 + δ) -2α |βB| |(1 + δ) -1 4 β z+1 B l J | exp 1 2 (ν (z,J) -µ) βB -(1 + δ) -1 4 β z+1 B l J -1 (ν (z,J) -µ) ≤ (1 + δ) -( 9α 8 + 1 4 ) exp 1 2δ(1 -2 -1 4 ) (ν (z,J) -µ) (β z+1 B l J ) -1 (ν (z,J) -µ)
using the concavity of the function δ → 1 -(1 + δ) -1 4 . The following inequalities (ν (z,J) -µ) (β z B u J ) -1 (ν (z,J) -µ) ≤ δ 2 4 (3α + 1) ( 9) and (ν (z,J) -µ) (β z+1 B l J ) -1 (ν (z,J) -µ) ≤ δ 2 (1 -2 -1 4 ) 9α + 2 8 (10)

are then sufficient to have l ≤ Φ ≤ u. We can check that condition (7) implies the two inequalities ( 9) and ( 10).

Third, we show that d H (u, l) ≤ ε. According to Proposition 7 of Appendix C in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF], we have that the quantity d 2 H (l, u) is equal to

(1 + δ) -2α + (1 + δ) 2α -2 α 2 +1 (1 + δ)β z B u J (1 + δ) -1 4 β z+1 B l J -1 4 (B u J ) -1 (1 + δ)β z + (1 + δ) 1 4 (B l J ) -1 β z+1 -1 2 = (1 + δ) -2α + (1 + δ) 2α -2
(1 + δ) Finally, we can construct an ε-bracket family (with respect to d H ) to cover F (K,α) . Let (Φ(.|µ 1 , βB 1 ), . . . , Φ(.|µ K , βB K )) be an element of F (K,α) . Let z ∈ {0, . . . , r } and J 1 , . . . , J K in {2, . . . , r} α-1 such that for all k ∈ {1, . . . , K} and all q ∈ {1, . . . , α},

β z+1 B l J k ,q ≤ β B k,qq ≤ β z B u J k ,q .
For all k, there exists a vector ν (z,J k ) such that condition (7) is satisfied for the mean vector µ k . For z, J k and ν (zJ k ) , the two associated functions defined by (8) are denoted u k and l k . Then we define L := (l 1 , . . . , l K ) and U := (u 1 , . . . , u K ). The set of all such brackets [L, U ] covers the family F (K,α) and is denoted R(ε, K, α). An upper bound of the bracketing entropy of F (K,α) is thus determined by the computation of the cardinal of R(ε, K, α).

If 2a ≥ √ β M c 1 ( √ λ M ∨ λ M 1-α ) then N [.] (ε, F (K,α) , d H ) ≤ card R(K, ε, α) ≤ r z=0 J α-1 q=1 2 a (1 + δ) j(q)+1 4 + z 2 √ β M λ M c 1 δ 2 a (1 + δ) α+z 2 - S J 4 β M c 1 λ M 1-α δ K ≤ r z=0 (r -1) α-1 2a √ β M c 1 δ α (1 + δ) 3α-1 4 + z 2 K (11) ≤ (r -1) K(α-1) 2a √ β M c 1 δ Kα (1 + δ) K( 3α-1 4 + z 2 ) (1 + δ) Kαr 2 (r + 1) ≤ A 2 α ε Kα+K(α-1)+1 (12)
INRIA using 1 + δ ≤ 2 and the definition of r, r and δ. The constant A 2 depends on parameters a, β m , β M , λ M and λ m . If 2a ≤ √ β M c 1 ( √ λ M ∨ λ M 1-α ), the upper bound ( 12) is again satisfied even if the constant A 2 is modified.

B Bracketing entropy of the model collection [LC]

For [LC] mixture shape, the variance matrices in mixture are all equal to the same positive definite matrix. Thus, we need to find an upper bound of the number bracketing of the following set where the constant A 3 only depends on λ m , λ M , and a, and Kα + α(α+1) 2 is the dimension of F (K,α) . This result is obtained by considering the following bracket family. Its construction is inspired by the bracket family used in the study of the bracketing entropy of L (K,α) for the model [L k C k ] stated in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF].

ν q = √ λ m ε 3α s q ; s q ∈ Z ∩ - 3 a α √ λ m ε , 3 a α √ λ m ε ,
is an ε-bracket set over F (K,α) .

Finally the bracketing number of F (K,α) is upper bounded by

N [.] (ε, F (K,α) , d H ) ≤ card(R(β)) × card(X (ε, a, λ m , α)) K ≤ 2 λ M β α(α+1) 2 × 6aα √ λ m ε Kα ≤ 6 √ 3 λ M α 2 ε λ m α(α+1) 2 × 6aα √ λ m ε Kα ≤ A 3 α ε Kα+ α(α+1) 2
where the constant A 3 only depends on λ m , λ M , and a.

3.

  Model selection step: The minimizer D of the criterion D → γ n (ŝ D ) + 2 ĈD is determined. Then, among the initial model collection, we find the model ( K, v) fulfilling D( K, v) = D and ŝ( K,v) = ŝ D . Finally, the parameter estimators associated to ( K, v) provide a data clustering using the MAP rule.

Figure 1 :

 1 Figure 1: Boxplots of the first eleven variables (VAR1,...,VAR11) on the four mixture components (G1,G2,G3,G4).

Figure 2 :

 2 Figure 2: On the top graph, the function D → -γ n (ŝ D ) is plotted. The linear regression is made for D ≥ D 0 = 300. The associated residuals are drawn on the bottom graph.

Figure 3 :Figure 4 :

 34 Figure 3: Evolution of the clustering error rate in function of the chosen number α of significant variables. Each curve corresponds to a fixed number K of mixture components.

Figure 5 :

 5 Figure 5: Oil production profiles normalized by reserves of three fields located in the North Sea.

Figure 6 :

 6 Figure 6: Slope method applied to the [LB k ] collection for the wavelet coefficient curve data.

Figure 7 :

 7 Figure 7: Normalized production profile means for each cluster.

Figure 8 :

 8 Figure 8: Boxplots of the logarithm of the reserves variable for each cluster.

INRIA

  Number of differentially expressed genes per experiment.

  Number of genes per cluster.

Figure 9 :

 9 Figure 9: Penalty determination on the linear behavior of the function D → -γ n (ŝ D ).

F

  (K,α) = {(Φ(.|µ 1 , Σ), . . . , Φ(.|µ K , Σ); µ k ∈ [-a, a] α , Σ ∈ D + (α) (λ m , λ M )} where D + (α) (λ m , λ M )is the set of α×α symmetric positive definite matrices whose eigenvalues are in the interval [λ m , λ M ]. Proposition 3. For all ε ∈ (0, 1], N [.] F (K,α) ≤ A 3

Proposition 4 .

 4 For all ε ∈ (0, 1], let δ = ε √ 3 α and β = λm ε 3 √ 3α 2 . The following set ([l 1 , u 1 ], . . . , [l K , u K ]); u k (x) = (1 + 2δ) α Φ (x| ν k , (1 + δ)A) l k (x) = (1 + 2δ) -α Φ x| ν k , (1 + δ) -1 A ; A ∈ R(β), ν k ∈ X (ε, a, λ m , α)whereR(β) = A = (A ij ) 1≤i,j≤α ; A ij = a ij β; a ij = a ji ∈ Z ∩ -λ M β , λ M βand X (ε, a, λ m , α) = ν = (ν 1 , . . . , ν α );

Table 1 :

 1 For each criterion, number of times that a model (K, α) is selected among the 1000 simulations.

	INRIA

Table 2 :

 2 Contingency table for the clustering obtained with the slope heuristics.

		cluster 1 cluster 2 cluster 3 total
	group 1	1331	185	176	1692
	group 2	95	99	1459	1653
	group 3	65	1494	96	1655
	total	1491	1778	1731	5000

  ≤ 9α 2 δ 2 = ε 2 .

								11 8 + (1 + δ) -11 8	-α-1 2	(1 + δ)	2α+7 8	8 + (1 + δ) -2α+7	-1 2
									2	2
	= [2 cosh {2α ln(1 + δ)} -2] + 2 -2 cosh	8 11	ln(1 + δ)	-α-1 2
	+ 2 cosh	11 8	ln(1 + δ)	-α-1 2	1 -cosh	8 7 + 2α	ln(1 + δ)	-1 2
	≤ 4 sinh(1)α 2 δ 2 + 2	α -1 2	11 8	δ 2 + 2	7 + 2α 8	1 2	δ 2

RR n°6550

Data is available on the website of the Norwegian Petroleum Directorate : www.npd.no/engelsk/cwi/ pbl/en/index.htm, and the website of the English Department of Trade and Industry (DTI): www.og.dti. gov.uk/fields/fields index.htm.

The DTI does not provide estimations of reserves of their fields, consequently we use the IHS database http://energy.ihs.com for the english fields of the sample INRIA

Acknowledgements

The authors are grateful to Gilles Celeux (INRIA), Pascal Massart (Université Paris-Sud 11) and Marie-Laure Martin-Magniette (INRA) for helpful discussions and valuable comments.

Appendices

In this paper, the two model collections [LB k ] and [LC] used for the real dataset examples in Section 5 are not entering the theoretical study of [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF]. Nevertheless it can be proved that the penalty functions have the same form as in Theorem 1. We sketch the proof in this appendix. According to the work of [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF], it is only required to determine an upper bound of the bracketing entropy of the collection S (K,v) for the Hellinger distance. To determine this upper bound, we cannot use Proposition 1 of Appendix A in [START_REF] Maugis | A penalized criterion for Gaussian mixture model selection[END_REF] which allows to recast the problem as the study of the bracketing entropy of the associated mixture density family. Indeed, in these two collections, the variance matrices of mixtures have the same volume or are identical. Here we use the following result which can be proved along the lines of the proof of Theorem 2 in [START_REF] Genovese | Rates of convergence for the Gaussian mixture sieve[END_REF].

The notion of bracketing number, denoted N [.] , for a family of nonnegative integrable functions

If F (K,α) denotes the set of K-uples of Gaussian densities used in the mixtures of L (K,α) , Proposition 1 allows to give an upper bound of the bracketing number of L (K,α) from an upper bound of the bracketing entropy of the F (K,α) for the Hellinger distance.

A Bracketing entropy of the model collection

For [LB k ] mixture shape, the variance matrices of a mixture have the same volume and are all diagonal. These variance matrices can be decomposed into

is the common volume and B k is a diagonal matrix with determinant 1. Thus, the family

where ∆ 1 (α) (λ m , λ M ) is the set of α × α diagonal matrices with determinant 1 and which eigenvalues are in the interval [λ m , λ M ] where 0 < λ m < λ M . According to Proposition 1, the aim is to build a bracket family for F (K,α) and to determine its cardinal in order to obtain the following result:

where the constant A 2 only depends on λ m , λ M , β m , β M and a and K(2α -1) + 1 is the dimension of F (K,α) .

Proof. The proof for the unidimensional case (α = 1) is already available in [START_REF] Genovese | Rates of convergence for the Gaussian mixture sieve[END_REF]. Let ε ∈ (0, 1] and assume K ≥ 2 and α ≥ 2 fixed. Let δ = ε/(3α). For j ∈ {2, . . . , r}, we define b

h denotes the smallest integer greater than or equal to h. For z ∈ {0, . . . , r } we consider

For a vector J = (j(1), . . . , j(α -1)) ∈ {2, . . . , r} α-1 , the diagonal matrices B l J and B u J are defined by

α-1 q=1 j(q). The q th diagonal coefficients of these matrices are denoted B l J,q and B u J,q respectively. First, a function Φ(.|µ, βB) such that β ∈ [β m , β M ], µ ∈ [-a, a] α and B ∈ ∆ 1 (α) (λ m , λ M ) is considered. Let z ∈ {0, . . . , r } be the unique integer of {0, . . . , r } such that β z+1 < β ≤ β z and let J be the unique vector of {2, . . . , r} α-1 such that ∀q ∈ {1, . . . , α -1}, B l J,q ≤ B qq ≤ B u J,q . Hence for all q ∈ {1, . . . , α},

For a couple (z, J), we also consider a regular lattice of mean vector ν (z,J) = ν (z,J) 1 , . . . , ν (z,J) α ∈ [-a, a] α such that for all q ∈ {1, . . . , α -1}, ν (z,J) q = (1 + δ) -j(q)+1 4 -z 2 λ M β M c 1 δs q , with s q ∈ {-N q , . . . , N q } where N q = a(1+δ)