¸ois Caron
email: francois.caron@stats.ox.ac.uk

Yee Whye

Thomas Brendan Murphy

BAYESIAN NONPARAMETRIC PLACKETT-LUCE MODELS FOR THE ANALYSIS OF PREFERENCES FOR COLLEGE DEGREE PROGRAMMES

In this paper we propose a Bayesian nonparametric model for clustering partial ranking data. We start by developing a Bayesian nonparametric extension of the popular Plackett-Luce choice model that can handle an infinite number of choice items. Our framework is based on the theory of random atomic measures, with prior specified by a completely random measure. We characterise the posterior distribution given data, and derive a simple and effective Gibbs sampler for posterior simulation. We then develop a Dirichlet process mixture extension of our model and apply it to investigate the clustering of preferences for college degree programmes amongst Irish secondary school graduates. The existence of clusters of applicants who have similar preferences for degree programmes is established and we determine that subject matter and geographical location of the third level institution characterise these clusters.

1. Introduction. In this paper we consider partial ranking data consisting of ordered lists of the top-m items among a set of objects. Data in the form of partial rankings arise in many contexts. For example, in this paper we shall consider data pertaining to the top ten preferences of Irish secondary school graduates who are applying to undergraduate degree programmes offered in Irish third level institutions. The third level institutions consist of universities, institutes of technologies and private colleges. This application is described in detail in Section 2.

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 a popular model for modeling such partial rankings of a finite collection of M items. It has found many applications, including choice modeling [START_REF] Luce | The choice axiom after twenty years[END_REF][START_REF] Chapman | Exploiting Rank Ordered Choice Set Data within the Stochastic Utility Model[END_REF], sport ranking [START_REF] Hunter | MM algorithms for generalized Bradley-Terry models[END_REF], and voting [START_REF] Gormley | Exploring voting blocs with the Irish electorate: a mixture modeling approach[END_REF]. Diaconis (1988, Chapter 9) provides detailed discussions on the statistical foundations of this model.

In the Plackett-Luce model, each item k ∈ [M ] = {1, . . . , M } is assigned a positive rating parameter w k , which represents the desirability or rating of a product in the case of choice modeling, or the skill of a player in sport rankings. The Plackett-Luce model assumes the following generative story for a top-m list ρ = (ρ 1 , . . . , ρ m ) of items ρ i ∈ [M ]: At each stage i = 1, . . . , m, an item is chosen to be the ith item in the list from among the items that have not yet been chosen, with the probability that ρ i is selected being proportional to its desirability w ρ i . The overall probability of a given partial ranking ρ is then (1)

P (ρ) = m i=1 w ρ i M k=1 w k - i-1 j=1 w ρ j ,
with the denominator in (1) being the sum over all items not yet selected at stage i.

In many situations the collection of available items can be very large and/or potentially unknown. In this case a nonparametric approach can be sensible, where the pool of items is assumed to be infinite and the model allows for the possibility of items not observed in previous top-m lists to appear in future ones. A naïve approach, building upon recent work on Bayesian inference for the (finite) Plackett-Luce model and its extensions [START_REF] Gormley | A grade of membership model for rank data[END_REF][START_REF] Guiver | Bayesian inference for Plackett-Luce ranking models[END_REF][START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF], is to first derive a Markov chain Monte Carlo sampler for the finite model, then to "take the infinite limit" of the sampler, where the number of available items becomes infinite, but such that all unobserved items are grouped together for computational tractability.

Such an approach, outlined in Section 3, is reminiscent of a number of previous approaches deriving the (Gibbs sampler for the) Dirichlet process mixture model as the infinite limit of (a Gibbs sampler for) finite mixture models [START_REF] Neal | Bayesian Mixture Modeling[END_REF][START_REF] Rasmussen | The Infinite Gaussian Mixture Model[END_REF][START_REF] Ishwaran | Exact and Approximate Sum-Representations for the Dirichlet Process[END_REF]. Although intuitively appealing, this is not a satisfying approach since it is not clear what the underlying nonparametric model actually is, as it is actually the algorithm whose infinite limit was taken. It also does not directly lead to more general and flexible nonparametric models with no obvious finite counterpart, nor does it lead to alternative perspectives and characterisations of the same model, or resultant alternative inference algorithms. [START_REF] Orbanz | Construction of Nonparametric Bayesian Models from Parametric Bayes Equations[END_REF] further investigates the approach of constructing nonparametric Bayesian models from finite dimensional parametric Bayesian models. [START_REF] Caron | Bayesian nonparametric models for ranked data[END_REF] recently proposed a Bayesian nonparametric Plackett-Luce model based on a natural representation of items along with their ratings as an atomic measure. Specifically, the model assumes the existence of an infinite pool of items {X k } ∞ k=1 , each with its own rating parameter, {w k } ∞ k=1 . The atomic measure then consists of an atom located at each X k with a mass of w k :

(

2) G = ∞ k=1 w k δ X k .
The probability of a top-m list of items, say (X ρ 1 , . . . , X ρm ), is then a direct extension of the finite case (1):

(3) P (X ρ 1 , . . . ,

X ρm |G) = m i=1 w ρ i ∞ k=1 w k - i-1 j=1 w ρ j .
Using this representation, note that the top item X ρ 1 in the list is simply a draw from the probability measure obtained by normalising G, while subsequent items in the top-m list are draws from probability measures obtained by first removing from G the atoms corresponding to previously picked items and normalising. Described this way, it is clear that the Plackett-Luce model is none other than a partial size-biased permutation of the atoms in G [START_REF] Patil | Diversity as a Concept and its Implications for Random Communities[END_REF], and the existing machinery of random measures and exchangeable random partitions [START_REF] Pitman | Combinatorial stochastic processes[END_REF]Lijoi and Prünster, 2010) can be brought to bear on our problem.

For example, we may use a variety of existing stochastic processes to specify a prior over the atomic measure G. [START_REF] Caron | Bayesian nonparametric models for ranked data[END_REF] considered the case, described in Section 4, where G is a gamma process. This is a completely random measure [START_REF] Kingman | Completely Random Measures[END_REF]Lijoi and Prünster, 2010) with gamma marginals, such that the corresponding normalised probability measure is a Dirichlet process [START_REF] Ferguson | A Bayesian Analysis of Some Nonparametric Problems[END_REF]. They showed that with the introduction of a suitable set of auxiliary variables, it is possible to characterise the posterior law of G given observations of topm lists distributed according to (3). A simple Gibbs sampler can then be derived to simulate from the posterior distribution which corresponds to the infinite limit of the Gibbs sampler for finite models. In Appendix, we show that the construction can be extended from gamma processes to general completely random measures, and we discuss extensions of the Gibbs sampler to this more general case.

In Section 5 we describe a Dirichlet process mixture model [START_REF] Ferguson | A Bayesian Analysis of Some Nonparametric Problems[END_REF][START_REF] Lo | On a class of Bayesian Nonparametric estimates: I. Density Estimates[END_REF] for heterogeneous partial ranking data, where each mixture component is a gamma process nonparametric Plackett-Luce model. As shown in Section 2, such a model is relevant for capturing heterogeneity in preferences for college degree programmes. As we will see, in this model it is important to allow the same atoms to appear across the different random measures of the mixture components, otherwise the model becomes degenerate with all observed items that ever appeared together in some partial ranking being assigned to the same mixture component. To allow for this, we use a tree-structured extension of the time varying model of [START_REF] Caron | Bayesian nonparametric models for ranked data[END_REF]. In Section 6 we apply this mixture model to the Irish college degree programme preferences data, showing that the model is able to recover clusters of students with similar and interpretable preferences.

Finally, we conclude in Section 7 with a discussion of the important contributions of this paper and proposals for future work.

2. Irish College Degree Programmes. Applications to college degree programmes in Ireland are handled by a centralised applications system called the College Application Office (CAO) (www.cao.ie); a degree programme involves studying a specific subject (broad or focussed) in a particular third level institution. The CAO handle applications for 35 different third level institutions including universities, institutes of technologies and private colleges. In autumn of each year, a list of all degree programmes for the subsequent year is made available to applicants. Quite often new degree programmes are added to the list of potential choices after the initial list has been published, thus meaning that the potential list of degree programme choices is evolving and not always completely known. Applications are completed early in the year in which the students plan to enter their college degree programme. The list of available degree programmes changes from year-to-year but has been generally growing in size year-on-year. Many degree programmes have a specific subject area, for example Mathematics, History or Computer Science, but others are more general, for example Science, Commerce or Arts. In the year 2000, which we are examining herein, there were 533 degree programmes available to be selected by the applicants. When students apply for degree programmes they rank up to ten degree programmes, in order of preference, from the list of all degree programmes that are being offered. Two examples of such applications for two different applicants are shown in Table 1.

