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

Christophe Biernacki
email: christophe.biernacki@math.univ-lille1.fr

Model-based clustering for multivariate partial ranking data

Keywords: Multivariate ranking, partial ranking, mixture model, Insertion Sort Rank, SEM algorithm, Gibbs sampling

This paper proposes the first model-based clustering algorithm dedicated to multivariate partial ranking data. This is an extension of the Insertion Sorting Rank (isr) model for ranking data, which is a meaningful and effective model obtained by modelling the ranking generating process assumed to be a sorting algorithm. The heterogeneity of the rank population is modelled by a mixture of isr, whereas conditional independence assumption allows the extension to multivariate ranking. Maximum likelihood estimation is performed through a SEM-Gibbs algorithm, and partial rankings are considered as missing data, what allows to simulate them during the estimation process. After having validated the estimation algorithm on simulations, three real datasets are studied: the 1980 American Psychological Association (APA) presidential election votes, the results of French students to a general knowledge test and the votes of the European countries to the Eurovision song contest. For each application, the proposed model shows relevant adequacy and leads to significant interpretation. In particular, regional alliances between European countries are exhibited in the Eurovision contest, which are often suspected but never proved.

Classification automatique de données de rang multivariées incomplètes

Résumé : Nous proposons le premier modèle de classification automatique pour données de rang multivariées potentiellement incomplètes. Ce modèle est une extension du modèle Insertion Sorting Rank (isr) pour données de rang, qui est un modèle efficace et signifiant obtenu en modélisant le processus de génération des données. L'hétérogénéité des données est traitée à l'aide d'un modèle de mélange, tandis qu'une hypothèse classique d'indépendance conditionnelle permet de prendre en compte les rangs multivariés. L'estimation des paramètres du modèle est réalisée par maximum de vraisemblance à l'aide d'un algorithme SEM-Gibbs. Les données incomplètes sont considérées comme des données manquantes, ce qui permet de les simuler durant le processus d'estimation. Après avoir validé la stratégie d'estimation sur données simulées, trois jeux de données ont été étudiés : les votes lors de l'élection du président de l'American Psychological Association de 1980, les résultats d'étudiants français lors d'un test de culture générale, et les votes des pays lors du concours de l'Eurovision. Pour chaque application, le modèle proposé a montré une très bonne qualité d'ajustement et à conduit à des interprétations intéressantes. Notamment, pour le concours de l'Eurovision, nous avons mis à jour des alliances géographiques entre pays voisins, ce qui a souvent été suspecté pour ce concours mais jamais prouvé.

Mots-clés : Données de rang multivariées, rangs partiels, modèle de mélange, tri par insertion, algorithme SEM, échantillonneur de Gibbs

Introduction

Ranking data occur when a number of subjects are asked to rank a list of objects according to their personal order of preference. These data are of great interest in human activities involving preferences, attitudes or choices like Politics, Economics, Biology, Psychology, Marketing, etc. For instance, the voting system single transferable vote occurring in Ireland, Australia and New Zeeland, is based on preferential voting [START_REF] Gormley | A mixture of experts model for rank data with applications in election studies[END_REF]. In Economics, it is sometimes more relevant to study the ranking of different economic actors according to some economical indicators rather than the only value of these indicators, since rank analysis focuses on comparisons between actors [START_REF] Schwab | The global competitiveness report 2012 -2013[END_REF]. Around the mid twentieth century, numerous probabilistic models for rank data have been proposed, based on different assumptions about the origin of a rank datum. For a survey, refer for instance to [START_REF] Marden | Analyzing and modeling rank data[END_REF]. Thurstone considers that a rank datum is the result of a ranking of latent continuous variables associated with each object to rank [START_REF] Thurstone | A law of comparative judgment[END_REF]. Paired comparison models [START_REF] Kendall | On the method of paired comparisons[END_REF][START_REF] Mallows | Non-null ranking models[END_REF] liken a rank to the result of a paired comparison process. Modelling parsimoniously each paired comparison leads to the famous Mallows Φ model [START_REF] Mallows | Non-null ranking models[END_REF] and its generalization to distance-based models [START_REF] Fligner | Distance based ranking models[END_REF]. Multistage models [START_REF] Fligner | Multistage ranking models[END_REF][START_REF] Plackett | The analysis of permutations[END_REF] assume that a rank is the result of an iterative process consisting in choosing at each step the best object among the remaining ones. Among this last class of models, 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 probably the most studied. More recently, [START_REF] Biernacki | A generative model for rank data based on sorting algorithm[END_REF] propose the Insertion Sorting Rank (isr) model as an effective and meaningful alternative for modelling ranking data. The isr model is set up by modelling the ranking generating process, assumed to be a sorting algorithm in which a stochastic event has been introduced at each comparison between two objects.

All these models consider homogeneous, full and univariate ranking data, what limits their scope. Indeed, in a lot of applications, the study of rank data discloses heterogeneity, due for instance to different political meanings, different economical strategies, different human preferences, etc. Heterogeneous rankings have thus attracted a great interest in the last decade: [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF] consider mixture of distance-based models and apply it to the modelling of the American Psychological Association's (APA) 1980 presidential election dataset [START_REF] Diaconis | A generalization of spectral analysis with applications to ranked data[END_REF], while [START_REF] Busse | cluster analysis of heterogeneous rank data[END_REF] adapt these models for tied and partial rankings. Mixture of multistage models [START_REF] Benter | Computer-based horse race handicapping and wagering systems: A report[END_REF] and Plackett-Luce models have also been successfully applied to the clustering of Irish election data and college admission data by [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF][START_REF] Gormley | A latent space model for rank 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]. If mixture of multistage models leads to interesting adequacy power, mixture of distance-based models have more meaningful parameters (and in a lower number), and moreover are simple to implement [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF]. On the other side, multivariate ranking data have been rarely studied, despite a strong interest in satisfaction surveys or polls. [START_REF] Bockenholt | Multivariate thurstonian models[END_REF] extends Thurstonian model to the multivariate case, but this extension suffers from numerical integration complexity.

Partial ranking is probably most frequent than full ranking: refer for instance to the 2002 General Election for the Irish House of Parliament dataset, studied in [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF], in which 96% of the electors did not rank all the 14 candidates, or the APA's 1980 presidential election which contains more than 60% of partial ranking. [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF]'s mixture model is extended to partial ranking by assuming a distribution on the missing entries according to a maximum entropy approach [START_REF] Busse | cluster analysis of heterogeneous rank data[END_REF]. [START_REF] Lebanon | Non-parametric modeling of partially ranked data[END_REF] propose a non-parametric estimator based on kernel smoothing for the estimation of the distribution of partial ranking data, and a visualisation technique based on multi-dimensional scaling in [START_REF] Kidwell | Visualizing incomplete and partially ranked data[END_REF].

Our contribution consists in defining a clustering algorithm for multivariate partial ranking data based on an extension of the isr model [START_REF] Biernacki | A generative model for rank data based on sorting algorithm[END_REF], initially devoted to univariate full ranking. For this, a mixture model will be considered, with a conditional independence assumption on the multivariate ranking components. The missing entries in the partial ranking will be considered as missing data and inferred in the estimation procedure. Thus, the proposed algorithm will be able to cluster ranking data sets with full and/or partial ranking, univariate or multivariate. To the best of our knowledge, this is the only clustering algorithm for ranking data with a so wide application scope.

The paper is organised as follows: Section 2 briefly reviews the isr model and extends this model for heterogeneous multivariate partial ranking data. Maximum likelihood estimation is considered in Section 3 by the mean of a SEM-Gibbs algorithm. Section 4 illustrates the relevance of the mixture of multivariate isr through simulation study and three real applications, and finally Section 5 concludes the paper.

2 The isr model for heterogeneous multivariate partial ranks

The univariate isr model

