Christophe Biernacki

Julien Jacques

A generative model for rank data based on an insertion sorting algorithm

Keywords: Full rank data, sorting process, insertion sort algorithm, EM algorithm, quiz data

An original and meaningful probabilistic generative model for full rank data modelling is proposed. Rank data arise from a sorting mechanism which is generally unobservable for statisticians. Assuming that this process relies on paired comparisons, the insertion sort algorithm is known as being the best candidate in order to minimize the number of potential paired misclassifications for a moderate number of objects to be ordered. Combining this optimality argument with a Bernoulli event during paired comparison step, a model that possesses desirable theoretical properties, among which are unimodality, symmetry and identifiability is obtained. Maximum likelihood estimation can also be performed easily through an EM or a SEM-Gibbs algorithm (depending on the number of objects to be ordered) by involving the latent initial presentation order of the objects. Finally, the practical relevance of the proposal is illustrated through its adequacy with several real data sets and a comparison with a standard rank data model.

Introduction

Ranking data is of great interest in human activities involving preferences, attitudes or choices like Web Page ranking, Sport, Politics, Economics, Educational Testing, Biology, Psychology, Sociology, Marketing, etc. Ranks are so meaningful that it is not unusual for them to be found as the result of a transformation of other kinds of data. In this paper only full rankings will be considered but possible extensions for partial, tied or incomplete rankings are being considered as future prospects of the proposed model.

The rank datum results from a ranking of m objects O 1 , . . . , O m by a judge (human or not). Two representations of these data are commonly used: ranking and ordering. The ranking representation x -1 = (x -1 1 , . . ., x -1 m ) contains the ranks assigned to the objects, and means that O i is in x -1 i th position (i = 1, . . ., m). A ranking is then an element of P m ,

The two most popular classes of models for rank data consist on modelling directly the hypothetical ranking process followed by the judge. For a complete review, refer to [START_REF] Marden | Analyzing and modeling rank data[END_REF]Chap. 5 to 10]. The first class is derived from a paired comparison process [START_REF] Kendall | On the method of paired comparisons[END_REF]: the judge constructs a rank by first comparing each pair of objects, and second ensuring the consistency of these paired comparisons (if O 1 is prefered to O 2 and O 2 to O 3 , O 1 must be prefered to O 3 ). It leads to the Babington Smith model [27, p. 116]. The number of parameters of this model is very large especially when m grows, therefore some simplifications have been considered. [START_REF] Bradley | Rank analysis of incomplete block designs. I. The method of paired comparisons[END_REF] associate to each object O j a score indicating an overall degree of preference of this object, and connect these scores to the Babington Smith model's parameters, which defines the Bradley-Terry-Mallows model [27, p. 117]. [START_REF] Mallows | Non-null ranking models[END_REF] goes forward into the simplification by imposing that the scores only depend on the sign of µ -1

x i -µ -1 x j , where µ is a "reference" rank. It leads to the famous Mallows Φ model:

p(x; µ, λ ) = C (λ ) -1 exp -λ K(x,µ) ,
where K is the Kendall distance between two ranks [START_REF] Kendall | A new measure of rank correlation[END_REF][START_REF] Critchlow | Metric methods for analyzing partially ranked data[END_REF] and where

C (λ ) = m-1 ∏ j=1 1 -exp(-(m -j + 1)λ ) 1 -exp(-λ )
is a normalization constant [START_REF] Fligner | Distance based ranking models[END_REF] with λ ∈ R a precision parameter. For instance, a high λ value leads to strong unimodality around µ. In the last decade, some extensions of this model have been proposed: [START_REF] Lebanon | Cranking: Combining rankings using conditional models on permutations[END_REF] use a generalization of Mallows Φ model to combine multiple input rankings, and [START_REF] Lebanon | Non-parametric modeling of partially ranked data[END_REF] model partially ranked data via semi-parametric model using Mallows kernel.

1

The second popular class of rank data models is multistage models, which consider the following iterative ranking process: the judge selects firstly the best object among the m ones, then the best among the m -1 remaining ones, and so on. The Plackett-Luce model [START_REF] Luce | Individual choice behavior: A theoretical analysis[END_REF][START_REF] Plackett | The analysis of permutations[END_REF] is then based on the probability that each O x i is ranked first among the m objects [27, p. 119]. It could be noticed that this model corresponds to Thurstonian model [START_REF] Thurstone | A law of comparative judgment[END_REF][START_REF] Böckenholt | Applications of Thurstonian models to ranking data[END_REF] with a Gumbel density. [START_REF] Benter | Computer-based horse race handicapping and wagering systems: A report[END_REF] proposed a variant by introducing a dampening parameter and [START_REF] Gormley | A latent space model for rank data[END_REF] combined the Plackett-Luce model with a latent space model in order to spatially model the ranked nature of the data. [START_REF] Fligner | Multistage ranking models[END_REF] introduce an alternative multistage model by considering another form of probability at each step of the ranking process. It leads to the Fligner and Verducci's strongly unimodal model [27, p. 120]. Assuming specific forms of this model could lead to Mallows Φ model or to a generalization of this latter named Φ component-model [27, p. 121].

Thus, the ranking processes that have motivated these two classes of rank data models can be interpreted as two different sorting processes, in which stochastic errors are introduced to define a probability distribution on the whole rank data space. The natural question involved by this interpretation is whether the used sorting algorithms are the most appropriate. For instance, in paired comparison models, it is not optimal to do so many comparisons since it leads to a sorting algorithm with excessively high computational complexity. Here, we propose a generative model for rank data based on (straight) insertion sort algorithm. This is one of the most powerful sorting algorithms among the usual ones when m ≤ 10 [21, Chap. 5], the situation expected to be the most frequent particularly in case of "human ranking". However, the proposed model can be applied in practice to any number of objects. This new kind of generative model enjoys good theoretical properties and has originality to involve the (potentially latent) initial presentation order of the objects.

The paper is organized as followed. Section 2 introduces and sets up the proposed model which is based on insertion sort algorithm, and its theoretical properties (unimodality, symmetry, identifiability) are also detailed. Maximum likelihood estimation is considered in Section 3 by means of an EM or a SEM-Gibbs algorithm (depending on the number of objects to be ordered) since a missing data interpretation of the proposed model can be exhibited. Numerical illustrations are presented in Section 4 to evaluate the relevance of the proposed model on real data sets both from a distributional adequacy and comparison to the usual Mallows Φ model point of view. Since this work sheds a new light on rank data modelling, numerous related perspectives are discussed in the last section (Section 5).

A generative model for rank data based on insertion sort

Motivation to adopt insertion sort algorithm

We assume it exists an ordering µ = (µ 1 , . . . , µ m ) on the m objects, so that a judge who perfectly sorts these objects returns this reference rank µ. Making also the natural as-Table 1: An example to illustrate the standard insertion sort process with µ = (1, 2, 3) and

y = (1, 3, 2) step j unsorted sorted start y = 1 3 2 - 1 3 2 1 2 2 3 ? ↔ 1 1 3 3 - 2 ? ↔ 1 3 1 2 ? ↔ 3 x = 1 2 3
sumption that a rank x = (x 1 , . . . , x m ) is the result of a sorting process relying on successive object paired comparisons, any difference between the final rank x and µ is necessarily attributed to some incorrect paired comparisons. As a consequence, reducing the gap between x and µ is strongly correlated to minimize the number of paired comparisons involved in the sorting process. For instance, for a moderate number of objects (m ≤ 10), an "optimal judge" should adopt the insertion sort algorithm which is optimal in this case [START_REF] Knuth | Sorting and Searching: Volume3. The art of Computer Programming[END_REF]Chap. 5]. Arguing that m ≤ 10 is a frequent situation in practice (in particular for "human rankings"), we propose to keep this sorting algorithm although in practice our model will not impose any upper bound on m.

We briefly sketch how the standard insertion sort algorithm performs. More details can be found in [START_REF] Knuth | Sorting and Searching: Volume3. The art of Computer Programming[END_REF]Chap. 5] for instance. We note y = (y 1 , . . . , y m ) the initial order of the objects presented to the judge. First, the current object in y to be sorted is placed before the already sorted objects, and is compared to the first object after. If the relative position of both objects in this pair is correct (according to µ), this pair order is unchanged and this process is restarted with the next object in y. Otherwise, the pair order is reversed and a new pair comparison is performed with the next object after (if it exists). And so forth, until obtaining the final ranking x. An example is also displayed in Table 1 with m = 3 to give step by step an overview of standard insertion sort algorithm.

Modelling a stochastic insertion sort algorithm

Our idea is to merge the previous deterministic process with the following stochastic paired comparisons process. We assume that each paired comparison is the result of a Bernoulli experiment whose outcome is a correct comparison (according to µ) with probability π ∈ [0, 1] and an incorrect comparison with probability 1-π. The parameter π models the reliability of the judge about the "true" rank µ. We also assume that each pair ranking operation is independent from others and that the probability π is constant along the sorting process. Thus, at the end of the proposed stochastic process, the final rank x can differ from the reference rank µ. From these simple ideas, we obtain a meaningful generative model for rank data that is presented at length now.