Places in these degree programmes are allocated on the basis of the applicants' performance in the Irish Leaving Certificate examination. Students typically take between seven and nine subjects in the Leaving Certificate examination. Points between zero and one hundred are awarded for each applicant's best six subjects in the Leaving Certificate examination and the points are totalled to give an overall points score. The allocation of applicants to most degree programmes is solely on the basis of the applicant's points score and applicants with a high points score are more likely to get their high preference choices. The minimum points score of all applicants accepted into a degree programme is publicly available and is called the points requirement. It is worth mentioning that even though degree programmes may have required Leaving Certificate subjects and grades as part of the minimum entry requirements, the subjects used in the applicant's points score calculation can be any six Leaving Certificate subjects.

The college applications system in Ireland is much debated in the educational sector and it receives much attention in the Irish media. The debate has two main parts: one part of the debate is whether the current system of allocating points to students on the basis of a single Leaving Certificate examination is a fair method, especially when the points can be gained from any Leaving Certificate subjects; the other part of the debate explores the choice behaviour of the applicants and whether students are choosing degree programmes in a coherent manner. We focus on the applicant's choices which are core to the second part of the debate.

Many people feel that students don't necessarily pick degree programmes on the basis of the courses offered but that they choose on other grounds, like the perceived prestige of the degree programme. However, other factors like geographical location of the third level institution may also have an impact on the applicant's choice behaviour. The two example applications in Table 1 illustrate that a number of factors influence applicants choices. The first applicant has selected degree programmes in medicine and other health sciences, so their choices appear to be largely based on the course material. However, the second application includes a wide variety of different degree programmes; the applicant's first choice degree programme leads to a career in Primary Teaching whereas the other degree programmes are in different areas. However, the institutions that have been chosen are geographically close (within 100 km).

In the year 1997, the Department of Education and Science, commissioned a review of the Irish college applications system. A report [START_REF] Hyland | Commission on the Points System: Final Report and Recommendations[END_REF]) reviewed the current system and made some recommendations concerning the future of the system. In addition, four research reports were published, one of which [START_REF] Tuohy | Demand for Third-Level Places[END_REF] examined the applicant's choices. [START_REF] Tuohy | Demand for Third-Level Places[END_REF] used a number of exploratory data analysis techniques to investigate the degree programmes selected, but without reference to the preference ordering, and he found that subject matter was an important factor in applicant choices. More recently, [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF] used a finite mixture of Plackett-Luce models to find clusters of applications with similar choice profiles. They fitted their model using maximum likelihood and chose the number of mixture components using the Bayesian Information Criterion (BIC). Their results also indicated that subject matter and geographical location were strong determinants of student choices. However, the model fitting paradigm used in their analysis could not find small clusters of applicants because of the manner that BIC penalises each additional mixture component. Further, McNicholas (2007) used association rule mining to further explore college applicant choices, but he restricted his attention to degree programme choice combinations that were selected by at least 0.5% of the applicants; thus, that analysis emphasised only high frequency choice behaviour.

O'Connell, [START_REF] O'connell | Who went to college in 2004? A national survey of new entrants to higher education[END_REF] conducted a survey of new college entrants (as opposed to applicants) in 2004 and found that the choice of college where they commenced their degree programme was influenced primarily by reputation and geographical location of the third level institution, and that the choice of degree programme was influenced by intrinsic interest in the subject matter and to a lesser extent future career prospects. Whilst that study only looks at students who entered college and the degree programme that they ultimately studied, it provides a further insight into the factors that influence choice of degree programme.

We investigate the complete degree programme choice data for the year 2000 cohort of applications to the College Application Office; these data correspond to top-10 rankings of college degree programmes for 53757 applicants. The model proposed herein has a number of appealing properties because it can account for choosing from the large number of degree programmes on offer, it allows for small differences in preference between degree programmes, it facilitates discovering large and small clusters of applicants with similar preferences and the fitting in the Bayesian paradigm facilitates a deep exploration of the clustering and co-clustering of applicants.

3. An extension of the Plackett-Luce model to countably infinite choice sets. We start this section with a review of a Bayesian approach to inference in finite Plackett-Luce models [START_REF] Gormley | A grade of membership model for rank data[END_REF][START_REF] Guiver | Bayesian inference for Plackett-Luce ranking models[END_REF][START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF], and taking the infinite limit to arrive at a nonparametric model. This will give good intuitions for how the model operates, before we rederive the same nonparametric model more formally in the next section using gamma processes.

Recall that we have M choice items indexed by [M ] = {1, . . . , M }, with item k ∈ [M ] having a positive desirability parameter w k . We will suppose that our data consists of L partial rankings of the M choice items, with the th ranking being denoted ρ = (ρ 1 , . . . , ρ m ), for = 1, . . . , L, where each ρ i ∈ [M ]. For notational simplicity we assume that all the partial rankings are of length m.

3.1.

Finite Plackett-Luce model with gamma prior. As noted in the introduction, the Plackett-Luce model constructs a partial ranking ρ = (ρ 1 , . . . , ρ m ) iteratively. At the ith stage, with i = 1, 2, . . . , m, we pick ρ i as the ith item from among those not yet picked with probability proportional to w ρ i . The probability of the partial ranking ρ is then as given in ( 1). An alternative Thurstonian interpretation, which will be important in the following, is as follows: For each item k let z k be exponentially distributed with rate w k :

z k ∼ Exp(w k )
Thinking of z k as the arrival time of item k in a race, let ρ i be the index of the ith item to arrive (the index of the ith smallest value among (z k ) M k=1 ). The resulting probability of the first m items to arrive being ρ can be shown to be the probability (1) from before. In this interpretation (z k ) can be understood as latent variables, and the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] can be applied to derive an algorithm to find a ML setting for the parameters (w k ) M k=1 given multiple partial rankings. Unfortunately the posterior distribution of (z k ) M k=1 given ρ is difficult to compute, so we can instead consider an alternative parameterisation: Let Z i be the waiting time for the ith item to arrive after the i -1th item. That is,

Z i = z ρ i -z ρ i-1
with z ρ 0 defined to be 0. Then it is easily seen that the joint probability of the observed partial rankings, along with the alternative latent variables (Z i ), is:

P ((ρ ) L =1 , ((Z i ) m i=1 ) L =1 |(w k ) M k=1 ) = L =1 m i=1 w ρ i exp   -Z i   M k=1 w k - i-1 j=1 w ρ j     (4)
In particular, the posterior of (Z i ) m i=1 is simply factorised, with

Z i |(ρ ) L =1 , (w k ) M k=1 ∼ Exp   M k=1 w k - i-1 j=1 w ρ j  
being exponentially distributed. The M step of the EM algorithm can be easily derived as well. The resulting algorithm was first proposed by [START_REF] Hunter | MM algorithms for generalized Bradley-Terry models[END_REF] as an instance of the MM (majorisationmaximisation) algorithm [START_REF] Lange | Optimization transfer using surrogate objective functions (with discussion)[END_REF] and its re-interpretation as an EM algorithm was recently given by [START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF].

Taking a further step, we note that the joint probability ( 4) is conjugate to a factorised gamma prior over the parameters, say w k ∼ Gamma( α M , τ ) with hyperparameters α, τ > 0. Now Bayesian inference can be carried out, for example, using a variational Bayesian EM algorithm, or a Gibbs sampler. In this paper we shall consider only Gibbs sampling algorithms. By regrouping the terms in the exponential in (4), the parameter updates are derived to be [START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF]:

w k |ρ, (Z i ), (w k ) k =k ∼ Gamma α M + n k , τ + L =1 m i=1 δ ik Z i (5)
where n k is the number of occurrences of item k among the observed partial rankings, and

δ ik = 0 if there is a j < i with ρ j = k, 1 otherwise.
Note that the definitions of n k and δ ik slightly differ from those in [START_REF] Hunter | MM algorithms for generalized Bradley-Terry models[END_REF] and [START_REF] Caron | Efficient Bayesian inference for generalized Bradley-Terry models[END_REF]. In these articles, the authors consider full m-rankings of subsets of [M ] whereas we consider here partial top-m rankings of all M items. 3.2. Taking the infinite limit. A Gibbs sampler for a nonparametric Plackett-Luce model can now be easily derived by taking the limit as the number of choice items M → ∞. If item k has appeared among the observed partial rankings, the limiting conditional distribution ( 5) is well defined since n k > 0. For items that did not appear in the observations, ( 5) becomes degenerate at 0. Instead we can define w * = k:n k =0 w k to be the total desirability among all the infinitely many unobserved items. Making use of the fact that sums of independent gammas with the same scale parameter is a gamma with shape parameter given by the sum of the shape parameters,

w * |ρ, (Z i ), (w k ) k:n k >0 ∼ Gamma α, τ + L =1 m i=1 Z i
The resulting Gibbs sampler alternates between updating the latent variables (Z i ), and updating the desirabilities of the observed items (w k ) k:n k >0 and of the unobserved ones w * .

This nonparametric model allows us to estimate the probability of seeing new items appearing in future partial rankings in a coherent manner. While intuitive, the derivation is ad hoc in the sense that it arises as the infinite limit of the Gibbs sampler for finite Plackett-Luce models, and is unsatisfying as it did not directly capture the structure of the underlying infinite dimensional object, which we will show in the next section to be a gamma process.