Rank data arise when judges or subjects are asked to rank several objects O 1 , . . . , O m according to a given order of preference. The resulting ranking can be designed by its ordering representation x = (x 1 , . . . , x m ) ∈ P m which signifies that Object O x h is the hth (h = 1, . . . , m), where P m is the set of the permutations of the first m integers. Based on the assumption that a rank datum is the result of a sorting algorithm based on paired comparisons, and that the judge who ranks the objects uses the insertion sort because of its optimality properties, [START_REF] Biernacki | A generative model for rank data based on sorting algorithm[END_REF] state the following isr model (see Appendix A.1 for details):

p(x; µ, π) = 1 m! y∈Pm p(x|y; µ, π) = 1 m! y∈Pm π G(x,y,µ) (1 -π) A(x,y)-G(x,y,µ) , (2.1)
where µ ∈ P m is the modal ranking, also named the reference or central ranking, π ∈ [ 1 2 , 1] is the probability of good paired comparison (in the sort algorithm) according to µ, and the sum over y ∈ P m corresponds to all the possible initial presentation orders of the objects to rank (with identical prior probabilities equal to 1/m!). The term G(x, y, µ) is equal to the number of good paired comparisons during the sorting process leading to return x when the presentation order is y, whereas A(x, y) corresponds to the total number of paired comparisons (good or wrong). Their precise definitions are given in Appendix A.2 and proofs can be found in [START_REF] Biernacki | A generative model for rank data based on sorting algorithm[END_REF]. This model has several interesting properties: the distribution is uniform when π = 1 2 ; µ is the mode of the distribution whereas its opposite ranking μ (μ = µ • ē with ē = (m, . . . , 1) the permutation of total inversion) is the rank of smallest probability; the mode of the distribution is uniformly more pronounced when π grows; identifiability of the parameters occurs once π > 1 2 .

Mixture of multivariate isr

Let now redefine x = (x 1 , . . . , x p ) as a multivariate rank, in which x j = (x j1 , . . . , x jmj ), for 1 ≤ j ≤ p, is a rank of m j objects. The population of multivariate ranks is assumed to be composed of K groups in proportions p k (p k ∈ [0, 1] and K k=1 p k = 1). Given a group k, the components x 1 , . . . , x p of the multivariate rank datum x are assumed to be sampled from independent isr distributions with reference ranking µ 1 k , . . . , µ p k ∈ P mj and good paired comparison probability π 1 k , . . . , π p k ∈ [ 1 2 , 1]. This conditional independence assumption, classical for categorical data and called latent class model [START_REF] Everitt | An introduction to latent variable models[END_REF][START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF], is considered since ranking can be viewed as specific categorical data and also because it is straightforward to implement. In addition, we will see in the study of the Eurovision contest (Section 4.4) that the model seems relatively robust to this assumption.

Inria

The unconditional probability of a rank x is then

p(x; θ) = K k=1 p k p j=1 1 m j ! y∈Pm j p(x j |y; µ j k , π j k ), (2.2)
where θ = (π j k , µ j k , p k ) k=1,...,K ,j=1,...,p and p(x j |y; µ j k , π j k ) is defined in (2.1).

Partial ranking

Each component x j of x can be full or partial. Let xj be the rank x j in which the unobserved positions are replaced by 0, and xj be x j with 0 in the place of the observed positions. With these notations we have x j = xj + xj . Let Ǐj ⊂ {1, . . . , m j } be the set of indices corresponding to observed ranking positions in x j , and Îj ⊂ {1, . . . , m j } the set of unobserved positions indices ( Ǐj ∪ Îj = {1, . . . , m j }). Let also define x = (x 1 , . . . , xp ) and x = (x 1 , . . . , xp ). In order to illustrate these notations, we consider the following example in dimension p = 1 with m 1 = 5: let assume that the judge does not rank objects in third and forth positions and returns x = (2, 5, 0, 0, 3). The objects O 1 and O 4 have then not been ranked. Possible x can be either (0, 0, 1, 4, 0) or (0, 0, 4, 1, 0), and then x can be either (2, 5, 1, 4, 3) or (2, 5, 4, 1, 3). Frequently, the objects in the top positions will be ranked and the missing ones will be at the end of the ranking, but our model does not impose such situation and is able to work with partial ranking whatever are the positions of the missing data.

3 Estimation

Likelihood expression

Let introduce the latent variable z which records the group membership of the observations. The latent variable z = (z 1 , . . . , z K ) is defined such that z k = 1 if the observation belongs to group k and z k = 0 otherwise. Let x = {x 1 , . . . , x n } be a sample of n multivariate rankings and z = {z 1 , . . . , z n } the corresponding latent variables. Let Ǐj i and Îj i be respectively the sets of indices of observed and unobserved positions in the jth component x j i of the ith observation x i . Similarly, let xj i and xj i correspond to the previous notations for the jth component of the ith observation, xi = {x 1 i , . . . , xp i } and xi = {x 1 i , . . . , xp i }. Let also define x = {x i ; i = 1, . . . , n} and x = {x i ; i = 1, . . . , n}. Let y i = (y 1 i , . . . , y p i ) ∈ P m1 × . . . × P mp be the presentation orders of the objects for the ith observation and y = {y 1 , . . . , y n }. Assuming that triplets (x i , y i , z i ) arise independently (i = 1, . . . , n), the observed-data loglikelihood of model (2.2) is:

l(θ; x) = n i=1 ln   K k=1 p k p j=1 1 m j ! y∈Pm j x∈X j i p(x|y; µ j k , π j k )   ,
where X j i = {x ∈ P mj : x h = xjh i , ∀h ∈ Ǐj i } is the set of all the rankings compatible with the observed part xj i of x j i . The maximization of this likelihood is not straightforward since several missing data occurs: the group membership z i , the presentation order y i and the unobserved ranking position xi . In such a situation, a convenient way to maximize the likelihood is to consider an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], which relies on the following completed-data log-likelihood:

l c (θ; x, y, z) = n i=1 K k=1 z k i p j=1 ln p k m j ! p(x j i |y j i ; µ j k , π j k ) .
Since this completed log-likelihood is not linear for all three type of missing data, the E step of the EM algorithm is intractable. In this work a SEM-Gibbs approach is proposed to overcome this difficulty.

SEM-Gibbs algorithm

The fundamental idea of this algorithm is to reduce the computational complexity that is present in both E and M steps of EM by removing all explicit and extensive use of the conditional expectations of any product of missing data. 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 , z i and xi at a so-called stochastic step (S step) from the conditional probabilities computed at the E step. Then these latent variables are directly used in the M step. However, the advantage with SEM-Gibbs algorithm relies on the fact that the latent variables are generated without calculating conditional probabilities at the E step, thanks to a Gibbs algorithm. The proposed SEM-Gibbs algorithm proceeds in the following two steps (SE-Gibbs step and M step), after starting from initial values x{0} = (x

j{0} i ) 1≤i≤n,1≤j≤p , y {0} = (y j{0} i ) 1≤i≤n,1≤j≤p , z {0} = (z k{0} i
) 1≤i≤n,1≤k≤K and θ {0} , and for Q SEM iterations.

The SE-Gibbs step It consists of three sub-steps that are now described at the qth iteration:

Generate

y j{q} i |{z {q-1} i , x j{q-1} i , (µ k , π k ) {q-1} }, where x j{q-1} i = {x jh{q-1} i ; 1 ≤ h ≤ m j } with x jh{q-1} i = xjh i if h ∈ Ǐj i and x jh{q-1} i = xjh{q-1} i otherwise.
For this, we consider a Gibbs sampler generating a chain y

p y jh i , y jh+1 i |y j1→h-1{q,r} i , y jh+2→mj {q,r-1} i , x j{q-1} i ; (µ k , π k ) {q-1}
for h ∈ {1, . . . , m j -2}, and (y

jmj -1{q,r} i , y jmj {q,r} i ) according to p y jmj -1 i , y jmj i |y j1→mj -2{q,r} i , x j{q-1} i ; (µ k , π k ) {q-1} ,
where k is such that z k{q-1} i = 1. The notation y ja→b i means y ja i , . . . , y jb i .

Generate

z {q} i |{y {q} i , x {q-1} i ; θ {q-1} } ∼ M(t 1{q} i , . . . , t K{q} i
) with, for k = 1, . . . , K:

t k{q} i ∝ p {q-1} k p j=1 p(x j{q-1} i |y j{q} i ; µ j{q-1} k , π j{q-1} k ) where x {q-1} i = {x j{q-1} i ; 1 ≤ j ≤
j Îj i (h){q,r} i , •) is generated according to p " xj Îj i (h) i , xj Îj i (h+1) i |x j Îj i (1→h-1){q,r} i , xj Îj i (h+2→ mj i ){q,r-1} i , xj i , y j{q} i ; (µ k , π k ) {q-1} « for h ∈ {1, . . . , mj i -2} and (x j Îj i ( mj i -1){q,r} i , xj Îj i ( mj i ){q,r} i ) according to p " xj Îj i ( mj i -1) i , xj Îj i ( mj i ) i |x j Îj i (1→ mj i -2){q,r} i , xj i , y j{q} i ; (µ k , π k ) {q-1} « .
where k is such that

z k{q} i = 1. The last value xj{q,Rj} i of the Gibbs chain is retained for xj{q+1} i .
The M step The M step consists in computing the parameter value θ {q} which maximizes the completed log-likelihood computed at the previous step:

θ {q} = argmax θ∈Θ l c (θ; {x, x{q} }, y {q} , z {q} ), where Θ = [ 1 2 , 1] p ⊗ p j=1 P mj × ∆ K with ∆ K = {p k : 1 ≤ k ≤ K; 0 ≤ p k ≤ 1, K k=1 p k = 1}
is the unit K-simplex, and where x{q} , y {q} and z {q} are simulated in the E step. The maximum for the mixing proportions are:

p {q} k = 1 n n i=1 z k{q} i
.

Given that the whole exploration of P m is intractable to estimate each modal rank µ j{q+1} k , a Gibbs sampling is used. The justification of the use of such an algorithm can be found in a Bayesian setting, in which the maximum a posteriori of the distribution p(µ j k |x, x{q} , y {q} , z {q} ; π j k ) is equivalent to the argument maximizing the completed log-likelihood when a uniform prior on µ j k is considered, since: p(µ|x, y, z; π, p) ∝ exp(l c (θ; x, y, z)

where θ = (µ, π, p) with µ = {µ j k } k=1,...,K,j=1,...,p , π = {π j k } k=1,...,K,j=1,...,p and p = {p k } k=1,...,K . The Gibbs sampling proceeds as follows. Starting from initial values π are estimated according to a Gibbs sampling alternating on R j iterations the two following steps:

(i) µ j{q,r+1} k |{x, x{q} , y {q} , z {q} ; π j{q,r} k
} is simulated following a similar scheme as for y

j{q} i : (µ jh{q,r+1} k , •) is simulated according to p µ jh k , µ jh+1 k |(µ j1→h-1{q,r+1} k , µ jh+2→mj
= P n i=1 z k i G(x j i ,x j{q} i ,y j{q} i ,µ j{q,r+1} k ) P n i=1 z k i A(x j i ,x j{q} i ,y j{q} i ) .
The retained value for π is the couple in the sequence of size R j maximizing the completed log-likelihood.

Choice of θ {q+1}

The SEM-Gibbs algorithm is stopped after a given number Q SEM of iterations. After removing a burn-in period of B SEM iterations, the estimation θ of θ is obtained as follows: for each distinct value of {µ j k ; 1

≤ j ≤ p, 1 ≤ k ≤ K} in the sequence µ j{q} k (B SEM ≤ q ≤ Q SEM ), take the mean πj k of the π j{q} k
and pk of the p {q} k , and retain the parameters (p k , µ j k , πj k ) 1≤j≤p,1≤k≤K leading to the highest log-likelihood. Since the log-likelihood is computationally intractable, it has to be approximated (Section 3.3). In such a situation, label switching can traditionally occur and could need to use some procedures as described in [START_REF] Celeux | Computational and inferential difficulties with mixture posterior distributions[END_REF] and [START_REF] Stephens | Dealing with label switching in mixture models[END_REF]. However, in our context, we can expect for no label switching when the clusters are well separated, while the case of non-well separated clusters will be avoided because not retained by the model selection criterion defined in Section 3.5.

Likelihood approximation

Since the observed-data log-likelihood is not tractable, it is approximated by:

l(θ; x) ≈ - n i=1 ln   1 Q l Q l q=B l 1 p(x i |(x i , y i , z i ) {q} ; θ)   (3.1)
where (x i , y i , z i ) {q} arise independently from p(x i , y i , z i |x i ; θ). The simulation of these triplets is carried out sequentially as in the SE-Gibbs step with Q L iterations.

Choice

of R j , Q SEM , B SEM , Q l , B l
The size R j of the Gibbs sampling, used in both SE-Gibbs and M steps, will be chosen greater than

mj (mj -1) 2
which is the maximum Kendall distance between two ranks of size m j , so that any rank of P mj can be reached with non-null probability from any arbitrary initialisation y j{q,0} i for instance. The number Q SEM of iterations of the SEM-Gibbs algorithm and the size B SEM of the burn-in period, as the sizes Q l and B l for the likelihood approximation, will be tuned empirically thanks to simulation study. We will see in Section 4.1 than these numbers have not to be very high to produce good estimations.

Model selection

In order to select the number K of components in the mixture, the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] is used:

BIC = -2l( θ; x) + (Kp + K -1) log(n),
where l( θ; x) is the maximum log-likelihood (in practice it will be approximated by (3.1)) and Kp + K -1 the number of continuous parameters (proportions p k and probabilities of good paired comparisons π j k ).

Inria

Estimation of the missing ranking positions

For each partial ranking, i.e. when xi (1 ≤ i ≤ n) is not full, the corresponding full ranking x i can be estimated using the Gibbs chain simulated in Section 3.3. The estimation of x i can be unconditional, using the mode of the empirical distribution generated by the chain (x {q} i ) q∈{BL,...,Q l } , or conditional to the cluster k when considering only the x{q} i 's of this chain such that z k{q} i = 1.

Numerical experiments

In this section, after having demonstrated the efficiency of the SEM-Gibbs estimation algorithm through a simulation study, the interest of the proposed mixture of multivariate isr is illustrated on three real datasets: the APA election dataset which contains full and partial rankings; the results of students to a general knowledge test, which consists of multivariate full rankings; the Eurovision dataset containing multivariate partial rankings.

In each application, the SEM-Gibbs algorithm is used with the following iteration numbers: R j = 10 (number of iterations of the Gibbs sampling, 1 ≤ j ≤ p), Q SEM = 100 (number of SEM iterations), B SEM = 10 (size of the burn-in period in the SEM-Gibbs algorithm) and Q L = 100 (number of iterations used for the likelihood approximation) and B l = 10 (size of the burnin period for the likelihood approximation). For each estimation, the SEM-Gibbs algorithm is launched 20 times, and the solution corresponding to the best approximated likelihood is retained.

Evaluation of the SEM-Gibbs algorithm on simulation

The goal of this experiment is to validate the proposed estimation algorithm by a simulation study, in particular in presence of partial rankings. For this, ranking data are simulated according to a mixture of two bivariate isr distributions (with equal proportions p 1 = p 2 = 0.5). The first component of the mixture is parametrized by µ 1 1 = (1, 2, 3, 4, 5), µ 2 1 = (3, 4, 1, 5, 2), π 1 1 = 0.8 and π 2 1 = 0.9, and the second by µ 1 2 = (5, 4, 3, 2, 1), µ 2 2 = (2, 5, 4, 1, 3), π 1 2 = 0.7 and π 2 2 = 0.95. Note that the components of the mixture are relatively well separated, and good estimation even for small sample size can be expected. In order to evaluate the proposed estimation strategy, two samples of size n = 200 and n = 4 000 are simulated and a mixture of two bivariate isr distributions is estimated with the SEM-Gibbs strategy described in Section 3. As to evaluate the estimation strategy in presence of partial ranking, some partial rankings are introduced as follows: for each dimension, in d% of the rankings two values are deleted (for instance a simulated rank (3, 4, 1, 5, 2) is replaced by (3, 0, 0, 5, 2)), in d 2 % of the rankings three values are deleted and in d 4 % of it four values are deleted. The efficiency of the estimation of the modal rankings µ j k is illustrated by computing the averaged, minimum and maximum Kendall distance between the true modal rank and the estimated one over 20 simulations. Similarly, the efficiency of the estimation of the probabilities π j k of good paired comparison is illustrated by the averaged, minimum and maximum absolute value of the difference between the true probability and its estimation. Results are given in Table 1, for four proportions of partial rankings (d = 0, 5, 10, 20).

The modal ranking estimation is perfect with both sample size until d = 5 (8.75% of missing values in the sampled ranks). For higher proportions of partial rankings, the estimations are still very satisfying. Similarly, estimation of the probabilities of good paired comparisons is very satisfying when d ≤ 5, and relatively correct for d ≥ 10. Comparing the two sample sizes, we note that the larger sample size leads to better estimations only in the case of full rankings. Indeed, for a given proportion of partial ranking, both sample sizes lead to similar results. This can be explained by the fact that the proportion of partial ranking has been fixed and then growing the sample size leads to grows equivalently the number of incomplete rankings. We conclude from these experiments that the SEM-Gibbs algorithm is an efficient algorithm to estimate a mixture of multivariate isr, even in presence of a moderate proportion of partial rankings, and we can use it in the following real applications.

The APA election

This dataset is famous in the ranking data literature [START_REF] Diaconis | A generalization of spectral analysis with applications to ranked data[END_REF][START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF][START_REF] Busse | cluster analysis of heterogeneous rank data[END_REF]. It consists of votes for the 1980 American Psychological Association presidential election. Five candidates were present, and the votes consist of ranking this five candidates in order of preference. Among these candidates, candidates noted 1 and 2 are research psychologists, 4 and 5 are clinical psychologists whereas candidate 3 is a community psychologist. A total of 15 449 votes were cast in the election among which 5 738 ranked all five candidates. In this application, we firstly compare the mixture of Mallows models using both Cayley and Kendall distances, already successfully applied on this dataset in [START_REF] Murphy | Mixtures of distance-based models for ranking data[END_REF], with our mixture of isr on the 5 738 full rankings. Mixture of isr is estimated by the SEM-Gibbs algorithm. In order to compare isr and Mallows mixtures, the BIC criterion is used: Figure 1 (left) shows its values for 1 to 10 groups. According to the BIC criterion, the mixture of isr distributions provides a better modelling of the APA dataset than the mixture of Mallows Φ model. Using now the BIC criterion to select the number of groups, we decide to select 4 groups, although BIC can be slightly lower for 10 groups but a plateau appears after 4 groups. In addition, it allows to simplify the interpretation of the parameter estimation (given in Table 2). Examining these estimated parameters, we remark that the two first groups, representing about 45% of the population of voters, ranks first the community psychologist, then the research

Inria

Full ranking

Full and partial ranking Value of the BIC criterion on the APA dataset. On the left, for 1 to 10 groups, with mixture of Mallows models (using Kendall and Cayley distance) and mixture of isr. On the right, for 1 to 5 groups, with mixture of isr taking into account partial ranking. . psychologists and finally the clinical psychologists. The difference between the first and second groups are the order of the two research psychologists and of the two clinical psychologists. The two others groups, representing about 55% of the voters, do opposite ranking, but with more dispersion since the probability of good paired comparison are slightly lower than for the two first groups. They rank first the clinical psychologists, then the research psychologists and finally the community psychologist. As previously, the difference between these two groups consists in the different rankings within each sub-category of psychologists. Our SEM-Gibbs algorithm allows moreover to take into account the partial rankings present in these election, which represent about 63% of the total votes. Using these partial rankings, the BIC criterion (Figure 1, right) leads to an hesitation between an homogeneous isr model, with parameters µ = (3, 1, 5, 4, 2) and π = 0.527, and a mixture with two components, in which a majority component (p 1 = 0.991) has the same parameters as the homogeneous isr, and a minority component (p 1 = 0.009) is parametrized by µ 2 = (5, 1, 3, 2, 4) and π 2 = 0.601.

We conclude this experiment with some words about the results of this election [START_REF] Diaconis | A generalization of spectral analysis with applications to ranked data[END_REF]. The elected candidate, using the system of proportional voting system, was candidate 1, which appears second in the modal ranking(s) of the isr model using partial and full rankings (with 1 or 2 groups). Examining the 4-components mixture estimated on only the full ranking, we see that the elected candidate has always been ranked in intermediary positions, never in the top or the last position of any cluster modal ranking. On the other side, the community candidate (3), which appears on the top position of two groups (and also when using partial rankings) has certainly not been elected since he was also ranked in the last position by a lot of people, as indicate the modal ranks of the third and fourth groups. Contrary to other voting system, in which we vote only for one candidate and then contribute only to its election, this election system based on proportional voting allows also to contribute to the non-election of other candidates.

The fact that the cluster structure is less evident when taking into account both partial and full ranking than only full ranking can be explained as follows: voters who have strong convictions on which candidate they want to elect and which one they want not, rank all the five candidates with their favourite candidate in the top position and the candidate they do not like in the last position. The resulting cluster structure is then very clear. On the other side, partial rankings are due to voters who rank only their favourite candidate(s) and then, the difference between voters and consequently the cluster structure appears to be less evident. Nevertheless, the elected candidate has always a stable position (second) when considering one or two clusters.

General knowledge test

In this application, the results of a general knowledge test given to two groups of students are studied. These students were in 2010 in third year (40 students) and fourth year (30 students) of the Statistics and Computer Science department of Polytech'Lille Engineering School (France). This test contains four questions about literature, sports, basic mathematics and cinema. Additionally to the students answers, their gender (girl or boy) and their year of study are registered. The four questions of the test were the following:

Literature. Rank these four French writers according to their date of birth:

O 1 = Hugo , O 2 = Molière, O 3 = Camus, O 4 = Rousseau. The correct answer is µ * = (2, 4, 1, 3).
Sport. Rank these 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).
Mathematics. Rank these four number in increasing order:

O 1 = π/3 , O 2 = ln 1, O 3 = e 2 , O 4 = (1 + √ 5)/2
. The correct answer is µ * = (2, 1, 4, 3). Figure 2 plots the test results for each question on polytopes. The size of the points is proportional to the number of observations. A mixture of multivariate isr distributions is estimated on these data, with a number of groups from 1 to 4. Using the BIC criterion (Figure 3), three groups are selected. The corresponding parameters are given by Table 3. Examining the estimated parameters leads to say that the Literature and Mathematics questions are easy for all students, since the modal rankings of the three groups correspond to the right RR n 8113

Literature Sport Mathematics Cinema k p k µ 1 k π 1 k µ 2 k π 2 k µ 3 k π 3 k µ 4 k π 4 k 1 0.4 (2,4,1,3) 0.839 (1,2,4,3) 1 (2,1,4,3
answers. The sport question is maybe more difficult since only the first group, representing 40% of the students, has a modal rank corresponding to the true ranking. The Cinema question is the more difficult since none group has the right answer as modal rank.

The first cluster, composed of the majority of students (40%), represents students knowing often the right answers (for 3 questions over 4) and relatively self-confident. The second cluster, of smaller size (27%), contains students with lower knowledge in Sport, but equally self-confident. Finally, the last cluster, of intermediary size (33%), corresponds to students rather low in almost all areas (except Mathematics) and with little confidence in them. We note here that the confidence is quite correlated with the knowledge of the true answer. Students are somewhere aware of their potential weakness.

Table 4 gives the confusion matrices between the estimated partition and the girl/boy and 3rd/4th year repartitions. This table exhibits that girls are more present in the third group than in the other groups. Since girls are generally known for being less interested in football than boys, this can explain why this third group have a modal ranking for the sport question so far from the true ranking.

The Eurovision Song Contest

The Eurovision Song Contest is an annual competition held among active member countries of the European Broadcasting Union. Each member country submits a song to be performed on live television and then casts votes for the other countries' songs to determine the most popular song in the competition. The vote consists in ranking ten preferred song in order of preference. We consider in these experiments the votes of the n = 34 countries who participate to the competitions from 2007 to 2012. During this six years, only 8 countries have participated to the six finals of the competition: 1: France, 2: Germany, 3: Greece, 4: Romania, 5: Russia, 6: Spain, 7: Ukraine and 8: United Kingdom. The studied dataset is then composed of multivariate ranking (p = 6 corresponding to the six contests between 2007 and 2012), each rank being of size m = 8 (only the votes for the 8 countries which participated to each final are considered) and all ranking being partial. The absence of full ranking signifies that none country participating to the votes has ranked all of the 8 considered countries in its 10 preferences. This dataset is challenging since the number of observations (n = 34) is small compared to the size of the rank (m = 8) and the presence of partial rankings (precisely, 57.7% of the rankings elements are missing).

A mixture of multivariate isr distributions is estimated on this dataset, with a number of groups from 1 to 6. The BIC criterion (Figure 4) leads to select 5 groups, whose parameters are given by Table 5. Interesting analysis of these data can be deduced from these parameters, especially using the modal rankings. Indeed, in 2007, all the groups seem to have globally voted in the same way. Indeed, among the 8 considered countries the best were Ukraine [START_REF] Celeux | Computational and inferential difficulties with mixture posterior distributions[END_REF], Russia [START_REF] Celeux | The SEM algorithm: A probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF] and Greece (3) who finished respectively second, third and seventh. These three countries appear in the first three positions of the modal rank of each group, except for the group 5, in which Romania if preferred to Greece. Another remark concerns the votes of the group 5. In each six years, Ukraine and Russia (7 and 5) have been well sorted by this group of countries, since they always appear in a good position in the modal rankings of this group. This suggests the existence of geographical voting alliances often suspected in this contest. In order to verify this assumption, we plot on Figure 5 the estimated countries clustering in 5 groups (group 1: red, group 2: blue, group 3: yellow, group 4: green and group 5: gray). This clustering has been obtained by simulating the latent variable z. For this, the SE-Gibbs step has been repeated 1 000 times, using the estimated parameters given in Table 5. This map confirms the existence of geographical alliances. Indeed, group 1 (red) contains essentially West European countries, group 2 (blue, which is the smallest group) contains some Northern countries, group 3 (yellow) contains Mediterranean countries, group 4 (green) is maybe more dispersed and finally group 5 (gray) contains essentially East European countries. Note that this last group is the one we previously detect to be used to rank Ukraine and Russia in the first positions. This clustering in 5 groups means that countries of a same cluster tend to have similar votes. The regional proximity of the country within each cluster confirms the assumption A last remark on this dataset concerns the independence assumption we made in order to work with multidimensional rankings. The geographical alliances we just exhibited tend to suggest that this assumption is not satisfied. Nevertheless, the interesting obtained results from an interpretation point of view allow to think that the proposed model is relatively robust to this assumption.