Based on this stochastic modelling of the insertion sort algorithm, the probability p(x|y; µ, π) to obtain a rank x from an initial presentation order y, with a reference rank µ, is given by: p(x|y; µ, π) = π G(x,y,µ) (1 -π) A(x,y)-G(x,y,µ) , (

where π is the probability of good paired comparison (according to µ), G(x, y, µ) is the number of good paired comparisons in the sorting process, and A(x, y) is the overall number of paired comparisons. In addition, denoting by p(y) the presentation order distribution, the marginal distribution of x is given by: p(x; µ, π) = ∑ y∈P m p(x|y; µ, π)p(y).

(2.2)

In this paper, we assume the presentation orders are unknown and uniformly distributed, thus p(y) = m! -1 for all y ∈ P m . In the sequel the rank data model defined by distribution (2.2) will be quoted as ISR(µ, π) for Insertion Sort Rank data model associated to parameters (µ, π).

The proof of (2.1) is given in B. At the beginning of this appendix, notations A(x, y) and G(x, y, µ) are also mathematically defined and illustrated through an example in Table 5. Note that Tables 1 (deterministic insertion algorithm) and 5 (stochastic insertion algorithm) are different because they lead to a different x value.

Remark

The fact that π remains constant along the sorting process makes sense as the judge's knowledge does not change on µ as well as the tiredness of the judge is negligible. However, this hypothesis could be weakened. This issue is discussed in Section 5.

Properties of ISR model

In this section, the main properties of the ISR model are stated. Precise statement of each property as a mathematical proposition and related proofs are available in C at the end of this paper. Proofs rely on applying permutation properties on both ranking and ordering notations on P m .

• Uniformity for π = 1 2 . In the case where paired comparison is performed at random (π = 1 2 ), the ISR is the uniform distribution on P m . In this case, the reference rank µ can be arbitrarily chosen. See Proposition 3 in C.

• Mode µ. One of the most important properties which can be expected from the ISR distribution is that the reference rank µ is the unique mode of the distribution if π > 1 2 . See Proposition 4 in C.

• Anti-mode μ. Let μ be defined by μ = µ • ē where ē = (m, . . ., 1) is the permutation of total inversion. This rank μ is the farthest from µ for the Kendall distance.

Then the unique anti-mode (the rank of smallest probability) is μ if π > 1 2 . See Corollary 1 in C.

• Link between µ and π. The mode µ is also uniformly more pronounced when π grows. See Proposition 5 in C.

• Symmetry. Distributions ISR(µ, π) and ISR( μ, 1 -π) are equivalent. This property will be especially useful to exhibit the identifiability conditions of the ISR distribution below. See Proposition 6 in C.

• Identifiability. The uniformity for π = 1 2 of the ISR distribution and its symmetry lead to impose π > 1 2 as a necessary condition for identifiability. In fact this condition is also sufficient to have identifiability. See Proposition 7 in C.

Estimation of the model parameters

The ISR model for rank data has two parameters: the probability π ∈ [ 1 2 , 1] and the reference rank, or modal rank, µ, which can take its values in P m . Note that the case π = 1 2 is kept, although this is a non-identifiability situation because it leads to the uniformity of the ISR distribution, that can be of interest for practical applications. Considering (x 1 , . . . , x n ) as an independent sample of n ranks from ISR(µ, π), we present in this section estimation of (µ, π) by maximizing the log-likelihood of the ISR model given by

l(µ, π; x 1 , . . . , x n ) = n ∑ i=1 ln 1 m! ∑ y∈P m p(x i |y; µ, π) .
We assume in the following that pairs (x i , y i ) arise independently (i = 1, . . . , n), where (y 1 , . . . , y n ) denote all the latent presentation orders.

Using an EM algorithm for a small number of objects (m ≤ 7)

As (y 1 , . . ., y n ) are unknown, we use an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] to maximize the observed data log-likelihood. Denoting by (µ, π) {0} , the parameter starting values for the EM algorithm and by (µ, π) {q} the value of the parameters at the step q (q ∈ N), the two steps (E and M) of this algorithm are described as follows. This algorithm is computationally feasible for a small number of objects, typically m ≤ 7. However, it could lead to computational difficulties for larger m values because of the sum over P m involved in both E and M steps (see discussion at the beginning of Section 3.3). Note that since the conditional probability (2.1) is invariant to an inversion of the first two elements of the presentation order (Lemma 1 of D), the number m! of presentation orders y to be considered in the calculation of the probability (2.2) may be reduced by half. This remark can be used in order to accelerate the EM algorithm.

The E step

The complete-data log-likelihood is given by

l c (µ, π; x 1 , . . . , x n , y 1 , . . ., y n ) = n ∑ i=1 ∑ y∈P m 1{y = y i } ln 1 m! p(x i |y; µ, π) .
The E step consists of computing the conditional expectation Q of l c expressed by:

Q((µ, π), (µ, π) {q} ; x 1 , . . . , x n ) = n ∑ i=1 ∑ y∈P m t {q} iy ln 1 m! p(x i |y; µ, π) ,
where the conditional probability t {q} iy that y i = y is given by:

t {q} iy = p(x i |y; (µ, π) {q} ) ∑ τ∈P m p(x i |τ; (µ, π) {q} )
.

The M step The M step consists of choosing the value (µ, π) {q+1} which maximizes the conditional expectation Q computed at the E step:

(µ, π) {q+1} = argmax (µ,π)∈P m ×[ 1 2 ,1] Q((µ, π), (µ, π) {q} ; x 1 , . . ., x n ).
For the modal rank µ, it is however numerically expensive to explore all the discrete space P m even for relatively small values of m. To overcome this difficulty, a specific strategy will be proposed in Section 3.4. The value of probability π, maximizing Q leads to following update:

π {q+1} = ∑ n i=1 ∑ y∈P m t {q} iy G(x i , y, µ {q} ) ∑ n i=1 ∑ y∈P m t {q} iy A(x i , y)
.

Note that this value of π {q+1} can be interpreted as the proportion of good manipulations (switching to right or stop) in the insertion sort algorithm.

The algorithm stops when the difference of the log-likelihood between two successive iterations is less than a given threshold. We discuss now how to efficiently start the algorithm.

Initialization strategy for π in EM

We propose here a straightforward asymptotic bound on π to initialize EM. We will discuss in Section 3.4 how to also restrict the possible values for µ. (3.1)

Proof. Using Lemma 6 (see D) and also the fact that, for any µ and y, p(µ|y; µ, π) = π A(µ,y) (see the proof in Lemma 5), it leads to the following bounds for the probability of µ:

π m(m-1)/2 ≤ p(µ; µ, π) ≤ π m-1 .
Since f 0 is a consistent estimator of p(µ; µ, π), it ends the proof.

As soon as πand π+ are greater than 1 2 , this result is useful to initialize π in EM by choosing uniformly π {0} in the interval given by (3.1). If only π+ ≥ 1 2 , the interval becomes [ 1 2 , π+ ]. If both bounds are lower than 1 2 , then the interval [ 1 2 , 1] must be used. In Table 4 of Section 4, bounds associated with all the data sets are greater than 1 2 and the retained intervals are quite narrow in comparison to [ 1 2 , 1], so the strategy seems to be efficient.

Using a SEM-Gibbs algorithm for a large number of objects (m ≥ 8)

As discussed at the beginning of Section 3.1, the EM algorithm is not appropriate as soon as m ≥ 8 because of the factorial term present in both E and M steps: for m = 8 or m = 9, EM could be still feasible although it would be extremely slow (respectively about 4.10 4 and 4.10 5 sums involve in E and M steps); for m ≥ 10, this computational times would become definitively inaccessible (more than 3.10 6 sums involve in E and M steps).

A so-called SEM-Gibbs algorithm may provide an efficient and elegant solution. The fundamental idea of this algorithm is to reduce the computational complexity that was present in both E and M steps of EM by removing all explicit and extensive use of the conditional probabilities t {q} iy . It relies on the SEM algorithm [START_REF] Geman | Stochastic relaxation, gibbs distributions and the bayesian restoration of images[END_REF][START_REF] Celeux | The SEM algorithm: A probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF] which generates the latent variables y i at a so-called stochastic step (S step) from the conditional probabilities t {q} iy computed at the E step. Then these latent variables are directly used during the M step. However, the advantage with SEM-Gibbs algorithm relies on the fact that the latent variables are generated without calculating conditional probabilities thanks to Gibbs algorithm. The proposed SEM-Gibbs algorithm proceeds in the following steps:

The SE-Gibbs step It consists of generating a sample y i{q} from t {q} iy (for all i ∈ {1, . . ., n}) like in a SE step of the standard SEM algorithm but without any calculation of conditional probabilities t {q} iy since it invokes instead the following Gibbs algorithm. Starting from an arbitrary sample y i{q,0} , generate r ∈ {1, . . ., R} sequences y i{q,r} (R being a given number) where (y i{q,r+1} j , •) ∼ p y i j , y i j+1 |(y 1 , . . . , y j-1 ) i{q,r+1} , (y j+2 , . . ., y m ) i{q,r} , x 1 , . . . , x n ; (µ, π) {q} for j ∈ {1, . . ., m -2} and where

(y m-1 , y m ) i{q,r+1} ∼ p y i m-1 , y i m |(y 1 , . . . , y m-2 ) i{q,r+1} , x 1 , . . . , x n ; (µ, π) {q} .
Note that both previous expressions do not involve any combinatorial calculation. If R is large enough, y i{q,R} arises from t {q} iy , thus we retain y i{q} = y i{q,R} (i = 1, . . ., n).

The M step

The M step consists of choosing the value (µ, p) {q+1} which maximizes the completed log-likelihood computed at the SE-Gibbs step:

(µ, π) {q+1} = argmax (µ,π)∈P m ×[ 1 2 ,1]
l c (µ, π; x 1 , . . . , x n , y 1{q} , . . . , y n{q} ).

For the modal rank µ, it is however numerically expensive to explore all the discrete space P m even for relatively small values of m. To overcome this difficulty, a specific strategy will be proposed in Section 3.4. For the probability π, maximizing l c (µ, π; x 1 , . . . , x n , y 1{q} , . . . , y n{q} ) leads to the following update, removing any combinatorial difficulty:

π {q+1} = ∑ n i=1 G(x i , y i{q} , µ {q} ) ∑ n i=1 A(x i , y i{q} )
.

Unlike EM, the random sequence of parameters (µ, π) {q} generated by SEM-Gibbs does not converge pointwise. Consequently its stopping rule can not rely on difference of the likelihood between two successive iterations. The simplest alternative solution is to stop the SEM-Gibbs algorithm after a given number Q of iterations. After removing a burn in period corresponding to the first B iterations, we retain a point estimate of (µ, π) in the following manner: for each distinct µ values in the sequence µ {q} (q = B, . . . , Q), take the mean πµ of the associated π {q} values and then retain the couple (µ, πµ ) leading to the highest log-likelihood l(µ, πµ ; x 1 , . . . , x n ). Since the log-likelihood calculation suffers from combinatorial issues, we use in addition the following approximation:

l(µ, π; x 1 , . . ., x n ) = - n ∑ i=1 ln ∑ y∈P m 1 p(x i |y; µ, π) p(y|x i ; µ, π) ≈ - n ∑ i=1 ln 1 S S ∑ s=1 1 p(x i |y i,s ; µ, π)
, where all y i,s arise independently from p(y|x i ; µ, π) (s = 1, . . ., S). Section 3.5 will validate on some real and artificial data sets that our SEM-Gibbs strategy leads to very accurate estimates despite its lower computational cost compared to the standard EM algorithm.

Reducing the number of reference ranks µ for EM and SEM-Gibbs

We propose a strategy to reduce (often drastically) the number of possible values for µ in the step M of both EM and SEM-Gibbs. This result relies on Proposition 1 and also on the following proposition.

Proposition 2. Let N x be the number of individuals equal to x ∈ P m among a random sample from ISR(µ, π) of size n. Denoting by h α (π) = #{x : p(N x ≥ N µ ; µ, π) ≥ α} the number of ranks for which the empirical frequency can be greater than or equal (with probability at least α ∈ [0, 1]) to the empirical frequency associated with the theoretical modal rank µ, the following inequality asymptotically holds for any µ ∈ P m and π ∈ [ 1 2 , 1]:

h α (π) ≤ h α π-.
Proof. We know from Proposition 1 that asymptotically π-≤ π. Proposition 5 in C allows to conclude the proof.

The idea is to browse the empirical modal rank in association with some other ranks having high empirical relative frequency. The following strategy can be used only if π-≥ 1 2 . Firstly, h α ( π-) is estimated with a parametric bootstrap [START_REF] Efron | An introduction to the bootstrap, volume 57 of Monographs on Statistics and Applied Probability[END_REF] of M replications from ISR(µ, π-). The key point is that it is independent from µ, so any µ ∈ P m can be used. Then the h α ( π-) most frequent distinct ranks in the sample (x 1 , . . . , x n ) are retained as possible µ values among the potential m!/2 possibilities and are used both as potential initial values µ {0} and also as values to browse during the M step.

The proposed strategy is aimed to significantly reduce the number of candidates for µ. The number of candidates for µ reduces when the size of the observed sample n grows since h α ( π-) p → 1 when n → ∞. So, the browsed ranks are asymptotically reduced to the empirical modal rank which is known to be a consistent estimator of µ.

Table 4 (Column "#µ") of Section 4 illustrates through numerical examples that this procedure effectively reduces the number of possible ranks for µ in comparison to the m!/2 possible values.