4. A Bayesian nonparametric Plackett-Luce model based on the gamma process. Let X be a measurable space of choice items. In the case of college applications, the space X is the space of all possible Irish programme courses. A gamma process is a completely random measure over X with gamma marginals. Specifically, it is a random atomic measure of the form (2), such that for each measurable subset A, the (random) mass G(A) is gamma distributed. Assuming that G has no fixed atoms (that is, for each element x ∈ X we have G({x}) = 0 with probability one) and that the atom locations {X k } are independent of their masses {w k } (that is, the gamma process is homogeneous), it can be shown that such a random measure can be constructed as follows (Kingman, 1967, Chapter 9): each X k is iid according to a base distribution H (which we assume is non-atomic with density h(x)), while the set of masses {w k } is distributed according to a Poisson process over R + with mean intensity where α > 0 is the concentration parameter and τ > 0 the inverse scale. We write this as G ∼ Γ(α, τ, H). Under this parametrisation, we have that G(A) ∼ Gamma(αH(A), τ ). λ(w)h(x) is known as the Lévy intensity of the homogeneous CRM G. The jump part λ(w) of the Lévy intensity verifies the necessary condition

λ(w) = αw -1 e -wτ
(6) ∞ 0 (1 -exp(-w))λ(w)dw < ∞
and plays a significant role in characterising the properties of the gamma process.

We shall interpret each atom X k as a choice item, with its mass w k > 0 corresponding to the desirability parameter. The Thurstonian view described in the finite model can be easily extended to the nonparametric one, where a partial ranking (X ρ 1 , . . . , X ρm ) can be generated as the first m items to arrive in a race. In particular, for each atom X k let z k ∼ Exp(w k ) be the time of arrival of X k and X ρ i the ith item to arrive. The first m items to arrive (X ρ 1 , . . . , X ρm ) then constitutes our partial ranking, with probability as given in (3). This construction is depicted on Figure 1. The top row of Figure 2 visualises some top-5 rankings generated from the model, with τ = 1 and different values of α. Figure 3 shows the mean number of items appearing in L top-m rankings. For m = 1, one recovers the well-known result on the number of clusters for a Dirichlet process model.

Again reparametrising using inter-arrival durations, let imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014 associated latent variables can be derived to be:

Z i = z ρ i -z ρ i-1 for i = 1, 2, . . . (with z ρ 0 = 0).
P ((X ρ 1 , . . . , X ρm ), (Z 1 , . . . , Z m )|G) =P ((z ρ 1 , . . . , z ρm ), and z k > z ρm for all k ∈ {ρ 1 , . . . , ρ m }) (7) = m i=1 w ρ i exp(-w ρ i z ρ i ) k ∈{ρ 1 ,...,ρm} exp(-w k z ρm ) = m i=1 w ρ i exp -Z i ∞ k=1 w k - i-1 j=1 w ρ j
Marginalising out (Z 1 , . . . , Z m ) gives the probability of (X ρ 1 , . . . , X ρm ) as in (3). Further, conditional on ρ = (ρ i ) m i=1 it is seen that the inter-arrival durations Z 1 . . . Z m are mutually independent, with

Z i |(X ρ 1 , . . . , X ρm ), G ∼ Exp ∞ k=1 w k - i-1 j=1 w ρ j
In the next section we shall characterise the posterior distribution over G given observed partial rankings and their associated latent variables. We end this subsection with two observations. Firstly, note that the jump part λ(w) of the Lévy intensity of the gamma process satisfies the following property:

∞ 0 λ(w)dw = ∞, (8)
This property is equivalent (via Campbell's Theorem) to the fact that there are an infinite number of atoms in G with probability one. In other words that we are dealing with a nonparametric model with an infinite number of choice items. It is also a necessary and sufficient condition for the homogeneous CRM G to have finite and strictly positive total mass 0 < G(X) < ∞ (Regazzini, Lijoi and Prünster, 2003). It therefore ensures that the generative Plackett-Luce probability (3) is well defined.

The second observation is with regard to a subtle but important difference between the atomic measure approach described in this section and the finite Plackett-Luce model of the previous section. In particular, here we specified the choice items X k as locations in a space X with a prior given by the base distribution H, while in the finite Plackett-Luce model we simply index the M choice items using 1, . . . , M . One may wonder if it is possible to simply index the infinitely many choice items using the natural numbers, and dispense with the atom locations {X k } altogether. This turns out to be impossible, if we were to make the following reasonable assumptions: That item desirabilities are a priori mutually independent, that they are positive with probability one, and that item desirabilities do not depend on the index of their corresponding items. With these assumptions, along with an infinite number of choice items, it is easy to see that the sum of all item desirabilities will be infinite with probability one, so that the Plackett-Luce generative model becomes ill-defined. Using the atomic measure approach, it is possible to satisfy all assumptions while making sure the Plackett-Luce generative model is well-defined. Note that the atoms locations X k are just used for modelling purposes. When considering inference, they are assumed to be known and need not to be defined explicitly so that to make inference on the item desirabilities. 4.1. Posterior characterisation. In this section we develop a characterisation of the posterior law of G under a gamma process prior and given Plackett-Luce observations consisting of L partial rankings. Posterior characterisation for our model is a variation of posterior characterisation for normalised random measures in density estimation (Prünster, 2002;[START_REF] James | Poisson process partition calculus with applications to exchangeable models and Bayesian nonparametrics[END_REF]James, Lijoi and Prünster, 2009;Lijoi and Prünster, 2010). We shall denote the th partial ranking as Y = (Y 1 , . . . , Y m ), where each Y i ∈ X. Note that previously our partial rankings (X ρ 1 , . . . , X ρm ) were denoted as ordered lists of the atoms in G. Since G is unobserved here, this is no longer possible, so we instead simply use a list of observed choice items (Y 1 , . . . , Y m ). Re-expressing the conditional distribution (3) of Y given G, we have:

P (Y |G) = m i=1 G({Y i }) G(X\{Y 1 , . . . , Y i-1 })
In addition, for each , we will also introduce a set of auxiliary variables Z = (Z 1 , . . . , Z m ) (the inter-arrival times) that are conditionally mutually independent given G and Y , with:

Z i |Y , G ∼ Exp(G(X\{Y 1 , . . . , Y i-1 })) (9)
The joint probability of the item lists and auxiliary variables is then (c.f. ( 7)):

P ((Y , Z ) L =1 |G) = L =1 m i=1 G({Y i }) exp(-Z i G(X\{Y 1 , . . . , Y i-1 }))
Note that under the generative process described in Section 4, there is positive probability that an item appearing in a list Y appears in another list Y with = . Denote the unique items among all L lists by X * 1 , . . . , X * K , and for each k = 1, . . . , K let n k be the number of occurrences of X * k among the item lists. Finally define occurrence indicators

δ ik = 0 if ∃j < i with Y j = X * k ; 1 otherwise. ( 10
)
Then the joint probability under the nonparametric Plackett-Luce model is:

P ((Y , Z ) L =1 |G) = K k=1 G({X * k }) n k × L =1 m i=1 exp(-Z i G(X\{Y 1 , . . . , Y i-1 })) = exp -G(X) i Z i K k=1 G({X * k }) n k exp -G({X * k }) i (δ ik -1)Z i (11)
Taking expectation of ( 11) with respect to G gives:

Theorem 1 The marginal probability of the L partial rankings and latent variables is:

P ((Y , Z ) L =1 ) = e -ψ( i Z i ) K k=1 h(X * k )κ n k , i δ ik Z i (12)
where ψ(z) is the Laplace transform of λ(w),

ψ(z) = -log E e -zG(X) = ∞ 0 (1 -e -zw )λ(w)dw = α log 1 + z τ
and κ(n, z) is the nth moment of the exponentially tilted intensity λ(w)e -zw :

κ(n, z) = ∞ 0 w n e -zw λ(w)dw = α (z + τ ) n Γ(n)
imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014

The proof, using the Poisson process characterisation of completely random measures and the Palm formula (James, Lijoi and Prünster, 2009), is given in the appendix.

Another application of the Palm formula (James, Lijoi and Prünster, 2009) now allows us to derive a posterior characterisation of G. The posterior CRM can be decomposed as the sum of a CRM with fixed atoms and a CRM whose jump part of the Lévy intensity is updated to λ * (w) in a conjugate fashion, similar to deriving a conjugate posterior for a parametric distribution.

Theorem 2 Given the observations and associated latent variables (Y , Z ) L =1 , the posterior law of G is also a gamma process, but with atoms with both fixed and random locations. Specifically,

G|(Y , Z ) L =1 = G * + K k=1 w * k δ X * k ( 13
)
where G * and w * 1 , . . . , w * K are mutually independent. The law of G * is still a gamma process,

G * |(X , Z ) L =1 ∼ Γ(α, τ , H), τ * = τ + i Z i
while the masses have distributions,

w * k |(Y , Z ) L =1 ∼ Gamma n k , τ + i δ ik Z i
Proof. Let f : X → R be measurable with respect to H. Then the characteristic functional of the posterior G is given by:

E[e -f (x)G(dx) |(Y , Z ) L =1 ] = E[e -f (x)G(dx) P ((Y , Z ) L =1 |G)] E[P ((Y , Z ) L =1 |G)] (14)
The denominator is as given in Theorem 1, while the numerator is obtained using the same Palm formula technique as Theorem 1, with the inclusion of the term e -f (x)G(dx) . Some algebra then shows that the resulting characteristic functional of the posterior G coincides with that of (13).

The proof details are given in the appendix.

4.2. Gibbs sampling. Given the results of the previous section, a simple Gibbs sampler can now be derived, where all the conditionals are of known analytic form. In particular, we will integrate out all of G * except for its total mass w * * = G * (X). This leaves the latent variables to consist of the masses w * * , (w * k ) K k=1 and the latent variables

((Z i ) m i=1 ) L =1 .
The update for Z i is given by (9), while those for the masses are given in Theorem 2:

Gibbs update for Z i :

Z i |rest ∼ Exp w * * + k δ ik w * k Gibbs update for w * k : w * k |rest ∼ Gamma n k , τ + i δ ik Z i Gibbs update for w * * : w * * |rest ∼ Gamma α, τ + i Z i (15)
Note that the latent variables are conditionally independent given the masses and vice versa. Hyperparameters of the gamma process can be simply derived from the joint distribution in Theorem 1. Since the marginal probability of the partial rankings is invariant to rescaling of the masses, it is sufficient to keep τ fixed at 1. As for α, if a Gamma(a, b) prior is placed on it, its conditional distribution is still gamma: Note that this update was derived with w * * marginalised out, so after an update to α it is necessary to immediately update w * * via (15) before proceeding to update other variables.

Gibbs update for α: α|rest ∼ Gamma a + K, b + log 1 + i Z i τ imsart-
In Section C in Appendix, we show that the construction can be extended from gamma processes to general completely random measures, and we discuss extensions of the Gibbs sampler to this more general case. In particular, we show that a simple Gibbs sampler can still be derived for the generalised gamma class of completely random measures. Plackett-Luce Components. In this section we propose a mixture model for heterogeneous ranking data consisting of nonparametric Plackett-Luce components. Using the same data augmentation scheme, we show that an efficient Gibbs sampler can be derived, and apply the model to a dataset of preferences for Irish degree programmes by high school graduates.

Mixtures of Nonparametric

5.1. Statistical model. Assume that we have a set of L rankings (Y ) for ∈ [L] of top-m preferred items, and our objective is to partition these rankings into clusters of similar preferences. We consider the following Dirichlet process (DP) mixture model:

π ∼ GEM(γ) (16) c |π ∼ Discrete(π) for = 1, . . . , L, Y |c , G c ∼ PL(G c )
where GEM(γ) denotes the Griffiths-Engen-McCloskey (GEM) distribution [START_REF] Pitman | Combinatorial stochastic processes[END_REF] with concentration parameter γ (also known as the stick-breaking construction) and PL(G) denotes the nonparametric Plackett-Luce model parameterised by the atomic measure G described in Section 4. The jth cluster in the mixture model is parameterised by an atomic measure G j and has mixing proportion π j .

To complete the model, we have to specify the prior on the component atomic measures G j . An obvious choice would be to use independent draws from a gamma process Γ(α, τ, H) for each G j . This unfortunately does not work. The reason is because if H is smooth then different atomic measures will never share the same atoms. On the other hand, notice that all items appearing in some observed partial ranking has to come from the same Plackett-Luce model, thus has to appear as atoms in the corresponding atomic measure. Putting these two observations together, the result is that any observed pair of partial rankings that share a common item will have to be assigned to the same component, and the mixture model will degenerate to using a few very larger components only. In consequence the model will not capture the fine-scale preference structure that may be present in the partial rankings. This is a similar problem that motivated the hierarchical DP [START_REF] Teh | Hierarchical Dirichlet Processes[END_REF], and the solution there as in here is to allow different atomic measures to share the same set of atoms, but to allow different atom masses.

Our solution, which is different from [START_REF] Teh | Hierarchical Dirichlet Processes[END_REF], is to make use of the Pitt-Walker [START_REF] Pitt | Constructing stationary time series models using auxiliary variables with imsart-aoas ver[END_REF] dependence model for gamma processes. Consider a tree-structured model where there is a single root G 0 and each component atomic measure G j is a leaf which connects directly to G 0 . The Pitt-Walker model allows us to construct the dependence structure between the root G 0 and the leaves (G j ) such that each G j marginally follows a gamma process Γ(α, τ, H). At the root, G 0 is first given a gamma process prior:

G 0 ∼ Γ(α, τ, H)
Since G 0 is atomic, we can write it in the form: Now for each j, define a random measure U j with conditional law:

G 0 = ∞ k=1 w 0k δ X k imsart-
U j |G 0 = ∞ k=1 u jk δ X k u jk |G 0 ∼ Poisson(φw 0k ) (17)
where φ > 0 is a parameter which, as we shall see, governs the strength of dependence between G 0 and each G j . Note that since G 0 has finite total mass, U j consists only of a finite number of atoms with positive masses; the other atoms all have masses equal to zero. Using the same Palm formula method as Section 4.1, we can show the following proposition:

Proposition 3 Suppose the prior law of G 0 is Γ(α, τ, H) and U j has conditional law given by (17). The posterior law of G 0 given U j is then:

G 0 = G * 0 + ∞ k=1 w * 0k δ X k
where G * 0 and (w * 0k ) ∞ k=1 are all mutually independent. The law of G * 0 is given by a gamma process while the masses are conditionally gamma,

G * 0 |U j ∼ Γ(α, τ + φ, H) w * 0k |U j ∼ Gamma(u jk , τ + φ)
Note that if u jk = 0, we define w * 0k to be degenerate at 0, thus the posterior of G 0 consists of a finite number of atoms in common with U j , along with an infinite number of atoms (those in G * 0 ) not in common. The total mass of G * 0 has distribution Gamma(α, τ + φ). The idea, inspired by [START_REF] Pitt | Constructing stationary time series models using auxiliary variables with imsart-aoas ver[END_REF], is to define the conditional law of G j given G 0 and U j to be independent of G 0 and to coincide with the conditional law of G 0 given U j as in Proposition 3. In other words, define

G j = G * j + ∞ k=1 w * jk δ X k (18)
where G * j ∼ Γ(α, τ + φ, H) and w * jk ∼ Gamma(u jk , τ + φ) are mutually independent. Note that if u jk = 0, the conditional distribution of w * jk will be degenerate at 0. Hence G j has an atom at X k if and only if U j has an atom at X k , that is, if u jk > 0. In addition, it also has an infinite number of atoms (those in G * j ) which are in neither U j nor G 0 . Since the conditional laws of G j and G 0 given U j coincide, and G 0 has prior Γ(α, τ, H), it can be seen that G j will marginally follow the same law Γ(α, τ, H) as well. More compactly, we can write the dependence model as:

U j |G 0 ∼ Poisson(φG 0 ) (19) G j |U j ∼ Γ α + U j (X), τ + φ, αH + U j α + U j (X)
As a final observation, the parameter φ can be interpreted as controlling the strength of dependence between G 0 and each G j . Indeed it can be shown that

E[G j |G 0 ] = φ φ + τ G 0 + τ φ + τ H
so that larger φ corresponds to each G j being more similar to G 0 . Larger φ may also favour a larger number of clusters as similar partial rankings are more likely to be clustered in different groups.

Our construction to inducing sharing of atoms has a number of qualitative differences from that of the hierarchical DP [START_REF] Teh | Hierarchical Dirichlet Processes[END_REF]. Firstly, the marginal law of each G j is known: it is marginally a gamma process. For the hierarchical DP the marginal laws of the individual random measures are not of simple analytical forms. Since normalising a gamma process gives a DP, our construction can be used as an alternative method to induce sharing of atoms across multiple random measures, each of which still has marginal DP law. Secondly, in our construction only a finite number of atoms will be shared across random measures (though the number shared can be controlled by the dependence parameter φ), while in the hierarchical DP all infinitely many atoms are shared. In [START_REF] Caron | Bayesian nonparametric models for ranked data[END_REF] we used the Pitt-Walker construction for a different purpose: we constructed a dynamical nonparametric Plackett-Luce model, where at each time t, G t is a gamma process, with the Pitt-Walker construction used to define a Markov dependence structure for the sequence of random measures (G t ).

The structure of ( 16), with a DP mixture with each component specified by a random atomic measure, is reminiscent of the nested DP of [START_REF] Rodríguez | The nested Dirichlet process[END_REF] as well, though our model has an additional hierarchical structure allowing the sharing of atoms among different component measures. In this respect, it also shares similarities with the hierarchical Dirichlet process model of Müller, [START_REF] Üller | A method for combining inference across related nonparametric Bayesian models[END_REF].

We focused here on a DP mixture for its simplicity, with a single parameter γ tuning the clustering structure. The model can be generalised to more flexible random measures, such as Pitman-Yor processes [START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF] or normalised random measures (Regazzini, Lijoi and Prünster, 2003;Lijoi, Mena and Prünster, 2007). 5.2. Posterior characterisation and Gibbs sampling. Assume for simplicity we have observed L top-m partial ranking Y = (Y 1 , . . . , Y m ) (the following will trivially extend to partial rankings of differing sizes). We extend the results of Section 4 in characterising the posterior and developing a Gibbs sampler for the mixture model.

Let X * = (X * k ) K k=1 be the set of unique items observed among Y 1 , . . . , Y L . For each cluster index j, let n jk be the number of occurrences of item X * k among the set of item lists Y in cluster j, that is, where c = j. Let ρ = (ρ i ) m i=1 be defined such as Y = (X * ρ 1 , . . . , X * ρ m ), and δ ik be occurrence indicators similar to (10).

As in Section 4, the observed items X * will contain the set of fixed atoms in the posterior law of the atomic measures G 0 , (G j ). We write the masses of the fixed atoms as

w 0k = G 0 ({X * k }), w jk = G j ({X * k })
, while the total masses of all other random atoms are denoted w 0 * = G 0 (X\X * ) and w j * = G j (X\X * ). We also write u jk = U j ({X * k }) and u j * = U j (X\X * ). As before, we will introduce latent variables for each = 1, . . . , L and i = 1, . . . , m:

Z i |Y , c , G c ∼ Exp w c * + K k=1 δ ik w c k (20)
The overall graphical model is described in Figure 4.

Proposition 4 Given the partial rankings (Y ) and associated latent variables (Z i ), (u jk ), (u j * ), and cluster indicators (c ), the posterior law of G j is a gamma process with atoms with both fixed and random locations. Specifically,

G j |(Y ), (Z i ), (u jk ), (u j * ), (c ) = G * j + K k=1 w jk δ X * k
where G * j and w j1 , . . . , w jK are mutually independent. The law of G * j is a gamma process,

G * j |(Y ), (Z i ), (u jk ), (u j * ), (c ) ∼ Γ   α + u j * , τ + φ + |c =j m i=1 Z i , H   , ( 21
)
while the masses have distributions,

w jk |(Y ), (Z i ), (u jk ), (u j * ), , (c ) ∼ Gamma n jk + u jk , τ + φ + |c =j m i=1 δ ik Z i (22)
Note that if n jk + u jk = 0, then w jk = 0 and G j will not have a fixed atom at X * k . To complete the posterior characterisation, note that conditioned on G 0 and G j the variables u j1 , . . . , u jK and u j * are independent, with u jk dependent only on w 0k and w jk and similarly for u j * . The conditional probabilities are:

p(u jk |w 0k , w jk ) ∝ f Gamma (w jk ; u jk , τ + φ)f Poisson (u jk ; φw 0k ) (23) p(u j * |w 0 * , w j * ) ∝ f Gamma (w j * ; α + u j * , τ + φ)f Poisson (u j * ; φw 0 * ) (24)
where f Gamma is the density of a Gamma distribution and f Poisson is the probability mass function for a Poisson distribution. The normalising constants are available in closed form [START_REF] Mena | On a construction of Markov models in continuous time[END_REF]:

p(w jk |w 0k ) = exp(-φw 0k )1 w jk ,0 + I -1 2 w jk φw 0k (τ + φ) φ(τ + φ)w 0k w jk 1/2 exp (-φ(w jk + w 0k ) -τ w jk ) (25) p(w j * |w 0 * ) =I α-1 2 w j * φw 0 * (τ + φ) (τ + φ) α+1 2 w j * φw 0 * α-1 2 exp(-φ(w j * + w 0 * ) -τ w j * ) (26)
where 1 a,b = 1 if a = b, 0 otherwise, and I is the modified Bessel function of the first kind. It is therefore possible to sample exactly from the discrete distributions ( 23) and ( 24) using standard retrospective sampling for discrete distributions, see for example [START_REF] Papaspiliopoulos | Retrospective Markov chain Monte Carlo methods for Dirichlet process hierarchical models[END_REF]. Alternatively, we describe in the appendix a Metropolis-Hastings procedure that worked well in the applications. Armed with the posterior characterisation, a Gibbs sampler can now be derived. Each iteration of the Gibbs sampler proceeds in the following order (details are in appendix):

1. First note that the total masses G j (X) are not likelihood identifiable, so we introduce a step to improve mixing. We simply sample them from the prior:

G 0 (X) ∼ Gamma(α, τ ) U j (X)|G 0 (X) ∼ Poisson(φG 0 (X)) G j (X)|U j (X) ∼ Gamma(α + U j (X), τ + φ)
The individual atom masses (w jk , w j * ) are scaled along with the update to the total masses. Then the Poisson masses (u jk ), (u j * ) are updated using ( 23) and ( 24). 2. The concentration parameter α and the masses w 0 * , (w j * ) and (u j * ) associated with other unobserved items are updated efficiently using a forward-backward recursion detailed in the appendix. 3. The masses (w 0k ) and w 0 * of the atoms in G 0 are updated via an extension of Proposition 3.

In particular, for each item k = 1, . . . , K, the masses are conditionally independent with distributions:

w 0k |u 1:J,k , φ ∼ Gamma J j=1 u jk , Jφ + τ
while the total mass of the remaining atoms have conditional distribution:

w 0 * |u 1:J * , φ ∼ Gamma α + J j=1 u j * , Jφ + τ
4. The latent variables (Z i ) are updated as in ( 20). 5. Conditioned on (Z i ), (u jk ) and (u j * ), the masses (w jk ) are updated via ( 22), while the total mass of the unobserved atoms is w j * ∼ Gamma(α * j , τ * j ) from ( 21). 6. The mixture weights π and the allocation variables c are updated using a slice sampler for mixture models [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF][START_REF] Kalli | Slice sampling mixture models[END_REF]. 7. Finally, the scale parameter γ of the Dirichlet process is updated using [START_REF] West | Hyperparameter estimation in Dirichlet process mixture models Technical Report[END_REF] and the dependence parameter φ is updated by a Metropolis-Hastings step using ( 25) and ( 26) with the latent (u jk ) and (u j * ) marginalised out.

The resulting algorithm is a valid partially collapsed Gibbs sampler [START_REF] Van Dyk | Partially collapsed Gibbs samplers[END_REF]. Note, however, that permutations of the above steps could result in an invalid sampler. The computational cost scales as O(K ×J ×m×L) where J is the average number of clusters. However, it is possible to parallelise over the different items in the algorithm to obtain an algorithm that scales as O(J × m × L).

6. Application: Irish College Degree Programmes. We now consider the application of the proposed model to study the choices made by the 53757 degree programme applicants to the College Application Office (CAO) in the year 2000. 6.1. Model Set-Up & Implementation Details . The following flat priors are used for the hyperparameters

p(α) ∝ 1/α p(φ) ∝ 1/φ p(γ) ∝ 1/γ
We run the Gibbs sampler with N = 20000 iterations. In order to obtain a point estimate of the partition from the posterior distribution, we use the approach proposed by [START_REF] Dahl | Model-based clustering for expression data via a Dirichlet process mixture model[END_REF]. Let c (i) , i = 1, . . . , N be the Monte Carlo samples. The point estimate c is obtained by

c = arg min c (i) ∈{c (1) ,...,c (N ) } k (δ c (i) k c (i) -ζ k ) 2 imsart-
ζ k = 1 N N i=1 δ c (i) k c (i)
and δ k = 1 if k = , 0 otherwise. Given this partition c, we run a Gibbs sampler with 2000 iterations to obtain the posterior mean Plackett-Luce parameters for each cluster. Clusters are then reordered by decreasing size. Table 2 shows the sizes of the 26 clusters which have a size larger than 10. In addition, a coclustering matrix was computed based on the first MCMC run which records for each pair of students the probability of them belonging to the same cluster. Figure 5 shows the coclustering matrix to summarise the clustering of the 53757 students, where students are rearranged by their cluster membership (members of the first cluster first, then members of the second cluster, etc.).

6.2.

Results. An examination of the Plackett-Luce parameter for each cluster reveals that the subject matter of the degree programme is a strong determinant of the clustering of students (Table 2). For example, clusters 6, 18 and 23 are characterised as construction, medicine and law, respectively. Besides the type of degree, geographical location is a strong determinant of degree programme choice. Clusters 12, 21 and 22 are respectively concerned with applications to college degree programmes in Cork, Galway and Limerick. There is a lot of heterogeneity in the subject area of the college degree programmes for these clusters, as can be seen for example for the Cork cluster 12 in Table 5. A number of clusters are also defined by a combination of both subject area and location, for example, for clusters 7 and 8 in Tables 3 and4, which correspond to computer science respectively outside and inside Dublin.

As mentioned in Section 2, there is a common perception in the Irish society and media that students pick degree programme based on prestige rather than subject area. Another perception is that the points requirement for a degree programme is a measure of prestige; in fact the points requirement is determined by a number of factors including the number of available places, the number of applicants who list the degree programme in their top-10 preferences and the quality of the applicants who apply for the degree programme. Such a selection-by-prestige phenomenon should be evidenced by a cluster of students picking degree programmes in medicine and law, both of which have very high points requirements, but no such cluster was found. In fact, medicine and law applicants are clustered separately into clusters 18 and 23, respectively. Therefore, the clustering suggests that students are primarily picking degree programmes on the basis of subject area FIG 5. Coclustering matrix of the 53757 college applicants for the CAO data. The posterior probability that two applicants belong to the same cluster is indicated by a color between blue (0) and red (1). Applicants are arranged by their cluster membership, and the clusters are ordered by size. The clusters are described in Table 2. and geographical considerations; this finding is in agreement with the results found in [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF][START_REF] Mcnicholas | Association rule analysis of CAO data[END_REF].

It is also of interest to look at the variability of the student choices within each cluster. This can be quantified by the normalised entropy, which takes its values between 0 and 1, and defined for each cluster j by -K k=1 ( w jk log w jk ) -w j * log w j * log(K + 1)

where w jk are the averaged normalised weights of item k in cluster j obtained from the second MCMC run; the normalised entropy values for each cluster are reported in Table 2. A low value indicates low variability in the choices within a cluster, whereas a large value indicates a lot of variability. Interestingly, cluster 15 has very low normalised entropy, where 56% of the students in that cluster are likely to take one of the three most popular degree programmes of that cluster (Drumcondra, Froebel or Marina) as their first choice; these degree programmes are the main primary teacher education degree programmes in Dublin and thus many members of this cluster have a strong interest in teacher education as a degree choice. Further, there is much more variability in cluster 7, where students choices are spread across various computing degree programmes, and only 23% of the students are likely to take one of the three most popular degree programmes as their first choice. The coclustering matrix reveals some interesting connections between clusters, which have not been explored in previous analyses of the CAO data. For example, the plot reveals that a number of applicants have high probability of belonging to clusters 4 and 19 which are both in the arts. Cluster 4 is characterised by arts degrees which do not require the applicants to select their major in advance, whereas cluster 19 is characterised by arts degrees where the student needs to specify their major in advance. It is worth observing that the clusters are fairly well separated, and very few clusters exhibit the phenomenon of sharing applicants, which is further evidence that the applicants are only selecting degree programmes of a particular type (as described by the cluster names in Table 2).

Marginal Posterior distributions of the hyperparameters α, γ and φ are respectively in the ranges [3, 8], [2, 5] and[ 100,200]. Correlation parameter φ is rather high. This is due to the fact that some degree programmes, such as Arts in University College Dublin or Cork often appear in the top-ten list of applicants, whatever their main subject matter is. Parameter γ is associated to the number of clusters, which is around 35. Parameter α relates to the variability of the weights within clusters (and thus to the entropy of the clusters).

7. Discussion. We have proposed a Bayesian nonparametric Plackett-Luce model for ranked data. Our approach is based on the theory of completely random measures, where we showed that the Plackett-Luce generative model corresponds exactly to a size-biased permutation of the atoms in the random measure. We characterised the posterior distribution, and derived a simple MCMC sampling algorithm for posterior simulation. Our approach can be seen as a multi-stage generalisation of posterior inference in normalised random measures (Regazzini, Lijoi and Prünster, 2003;James, Lijoi and Prünster, 2009;[START_REF] Griffin | Posterior simulation of normalized random measure mixtures[END_REF][START_REF] Favaro | MCMC for normalized random measure mixture models[END_REF].

We also developed a nonparametric mixture model consisting of nonparametric Plackett-Luce components to model heterogeneity in partial ranking data. In order to allow atoms to be shared across components, we made use of the Pitt-Walker construction, which was previously only used to define Markov dynamical models. Applying our model to a dataset of preferences for Irish college degree programmes, we find interesting clustering structure supporting the observation that students were choosing programmes mainly based on subject area and geographical considerations.

It is worthwhile comparing our mixture model to another nonparametric mixture model, DPM-GM, where each component is a generalised Mallows model [START_REF] Busse | Cluster analysis of heterogeneous rank data[END_REF]Meilȃ and Bao, 2008;Meilȃ and Chen, 2010). In the generalised Mallows model the component distributions are characterised by a (discrete) permutation parameter whereas in the Plackett-Luce model the component distributions are characterised by a continuous rating parameter. Thus the Plackett-Luce model offers greater modelling flexibility to capture the strength of preferences for each item. On the other hand, the scale parameters in the generalised Mallows model can accommodate varying precision in the ranking. Additionally, inference for the generalised Mallows models can be difficult.

The mixture model established the existence of clusters of applicants with similar degree programme preferences and characterises these clusters and their coherence in terms of choices. The results support the previous hypotheses that subject matter and geographical location are the primary drivers of degree programme choice [START_REF] Gormley | Analysis of Irish third-level college applications data[END_REF][START_REF] Mcnicholas | Association rule analysis of CAO data[END_REF]. These factors are important because they reflect the intrinsic interest in the subject matter of the degree programmes and the economic and practical aspects of choosing a third level institution for study. The geographical location influence is further supported by results on acceptances to degree programmes (O'Connell, [START_REF] O'connell | Who went to college in 2004? A national survey of new entrants to higher education[END_REF] and studies on how students fund their education which found that 45% of Irish university students live in their family home [START_REF] Clancy | Financing Third-Level Students in Ireland[END_REF] and thus attend an institution that is geographically close by.

An interesting extension of the proposed model would be to consider inhomogeneous completely random measures, where the preferences would depend on a set of covariates (e.g. location).

Curie Intra-European Fellowship Programme1 .

APPENDIX A: PROOF OF THEOREM 1

The marginal probability ( 12) is obtained by taking the expectation of ( 11) with respect to G. Note however that ( 11) is a density, so to be totally precise here we need to work with the probability of infinitesimal neighborhoods around the observations instead, which introduces significant notational complexity. To keep the notation simple, we will work with densities, leaving it to the careful reader to verify that the calculations indeed carry over to the case of probabilities.

P ((Y , Z ) L =1 ) =E P ((Y , Z ) L =1 |G) =E e -G(X) i Z i K k=1 G({X * k }) n k e -G({X * k }) i (δ ik -1)Z i
The gamma prior on G = ∞ j=1 w j δ X j is equivalent to a Poisson process prior on N = ∞ j=1 δ (w j ,X j ) defined over the space R + × X with mean intensity λ(w)h(x). Then,

=E   e -wN (dw,dx) i Z i K k=1 ∞ j=1 w n k j 1(X j = X * k )e -w j i (δ ik -1)Z i   (27)
We now recall the Palm formula (see e.g. Bertoin (2006, Lemma 2.3)).

Proposition 5 Palm Formula. Let N be a Poisson process on S with mean measure ν. Let S p denote the set of point measures on S, f : S → [0, +∞[ and G : S × S p → [0, +∞[ be some measurable functional. Then we have the so-called Palm formula

(28) E S f (x)G(x, N )N (dx) = S E[G(x, N + dx)]f (x)ν(dx)
where the expectation is with respect to N .

Applying the Palm formula for Poisson processes to pull the k = 1 term out of the expectation,

= E   e -w(N +δ w * 1 ,x * 1 )(dw,dx) i Z i K k=2 ∞ j=1 w n k j 1(X j = X * k )e -w j i (δ ik -1)Z i   × (w * 1 ) n 1 h(X * 1 )e -w * 1 i (δ i1 -1)Z i λ(w * 1 )dw * 1 =E   e -wN (dw,dx) i Z i K k=2 ∞ j=1 w n k j 1(X j = X * k )e -w j i (δ ik -1)Z i   × h(X * 1 ) (w * 1 ) n 1 e -w * 1 i δ i1 Z i λ(w * 1 )dw *
Now iteratively pull out terms k = 2, . . . , K using the same idea, and we get:

=E e -G(X) i Z i K k=1 h(X * k ) (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k =e -ψ( i Z i) K k=1 h(X * k )κ n k , i δ ik Z i (29)
This completes the proof of Theorem 1.

APPENDIX B: PROOF OF THEOREM 2

The proof is essentially obtained by calculating the numerator and denominator of ( 14). The denominator is already given in Theorem 1. The numerator is obtained using the same technique with the inclusion of the term e f (x)G(dx) , which gives:

E e -f (x)G(dx) P ((Y , Z ) L =1 |G) =E e -(f (x)+ i Z i )G(dx) K k=1 h(X * k ) (w * k ) n k e -w * k (f (X * k )+ i δ ik Z i ) λ(w * k )dw * k
By the Lévy-Khintchine Theorem (using the fact that G has a Poisson process representation N ),

= exp -(1 -e -w(f (x)+ i Z i ) )λ(w)h(x)dwdx × K k=1 h(X * k ) (w * k ) n k e -w * k (f (X * k )+ i δ ik Z i ) λ(w * k )dw * k (30)
Dividing the numerator (30) by the denominator (29), the characteristic functional of the posterior G is:

E e -f (x)G(dx) |(Y , Z ) L =1 = exp -(1 -e -wf (x) )e -i Z i λ(w)h(x)dwdx × K k=1 h(X * k ) e -f (X * k ) (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k (w * k ) n k e -w * k i δ ik Z i λ(w * k )dw * k
Since the characteristic functional is the product of K + 1 terms, we see that the posterior G consists of K + 1 independent components, one corresponding to the first term above (G * ), and the others corresponding to the K terms in the product over k. Substituting the Lévy measure λ(w) for a gamma process, we note that the first term shows that G * is a gamma process with updated inverse scale τ * . The kth term in the product shows that the corresponding component is an atom located at X * k with density (w * k ) n k e -w * k i δ ik Z i λ(w * k ); this is the density of the gamma distribution over w * k in Theorem 2. This completes the proof.

APPENDIX C: GENERALISATION TO COMPLETELY RANDOM MEASURES

The posterior characterisation we have developed along with the Gibbs sampler can be easily extended to completely random measures (CRM) [START_REF] Kingman | Completely Random Measures[END_REF]Regazzini, Lijoi and Prünster, 2003;Lijoi and Prünster, 2010). To keep the exposition simple, we shall consider homogeneous imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014

CRMs without fixed atoms. These can be described, as for the gamma process before, with atom locations {X k } iid according to a non-atomic base distribution H, and with atom masses {w k } being distributed according to a Poisson process over R + with a general Lévy measure λ(w) which satisfies the constraints (8) leading to a normalisable measure G with infinitely many atoms. We will write G ∼ CRM(λ, H) if G follows the law of a homogeneous CRM with Lévy intensity λ(w) and base distribution H.

Both Theorems 1 and 2 generalise naturally to homogeneous CRMs. In fact the statements and the proofs in the appendix still hold with the more general Lévy intensity, along with its Laplace transform ψ(z) and moment function κ(n, z):

Theorem 1' The marginal probability of the L partial rankings and latent variables is:

P ((Y , Z ) L =1 ) = e -ψ( i Z i ) K k=1 h(X * k )κ n k , i δ ik Z i where ψ(z) is the Laplace transform of λ(w), ψ(z) = -log E e -zG(X) = ∞ 0
(1 -e -zw )λ(w)dw and κ(n, z) is the nth moment of the exponentially tilted Lévy intensity λ(w)e -zw :

κ(n, z) = ∞ 0
w n e -zw λ(w)dw

Theorem 2' Given the observations and associated latent variables (Y , Z ) L =1 , the posterior law of G is also a homogeneous CRM, but with atoms with both fixed and random locations. Specifically,

G|(Y , Z ) L =1 = G * + K k=1 w * k δ X * k
where G * and w * 1 , . . . , w * K are mutually independent. The law of G * is a homogeneous CRM with an exponentially tilted Lévy intensity:

G * |(X , Z ) L =1 ∼ CRM(λ , H) λ * (w) = λ(w)e -w i Z i
while the masses have densities:

P (w * k |(Y , Z ) L =1 ) = (w * k ) n k e -w * k i Z i λ(w * k ) κ(n k , i Z i )
.

Examples of CRMs that have been explored in the literature for Bayesian nonparametric modelling include the stable process [START_REF] Kingman | Random discrete distributions[END_REF], the inverse Gaussian process (Lijoi, Mena and Prünster, 2005), the generalised gamma process [START_REF] Brix | Generalized Gamma Measures and Shot-noise Cox Processes[END_REF], and the beta process [START_REF] Hjort | Nonparametric Bayes Estimators Based on Beta Processes in Models for Life history Data[END_REF]. The generalised gamma process forms the largest known simple and tractable family of CRMs, with the gamma, stable and inverse Gaussian processes included as subfamilies. It has a Lévy intensity of the form:

λ(w) = α Γ(1 -σ) w -1-σ e -τ w
where the concentration parameter is α > 0, the inverse scale is τ ≥ 0, and the index is 0 ≤ σ < 1.

The gamma process is recovered when σ = 0, the stable when τ = 0, and the inverse Gaussian when σ = 1/2. The Laplace transform and the moment function of the generalised gamma process are:

ψ(z) = α σ ((τ + z) σ -τ σ ) κ(n, z) = α (τ + z) n-σ Γ(n -σ) Γ(1 -σ) .
The Gibbs sampler developed for the gamma process can be generalised to homogeneous CRMs as well. Recall that given the observed partial rankings, the parameters consist of the ratings (w * k ) K k=1 of the observed items and the total ratings w * * of the unobserved ones, while the latent variables are (Z i ). A corollary of Theorems 1' and 2' which will prove useful is the joint probability of these along with the observed partial rankings:

P ((Y i , Z i ), (w * k ), w * * ) = e -w * * ( i Z i )) f (w * * ) K k=1 h(X * k )(w * k ) n k e -w * k ( i δ ik Z i ) λ(w * k ) (31)
where f (w) is the density (assumed to exist) of the total mass w * * under a CRM with the prior Lévy intensity λ(w). Note that integrating out the parameters (w * k ), w * * from (31) gives the marginal probability in Theorem 1'. From the joint probability (31), the Gibbs sampler can now be derived: Gibbs update for Z i :

Z i |rest ∼ Exp w * * + k δ ik w * k Gibbs update for w * k : P (w * k |rest) ∝ (w * k ) n k e -w * k i Z i λ(w * k ) Gibbs update for w * * : P (w * * |rest) ∝ e -w * * ( i Z i )) f (w * * )
To be concrete, consider the updates for a generalised gamma process. The conditional distribution for w * k can be seen to be Gamma(n k -σ, τ + i Z i ), while the conditional distribution for w * * can be seen to be an exponentially tilted stable distribution. This is not a standard distribution (nor does it have known analytic forms for its density), but can be effectively sampled using recent techniques [START_REF] Devroye | Random variate generation for exponentially and polynomially tilted stable distributions[END_REF]. Another approach is to marginalise out w * * first:

P ((Y i , Z i ), (w * k )) = e -ψ( i Z i ) K k=1 h(X * k )(w * k ) n k e -w * k ( i δ ik Z i ) λ(w * k )
The MCMC algorithm then consists of sampling the ratings (w * k ) and auxiliary variables (Z i ). Marginalising out w * * introduces additional dependencies among the latent variables Z i . Fortunately, since the Laplace transform for a generalised gamma process is of simple form, it is possible to update the latent variables (Z i ) using a variety of standard techniques, including Metropolis-Hastings, Hamiltonian Monte Carlo, or adaptive rejection sampling. For these techniques to work well we suggest reparametrising each Z i using its logarithm log Z i instead.

APPENDIX D: GIBBS SAMPLER FOR THE MIXTURE OF NONPARAMETRIC PLACKETT-LUCE COMPONENTS

Let J be the number of different values taken by c (number of clusters). Please note that the number of clusters is not set in advance and its value may change at each iteration. The Gibbs sampler proceeds with each of the following updates in turn:

1. a. Update G 0 (X) given α, then for j = 1, . . . , J, update G j (X) given (G 0 (X), α, φ, c) b. For j = 1, . . . , J, update (u j , u j * ) given (w 0 , w 0 * , w j , w j * , φ, α, c) 2. a. Update α given (Z, φ, c) b. Update w 0 * given (Z, φ, c, α) c. For j = 1, . . . , J, update u j * given (Z, φ, c, α, w 0 * ) d. For j = 1, . . . , J, update w j * given (Z, α, u j * , φ, c) where x 0 is defined above.

2.c) For j = 1, . . . , J, update u j * given (Z, φ, c, α, w 0 * ) We can sample from the full conditional which is given, for j = 1, . . . , J by

u j * |(Z, φ, c, α, w 0 * ) ∼ Poisson 1 + φ 1 + φ + Z j φw 0 *
where Z j is defined above.

2.d) For j = 1, . . . , J, update w j * given (Z, α, u j * , φ, c) We can sample from the full conditional which is given, for j = 1, . . . , J by w j * |u j * , Z, c, α ∼ Gamma α + u j * , τ + φ + Z j where Z j is defined above.

3) Update (w 0k ), w 0 * given (U 1:J , α) The allocation variables (c 1 , . . . , c L ) are updated using the slice sampling technique described in [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF][START_REF] Kalli | Slice sampling mixture models[END_REF][START_REF] Fall | Gibbs sampling methods for Pitman-Yor mixture models Technical Report[END_REF]. It builds on the introduction of additional latent slice variables, and does not require to set any truncation. For completeness, we briefly recall here the details of the sampler. From Eq. ( 16), we have

• Set k = k + 1 • Sample v k ∼ Beta(1, γ) • Set π k = π * v k k-1 j=J+1 (1 -v j )
• Sample G k given G 0 using Eq. ( 19) 2. For = 1, . . . , L, sample c from

p(c = k) ∝ 1(π k > ω ) PL(Y ; G c ) 7) Update γ given c
The scale parameter γ of the Dirichlet process is updated using the data augmentation technique of [START_REF] West | Hyperparameter estimation in Dirichlet process mixture models Technical Report[END_REF]. 8) Update φ given w 0 , w 0 * , w 1:J , w 1:J * , α, φ