Conclusion

Mixture of isr models provides an efficient tool to model heterogeneous population of ranks. The estimation strategy, based on a SEM-Gibbs algorithm, is moreover an elegant way to take into account partial ranking data, which are very frequent in real applications. Using this model on the challenging Eurovision dataset has allowed to exhibit and confirm regional alliances between countries participating this contest, which has been often suspected and criticized. In [START_REF] Gatherer | Comparison of eurovision song contest simulation with actual results reveals shifting patterns of collusive voting alliances[END_REF], a study of the Eurovision contest votes from 1975 to 2000, based on Monte Carlo simulations, has also concluded to regional alliances. A more precise examination of the structure of the geographical blocs of alliances allowed them to exhibit that centrally placed countries within these blocs have a higher probability of being future winners.

In this work, the conditional independence assumption has been used to take into account multivariate ranking. This is a first approach that has the advantage of simplicity. The introduction of correlation between rankings is a great challenge which has to be raised in order to define even more relevant models for multivariate ranking.

Computer code

An R package for clustering of multivariate partial ranking data using our methodology is available on request from the authors and will be soon available on the CRAN website1 .

A Definitions and notations for isr A.1 Definitions and notations

The sorting algorithm modelled by isr operates as follows. First, the current object in the presentation order y is placed before (on the left of) 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 with probability π and this process is restarted with the next object in y. Otherwise, the pair order is reversed (with probability π) and a new pair comparison is performed with the next object after (if it exists). And so forth, until obtaining the final ranking x. Table 6 illustrates this algorithm on an example.

A.2 Notations

Let x -1 = (x -1 1 , . . . , x -1 m ) ∈ P m be the ranking representation of a rank data, which contains the ranks given to the objects and means that O i is in the x -1 i th position.

A(x, y) = m j=1 A - j (x, y) + A + j (x, y) is the total number of all paired comparisons for the whole sorting process, where -A - j (x, y), the number of all comparisons of the current object with the objects already ranked (according to x) on its left (if they exist), is the cardinal of A - j (x, y) = {i : x -1 yi < x -1 yj , 1 ≤ i < j}, the set of the indices of the presentation order y for which the already sorted objects O y1 , . . . , O yj-1 are ranked in x before the current object O yj , and consequently on its left.

-A + j (x, y), which indicates if the current object O yj is compared, at the j step of the sorting, with the object ranked in x just on its right, is the cardinal of A + j (x, y) = {i : i = arg min 1≤i ′ <j {i ′ : x -1 y i ′ > x -1 yj }} containing the index of the rank y designating the object sorted in x just after (so on the right of) O yj among the already sorted objects O y1 , . . . , O yj-1 , if it exists. Inria G(x, y, µ) = m j=1 i∈A - j (x,y) δ yiyj (µ) + i∈A + j (x,y) δ yjyi (µ) is the total number of good paired comparisons for the whole sorting process, with δ ii ′ (µ) = 1{µ -1 i < µ -1 i ′ } is equal to 1 if O i is correctly ranked before O i ′ (according to µ), 0 otherwise (i, i ′ = 1, . . . , m, i = i ′ ).

RR n 8113

  } i, in which the last value y j{q,Rj } i (for R j large enough, see Section 3.4) is retained for y sampler generates r ∈ {1, . . . , R j } sequences y j{q,r} i where (y jh{q,r} i , •) is generated according to the following distribution:

  Figure1: Value of the BIC criterion on the APA dataset. On the left, for 1 to 10 groups, with mixture of Mallows models (using Kendall and Cayley distance) and mixture of isr. On the right, for 1 to 5 groups, with mixture of isr taking into account partial ranking. .

Table 2 :

 2 Results of estimation of a mixture of isr with 4 groups on the APA full ranking dataset.

Cinema.

  Rank 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).

Table 3 :

 3 Parameters of a mixture of multivariate isr distribution with 3 groups on the general knowledge test data.

Figure 2 :

 2 Figure 2: Empirical distribution of the 70 student's answers to the general knowledge test.

Figure 3 :

 3 Figure 3: Value of the BIC criterion with mixture of isr for the general knowledge test.

Inria

Figure 4 :

 4 Figure 4: Value of the BIC criterion with mixture of isr for the Eurovision datatest.

Figure 5 :

 5 Figure 5: Classification of European countries according to their vote to the Eurovision competition.

Table 6 :

 6 An example to illustrate the stochastic insertion sort algorithm considered by the isr model, with µ = (1, 2, 3), y = (1, 3, 2) and x = (3, 1, 2). The notation x (j) means the ranking of the j first objects in y in the order imposed by xWe now define the notations used in Section 2.1.