Remarks

• Proposition 2 gives asymptotic guaranties on the proposed strategy to reduce the number of candidates for mode µ of the distribution. The practical efficiency of this strategy is illustrated in the next section through simulation studies.

• On the real data sets used in Section 4, we noticed that the proposed strategy to preselect µ is less time-consuming since it takes less than 1% of the the entire process composed of preselection and of estimation (whatever be the estimation algorithm: EM or SEM-Gibbs).

• Note that the selection of the possible ranks should be carried out only once before the start of the EM or the SEM-Gibbs algorithms. This will be the strategy we follow in experiments throughout this paper.

Validation of the SEM-Gibbs algorithm based on simulations

In order to illustrate the estimation accuracy of the SEM-Gibbs approach, we propose to evaluate it both on simulated and real data sets. We introduce here the normalized Kendall distance K(•,

•) = K(•, •)2[m(m -1)
] -1 because its value is between 0 and 1, indicating respectively minimum and maximum disagreement.

• Simulated data. 20 samples of size n = 100 are simulated according to ISR with three values of π (π ∈ {0.6, 0.75, 0.9}) and two values of m (m ∈ {5, 10}). These parameters are then estimated by SEM-Gibbs approach with the following settings: Q = 100 SEM iterations, B = 10 burn in iterations and R = 10 Gibbs iterations inside the SE step. Results are displayed in Table 2. We tried with higher values for Q, B and R, but the results were essentially unchanged. Despite the small number of both Gibbs and SEM iterations, it appears that the accuracy of estimated parameters is very satisfactory, both with low and larger value of m. Note also that estimation of µ is a harder task when π is low, as expected (recall that information on µ is weak for small π).