We sample φ using a MH step. Propose φ * = φ exp(σε) where σ > 0 and ε ∼ N (0, 1). And accept it with probability

FIG 1 .

 1 FIG 1. Bayesian nonparametric Plackett-Luce model. Left: an instantiation of the atomic measure G encapsulating both the items and their ratings. Right: Arrival times z k and latent variables Z k = zρ k -zρ k-1 . The top 5 items are (ρ1, ρ2, . . . , ρ5).

FIG 2 .FIG 3 .

 23 FIG 2. Visualisation of top-5 rankings with rows corresponding to different rankings and columns to items sorted by size biased order. A lighter shade corresponds to a higher rank. Results are shown for a gamma process with λ(w) = αw -1 exp(-τ w) with τ = 1 and different values of α. The parameter α tunes the variability in the partial rankings. The larger α the higher the variability. As the probability of partial rankings 3 is invariant to rescaling of the weights, the scaling parameter τ has no effect on the partial rankings.

FIG 4 .

 4 FIG 4. Graphical model of the Dirichlet process mixture of nonparametric Plackett-Luce components. The variables at the top are hyperparameters, (ρ ) are the observed partial rankings, while the other variables are unobserved variables.

  We now consider sampling of u j * , j = 1, . . . , J. We can use a MH step. Sample w * j * ∼ Poisson(φw 0 * ) and accept with probabilitymin 1, f Gamma (u j * ; α + u * j * , τ + φ) f Gamma (u j * ; α + u * j * , τ + φ) 2.a) Update α given (Z, φ, c)We can sample from the full conditional which is given by α|(Z, γ, φ, c) ∼ Gamma (a + K, b + y 0 + log(1 + x 0 )) Update w 0 * given (Z, φ, c, α) We can sample from the full conditional which is given by w 0 * |(Z, φ, c, α) ∼ Gamma (α, τ + x 0 )