2 2 . 1 2 . 2 2 . 3 3 5 3 . 1 3 . 3 3 . 4 3 . 5 3 . 6 9 4 5 16 A

 212223331333435369516 The isr model for heterogeneous multivariate partial ranks 4 The univariate isr model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Mixture of multivariate isr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Partial ranking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Estimation Likelihood expression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2 SEM-Gibbs algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Likelihood approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Choice of R j , Q SEM , B SEM , Q l , B l . . . . . . . . . . . . . . . . . . . . . . . . . 8 Model selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Estimation of the missing ranking positions . . . . . . . . . . . . . . . . . . . . . Numerical experiments 9 4.1 Evaluation of the SEM-Gibbs algorithm on simulation . . . . . . . . . . . . . . . 9 4.2 The APA election . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.3 General knowledge test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.4 The Eurovision Song Contest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Conclusion Definitions and notations for isr18 A.1 Definitions and notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 A.2 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Inria RESEARCH CENTRE LILLE -NORD EUROPE Parc scientifique de la Haute-Borne 40 avenue Halley -Bât A -Park Plaza 59650 Villeneuve d'Ascq Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

Table 1 :

 1 Results of simulation: distance between the true parameter and the estimated one for different proportion d of partial ranking.

	n→	200	K(µ, μSEM-gibbs ) 4 000 200 4 000 200 4 000	200	4 000	|π -πSEM-gibbs | 200 4 000	200	4 000
	d j k	mean	best		worst	mean	best	worst
	0 1 1	0	0	0	0	0	0	0.025 0.005 0.023 0.004 0.027 0.006
	0 1 2	0	0	0	0	0	0	0.013 0.006 0.011 0.005 0.017 0.007
	0 2 1	0	0	0	0	0	0	0.002 0.002 0.001 0.001 0.004 0.002
	0 2 2	0	0	0	0	0	0	0.022 0.002 0.021 00.02 0.023 0.003
	5 1 1	0	0	0	0	0	0	0.014 0.071 0.011 0.070 0.016 0.072
	5 1 2	0	0	0	0	0	0	0.044 0.070 0.038 0.070 0.048 0.071
	5 2 1	0	0	0	0	0	0	0.017 0.041 0.015 0.041 0.020 0.042
	5 2 2	0	0	0	0	0	0	0.030 0.022 0.028 0.021 0.032 0.023
	10 1 1	0	0.025	0	0	0	0.5	0.021 0.101 0.014 0.089 0.039 0.299
	10 1 2 0.025	0	0	0	0.1	0	0.088 0.084 0.083 0.057 0.109 0.087
	10 2 1	0	0	0	0	0	0	0.028 0.065 0.023 0.046 0.049 0.039
	10 2 2	0	0	0	0	0	0	0.051 0.030 0.041 0.013 0.066 0.031
	20 1 1	0	0	0	0	0	0	0.110 0.128 0.090 0.067 0.130 0.136
	20 1 2	0.41	0.05	0.1	0	0.6	0.5	0.164 0.140 0.144 0.133 0.181 0.197
	20 2 1	0	0	0	0	0	0	0.050 0.066 0.034 0.004 0.066 0.074
	20 2 2	0	0.045	0	0	0	0.5	0.067 0.099 0.049 0.061 0.085 0.431

Table 4 :

 4 Repartition (in %) of girls/boys and 3rd/4th years student in the three class estimated by the mixture of isr.

		sexe			year of study
	class girl boys		class 3rd	4th
	1	5.7	35.7 41.4	1	28.6	12.9	41.5
	2	8.6	18.6 27.2	2	17.1	10	27.1
	3	21.4	10	31.4	3	11.4	20	31.4
		35.7 64.3 100		57.1	42.9	100

Table 5 :

 5 Parameters of a mixture of multivariate isr distribution with 5 groups on the Eurovision dataset.

RR n 8113

Inria

http://cran.r-project.org/ RR n 8113

Acknowledgement

We thank Quentin Grimonprez for developing a C++ implementation of the SEM-Gibbs algorithm, and Pr. Brendan Murphy for providing the idea to work on the Eurovision contest.


Procópio Stein
email: procopiostein@gmail.com

Anne Spalanzani
email: anne.spalanzani@inria.fr

Vítor Santos

Christian Laugier
email: christian.laugier@inria.fr

Robot Navigation Taking Advantage of Moving Agents

A crucial requirement for service robots is to be able to move in dynamic environments shared with humans as well as interact with them. Navigation in such environments is a challenging task, as the environment is constantly changing, future states have to be predicted and planning and execution must be carried on-line.

However, even in very complex situations, humans can easily find a path that avoid both dynamic agents and static obstacles. This paper proposes a technique to take advantage of the human movement in such populated environments, selecting a leader to be followed in a probabilistic fashion, according to the robot's desired destination.

In this way, the robot can take advantage of the paths traveled by humans, effortlessly avoiding dynamic and static features as the human leader does, relieving the robot from the burden of having to generate its own path.

I. INTRODUCTION

With advances in mobile robotics and lowering costs of computers, it is becoming more and more common for us to find robots among groups of people. Service robots (home care, hospital, museum guides) are real example cases where robots have to be able to move and interact with humans in an ever changing environment. The success of interactions and human acceptance of service robots is directly related to the way they behave and approach others, as well as their capability to adapt to the environment.

However, navigation in dynamic environments is still an open and challenging issue for the robotic community. In such environments, sensor's measurements are prone to noise and the measurements have a short lifespan, being valid only for small time periods. Static features, that could guide a navigation algorithm, might not be detected. There are also limitations in the time spent by the navigation algorithm to provide solutions, and optimal approaches are unsuited for this task, as generated paths might be valid only for few time-steps.

Despite these difficulties, humans can easily navigate in dynamic environments. Moreover, their movement can provide indirect information of the environment, as they usually do not move at random. Instead of that, they move according to typical patterns, and their movements are related to features that they are interested in, such as doors, elevators, stairs or other people [START_REF] Mombaur | From human to humanoid locomotionan inverse optimal control approach[END_REF].

Taking advantage of that, several works propose to model those typical paths with probabilistic approaches as Gaussian Processes [START_REF] Tay | Modelling smooth paths using gaussian processes[END_REF] [3] [START_REF] Ellis | Modelling pedestrian trajectory patterns with gaussian processes[END_REF] or Hidden Markov Models [START_REF] Vasquez Govea | Growing hidden markov models: An incremental tool for learning and predicting human and vehicle motion[END_REF], to address the problem of navigating in populated environments [START_REF] Bennewitz | Learning motion patterns of people for compliant robot motion[END_REF], [START_REF] Fulgenzi | Probabilistic motion planning among moving obstacles following typical motion patterns[END_REF]. Correctly detecting and learning a likely path permit the robot to avoid trajectories that have a risk of future collision with a pedestrian, as well as avoiding social disturbances [START_REF] O'callaghan | Learning navigational maps by observing human motion patterns[END_REF], [START_REF] Sehestedt | Robot path planning in a social context[END_REF].

Those approaches, however, do not take into account changes people perform in their typical paths to avoid and adapt to other moving people [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF]. These conditions allied to excessive future uncertainty may lead to situations where every generated path lead to colisions or frozen situations, as shown by [START_REF] Trautman | Unfreezing the robot: Navigation in dense, interacting crowds[END_REF].

The key insight of this paper is to present a navigation technique for dynamic environments that takes advantage of the typical movement of humans. This approach relies on the fact that people try to guarantee their safety and the safety of other people in the same environment, avoiding obstacles and avoiding hitting other persons.

In other words, people walking in populated scenarios can provide rich indirect information about their surroundings, as they are constantly dealing with large amount of high-level information and reacting to it, while following a goal.

A robot moving to a certain destination could identify behaviors of humans and detect a leader, someone moving along a typical pattern that would pass close to the destination point. After identification, the robot could follow people along that path, as "moving with the flow", relying more and more on people in front of it, as leaders or as part of a swarm.

A similar approach has been developed by [START_REF] Müller | Socially inspired motion planning for mobile robots in populated environments[END_REF]. With the difference that the main goal of that work is to implement a human-like motion behavior, and the choice of a leader is deterministic, based on the motion direction of the subject regarding the trajectory planned by the robot. Although such technique can be used for the same objectives that are proposed in this work, in some environments the algorithm may fail to find a leader as the extrapolation of the initial movement of a candidate may not match his/her actual goal, as illustrated in Fig. 1.