• Real data. For all data sets considered below in Section 4 we have run EM and SEM-Gibbs, except the Election data set where m = 14 (EM is not numerically available in this case). Tuning parameters of EM and SEM-Gibbs are also the same as in Section 4. However, SEM-Gibbs is run here 10 times at random to evaluate stability of its results for µ, π and l (here l is exact likelihood, not the approximated [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF][START_REF] Lebanon | Non-parametric modeling of partially ranked data[END_REF], it is also linked with some multistage models (see Section 1).

• Football quiz. This quiz consists of ranking four national football teams according to increasing number of wins in the football World Cup:

O 1 = France, O 2 = Germany, O 3 = Brasil, O 4 = Italy
. The correct answer is µ * = (1, 2, 4, 3).

• Cinema quiz. This quiz consists of ranking chronologically the following Quentin Tarantino movies:

O 1 = Inglourious Basterds, O 2 = Pulp Fiction, O 3 = Reservoir Dogs, O 4 = Jackie Brown
. The correct answer is µ * = (3, 2, 4, 1).

• Rugby. This data set is the result of four nations rugby league, from 1883 to 1909 (except years 1888 and 1889 because only three nations were in the tournament, and except years 1886, 1890, 1897, 1898 and 1906 due to tie), which includes

O 1 = England, O 2 = Scotland, O 3 = Ireland and O 4 = Wales.
• Words. [START_REF] Fligner | Distance based ranking models[END_REF] examined the data collected under the auspices of the Graduate Record Examination Board. A sample of 98 college students were asked to rank five words according to strength of association (least to most associated) with the target word "Idea": • Sports. This data set is due to Louis Roussos [START_REF] Marden | Analyzing and modeling rank data[END_REF] who asked 130 students at the University of Illinois to rank seven sports according to their preference in participating:

O 1 = Thought, O 2 =
O 1 = Baseball, O 2 = Football, O 3 = Basketball, O 4 = Tennis, O 5 = Cycling, O 6 = Swimming, O 7 = Jogging.
• Election. The general election system for the Irish House of Parliament is based on ranking of the candidates in order of preference. [START_REF] Gormley | A latent space model for rank data[END_REF] presented and studied the 2002 General Election data set, consisting of the 64,081 rankings of the 14 candidates among which 2,490 are full rankings. Only these 2,490 full rankings are used here. Refer to [START_REF] Gormley | A latent space model for rank data[END_REF] for more details on this data set.

The empirical distribution of the first three data sets (for which the number of objects to rank is 4) is graphically displayed on the left column of Figure 1 in the ranking space (orderings are displayed on each node).

Estimation results

For each data set, the ISR distribution and the Mallows Φ model are estimated. An R package is available on the authors website 1 . For the ISR model, the estimation is carried out using EM algorithm when m ≤ 7, and SEM-Gibbs algorithm is used when m > 7 (Election data set). For ISR the convergence threshold for the growth of log-likelihood in EM algorithm was fixed to 10 -6 and only one initialization of π in [ π-, π+ ] has been used (no change on the results have been observed with several initializations), and B = 30, Q = 100, R = 10 for SEM-Gibbs algorithm is used. For Mallows Φ model, the numerical optimization has been carried out with a quasi-Newton method and the convergence threshold of ISR (10 -6 ). The ISR distribution of the first three data sets is graphically displayed on the right column of Figure 1 for a visual comparison with the empirical distribution. In addition, a χ 2 adequacy test, where the distribution under the null assumption is estimated by boot- strap [START_REF] Efron | An introduction to the bootstrap, volume 57 of Monographs on Statistics and Applied Probability[END_REF] based on 1,000 replications, is performed for both models and for all data sets and the results are displayed in Table 4 (Column " p-value"). Depending on selected threshold on the p-value (typically models with p-value greater than 0.05 are not rejected) we notice that both models can be suitable for some data sets but not for all of them. Moreover, when comparing maximum log-likelihood values (Column "l"; the highest likelihoods are in bold), ISR leads to a greater maximum likelihood than Mallows Φ model for 4 data sets among 6. Consequently, ISR could be a natural competitor to be considered beside other classical models in any rank data analysis. Table 4: ISR and Mallows Φ models estimation results: estimate parameters μ, π (ISR) and λ (Mallows), maximum log-likelihood l, estimated p-value of the χ 2 adequacy test, number of possible µ explored (#µ; For ISR it corresponds to ĥα ( π-) with α = 0.05 and M = 100 replications), lower and upper bounds π-and π+ for π We note that the strategy selecting the number of possible reference ranks to explore (Section 3.4) is effective. Indeed, only one candidate for µ has been selected by this strategy for the three data sets Football, Words and Sports (Column "#µ") and for other data sets the number of candidates is relatively small in comparison with m!. Concerning Mallows Φ model, the estimation of µ is carried out by an empirical iterative local search (in the sense of the Kendall distance) around the modal rank [START_REF] Fligner | Multistage ranking models[END_REF] which appears to be effective.

We discuss now the meaningful interpretation of ISR parameters by further analysing Table 4.

• Football quiz. The estimation of the reference rank µ coincides with the real rank.

In addition, the accuracy level of students knowledge in football is quite high since it is reflected by a high probability π (0.834) of well paired comparison. This underlines that the right answer is, on the whole, known with high level of confidence by this population of students.

• Cinema quiz. The estimation of the reference rank µ does not coincide with the real rank even if the chronological order is correct for three films of four. However, the accuracy level of students knowledge in cinema is equal to π = 0.723. The students seem to have better knowledge in football than in cinema.

• Rugby. The ISR model estimation on the Rugby data set suggests a ranking between these four nations: during this time Scotland were the best, then Wales, England and finally Ireland. But the low value of the probability π (0.681) means that the confidence in this ranking is not very high.

• Words. The high value of π (0.879) for the Fligner and Verducci's Word data set shows that the questioned students overall share the same opinion for the association with the word Idea: Thought is the most associated followed by Theory, Dream, Attention and finally Play which is the least associated.

• Sports. The reference rank (1, 3, 2, 4, 5, 7, 6) estimated for the ISR model reflects a preference of the students at the University of Illinois for collective sports: Baseball, Basketball and Football are at the top three places while individual sports are at the end of the ranking: Cycling, Jogging and Swimming. Tennis, which is intermediate between a collective sport and an individual sport, is rationally ranked between these two groups.

• Election. This last data set calls for an interesting remark: except the first two candidates, it suggests that all other candidates tends to be ordered similar to their initial order in the list. Note that we considered only 4% of the voters in our study, those who "bothered" to rank every one of the 14 candidates. Probably most of them have the common trait that they did not care who is ranked after the first two candidates since, in the Irish voting system, two-thirds of voters see their first choice elected. They may express such an indifference beyond the second rank by choosing the initial order in the list of candidates whereas the remaining 96% of voters prefer not to complete the rest of the list. Notice also that about 45% of these remaining 96% of voters select only one or two candidates. A more comprehensive study including partial ranks would be obviously required to analyse more precisely the results of these elections.

From the Mallows Φ parameters point of view, most results are highly consistent with ISR: main modal ranks are identical and the dispersion parameter λ is also well correlated with π, though λ is more abstract and could be less easy to understand by a practitioner. Only the modal rank of the last two data sets (Sports and Election) differs: for the Sports data set, Mallows Φ model classifies Tennis inside the collective sports collection instead of being put at the borderline of collective and individual sports, as ISR does; For the Election data set, Mallows Φ model provides the same first two candidates as ISR but in the reverse order and, in addition, all other candidates are very different between both methods.

Discussion

In this paper we propose that rank data could be considered as the result of a paired comparisons sort algorithm, where the possibility of wrong comparisons exists and occurs randomly according to a Bernoulli model. It opens a new way to propose many distributions on rankings, all of them benefiting from very meaningful parameters (the reference rank µ and the probability π of good paired comparison). In order to minimize the number of paired comparisons, and consequently number of potential wrong comparisons, the insertion sort algorithm has been kept in this paper for its optimality when m ≤ 10 since it is expected to be a frequent limit in "human rankings". The resulting distribution, the so-called ISR, has been established and many desirable properties were being pointed out. In addition, the latent variable interpretation of model allows the derivation of specific EM and SEM-Gibbs algorithms which can be easily accelerated by drastically reducing the number of potential reference ranks µ to consider. Although our approach is not able to deal with very large data sets like web pages rankings or some kinds of biological data, these computational gains allow to deal with usual data sets of moderate size typically provided by a "human ranking" process.

In fact, the ISR model can be considered as the precursor to a wider family of ranking models. The insertion sort algorithm was kept in this work from some optimality arguments. However, it is possible to easily set up a new model by changing it. For instance, we can use a selection sort algorithm instead, or any other standard algorithms [?, see]]Knu1973. Obviously, properties of the new model obtained in this way would need to be established again. In each case, it could also be possible to weaken the assumption that π remains constant all along the sorting process, for instance to model the tiredness of the judge during the sorting process.

Another interesting prospect initialized by the present work is the possibility to include some information about the initial ranking y in the model and its corresponding estimation. Indeed, in questionnaires this initial order is often known and it is useful information which can be naturally used by our class of models. It is also possible to consider some more diffuse information about y, for instance to ignore the exact y value but to know that all y are the same for all questionnaires (realistic situation for many ranks coming from quiz studies), or other realistic variants [START_REF] Biernacki | Modèles génératifs de rangs relatifs à un algorithme de tri par insertion[END_REF]. In the same spirit, the model is also flexible enough to take into account the usual behaviour that some individuals may rank the items from the beginning to the end and others to do the opposite.

Although the ISR is unimodal (as many other distributions for ranks), multimodality can be easily taken into account through mixture of ISR distributions after leading a specific identifiability study. For instance, we can think that in our football quiz, girls and boys responses will probably not follow the same distribution, as it is suggested by very low estimated p-value [START_REF] Jacques | model-based clustering for rank data based on an insertion sorting algorithm[END_REF]. This extension is natural since several mixtures of rank data models have already been considered with success to treat heterogeneity of rank population: mixture of Mallows Φ models [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF], mixture of Plackett-Luce or Benter models [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF][START_REF] Gormley | Exploring voting blocs within the irish electorate: A mixture modeling approach[END_REF][START_REF] Gormley | A mixture of experts model for rank data with applications in election studies[END_REF] and more recently mixtures of weighted distance-based models [START_REF] Lee | Mixtures of weighted distance-based models for ranking data with applications in political studies[END_REF].

Finally, there is also a need to adapt our models to other situations than full rank data. This approach needs to be extended to other types of ranks, frequently encountered in practice, as partially or incomplete ranked data [START_REF] Lebanon | Cranking: Combining rankings using conditional models on permutations[END_REF][START_REF] Kidwell | Visualizing incomplete and partially ranked data[END_REF][START_REF] Lebanon | Non-parametric modeling of partially ranked data[END_REF], which would be very useful to further analyse the Irish election data set, tied data or even ranks resulting from multiple preference responses.

We conclude the proof by noticing that p(x ( j) |y; µ, π) = p(x ( j) |x ( j-1) , y; µ, π)p(x ( j-1) |y; µ, π), from the following implied relationship between events: x ( j) ⇒ x ( j-1) .

C. Mathematical statement of the ISR properties and related proofs

In the following, composition τ • x will be noted shortly τx for any τ and x in P m .

Proposition 3. (Uniformity for π = 1 2 .) For all x, µ ∈ P m , p(x; µ, 1 2 ) = m! -1 . Proof. Let e be the identity permutation of P m . Firstly using Lemma 3 of D and then using the fact that p(.|e; µ, 1 2 ) is a probability distribution on P m , we have p(x; µ, The last equality comes from the fact that A(µ, y ′ ) = G(µ, y ′ , µ).

Corollary 1. (Anti-mode μ.)For all x = μ ∈ P m and π > 1 2 , p( μ; µ, π) < p(x; µ, π). The proof is symmetrical to that of Proposition 4. Proposition 5. (Link between µ and π.) For all x, µ ∈ P m , p(µ; µ, π)p(x; µ, π) is an increasing function of π ≥ 1 2 . Proof. Noting ∆(π) = p(µ; µ, π)p(x; µ, π), ∂ ∆(π)/∂ π, we can written

∂ ∆(π) ∂ π = 1 m! ∑ y∈P m
A(µ, y)π A(µ,y)-1 -G(x, y, µ)π G(x,y,µ)-1 (1 -π) A(x,y)-G(x,y,µ) + c, where c is a non-negative term independent from π. Since π ≥ 1 2 , we deduce that G(x, y, µ)π G(x,y,µ)-1 (1 -π) A(x,y)-G(x,y,µ) ≤ G(x, y, µ)π A(x,y)-1 .

Using the fact that A(µ, y) ≥ G(x, y, µ), we deduce that ∂ ∆(π)/∂ π ≥ 0. Proposition 6. (Symmetry.) For all x, µ ∈ P m and all π ∈ [0, 1], p(x; μ, 1-π) = p(x; µ, π).

Proposition 1 . 1 0

 11 Denoting by f 0 the empirical modal relative frequency, the interval [ π-, π+ ] asymptotically contains π where π-= f 1 mand π+ = f 2 m(m-1) 0 .

  Play, O 3 = Theory, O 4 = Dream and O 5 = Attention.

Figure 1 :

 1 Figure 1: Empirical (left) and estimate ISR (right) distributions for the three data sets where m = 4. The area of the dots is proportional to the corresponding probability.

1 2 2 ) = 1 . 4 .

 1214 ) ∝ ∑ y∈P m p(x|y; µ, 1 2 ) = ∑ y∈P m p(y -1 x|y -1 y; µ, 1 2 ) = ∑ y∈P m p(y -1 x|e; µ, 1 Proposition (Mode µ.) For all x = µ ∈ P m and π > 1 2 , p(µ; µ, π) > p(x; µ, π).Proof. Using the fact that {π > 1 2 ⇔ π > 1 -π}, x = µ and then Lemma 2, we obtain:m! p(x; µ, π) < ∑ y∈P m π A(x,y) = ∑ y∈P mπ A((µx -1 )x,(µx -1 )y) = ∑ y ′ ∈P m π A(µ,y ′ ) = m! p(µ; µ, π).

  α = 0.1 for the Sports data set to avoid to many µ due to the small π value ‡ Useless for SEM-Gibbs

				723	-111.94	0.042 14 0.630 0.794
		Mallows	(4, 3, 2, 1)	0.628	-112.12	0.029	2	-	-
	Rugby	ISR	(2, 4, 1, 3)	0.681	-58.68	0.538 12 0.585 0.765
		Mallows	(2, 4, 1, 3)	0.528	-58.33	0.395	2	-	-
	Words	ISR	(2, 5, 4, 3, 1)	0.879	-275.43	0.001	1 0.762 0.897
		Mallows	(2, 5, 4, 3, 1)	1.431	-251.27	0.019	1	-	-
	Sports	ISR	(1, 3, 2, 4, 5, 7, 6)	0.564 -1102.12	0.999 2 † 0.534 0.836
		Mallows (1, 3, 4, 2, 5, 6, 7)	0.083 -1102.84	1	11	-	-
	Election ISR	(13, 4, 1, 2, 3, 5, 6	0.682 -48329.76 0.999	6	‡	‡
			7, 8, 9, 10, 11, 12, 14)
		Mallows (4, 13, 2, 5, 1, 14, 7	0.164 -60157.38 0.999 38	-	-
			6, 10, 8, 9, 12, 3, 11)

†

http://math.univ-lille1.fr/∼jacques/soft.html

Acknowledgment

We thank Pr. Brendan Murphy and Dr. Claire Gormley for providing the Irish election data set. We also thank the reviewers for their contributions to greatly improve this work.

Table 2: Simulated data to validate SEM-Gibbs algorithm (20 replicates of size n = 100): π and K( μ, µ) indicate respectively the mean of the π estimate (in parenthesis its standard deviation) and the mean of the standardized Kendall distance between the µ estimate and the true one (in parenthesis its standard deviation). m = 5 m = 10 π π K( μ, µ) π K( μ, µ) 0.6 0.572 0.210 0.566 0.180 (0.032) (0.141) (0.250) (0.063) 0.75 0.759 0.000 0.700 0.056 (0.021) (0.000) (0.022) (0.036) 0.9 0.895 0.000 0.884 0.002 (0.015) (0.000) (0.008) (0.007) Table 3: Real data to validate SEM-Gibbs algorithm (10 independent runs of SEM-Gibbs): μEM , μSEM-gibbs , πEM , πSEM-gibbs , l EM , l SEM-Gibbs denote respectively µ, π and l estimated by EM and by SEM-Gibbs. Column "mean" displays the mean of each statistics on µ, π and l given at top of the table, over 10 runs of SEM-Gibbs. Columns "best" and "worst" give the value of these statistics with parameters obtained respectively with the best and the worst likelihood over the 10 runs of SEM-Gibbs. 3: in column "best", the best SEM-Gibbs (according to l) and EM coincide; in columns "means" and "worst" SEM-Gibbs is very close to EM. Thus, it confirms the good behaviour of SEM-Gibbs we validated previously with simulated data.

Numerical illustration

Presentation of the six real data sets

The ISR distribution is now compared to Mallows Φ model on six real data sets: two general knowledge quizzes (the answers of the 40 students being questioned are in E), four nations rugby league rankings, Fligner and Verducci's words associations rankings [START_REF] Fligner | Distance based ranking models[END_REF], Louis Roussos's sports rankings [START_REF] Marden | Analyzing and modeling rank data[END_REF] and 2002 General Irish Election data set [START_REF] Gormley | A latent space model for rank data[END_REF].

Mallows Φ model has been chosen as the reference model since, in addition to being one of the model based on paired comparison which has been the most studied [START_REF] Mallows | Non-null ranking models[END_REF][START_REF] Fligner | Distance based ranking models[END_REF][START_REF] Fligner | Multistage ranking models[END_REF],

A. Additional notations useful for the proofs

Additional notations are needed in the proofs displayed in the appendices below. All these notations are illustrated in detail through an example in Table 5. In addition, to clarify them, the calculation for step j = 2 (for instance) is detailed at the end of the present appendix.

In the following, j = 1, . . . , m denotes the step in the sorting algorithm consisting in ranking the object O y j .

A(x, y) = ∑ m j=1 A j (x, y) and A j (x, y) = A - j (x, y) + A + j (x, y) designate the total number of all paired comparisons respectively for the whole process and for the step j. In this definition, A - j (x, y) is the number of all comparisons of the current object O y j with the objects already ranked (according to x) before it (if they exist) and A + j (x, y) indicates if the current object O y j is compared, at the j step of the sorting, with the object ranked in x just after it. Formally, A - j (x, y) corresponds to the cardinal of A - j (x, y) = {i : x -1 y i < x -1 y j , 1 ≤ i < j} which is the set of the indices of the presentation order y for which the already sorted objects O y 1 , . . . , O y j-1 are ranked in x before the current object O y j . In a similar way, A + j (x, y) is the cardinal of the set A + j (x, y) = {i : i = arg min 1≤i ′ < j {i ′ : x -1 y i ′ > x -1 y j }} which corresponds to the index of the rank y designating the object sorted in x just after O y j among the already sorted objects O y 1 , . . . , O y j-1 , if it exists. This set has at most one element.

G(x, y, µ) = ∑ m j=1 G j (x, y, µ) and G j (x, y, µ) = G - j (x, y, µ) + G + j (x, y, µ) are the total number of good paired comparisons respectively for the whole process and for the step j. Formally, G - j (x, y, µ) = ∑ i∈A - j (x,y) δ y i y j (µ) is the number of good comparisons (according to µ) of the current object O y j with the objects already ranked before it (if they exist), where

is the indicator of good comparison (according to µ) of the current object O y j with the object already ranked just after it (if it exists).

Detail of step j = 2 to be read in conjunction with Table 5

• A - 2 = {}: the new object 3 (to be sorted) has been compared to the already sorted object 1 but 1 is not sorted before 3 in the final ordering x; Table 5: An example to illustrate the notations with µ = (1, 2, 3), y = (1, 3, 2) and x = (3, 1, 2). The notation x ( j) , defined in B, means the ranking of the j first objects in y in the order imposed by x step j unsorted sorted

the new object 3 (to be sorted) has been compared to the already sorted object 1 and 1 is sorted after 3 in the final ordering x;

• A - 2 = 0 and A + 2 = 1: number of comparisons listed respectively in A - 2 and A + 2 ;

• A 2 = 1: total number of comparisons to sort the new object 3.

• G - 2 = 0 and G + 2 = 1: number of good comparisons listed respectively in A - 2 and A + 2 ;

• G 2 = 1: total number of good comparisons to sort the new object 3.

B. Building the ISR distribution

The goal of this appendix is to prove that (2.1) corresponds to the stochastic insertion sort algorithm with probability π of good paired comparison, and independence between the paired comparisons.

Proof. Let x ( j) be the ordering of the first j (1 ≤ j ≤ m) objects in y in the order imposed by x (so x (m) = x). An example of this notation is in Table 5. Thus, there exists following relationship between x ( j) and x ( j-1) :

Equation (2.1) is now proved by induction on j. It is true for j = 1 while there is only one object y 1 to sort: p(x (1) |y; µ, π) = 1. Since the result of the ranking x ( j) from x ( j-1) is the result of A j (x, y) independent Bernoulli experiments of parameter π, then, conditionally to x ( j-1) , the probability of x ( j) is p(x ( j) |x ( j-1) , y; µ, π) = π G j (x,y,µ) (1 -π) A j (x,y)-G j (x,y,µ) .

Proof. Using Lemma 4, we can write: p(x; μ, 1 -π) ∝ ∑ y∈P m π A(x,y)-(A(x,y)-G(x,y,µ)) (1 -π) A(x,y)-G(x,y,µ) ∝ p(x; µ, π).

Proposition 7. (Identifiability.) The ISR distribution is identifiable since π > 1 2 . Proof. The identifiability problem can concern parameters π and/or µ.

• First, there exists no couple (µ, µ ′ ) ∈ P 2 m with µ = µ ′ such that p(x; µ, π) = p(x; µ ′ , π) for any x ∈ P m and any π > 1 2 . Indeed, choosing x = µ, from Lemma 5 we have p(µ; µ, π) = p(µ; µ ′ , π).

• Second, for a given µ ∈ P m , assume there exists π = π ′ such that p(x; µ, π) = p(x; µ, π ′ ) for any x ∈ P m . In particular, for x = µ, in the proof of Lemma 5 we obtained that G(x, y, x) = A(x, y), thus ∑ y∈P m π A(µ,y) = ∑ y∈P m π ′ A(µ,y) . The strictly increasing function p → π n on the interval [ 1 2 , 1] for all n ∈ N * ensures that π = π ′ . • Assume finally there exists (µ, µ ′ ) ∈ P 2 m with µ = µ ′ and π < π ′ such that p(x; µ, π) = p(x; µ ′ , π ′ ) for any x ∈ P m . In the proof of Lemma 5, it is also obtained that G(x, y, µ) < A(x, y) when x = µ, thus p(x|y; µ, π) < π A(x,y) < π ′ A(x,y) = p(x|y; x, π ′ ), and then by averaging over all y in P m gives p(x; µ, π) < p(x; x, π ′ ). Choosing x = µ ′ ensures the identifiability of the ISR model.

D. Lemmas

Lemma 1. Let ẽ = (2, 1, 3, . . ., m) be the permutation inverting the first two elements. For all x, y, µ ∈ P m and π ∈ [0, 1], p(x|y; µ, π) = p(x|y ẽ; µ, π).

Proof. We use notations x ( j) that have already been introduced in B. The key point of the proof is to notice that the first two objects in y lead to the same paired comparison at the second step of the sorting process whatever is their order in y, so p(x (2) |y ẽ, π) = p(x (2) |y, π). Combining this result with the fact that p(x|x (2) , y ẽ, π) = p(x|x (2) , y, π), since ẽ only affects the first two objects, this concludes the proof. Lemma 2. For all x, y, τ ∈ P m , A(x, y) = A(τx, τy).

Proof. First we prove that A - j (x, y) = A - j (τx, τy). For any j = 1, . . ., m, we have (notice that i is always such that 1 ≤ i < j)

Using the fact that A + j (x, y) = 1{A - j (x, y) + 1 ≤ j -1} we also deduce that A + j (x, y) = A + j (τx, τy). Consequently, A j (x, y) = A j (τx, τy) and, so, A(x, y) = A(τx, τy).

Lemma 3. For all x, y, µ, τ ∈ P m , p(x|y; µ, 1 2 ) = p(τx|τy; µ, 1 2 ).

Proof. When π = 1 2 , we obtain by using Lemma 2 p(τx|τy; µ, 1 2 ) = 1 2 A(τx,τy) = 1 2 A(x,y) = p(x|y; µ, 1 2 ).

Lemma 4. For all x, y, µ ∈ P m G(x, y, μ) = A(x, y) -G(x, y, µ).

Proof. Let ē be the permutation of total inversion previously introduced in Section 2.3 and i, i ′ = 1, . . . , m, i = i ′ . We first prove that G

and then

In a similar manner, we can prove that G + j (x, y, μ) = A + j (x, y) -G + j (x, y, µ). The proof follows immediately from these two results. Lemma 5. For all x, µ ∈ P m , x = µ and π > 1 2 , p(x; µ, π) < p(x; x, π).

Proof. First note that G(x, y, µ) < A(x, y) for µ = x. Since {π > 1 2 ⇔ 1 -π < π}, we deduce for µ = x that p(x|y; µ, π) < π A(x,y) . Also note that G(x, y, x) = A(x, y), thus p(x|y; x, π) = π A(x,y) . Consequently, we have p(x|y; µ, π) < p(x|y; x, π) and the proof is concluded by averaging over all possible presentation orders y in P m . Lemma 6. For all µ, y ∈ P m , m -1 ≤ A(µ, y) ≤ m(m -1)/2.

Proof. Left bound: there is no comparison when the first element arises and at least one comparison for each of the m-1 other elements. Right bound: there is still no comparison when the first element arises and at most j -1 comparisons when for the jth new object to rank, so A(µ, y) ≤ ∑ m j=1 ( j -1) = m(m -1)/2.

E. Quiz data sets


Julien Jacques
email: julien.jacques@polytech-lille.fr

Cristian Preda
email: cristian.preda@polytech-lille.fr

Model-based clustering for multivariate functional data

Keywords: Multivariate functional data, density approximation, model-based clustering, multivariate functional principal component analysis, EM algorithm

Model-based clustering for multivariate functional data

Introduction

Functional data analysis or "data analysis with curves" is an active topic in statistics with a wide range of applications. New technologies allow to record data with accuracy and at high frequency (in time or other dimension), generating large volume of data. In medicine one has growth curves of children and patient's state evolution, in climatology one records weather parameters over decades, chemometric curves are analysed in chemistry and physics (spectroscopy) and special attention is paid to the evolution of indicators coming from economy and finance. See [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF] for more details.

The statistical model underlying data represented by curves is a stochastic process with continuous time, X = {X(t)} t∈[0,T ] . Most of the approaches dealing with functional data consider the univariate case, i.e. X(t) ∈ R, ∀t ∈ [0, T ], a path of X being represented by a single curve. Despite its evident interest, the multidimensional case, X = {X(t)} t∈[0,T ] with X(t) = (X 1 (t), . . . , X p (t)) ′ ∈ R p , p ≥ 2 is rarely considered in literature. In this case a path of X is represented by a set of p curves. The dependency between these p measures provides the structure of X. One finds in [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF] a brief example of bivariate functional data, X(t) = (X 1 (t), X 2 (t)) ′ ∈ R 2 , as a model for gait data (knee and hip measures) used in the context of functional principal component analysis (FPCA) as an extension of the univariate case. For a more theoretical framework, we must go back to the pioneer works of [START_REF] Besse | Etude descriptive d'un processus[END_REF] on random variables in a general Hilbert space. [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF] provides a complete analysis of multivariate functional data from the point of view of factorial methods (principal components and canonical analysis).

In this paper we consider the problem of clustering multivariate functional data. Cluster analysis aims to identify homogeneous groups of data without using any prior knowledge on the group labels of data. When clustering functional data, the main difficulty is due to infinite dimensional space that data belong to. Consequently, most of clustering algorithms for functional data consists of a first step of transforming the infinite dimensional problem into a finite dimensional one and a second step using a clustering method designed for finite dimensional data. Recently, several new approaches considers the k-means algorithm applied to a B-spline fitting [START_REF] Abraham | Unsupervised curve clustering using B-splines[END_REF], to defined principal points of curves [START_REF] Tarpey | Clustering functional data[END_REF] or to a truncation of the Karhunen-Loeve expansion [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF]. [START_REF] Sangalli | K-means alignment for curve clustering[END_REF] also use a k-means algorithm to cluster misaligned curves. As in the finite dimensional setting, where Gaussian model-based clustering generalizes the k-means algorithm, some other works introduce more sophisticated model-based techniques: [START_REF] James | Clustering for sparsely sampled functional data[END_REF] define an approach particularly effective for sparsely sampled functional data, [START_REF] Ray | Functional clustering by Bayesian wavelet methods[END_REF] propose a nonparametric Bayes wavelet model for curves clustering based on a mixture of Dirichlet processes, Frühwirth-Schnatter and Kaufmann (2008) build a specific clustering algorithm based on parametric time series models, [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] extend the high-dimensional data clustering algorithm (HDDC, [START_REF] Bouveyron | High Dimensional Data Clustering[END_REF]) to the functional case and [START_REF] Jacques | Model-based clustering of functional data[END_REF] build a model-based clustering based on an approximation of the notion of functional variable density.

The case of multivariate functional data is more rarely considered in literature: Singhal and Seborg (2005) and [START_REF] Ieva | ECG signal reconstruction, landmark registration and functional classification[END_REF] use a k-means algorithm based on specific distances between multivariate functional data, whereas [START_REF] Kayano | Functional Cluster Analysis via Orthonormalized Gaussian Basis Expansions and Its Application[END_REF] consider Self-Organizing Maps based on the coefficients of multivariate curves into orthonormalized Gaussian basis expansions. [START_REF] Tokushige | Crisp and fuzzy k-means clustering algorithms for multivariate functional data[END_REF] extend crisp and fuzzy k-means algorithms for multivariate functional data by considering specific distance between functions, but applied their algorithms only on univariate functional data.

In the finite dimensional setting, model-based clustering algorithms consider that data is sampled from a mixture of probability densities. This is not directly applicable to functional data since the notion of probability density generally does not exist for functional random variable [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. Consequently, model-based clustering algorithms assume a parametric distribution on some finite set of coefficients characterizing the curves. In [START_REF] Jacques | Model-based clustering of functional data[END_REF], the authors use the density surrogate defined in [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF] to build a model-based clustering for univariate functional data. This density surrogate, based on the truncation of the Karhunen-Loeve expansion, relies on the probability density of the principal components of the curves [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF], which is assumed to be Gaussian. In this paper we propose an extension of [START_REF] Jacques | Model-based clustering of functional data[END_REF] approach to multivariate functional data. For this, we firstly introduce principal component analysis for multivariate functional data and assume a cluster-specific Gaussian distribution for the principal component scores. The elements derived from FPCA are estimated using approximation of the multivariate curves into a finite dimensional functional space. The number of principal components used in the density surrogate as well as the computation of the principal component scores are cluster specific.

The main advantage of our model is its ability to take into account dependency between the p curves of the multidimensional data, thanks to the principal component analysis for multivariate functional data.

The paper is organized as follows. Section 2 introduces principal components analysis for multivariate functional data. Estimation and approximation details are provided and the task of normalizing the curves is discussed. Section 3 defines an approximation of the probability density for multivariate functional random variable. The model-based clustering approach and parameter estimation via an EM-like algorithm are presented in Section 4. Comparisons with existing methods on simulated and real datasets are presented in Section 5, and a discussion concludes the paper in Section 6.

Principal component analysis for multivariate functional data (MFPCA)

Principal component analysis for multivariate functional data has already been suggested in [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF] and [START_REF] Berrendero | Principal components for multivariate functional data[END_REF]. In [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF] the authors propose to concatenate the observations of the functions on a fine grid of points (or the coefficients in a suitable basis expansion) into a single vector and then to perform a standard principal component analysis (PCA) on these concatenated vectors. When a basis expansion is used, this method forces to consider only orthonormal basis since the metric induced by the scalar product between the basis functions is not taken into account. In [START_REF] Berrendero | Principal components for multivariate functional data[END_REF], the authors propose to summarize the curves with functional principal components instead of scalar ones as in usual FPCA. For this purpose, they carry out classical PCA for each value of the domain on which the functions are observed and suggest an interpolation method to build the functional principal components. Our approach is closely related to [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF] but in addition, we take into account the possible use of non orthonormal basis. In particular, our method allows to use different basis for each dimension of the multivariate curves. Let X 1 , . . . , X n be an i.i.d. sample of X. The observation of X 1 , . . . , X n provides a 3 set of n p-variate curves, called multivariate functional data. From this set of multivariate curves, one can be interested in optimal representation of curves in a function space of reduced dimension (principal component analysis), or in clustering, by determining an optimal partition of the observed curves with respect to some distance or homogeneity criterion. In order to address these two questions in a formal way, we need the hypothesis that considers X = (X 1 , . . . , X p ) ′ as a L 2 continuous stochastic process:

∀t ∈ [0, T ], lim h→0 E X(t + h) -X(t) 2 = lim h→0 T 0 p ℓ=1 E (X ℓ (t + h) -X ℓ (t)) 2 dt = 0.
Notice that L 2 -continuity of X implies L 2 -continuity of each component of X, i.e. X ℓ is a L 2 -continuous stochastic process for all ℓ = 1, . . . , p. The L 2 -continuity is a quite general hypothesis, as most of the real data applications satisfies this one. Let

µ ℓ = {µ ℓ (t) = E[X ℓ (t)]} t∈[0,T ] denotes the mean function of X ℓ (1 ≤ ℓ ≤ p) and µ = (µ 1 , . . . µ p ) ′ = E[X],
denotes the mean function of X. The covariance operator V of X:

V : L 2 ([0, T ]) p → L 2 ([0, T ]) p f V -→ Vf = T 0 V (•, t)f (t)dt,
is an integral operator with kernel V defined by:

V (s, t) = E [(X(s) -µ(s)) ⊗ (X(t) -µ(t))] , s, t ∈ [0, T ]
where ⊗ is the tensor product on R p . Thus, for any s, t ∈ [0, T ], V (s, t) is a p × p matrix with elements V (s, t)[j, ℓ] = Cov(X j (s), X ℓ (t)), j, ℓ = 1, . . . p.

Principal components analysis of X

Under the hypothesis of L 2 -continuity, V is an Hilbert-Schmidt operator, i.e compact, self-adjoint and such that j≥1 λ 2 j < +∞. The spectral analysis of V provides a countable set of positive eigenvalues {λ j } j≥1 associated to an orthonormal basis of eigenfunctions {f j } j≥1 , f j = (f 1 j , . . . , f p j ), called principal factors and are solutions of:

Vf j = λ j f j , (1)
with λ 1 ≥ λ 2 ≥ . . . and

T 0 p ℓ=1 f ℓ j (t)f ℓ j ′ (t)dt = 1 if j = j ′ and 0 otherwise.
The principal components C j of X are zero-mean random variables defined as the projections of X on the eigenfunctions of V:

C j = T 0 X(t) -µ(t), f j (t) R p dt = T 0 p ℓ=1 (X ℓ (t) -µ ℓ (t))f ℓ j (t)dt.
Similar to the univariate setting, the principal components {C j } j≥1 are zero-mean uncorrelated random variables with variance λ j , j ≥ 1. [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF] shows that the following Karhunen-Loeve expansion holds in multidimensional context:

X(t) = µ(t) + j≥1 C j f j (t), t ∈ [0, T ]. ( 2
)
Principal components and principal factors of MFPCA have the same interpretation as in the functional univariate case. The truncation of (2) at the first q terms provides a reduced dimensional space where classical tools (clustering, regression, ...) from multivariate analysis can be used to describe X.

Computational methods for MFPCA

Let x 1 , . . . , x n , with x i = (x 1 i , . . . , x p i ), be the observation of the sample X 1 , . . . , X n . The estimators for µ(t) and V (s, t), for s, t ∈ [0, T ], are:

μ(t) = 1 n n i=1 x i (t) and V (s, t) = 1 n -1 n i=1 (x i (s) -μ(s)) ⊗ (x i (t) -μ(t)).
In [START_REF] Deville | Méthodes statistiques et numériques de l'analyse harmonique[END_REF] it has been shown that μ and V converges to µ and V in L 2 -norm with convergences rate of O(n -1/2 ). Often in practice, data are observed at discrete time points and with some noise. In order to get the functional feature of data, smoothing and interpolation methods are used considering the true curve belongs to a finite dimensional space spanned by some basis of functions. This approximation also reduces the eigen-analysis problem (1) to the one in finite dimensional setting. Let assume that each curve

x ℓ i (1 ≤ i ≤ n, 1 ≤ ℓ ≤ p)
can be expressed as a linear combination of basis functions {φ j ℓ } j=1,m ℓ :

x ℓ i (t) = m ℓ j=1 a iℓj φ j ℓ (t), t ∈ [0, T ]. (3)
This can be written with the matrix formulation

x i (t) = Φ(t)a ′ i
with a i = (a i11 , . . . , a i1m 1 , a i21 , . . . , a i2m 2 , . . . , a ip1 , . . . , a ipmp ) being the vector of the basis expansion coefficients, and

Φ(t) =     φ 1 1 (t) . . . φ m 1 1 (t) 0 . . . 0 0 . . . 0 φ 1 2 (t) . . . φ m 2 2 (t) 0 . . . 0 . . . 0 . . . 0 φ 1 p (t) . . . φ mp p (t)     .
Let à be the n × p ℓ=1 m ℓ -matrix, whose rows are the vectors a i , and M(t) the n × p matrix with values x ℓ i (t) of functions

x ℓ i at times t ∈ [0, T ] (1 ≤ i ≤ n, 1 ≤ ℓ ≤ p)
. With these notations, we have

M(t) = ÃΦ ′ (t). ( 4
)
Under the basis expansion assumption (3), the estimator V of V , for all s, t ∈ [0, T ], is given by:

V (s, t) = 1 n -1 (M(s) -μ′ (s)) ′ (M(t) -μ′ (t)) = 1 n -1 Φ(s)A ′ AΦ ′ (t), (5)
where M(s) -μ′ (s) means that the row vector μ′ (s) is subtracted to each row of M(s), and A = (I n -1I n (1/n, . . . , 1/n)) Ã where I n and 1I n are respectively the identity n × nmatrix and the unit column vector of size n.

From ( 1) and ( 5), each principal factor f j belongs to the linear space spanned by the basis Φ:

f j (t) = Φ(t)b ′ j (6) with b j = (b j11 , . . . , b j1m 1 , b j21 , . . . , b j2m 2 , . . . , b jp1 , . . . , b jpmp ).
Using the estimation V of V , the eigen problem (1) becomes

T 0 V (s, t)f j (t)dt = λ j f j (s),
which, by replacing V (s, t) and f j (s) by their expressions given in ( 5) and ( 6), is equivalent to

T 0 1 n -1 Φ(t)A ′ AΦ ′ (t)f ′ j (t)ds = λ j Φ(s)b ′ j , (7)
⇔ 1 n -1 Φ(s)A ′ A T 0 Φ ′ (t)Φ(t)dt W b ′ j = λ j Φ(s)b ′ j , (8)
where W = T 0 Φ ′ (t)Φ(t)dt is defined as the symmetric block-diagonal p ℓ=1 m ℓ × p ℓ=1 m ℓ -matrix of the inner products between the basis functions. Since ( 8) is true for all s, we have:

1 n -1 A ′ AWb ′ j = λ j b ′ j .
By defining u j = b j W 1/2 , the multivariate functional principal component analysis is reduced to the usual PCA of the matrix

1 √ n-1 AW 1/2 : 1 n -1 W 1/2 ′ A ′ AW 1/2 u ′ j = λ j u ′ j .
The coefficient b j , j ≥ 1, of the principal factors f j are obtained by b j = (W 1/2 ′ ) -1 u ′ j , and the principal component scores, are given by

C j = AWb ′ j j ≥ 1.
Note that the principal components scores C j are also the solutions of the eigenvalues problem:

1 n -1 AWA ′ C j = λ j C j .

Normed principal component analysis

When the X ℓ 's components of X (1 ≤ ℓ ≤ p) are of different natures (different measure units for example), the question of normalizing data occurs naturally. It is well known that the principal components are defined as maximizing the variance with respect to the total variance T race(V) = j≥1 λ j . Since,

T race(V) = T 0 p j=1
Var(X j (t))dt, it is clear that components X j with large variances plays an important role in defining the principal components. This source of variability is in general not interesting and hides more interesting features of the data structure. For this reason, except in particular situations (same scale for all X j 's, for example), normalization is suitable. As in the classical framework of PCA, this is done by introducing some metrics. One way to introduce a metric in MFPCA is to consider it as a canonical analysis [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF], in which the principal components are defined as solutions of the following eigen problem:

T 0 P t (C j )dt = λ j C j , j ≥ 1 (9)
where P t is the orthogonal projection operator associated with X defined as

P t (C j ) =< X(t), [V (t, t)] -1 E[X(t)C j ] > R p . (10)
Combining ( 10) and ( 9) one obtains

C j = T 0 < X(t) -µ(t), f j (t) > R p dt (11)
where f j is the solution of the eigenvector problem

T 0 [V (s, s)] -1 V (s, t)f j (t)dt = λf (s). ( 12
)
Clearly, [V (s, s)] -1 must exist for each s ∈ [0, T ]. Under this hypothesis, the principal factors of the normalized MFPCA are the eigenfunctions of the integral operator with kernel [V (s, s)] -1 V (s, t) as in ( 12). The Karhunen-Loeve expansion of X becomes

X(t) = µ(t) + ∞ j=1 C j [V (t, t)]f j (t),
where the principal components C j , defined by ( 11), have zero mean and variance λ j .

Normalization in practice. Observe that if R(t, t) is the square root of the matrix V (t, t), i.e. V (t, t) = R(t, t)R(t, t) ′ , then the MFPCA of X with metric V (t, t) is equivalent to the MFPCA of Y defined by Y (t) = R(t, t) -1 X(t),
with identity metric as in Section 2.1. In practice, if X is observed at times t 1 , . . . , t r , r > 1, then Y is defined from X as

Y (t i ) = R(t i , t i ) -1 X(t i ), i = 1, . . . , r
and approximation and estimation methodology in Section 2.2 is applied to Y .

Approximation of the density for multivariate functional data

As the notion of probability density is not well defined for functional data [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF], we can use an approximation of the density based on the Karhunen-Loeve expansion (2). Considering the principal components indexed upon the descending order of the eigenvalues (λ 1 ≥ λ 2 ≥ . . .), and denoting X (q) as the approximation of X by truncating (2) at the q first terms, q ≥ 1, we have

X (q) (t) = µ(t) + q j=1 C j f j (t). ( 13
)
Then, X (q) is the best approximation of X, under the mean square criterion, among all the approximations of the same type (linear combination of q deterministic functions of t with random coefficients, [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF]).

Based on the approximation X (q) of X, [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF] shows that, in the case p = 1, the probability of X belonging to a ball of radius h centred at x ∈ L 2 ([0, T ]) p can be written as

log P ( X -x ≤ h) = q j=1 log f C j (c j (x)) + ξ(h, q(h)) + o(q(h)), (14)
where f C j is the probability density of C j and c j (x) is the jth principal component score of x,

c j (x) =< f j , x > L 2 ([O,T ]) p .
The functions q and ξ are such that q grows to infinity when h tends to zero and ξ depends only on h and q(h).

The equality ( 14) suggests the use of the multivariate probability density of the vector C (q) = (C 1 , . . . , C q ) as a surrogate f (q) X for the "density" of X:

f (q) X (x) = q j=1 f C j (c j (x)). ( 15
)
Jacques and Preda (2012) use successfully a similar density surrogate for the clustering of univariate functional data.

A model based-clustering for multivariate functional data

The aim of model-based clustering is to identify homogeneous groups of data sampled from a mixture density model. In this section, we build a mixture model based on the surrogate (15) for the density of X.

Let us consider that there exists a latent group variable Z, of K categories (K groups), Z = (Z 1 , . . . , Z K ) ∈ {0, 1} K such that Z g = 1 indicates that the multivariate curve X belongs to the cluster g, 1 ≤ g ≤ K, and Z g = 0 otherwise. For each i = 1, . . . , n, let Z i be the group indicator corresponding to X i .

In the following we suppose that X |Z g =1 is such that each X ℓ |Z g =1 is a zero-mean stochastic process (1 ≤ ℓ ≤ p). The number K of groups is assumed to be known. In the contrary case, an approximation of the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF], based on the approximated likelihood (17), could be used to select K.

The mixture model

Let assume that each couple (X i , Z i ), i = 1, . . . , n, is an independent realization of the random vector (X, Z) where X has a density surrogate depending on its group belonging:

f (qg) X |Z g =1 (x; Σ g ) = qg j=1 f C j |Z g =1 (c j,g (x); λ j,g )
where q g is the number of principal components retained in the approximation (15) for the group g, c j,g (x) is the jth principal component score of X |Z g =1 for X = x, f C j |Z g =1 its probability density and Σ g the diagonal matrix of the principal components variances (λ 1,g , . . . , λ qg,g ). Conditionally on the group g, the probability density f C j |Z g =1 of the jth principal component of X is assumed to be the univariate Gaussian density with zero mean (the principal component are centred) and variance λ j,g . This assumption is in particular satisfied when X |Z g =1 is a Gaussian process.

Remark (Data generation). For a given cluster g, 1 ≤ g ≤ K, provided that the q g eigenfunctions f j and eigenvalues λ j of the covariance operator of X |Z g =1 are known, then, generating an approximation X (qg)

|Z g =1 of X |Z g =1 reduces to generating a real random variables C j according to centred Gaussian distributions with variance λ j (1 ≤ j ≤ q g ). Of course, that does not generate the true X |Z g =1 . However, the main structure of clusters is assumed to be characterized by this type of approximations.

The vector Z = (Z 1 , . . . , Z K ) is assumed to have one-order multinomial distribution M 1 (π 1 , . . . , π K ), with π 1 , . . . , π K the mixing proportions ( K g=1 π g = 1). Under this model we can deduce a surrogate for the unconditional density of X:

f (q) X (x; θ) = K g=1 π g qg j=1 f C j |Z g =1 (c j,g (x); λ j,g ) (16)
where θ = {(π g , λ 1,g , . . . , λ qg,g ) 1≤g≤K } and q = (q 1 , . . . , q K ). From this density surrogate, we deduce the pseudo likelihood:

l (q) (θ; X) = n i=1 K g=1 π g qg j=1 1 2πλ j,g exp - 1 2 C 2 i,j,g (X i ) λ j,g (17)
where C i,j,g (X i ) is the jth principal score of the curve X i for the group g and X = (X 1 , . . . , X n ).

Remark (Identifiability). When the approximation orders are different, identifiability of the mixture model ( 16) is straightforward. When all the approximation orders are equal, the identifiability of model ( 16) deduces directly from the identifiability of mixture of multivariate Gaussian [START_REF] Titterington | Statistical analysis of finite mixture distributions[END_REF].

Parameter estimation

In the unsupervised context the estimation of the mixture model parameters is not so straightforward as in the supervised context since the groups indicators Z i are unknown. On one hand, we need to use an iterative algorithm which alternates the estimation of the group indicators, the estimation of the PCA scores for each group and then the estimation of the mixture model parameters. On the other hand, the parameter q = (q 1 , . . . , q g ) will be estimated by an empirical method, similar to those used to select the number of components in usual PCA.

A classical way to maximise a mixture model likelihood when data are missing (here the clusters indicators Z i ) is to use the iterative EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | Finite Mixture Models[END_REF]. In this work we use an EM-like algorithm including, between the standard E and M steps, a first step in which the principal components scores of each group are updated and a second one in which the approximation order q are selected. Our EM-like algorithm consists in maximizing the pseudo completed log-likelihood

L (q) c (θ; X, Z) = n i=1 K g=1 Z g i log π g + qg j=1 log f C j |Z g i =1 (C i,j,g (X i ); λ j,g ) ,
where Z = (Z 1 , . . . , Z n ). Let θ (h) contains the current values of the estimated parameters at step h, h ≥ 1.

E step.

As the group indicators Z g i 's are unknown, the E step consists in computing the conditional expectation of the pseudo completed log-likelihood:

Q(θ; θ (h) ) = E θ (h) [L (q) c (θ; X, Z)|X = x] = n i=1 K g=1 t i,g   log π g + qg j=1 log f C j |Z g i =1 (c i,j,g (x i ); λ j,g )  
where c i,j,g (x i ) is the value of the random variable C i,j,g (X i ) for X i = x i , t i,g is the probability for the multidimensional curve X i to belong to the group g conditionally to C i,j,g (X i ) = c i,j,g (x i ):

t i,g = E θ (h) [Z g i |X = x] ≃ π g qg j=1 f C j |Z g i =1 (c i,j,g (x i ); λ j,g ) K l=1 π l q l j=1 f C j |Z l i =1 (c i,j,l (x i ); λ j,l ) . ( 18
)
The approximation in ( 18) is due to the use of the surrogate for the density of X given by ( 15).

Principal score updating step. The computation of the principal component scores has been described in Section 2.2. Here, the principal component scores C i,j,g of the multidimensional curve X i in the group g is updated according to the current conditional probability t i,g computed in the previous E step. This computation is carried out by weighting the importance of each curve in the construction of the principal components with the t i,g 's. Consequently, the first step consists in centring the curve X i within the group g by subtracting the mean curve computed using the t i,g 's: the basis expansion coefficients matrix A becomes A g = (I n -1I n (t 1,g , . . . , t n,g )) Ã. The jth principal component scores C j,g is then the jth eigenvector of the matrix A g WA ′ g T g associated with the jth eigenvalue λ j,g , with T g = diag(t 1,g , . . . , t n,g ).

Group specific dimension q g estimation step. The estimation of the group specific dimension q g is an open problem. In this work we propose to use, once the group specific FPCA have been computed at the previous step, the scree-test of [START_REF] Cattell | The scree test for the number of factors[END_REF] in order to select each group specific dimension q g . The advantage of using this test is that one hyperparameter (the threshold of the Cattell scree-test) allows to estimate K approximation orders.

M step. The M step consists of computing the mixture model parameters θ (h+1) which maximizes Q(θ; θ (h) ) according to θ. It leads simply to the following estimators

π (h+1) g = 1 n n i=1
t i,g , and λ j,g (h+1) = λ j,g , 1 ≤ j ≤ q g , where λ j,g is the variance of the jth principal component of the cluster g computed in the principal score updating step.

Convergence and numerical considerations. The proposed estimation algorithm is not a proper EM algorithm, since the growth of the pseudo likelihood is not guaranteed between two steps. The main reason is the use of a pseudo likelihood built on an approximation of the notion of density. Indeed, since only a finite number of principal components are used, and since these principal components are computed with different weights at each iteration of the algorithm (the t i,g 's), the 'data' on which the pseudo likelihood is computed, i.e. the principal components scores, are not the same at each step. To avoid this phenomenon, all principal components should be used, which is generally not applicable for functional data since they form an infinite set. Thus, the values of the likelihood can not be directly compared between two iterations. The properties of the EM algorithm, which guarantee the convergence to a local maxima of the likelihood in the classical framework does not work any more. Moreover, the approximation orders q g , 1 ≤ g ≤ K, are updated at each iteration, and this can also induce an artificial increase or decrease of the pseudo likelihood.

In practice, we adopt the following empirical strategy which allows us to perform numerical applications:

the algorithm is launched several times with random initializations for a small number of iterations, and the best reached solution is used as the initialization point for a longer algorithm [START_REF] Biernacki | Initializing EM using the properties of its trajectories in Gaussian mixtures[END_REF]. Typically, 10 small runs with 10 iterations are used in the following experiments.

the values of q g are initialized to 1, and they are only allowed to increase.

the number S of iterations is set as follows: for a given S (200 for instance), the algorithm is executed 20 times, and S is considered to be large enough if the maximum of the pseudo likelihood has been achieved before 3S/4 iterations for the 20 executions.

Anticipating the application results, Figure 4 and Figure 6 illustrate the trajectories of the pseudo likelihood and the approximation orders on simulated and real datasets.

Numerical experiments

This section is devoted to compare our approach, which we will call Funclustas in univariate case [START_REF] Jacques | Model-based clustering of functional data[END_REF]-with other existing methods. The evaluation of a clustering algorithm is always a difficult and subjective task. Following [START_REF] Guyon | Clustering: Science or art[END_REF], three evaluations strategies are considered in this paper. First, Funclust is compared to other clustering methods for univariate functional data using three classification benchmark datasets. Second, a simulation study allows to compare Funclust with another clustering method for multivariate functional data based on kmeans. Third, a real clustering application on a climatology dataset is carried out. The clusters obtained by Funclust and the k-means based method are then compared from the interpretation point of view.

Remark (Data registration). When working with functional data, a curve registration step is often needed to remove the amplitude and phase variation of curves (Ramsay and Silverman, 2005, Chap. 7). In our opinion, in the clustering context, the amplitude and phase variability of curves are interesting elements to define clusters. For instance, in the Canadian weather example which will be analysed in the sequel, the geographical interpretation of the clusters is mainly due to amplitude variability. Similarly, for the Growth dataset, it is shown in [START_REF] Liu | Simultaneous curve registration and clustering for functional data[END_REF] that performing registration before or simultaneously with clustering failed in retrieving the gender of subjects, probably because the main gender difference is due to a time wrapping effect. For this reason, we do not perform data registration in this work before our clustering study.

The R code for Funclust is available on request from the authors.

Benchmark study in the case of univariate functional data

The data. Three real datasets are considered: the Kneading, Growth, and ECG datasets. These three datasets, already studied in [START_REF] Jacques | Model-based clustering of functional data[END_REF], are plotted on Figure 1. The first dataset (Kneading) comes from Danone Vitapole Paris Research Center and concerns the quality of cookies and the relationship with the flour kneading process. The kneading dataset is described in detail in [START_REF] Leveder | Discrimination de courbes de pétrissage[END_REF]. There are 115 different flours for which the dough resistance is measured during the kneading process for 480 seconds. One obtains 115 kneading curves observed at 241 equispaced instants of time in the interval [0, 480]. The 115 flours produce cookies of different quality: 50 of them have produced cookies of good quality, 25 produced medium quality and 40 low quality. This data has been already studied in a supervised classification context [START_REF] Leveder | Discrimination de courbes de pétrissage[END_REF][START_REF] Preda | PLS classification of functional data[END_REF]. This data is known to be hard to discriminate, even for supervised classifiers, partly because of the medium class. The second dataset (Growth) comes from the Berkeley growth study [START_REF] Tuddenham | Physical growth of california boys and girls from birth to eighteen years[END_REF] and is available in the fda package of the software R. In this dataset, the heights of 54 girls and 39 boys were measured at 31 stages, from 1 to 18 years. The goal is to cluster the growth curves and to determine whether the resulting clusters reflect gender differences. The last dataset (ECG) is taken from the UCR Time Series Classification and Clustering website1 . This dataset consists of 200 curves from 2 groups sampled at 96 time instants, and has already been studied in [START_REF] Olszewski | Generalized feature extraction for structural pattern recognition in time-series data[END_REF].

Experimental set-up. For each dataset, the labels indicating the group membership of observations are available. These labels have been provided by human experts (Kneading and ECG datasets) or simply by the nature of the individuals (gender for the Growth dataset). In order to compare the ability of the clustering methods to retrieve the class labels, we choose to use the correct classification rate (CCR) which measures the adequacy of the resulting clusters with the known partition. This measure varies between 0 and 1, and larger the CCR, better the correspondence between the clustering and the known partition. In order to deal with the labelling problem, all the possible permutations are tested to label the K groups, and the best CCR is retained.

In this benchmark study, Funclust is compared to three challengers dedicated to the clustering of univariate functional data: FunHDDC [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] and k-centres (kCFC, [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF]) which are the closest methods and fclust [START_REF] James | Clustering for sparsely sampled functional data[END_REF] which is known to be a good challenger. Note that since no code is available for kCFC (to the best of our knowledge), only the comparison on the Growth dataset is possible, thanks to the classification results presented in [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF].

The number of iterations and the initialization are set as explained in Section 4.2. The threshold of the Cattell scree test allowing to select the approximation order q k is fixed to 0.05. With this experimental set-up, Funclust estimation is obtained in about 30 seconds for each dataset, on a laptop (2.80GHz CPU) and with a code in R software.

Results. The estimated approximation orders q g for Funclust are the following: Kneading (q 1 = 2, q 2 = 1, q 3 = 3), Growth (q 1 = 2, q 2 = 3), ECG (q 1 = 9, q 2 = 4). These q g 's are relatively close (or slightly greater) to the numbers q * g 's of principal components retained by the Cattell scree test (with the same threshold of 0.05) when carrying out FPCA on the true classes: Kneading (q * 1 = 1, q * 2 = 1, q * 3 = 2), Growth (q * 1 = 1, q * 2 = 1), ECG (q * 1 = 4, q * 2 = 5). The correct classification rates (CCR) according to the known partitions are given in Table 1. Funclust performs better than the other methods on two datasets among three (Kneading and ECG). On the last dataset, the results are relatively poor (69.89% whereas some method are about 97%), but they can be greatly increased (95.70%) if the dimensions q g are fixed to 2. This dataset illustrates that the choice of the approximation order is a very important question, and that further works have to be carried out in this direction.

Simulation study in the case of multivariate functional data

The data. In this simulation study, the number of clusters is assumed to be known: K=2. A sample of n = 50 curves are simulated according to the following model inspired by [START_REF] Ferraty | Curves discrimination: a nonparametric approach[END_REF] and [START_REF] Preda | Regression models for functional data by reproducing kernel hilbert spaces methods[END_REF]: for t ∈ [1, 21],

Class 1 :

X 1 (t) = -5 + t/2 + U 2 h 3 (t) + U 3 h 2 (t) + √ 0.1ǫ(t), X 2 (t) = -5 + t/2 + U 1 h 1 (t) + U 2 h 2 (t) + U 3 h 3 (t) + √ 0.5ǫ(t), Class 2 : X 1 (t) = U 3 h 2 (t) + √ 10ǫ(t), X 2 (t) = U 1 h 1 (t) + U 3 h 3 (t) + √ 0.5ǫ(t),
where U 1 ∼ N (0.5, 1/12), U 2 ∼ N (0, 1/12) and U 3 ∼ N (0, 2/3) are independent Gaussian variables and ǫ(t) is a white noise, independent of U i 's and of unit variance. The function h 1 , h 2 and h 3 (plotted on Figure 2) are defined, for t ∈ [1, 21], by h 1 (t) = (6 -|t -11|) + where (•) + mean the positive part, h 2 (t) = (6 -|t -7|) + and h 3 (t) = (6 -|t -15|) + . The mixing proportions π i 's are chosen to be equal, and the curves are observed in 1001 equidistant points (t = 1, 1.02, . . . , 21). Figure 3 plots the simulated curves. Experimental set-up. Funclust is compared to the multivariate functional data clustering methods described in [START_REF] Ieva | ECG signal reconstruction, landmark registration and functional classification[END_REF], based on k-means method with the following distances:

d 1 (X, Y ) = 2 j=1 T 0 (X j (t) -Y j (t)) 2 dt and d 2 (X, Y ) = 2 j=1 T 0 (X ′ j (t) -Y ′ j (t)) 2 dt,
where X ′ j (t) is the derivative of X j (t). A discussion on these distances in the case of univariate functional data can be found in [START_REF] Ferraty | Curves discrimination: a nonparametric approach[END_REF]. Since no public implementation is available for this method, we built our own implementation in the R software. In particular, the distance d 2 was computed using the package fda. In addition to these methods, quoted with kmeans-d 1 and kmeans-d 2 in the sequel, we consider usual k-means applied on the values of the functions at each observation points t = 1, 1.02, . . . , 21 (quoted as kmeans-points) and on the coefficients resulting from a linear spline smoothing with 30 equidistant knots (kmeans-spline). Linear spline smoothing has also been used by Funclust, with initialization and iterations number fixed following Section 4.2, and with a Cattell scree test threshold fixed to 0.05. Since both components X 1 and X 2 have similar covariance structure, the curves have not been normalized. Results. The convergence of Funclust is illustrated by Figure 4.

Canadian temperature and precipitation

In this last application, the Canadian temperature and precipitation data (available in the R package fda and presented in detail in [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF]) are used to compared Funclust with [START_REF] Ieva | ECG signal reconstruction, landmark registration and functional classification[END_REF]'s method (kmeans-d 1 and kmeans-d 2 ). The dataset consists in the daily temperature and precipitation at 35 different locations in Canada averaged over 1960 to 1994. The goal is to provide a clustering into 4 groups, and to give a geographical interpretation of the resulting clusters. Since the units of both curves are different (Celsius degrees and millimetres), the data are normalized following methodology presented at the end of Section 2.3. Figure 5 plots original and reduced curves. The curve has been smoothed using Fourier basis with 65 knots, as in [START_REF] Ramsay | Functional data analysis, 2nd Edition[END_REF]. Funclust, kmeans-d 1 and kmeans-d 2 are applied on this dataset. For Funclust, the initialization and the iterations number has been chosen following Section 4.2 and the threshold of the Cattell scree test has been fixed to 0.2. The convergence of Funclust is illustrated by Figure 6. Figure 7 presents the clustering into 4 groups of the Canadian weather stations using Funclust. We can observe four distinct groups of stations. The green group is mostly made of northern continental stations, whereas Atlantic stations and southern continental stations are mostly gathered in the red group. The blue group mostly contains Pacific stations and the last group (black) contains only the northernmost station Resolute (N.W.T.). We recall that all these results have been obtained without using the geographical positions of the stations.

In comparison, Figure 8 shows the clustering with kmeans-d 1 and kmeans-d 2 methods, using the same normalized curves and the same basis approximation. This clustering seems less pertinent than Funclust clustering since the blue and red group contains both Atlantic and Pacific stations. Nevertheless, the black group mainly contains continental stations whereas in Funclust clustering, continental and Atlantic stations are gathered together. This last fact is probably due to the Resolute station which is so different from the others (the temperature and precipitation are the lowest) that Funclust clusters this station alone apart in a group.

Discussion

In this paper we propose a clustering procedure for multivariate functional data based on an approximation of the notion of density for multivariate random function. We

Figure 1 :

 1 Figure 1: Kneading, Growth and ECG datasets.

Figure 2 :Figure 3 :

 23 Figure2: Function h 1 (t) = (6 -|t-11|) + , h 2 (t) = (6 -|t-7|) + and h 3 (t) = (6 -|t-15|) + for t ∈ [1, 21].

Figure 4 :

 4 Figure4: Convergence of the pseudo EM algorithm (top: pseudo likelihood evolution, bottom: approximation orders evolutions). The red part in the pseudo likelihood stops when the maximum is achieved.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Temperature and precipitation curves for 35 locations in Canada, averaged over 1960 to 1994. The top figures are the original curves and the bottom figures are the reduced ones.

Table 1 :

 1 Correct classification rates (CCR) in percentage for Funclust, FunHDDC (best model according BIC), fclust and kCFC on the Kneading, Growth and ECG datasets.

	method	Kneading Growth ECG
	Funclust	66.96	69.89	84
	FunHDDC	62.61	96.77	75
	fclust	64	69.89	74.5
	kCFC	-	93.55	-

Table 2 :

 2 Table2presents the mean and standard deviation of the correct classification rates (CCR), for 100 simulations. The results confirm the good behaviour of Funclust that we have already noticed for univariate functional data. Mean and standard deviation (for 100 simulations) of correct classification rates (CCR) in percentage for Funclust, kmeans-d 1 , kmeans-d 2 and k-means applied on observations points and spline coefficients.

	method	Funclust kmeans-d 1 kmeans-d 2 kmeans-points kmeans-spline
	Mean CCR	86.80	86.32	85.76	80.60	86.14
	Std CCR	14.51	16.30	10.28	14.94	15.80

http://www.cs.ucr.edu/∼eamonn/time series data/

introduce the principal component analysis of multivariate functional data. Assuming normality for the principal components within clusters, we define a mixture model for multivariate functional data clustering. The obtained model is an extrapolation of the finite dimensional Gaussian mixture model to the infinite dimensional setting. With respect to other clustering techniques for multivariate functional data, our methodology has the advantage to take into account the dependency between the curves defining the multivariate data. An EM-like algorithm is proposed for parameter estimation. The results obtained on real and simulated data illustrate the efficiency of our methodology. Some questions still remain open and further research are to be undertaken to provide answers. First of all, as previously discussed, the selection of the approximation orders is a great challenge for which we actually use an empirical method. Moreover, since only an approximation of the likelihood is available, the question of the convergence of the estimation algorithm is currently without response. However, the heuristic strategy used in this paper provide interesting clusters.