For

  each item k = 1, . . . , K, sample w 0k |u 1:J,k , φ ∼ Gamma = 1, . . . , L, update Z given (w c , w c * , c ) For = 1, . . . , L and i = 1, . . . m, sampleZ i |c, w, w * ∼ Exp w c , * + K k=1 δ ik w c ,k5) For j = 1, . . . , J, update (w jk ), w j * given (Z, α, u j , u j * , φ, c)For each cluster j = 1, . . . , J• For each item k = 1, . . . , K, samplew jk |u jk , {ρ |c = j} ∼ Gamma   n jk + u jk , τ + φ + + n jk > 0,otherwise, set w jk = 0. • Sample the total mass w j * |u j * , {ρ |c = j} ∼ Gamma   α + u j * , τ + φ + = 1, . . . , L, update c and the weights π given w 1:J , w 1:J *

  L(Y ; G k )where the π k admit the following stick-breaking representation (33)π 1 = v 1 , π k = v k j<k (1 -v j )where the v k are i.i.d. from Beta(1, γ). For each observation Y , slice sampling introduces latent variable ω such that the joint distribution of Y , ω and c is given by (34)f (Y , ω , c |π, G) = 1(ω < π c ) PL(Y ; G c )For simplicity, assume that the c take values in {1, 2, . . . , J}. Let µ k be the number of allocation variables taking value k ∈ {1, . . . , J}. The sampler samples ω and v as a block given c, then c given v and ω. 1. (a) Sample (π 1 , . . . , π J , π * ) ∼ Dirichlet(µ 1 , . . . , µ J , γ) (b) For = 1, . . . , L, sample ω ∼ Unif([0, π c ]) (c) Set k = J. While k j=1 π k < (1 -min(ω 1 , . . . , ω L ))

  j * |φ * , w 0 * ) p(w j * |φ, w 0 * ) K k=1 p(w jk |φ * , w 0k ) p(w jk |φ, w 0k )  