In this paper, however, the choice of a leader is performed probabilistically, using the Growing Hidden Markov Models (GHMM) technique, an extension to the HMM capable of learning the models parameters and structure in an incremental fashion. This technique not only provides a prediction of the future states of a moving agent, but also its goal, which allows the robot to plan further ahead, with a probabilistic knowledge of the goals pursued by the moving persons.

The global path planning is implemented with the RiskRRT algorithm [START_REF] Fulgenzi | Probabilistic motion planning among moving obstacles following typical motion patterns[END_REF], which takes into account the risk of collision with dynamic obstacles while generating the tree. The use of this algorithm guarantees that the robot can find path solutions in dynamic environments when a leader is not found.

In section II the technique to choose a leader in dynamic environment is presented, followed by the explanation of a leader following technique, in section III. The experiments and their results are presented in section IV, and after that, the conclusions of this work are presented in section V.

II. CHOOSING THE LEADER

As the proposition of this work is to take advantage of the movements of a person, the method used to choose which person to follow plays a major role. Here, the choice of a leader is implemented based on the distance between the goal given to the robot and the predicted goal of the leader candidate.

The prediction of the motion and goal of a leader is not an easy task in dynamic environments. A simplistic approach may extrapolate the current orientation and speed of a moving person in an attempt to determine his/her likely goal. But due to environment or dynamic restrictions, the subjects can completely change their trajectory in the subsequent time steps, invalidating the extrapolation, as shown in Fig. 1. Here the dashed lines represent the goal and motion prediction using simple extrapolation or a Kalman Filter. Note however that due to the environment structure and interest points, the actual goal and paths, represented by solid lines, highly differ from the ones predicted using simplistic assumptions.

Knowing that the movement of humans is highly dependent on the environment structure and interest points, recent approaches take advantage of the typical paths in order to make predictions of motions and goals of humans. This approach can overcome the limitations posed by simple extrapolation techniques, as it allows to take into account the structure of the environment as well as the most common motion patterns.

In the current work, in order to predict the goal of moving agents, the Growing Hidden Markov Models (GHMM) algorithm is used [START_REF] Vasquez Govea | Growing hidden markov models: An incremental tool for learning and predicting human and vehicle motion[END_REF]. It implements an approach where the learning and prediction phases are on-line concurrent processes, resulting in a learn and predict paradigm. The structure of the GHMMs are the same as the regular HMMs, with the difference that as new observations sequences are incorporated into the model, the transition structure and the The GHMM algorithm consists of the use of the Growing Neural Gas (GNG) algorithm [START_REF] Fritzke | A growing neural gas network learns topologies[END_REF], used to estimate the model structure as well as the transition probabilities of a Hidden Markov Model (HMM). As the algorithm is adaptive, it is capable of creating or removing states to cope with new observations.

One very important aspect of the GHMM algorithm is that it is based on the hypothesis that moving agents always try to reach a goal in the environment. Therefore, each goal in the scenario has a different HMM associated to it. As the intention of the current work is to predict the most likely goal of a moving agent, the GHMM inherently provides a direct solution to that problem.

III. FOLLOWING THE LEADER

Once a moving agent has been detected as a leader, due to a similar goal with the robot's goal, there still remains the problem on how to follow the person. Simplistic approaches that try to use the current leader position as a subgoal may bring the robot to situations where undetected obstacles are present between the robot and the leader.

Following the leader path, generating subgoals along the tracked trajectory can overcome situations were obstacles that were not detected appear between the robot and the leader. However, if the robot is not close enough to the leader, it can lose track from the person or different agents that are not leaders may appear in the scene, blocking the tracked path.

Besides that, the motion algorithm has to be able to maintain a navigation solution even in situations where no leader is found, or after a leader is lost due to sensor occlusions or scene exit by the leader being followed. These constraints are addressed with the use of a variation of the RRT algorithm that efficiently explores the free space while at the same time takes into account the risk of collision with moving agents in a scenario.

A. Risk Rapid-exploring Random Tree

The Risk Rapid-exploring Random Tree (RiskRRT) is a variation of the classic RRT algorithm presented by [START_REF] Lavalle | Randomized kinodynamic planning[END_REF] developed for navigation in dynamic environments. It takes into account the risk of traveling along generated paths according to predicted objects' motion. It combines a part dedicated to perception (of static and moving obstacles) with another for planning trajectories. Navigation and planning are done in parallel.

The configuration-time space is searched randomly, and a tree T is grown from the initial configuration all over the configuration space. The algorithm chooses a point P in the configuration space and tries to extend the current search tree toward that point.

The points P are randomly sampled on the map, but at the beginning, and then once every 100 times, the goal itself is chosen; this bias, which has been empirically set, speeds up the exploration toward the goal. The node chosen for extension is the most promising node: all the nodes in T are weighted taking into account the risk of collision and the estimated length of the total path:

w(q N ) = L π (q N ) dist(q 0 , q N , P ) (1)
w(q N ) = w( q N ) q wq (2)
At numerator, the likelihood of π(q N ) is normalized with respect to the length of the path N ; at denominator, dist(.) is the sum between the length of the path from the root q 0 to the node q N (which is known) and an estimation of the length of the path to P .

The weights are normalized over the set of nodes in the tree [START_REF] Tay | Modelling smooth paths using gaussian processes[END_REF]. The node to grow next is then chosen taking the maximum over the weights or drawing a random node proportionally to the weight. The new node q + is obtained applying an admissible control from the chosen node q toward P . The weight of q + is computed. If w(q + ) ≥ w(q) the tree is grown again from q + toward P otherwise another point is sampled from the space.

The likelihood of each partial path can also be expressed as the multiplication of the independent probability of collision with static (P cs ) and dynamic components (P cd ):

L π (q N ) = N n=0 (1 -P cs (q N )) • N n=0 (1 -P cd (q N )) (3)
The prediction approach for forecasting the position of moving obstacles in the near future is done using the GHMM predictor. With the information of probable occupied positions in the future, the robot can anticipate the behavior of the agents. The selection of the best trajectories is done by taking into account the probability of collision for each path.

The probability of collision, or risk, can be seen in this case as a measure of the feasibility of a path, with the maximum accepted risk specified as a threshold. The RiskRRT algorithm also takes into account the interactions among Fig. 4. RiskRRT algorithm avoiding a social interaction zone humans so the robot can behave in a socially acceptable way. Therefore, the risk function must rely on safety but also in human friendly navigation.

The current version of the RiskRRT algorithm also takes into account the interactions among humans in order to behave in a socially acceptable way. The robot behavior should follow social conventions, respecting proximity constraints, avoiding people interacting or joining a group engaged in conversation without disturbing them. Therefore, the risk function must rely on safety but also in human friendly navigation.

After these extensions the "probability of success" calculated for every partial path is given by the probability of not encountering a collision along the path and not entering a social interaction zone. For more details about this method, refer to [START_REF] Rios-Martinez | Understanding human interaction for probabilistic autonomous navigation using Risk-RRT approach[END_REF]. Fig. 4 depicts a situation where the RiskRRT algorithm generates a path that avoids the social interaction between two persons.

B. Algorithm

The developed program to follow a leader is shown in III-B. The program starts after receiving a desired goal for the robot, which is used to initiate the RiskRRT algorithm. At the same time, the tracking program starts, in order to detect moving agents in the scenario. The detected agents are fed into the GHMM predictor, which outputs the prediction of a goal for each person.

After that, the Euclidean distance between the robot's current goal and each subject's predicted goal is computed. If this distance falls inside a threshold, the corresponding subject is chosen as a leader.

In the case a leader is found, the robot starts to track its path, but does not immediately starts to follow the leader. This only takes place once the elected leader is closer to the goal than the robot. The reason for this criteria is to avoid situations where the robot would have to move away from the goal in order to follow the path of a leader.

Once the criteria to start following a leader is satisfied, the robot calls an update routine in RiskRRT planning algorithm. This routine causes the RiskRRT algorithm to use the new subgoal instead of the original goal, in order to find a path.