TABLE 1

 1 Two samples from the CAO preference data. Each rank observation is an ordered list of up to ten degree programmes.

	Rank CAO code	College	Degree Programme
	1	DN002	University College Dublin	Medicine
	2	GY501	NUI -Galway	Medicine
	3	CK701	University College Cork	Medicine
	4	DN006	University College Dublin	Physiotherapy
	5	TR053	Trinity College Dublin	Physiotherapy
	6	DN004	University College Dublin	Radiotherapy
	7	TR007	Trinity College Dublin	Clinical Speech
	8	FT223	Dublin IT	Human Nutrition
	9	TR084	Trinity College Dublin	Social Work
	10	DN007	University College Dublin	Social Science
	Rank CAO code		College	Degree Programme
	1	MI005	Mary Immaculate Limerick Education -Primary Teaching
	2	CK301	University College Cork	Law
	3	CK105	University College Cork	European Studies
	4	CK107	University College Cork	Language -French
	5	CK101	University College Cork	Arts

TABLE 2

 2 Description of the different clusters. The size of the clusters, the entropy and a cluster description are provided.

	Cluster Size	Entropy Description	Cluster Size	Entropy Description
	1	3325 0.72	Social Science/Tourism	14	1918 0.71	Engineering
	2	3214 0.71	Science	15	1835 0.48	Teaching/Arts
	3	3183 0.64	Business/Commerce	16	1835 0.68	Art/Music -Dublin
	4	2994 0.58	Arts	17	1740 0.71	Engineering -Dublin
	5	2910 0.63	Business/Marketing -Dublin	18	1701 0.55	Medicine
	6	2879 0.68	Construction	19	1675 0.70	Arts/Religion/Theology
	7	2803 0.66	CS -outside Dublin	20	1631 0.76	Arts/History -Dublin
	8	2225 0.67	CS -Dublin	21	1627 0.66	Galway
	9	2303 0.67	Arts/Social -outside Dublin	22	1392 0.70	Limerick
	10	2263 0.63	Business/Finance -Dublin	23	1273 0.65	Law
	11	2198 0.65	Arts/Psychology -Dublin	24	1269 0.72	Business -Dublin
	12	2086 0.63	Cork	25	1225 0.79	Arts/Bus./Language -Dublin
	13	2029 0.64	Comm./Journalism -Dublin	26	47	0.96	Mixed
	where the coclustering matrix ζ is obtained with