As a result, the algorithm explores the open space and finds a path that poses the lesser risk to bring the robot to the chosen subgoal, which lies over the leader path. The getSubGoal routine decides when to pass a new subgoal of the path to the RiskRRT algorithm, based on the robot Leader Follower end while 26: end procedure distance to the previously passed subgoal. This sequence of steps makes the robot follow the tracked leader's path.

This sequence of steps makes the robot follow the tracked leader path. The use of the Risk-RRT algorithm to reach and follow a leader's path has two main advantages. Firstly, it provides a reliable method to navigate until the leader's path start, since in a dynamic environment the space between the robot and its first subgoal may be occupied by moving agents or static obstacles.

In second place, once the robot reaches and starts to follow the leader path, the algorithm is capable of reusing nodes of its exploration tree to efficiently generate new paths for each new subgoal received by the update routine. The reuse of previously generated nodes, reduce the computational load of the algorithm, while still taking into account the risk of navigation.

Finally, in the case that a leader is not found, or the current leader is lost, the update routine sends once again to the RiskRRT algorithm its final goal, as chosen at the beginning of the program.

IV. EXPERIMENTS

The experiments were performed using several independent modules using the Robot Operating System (ROS) [START_REF] Quigley | ROS: an open-source robot operating system[END_REF]. The modular architecture provided by ROS gives a series of advantages when implementing experiments, as modules can be added or removed, allowing the test of different techniques without the need of modification of the remaining modules. To conduct the experiments, the main hall of INRIA Rhône Alpes has been chosen (Fig. 5). It is an interesting choice as it has a large flow of people during different times of the day, entering and leaving the building during lunch hours and at the beginning and the end of a working day. These conditions allow an easy understanding of the typical paths present in the scenario.

A. Real Data Acquisition

The trajectories used in the experiments to test the robot capability to chose and follow a leader are real human trajectories, that were previously tracked and recorded. An overhanging camera with wide angle lens provides an overview of the test area. The implemented tracker is based on the work of [START_REF] Snape | Independent navigation of multiple mobile robots with hybrid reciprocal velocity obstacles[END_REF]. In the current work, fiducial markers were worn as hats by subjects in order to provide a robust and fast deployment tracker system, as shown in Fig. 6.

The GHMM is trained using a set of the real data acquired with the tracking system. Volunteers were asked to move naturally among interest points in the environment, as the entrance of the hall and the two doors. Fig. 7 shows a sample of the these trajectories.

B. Test Scenario

Two types of tests were conducted, one that evaluates the leader detection technique when several subjects move close to each other (Fig. 8), and another test that evaluates the advantage of the proposed technique to avoid moving agents (Fig. 9). The human perception and detection are performed by a second computer, using the techniques mentioned in the previous subsection.

The tests were conducted using a simulated robot, while the scenario agents represent real data recorded from the motion of humans. The scenario is simulated using The circles are the persons detected by the overhanging camera, and the triangles are their respective predicted goals. They have a letter associated to identify their colors (Red, Green and Blue). The robot goal is marked as an X, located at the lower left of the test area. Finally, the dots represent the RiskRRT tree nodes and the solid line is the best path found by the planning algorithm.

C. Leader Detection Test

In the first test, shown in Fig. 8, three humans start to move just in front of the robot, and pursue one different goal each. After some iterations, as the subjects start to move in the scenario, the prediction algorithm gives an estimation for two of them (red and green). Based on that estimations, the leader following algorithm makes the choice of following the red subject, as its predicted goal lies within a distance threshold from the robot's goal.

Once that the leader is closer to the goal than the robot, by an empiric factor, the planned trajectory is computed from the robot to the subgoal. This corresponds to the first position of the chosen leader's path, and the robot starts to move along the trajectory taken by the human, in the direction of its desired goal.

D. Dynamic Agents Avoidance Test

The objective of this test is to evaluate the benefits of following a leader in order to avoid other dynamic agents and is shown in Fig. 9. The way the robot selects and follow a leader occurs in the same fashion as in the previous test. The robot goal is again in the left bottom corner of the image, but here there are now two humans that move from the door to the stairs, in the opposite direction of the robot's desired trajectory.

After the leader is chosen, using the same approach as in the previous test, the robot starts to follow him/her. As the leader approaches the two humans moving in the opposite direction, they naturally give room for him/her to pass. As the robot is closely following the path taken by the leader, it is able to continue to move without the need to take evasive

E. Discussion

The tests assessed the capability of the system to predict the goal of real moving agents, as well as the ability of the designed algorithm to properly follow a chosen leader, while avoiding other dynamic agents.

Results show that the leader following algorithm makes a proper choice of a leader, based on a probabilistic approach for goal prediction, even when the initial movement and is not directed toward his/her goal. This is an important advantage of a probabilistic approach for goal detection, based on previous knowledge of the most common trajectories in the environment.

The navigation technique employed to follow the leader's path continuously explore the surrounding space for alternative trajectories. In the case the leader is lost of the path being followed becomes blocked, the branches that were generated in different directions can be used to find a new path, without the need of replanning from scratch. This is a very important characteristic while navigating in dynamic environments, specially due to time constraints.

The resulting behavior is that the robot is able to follow the leader while at the same time exploring the open spaces to accommodate new and unpredicted situations.

In the second test scenario, the advantages of following a leader in a dynamic environment becomes evident. Classical approaches that would attempt to plan a trajectory taking into account the predicted motion of the incoming humans would fail to find an optimal solution, as a straight line to the robot's goal would be blocked. However, as the robot follows a human that is able to correctly assume that the persons moving in the opposite direction will adapt their movement to avoid a collision, it is able to follow a straight trajectory to the goal. The result of this experiment clearly shows the benefit of the proposed technique, as the robot follows an optimal trajectory as a consequence of following a leader that has a better understanding on how to behave in such situations.

V. CONCLUSIONS

This work presents a method to take advantage of human motion in dynamic environments by selecting and following a leader. Its has two main contributions. The first is the methodology used to leader selection, which takes into account the typical paths in an environment and provides a probabilistic inference of subject's goal. The second is the modification of an algorithm designed for motion planning in dynamic environment, in order to adapt it to the task of leader following, while maintaining its original characteristics.

Tests used real data for the leader selection part, while the leader following algorithm was tested in a simulated and real environment. The results validated the proposed approach, with the robot being able to properly identify leaders among several subjects and follow him/her until its desired goal. It is our belief that the prediction algorithm can be further improved if it takes into account not only the position of the agents but also their speed and orientation, making it able to anticipate even more the future motion and goals of persons.

New ideas arose throughout this work, as the possibility to find leader that help the robot in portions of its path, which would required a more refined technique of leader selection, with the robot alternating between aided and unaided navigation.

Experiments will continue in different scenarios, with more tests in specific situations as leader obstruction/loss and also crowded environments. Currently a crowd simulator framework is in development, where extensive testing can take place and the impact of the robot behavior in the humans around can be taken into account.

Fig. 1 .

 1 Fig. 1. Comparison of goal prediction using simple extrapolation and the actual goals of two moving agents

Fig. 2 .Fig. 3 .

 23 Fig. 2. GHMM learning algorithm overview (adapted from [5])

  goalP red ← GHM M (agents)7:for i = 1 → agents.size() do 8: d ← Distance(goal, goalP red[i]) if f oundLeader = true then 17: path ← trackP ath(leader) 18: if leader is closer then 19:subgoal ← getSubGoal(path)

Fig. 5 .

 5 Fig. 5. INRIA Rhône Alpes entrance hall

Fig. 7 .

 7 Fig. 7. Real data used in the GHMM initial training

Fig. 8 .

 8 Fig. 8. Results of a typical test of leader detection and following

Fig. 9 .

 9 Fig. 9. Leader following allows the robot to avoid two incoming persons

ACKNOWLEDGMENT

We thank Dizan Vazques, Jorge Rios-Martinez and Arturo Escobedo for very useful insights and discussions.

This research has been funded by the INRIA PAL Project and by Fundac ¸ão para a Ciência e a Tecnologia SFRH/BD/46604/2008.