TABLE 3

 3 Cluster 7: Computer Science -outside Dublin

	Rank Aver. Norm. Weight	College	Degree Programme
	1	0.081	Cork IT	Computer Applications
	2	0.075	Limerick IT	Software Development
	3	0.072	University of Limerick	Computer Systems
	4	0.064	Waterford IT	Applied Computing
	5	0.061	Cork IT	Software Dev & Comp Net
	6	0.046	IT Carlow	Computer Networking
	7	0.038	Athlone IT Computer and Software Engineering
	8	0.036 University College Cork	Computer Science
	9	0.033	Dublin City University	Computer Applications
	10	0.033	University of Limerick	Information Technology

imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014

  * jk = 0 if u * jk = 0 otherwise sample w * jk ∼ Gamma(u jk , τ + φ). Accept (u * jk , w *

				jk )
	with probability
	min 1,	exp(-w * jk exp(-w jk	|c =j	m i=1 Z i )

imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014 then set w

imsart-aoas ver. 2012/08/31 file: aoas_final.tex date: January 14, 2014

 1 The contents reflect only the authors views and not the views of the European Commission.imsart-aoas ver.

2012/08/31 file: aoas_final.tex date: January 14, 2014

F. CARON ET AL.

Acknowledgements. The authors thank Igor Prünster for very helpful feedback on an earlier version of this work. F.C. acknowledges the support of the European Commission under the Marie

3. Update (w 0k ), w 0 * given (U 1:J , α) 4. For = 1, . . . , L, update Z given (w c , w c * , c ) 5. For j = 1, . . . , J, update (w j , w j * ) given (Z, α, u j , u j * , φ, c) 6. For = 1, . . . , L, update c and the mixture weights π given w 1:J , w 1:J * 7. Update γ given c 8. Update φ given w 0 , w 0 * , w 1:J , w 1:J * , α, φ

The step are now fully described.

1.a) Update G 0 (X) given α, then for j = 1, . . . , J, update G j (X) given (G 0 (X), α, φ, c)

We have G 0 (X)|α ∼ Gamma(α, τ )

and for j = 1, . . . , J G j (X) ∼ Gamma(α + M j , τ + φ)

where M j ∼ Poisson(φG 0 (X)).

1.b) For j = 1, . . . , J, update (u j , u j * ) given (w 0 , w 0 * , w j , w j * , φ, α, c)

Consider first the sampling of u j . We have, for j = 1, . . . , J and k = 1, . . If w jk = 0, we only have two possible moves: u jk = 0 or u jk = 1, given by the following probabilities

Note that the above Markov chain is not irreducible, as the probability is zero to go from a state (u jk > 0, w jk > 0) to a state (u jk = 0, w jk = 0), even though the posterior probability of this event is non zero in the case item k does not appear in cluster j. We can add such moves by jointly sampling (u jk , w jk ). For each k that does not appear in cluster j, sample u * jk ∼ Poisson(φw 0k )