Mame Diarra Fall
email: diarra.fall@univ-orleans.fr

Éric Barat
email: eric.barat@cea.fr

Gibbs Sampling Methods for Pitman-Yor Mixture Models

Keywords: Bayesian nonparametrics, Dirichlet process mixture model, Pitman-Yor process mixture model, Gibbs sampler, Slice sampling

We introduce a new sampling strategy for the two-parameter Poisson-Dirichlet process mixture model, also known as Pitman-Yor process mixture model (PYM). Our sampler is therefore applicable to the well-known Dirichlet process mixture model (DPM).

Inference in DPM and PYM is usually performed via Markov Chain Monte Carlo (MCMC) methods, specifically the Gibbs sampler. These sampling methods are usually divided in two classes: marginal and conditional algorithms. Each method has its merits and limitations. The aim of this paper is to propose a new sampler that combines the main advantages of each class. The key idea of the proposed sampler consists in replacing the standard posterior updating of the mixing measure based on the stick-breaking representation, with a posterior updating of [Pit96b] which represents the posterior law under a Pitman-Yor process as the sum of a jump part and a continuous one. We sample the continuous part in two ways, leading to two variants of the proposed sampler. We also propose a threshold to improve mixing in the first variant of our algorithm.

The two variants of our sampler are compared with a marginal method, that is the celebrated Algorithm 8 of [Nea00], and two conditional algorithms based on the stick-breaking representation, namely the efficient slice sampler of [KGW11] and the truncated blocked Gibbs sampler of [IJ01]. We also investigate effects of removing the proposed threshold in the first variant of our algorithm and introducing the threshold in the efficient slice sampler of [KGW11]. Results on real and simulated data sets illustrate that our algorithms outperform the other conditionals in terms of mixing properties.

Introduction

Bayesian nonparametrics have recently gained popularity in a great number of applications in statistics and machine learning (density estimation, clustering, image segmentation and reconstruction, language modeling etc.). Dirichlet process mixture models (DPM) ( [START_REF] Ferguson | Bayesian density estimation by mixtures of Normal distributions[END_REF], [START_REF] Lo | On a class of Bayesian Nonparametric estimates: I. density estimates[END_REF]) have become ubiquitous in Bayesian nonparametric modeling as reviewed by [START_REF] Müller | Nonparametric Bayesian data analysis[END_REF] and recently by [START_REF] Müller | Bayesian Nonparametric Inference-Why and How[END_REF]. This makes crucial the use of effective sampling strategies for DPM and their two-parameter generalization a.k.a as Pitman-Yor mixtures (PYM). In particular, the mixing properties of MCMC samplers appears as a key point in order to address high dimension applications using large datasets. In this paper, we investigate such a task.

A DPM assumes that the random density function can be written as

f (x) = p(x|θ)dH(θ) with H ∼ DP(α, G 0 ), (1)
where {p(•|θ) : θ ∈ Θ} is a family of non-negative (possibly multivariate) kernels defined on a complete and separable metric space X such that X p(x|θ)λ(dx) = 1 for all θ ∈ Θ and for some σ-finite measure λ. The prior over the mixing distribution H is a Dirichlet process (DP) [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] with parameters α > 0 and base distribution G 0 that is the prior guess at the shape of

H, E[H(•)] = G 0 (•).
As it is well-known, the DP selects an almost surely discrete random probability measure that can be represented as

H(•) = ∞ k=1 w k δ θ * k (•), (2)
where (w k ) ∞ k=1 are non-negative weights that sum to unity, and (θ * k ) ∞ k=1 a sequence of Θ-valued random locations. The random probability measure H can be constructed using the stick-breaking representation, where the weights are

w k = v k l<k (1 -v l )
for a sequence (v k ) k≥1 of independent random variables on (0, 1). One can replace the Dirichlet process in (1) with any almost surely discrete random probability measure, for example a Pitman-Yor process (PYP). The definition of the PYP is set forth in the next paragraph, but here we anticipate that the Pitman-Yor process mixture model is obtained as follows

f (x) = p(x|θ)dH(θ) with H ∼ PY(d, α, G 0 ), (3)
with d ∈ [0, 1), α > -d. The PYM is an interesting alternative to the DPM that allows more flexibility in the modeling. Alternatively, the model (3) can be expressed hierarchically as follows

x i |θ i ind ∼ p(x i |θ i ) θ i |H iid ∼ H (4)
H|α, G 0 ∼ PY(d, α, G 0 ).

By exploiting the discreteness of H, the PYM also provides a flexible model for clustering items of various kinds in a hierarchical setting without explicitly specifying the number of components. In a Bayesian context, one is interested in the posterior distribution of the random density f . However, this later does not exhibit any closed form and inference is necessarily simulation-based. After the first MCMC method introduced by [START_REF] Escobar | Estimating normal means with a Dirichlet process prior[END_REF] for DP mixture models, many authors contributed to its improvement ( [START_REF] Escobar | Bayesian density estimation and ence using mixtures[END_REF], [START_REF] Maceachern | Estimating mixture of Dirichlet process models[END_REF]) and to handle non-conjugate cases [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. These aforementioned techniques integrate with respect to the mixing measure H and, thus, remove the infinite dimensional aspect of the problem. This leads to the so-called marginal (or Pólya urn) approaches. Efficient versions of marginal samplers usually achieve good mixing performances [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. An alternative to marginal methods is given by the conditional algorithms which explicitly represent the mixing measure using, for instance, its stick-breaking series representation ([IJ01], [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF]). The challenge in conditional approaches is to deal with the countably infinite representation of H in Equation (2). In [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF], authors resort to an approximation and truncate the mixing measure at a deterministic value. An alternative which avoids hard truncations was proposed by [START_REF] Muliere | Approximating distributions of random functionals of Ferguson-Dirichlet priors[END_REF] who provided an approximation of the Dirichlet process by means of a random truncation of its stick-breaking representation. The idea of random truncation has also been developed in [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF] with a Metropolis-Hastings sampling scheme, and in [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF] using the slice sampling strategy. This later algorithm has been improved by [START_REF] Papaspiliopoulos | A note on posterior sampling from Dirichlet mixture models[END_REF] and [START_REF] Kalli | Slice sampling mixture models[END_REF].

The key advantage of conditional methods using the stick-breaking construction is in updating the mixing measure H as well as the other parameters in the model. This makes possible direct inference on H. Furthermore, components weights are explicitly represented and updated. This property makes these algorithms able to update blocks of parameters and easy to parallelize in order to take advantage of recent parallel computation hardware architectures, which is well suited particularly for large data sets. On the flip side, by integrating mixture components out of the model, marginal techniques are based on incremental updates which are prejudicial when working with huge data sets. Another drawback of marginalizing over the mixing measure is that computing posterior conditionals require additional sampling steps (see [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF]).

However, dealing with exchangeable prediction rules, marginal methods exhibit most of the time better mixing properties than conditionals, and our experimental comparison in Section 5 corroborates this assessment. Also, the random weights are collapsed by marginalization and this results in a crucial reduction of the parameters space dimension. A limitation of conditional methods based on the stick-breaking representation is that the sampler operates in the space of non-exchangeable cluster labels, as pointed out in the paper by [START_REF] Porteous | Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation[END_REF]. Indeed, in this representation, weights are explicitly defined by the prior and components are represented with a size-biased ordering over their labels. This means that components with lower labels have higher prior probabilities than components with higher labels. As a consequence, components are not interchangeable and cluster prior labelling contributes to the posterior sampling. In this situation, the sampler needs to mix efficiently over clusters labels to avoid any clustering bias. Then, [START_REF] Porteous | Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation[END_REF] recommend systematic use of two additional Metropolis-Hastings moves ("label-swap" and "label-permute") in order to improve mixing over clusters. When working with non-exchangeable clusters labels, this additional step seems to be the only way to improve the mixing over clusters (see also [START_REF] Papaspiliopoulos | A note on posterior sampling from Dirichlet mixture models[END_REF]). In contrast, in marginal methods using Pólya urn representation, the sampling occurs in the space of equivalence classes over exchangeable clusters labels where clusters identities are arbitrary and insignificant. This is the adequate space to live for the sampler because cluster labels are irrelevant.

In this context, we introduce a new conditional sampling scheme which is formulated in the space of equivalence classes over clusters labels where clusters identities are irrelevant. Our approach can be seen as marginalizing the clusters labels ordering. Before going further, we point out that instead of using the stick-breaking representation for the underlying mixing measure, recent conditional samplers exploit other constructive representations. For instance, the use of the so-called Ferguson and Klass representation [START_REF] Ferguson | A representation of independent increment processes without Gaussian components[END_REF] of independent increment processes has been recently considered in the literature. See, e.g., [START_REF] Griffin | Posterior simulation of Normalized Random Measure mixtures[END_REF] and [START_REF] Nieto-Barajas | A sensitivity analysis for Bayesian nonparametric density estimators[END_REF], and references therein for some recent contributions in this direction. Such approaches are interesting in the sense that they allow to consider classes of priors, in general, wider than the PYP and the DP. However, they become non-trivial to implement, even when applied to the DP model. Given the importance of the DPM and PYM, the dominant priors in Bayesian nonparametrics, it seems important to devote attention to the development of alternative, simple and efficient algorithms. It is the purpose of this paper to provide a simple and fast way to infer DPM and PYM. Also, it is worth mentioning that since this work was presented in a technical report, it has been successfully applied, for example in [CTM].

The remainder of this paper is structured as follows. In Section 2, we recall some preliminaries about Pitman-Yor processes. Section 3 briefly recalls the basis of the algorithms that are used to compare our samplers and discuss the advantages and limitations of marginal and conditional methods. Afterwards, in Section 4, we present the two variants of the proposed MCMC. We evaluate performance of the algorithms through application to real and simulated data sets in Section 5. Finally, we conclude the paper in Section 6 with discussions and extensions for further work.

Preliminaries on Pitman-Yor processes

In this section, we provide a succinct description of the Pitman-Yor process. We refer interested readers to [START_REF] Pitman | Some developments of the Blackwell-MacQueen urn scheme[END_REF], [START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF], [START_REF] Pitman | Combinatorial stochastic processes[END_REF], [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF], [START_REF] Ishwaran | Some further developments for stickbreaking priors: Finite and infinite clustering and classification[END_REF] and references therein for more details on Pitman-Yor processes.

Definition 1 (Two-parameter GEM distribution) Let d ∈ [0, 1), α > -d, (v j ) j≥1 a sequence of independent random variables such that for all j, v j ∼ Beta(1 -d, α + jd). Define the sequence of weights (w j ) by the stick-breaking scheme as follows

w 1 = v 1 , w 2 = v 2 (1 -v 1 ), . . . , w j = v j j-1 i=1 (1 -v i ).
The sequence w = (w j ) j≥1 is said to follow a two-parameter GEM distribution (after Griffiths, Engen and McCloskey), with parameters d and α, and denoted by w ∼ GEM(d, α).

Definition 2 (Pitman-Yor process) Let w ∼ GEM(d, α) and G 0 a diffuse (non-atomic) probability measure on a measurable space (Θ, B) (i.e G 0 ({θ}) = 0 for each θ ∈ Θ). Let Θ * = (θ * j ) j≥1 be iid G 0 , independently of w. Then

H(•) d = ∞ k=1 w k δ θ * k (•), (5)
where δ θ * (•) denotes the Dirac measure giving mass 1 at θ * , is distributed according to a two-parameter Poisson-Dirichlet distribution on (Θ, B). We refer to H as a Pitman-Yor process with parameters d and α and base measure G 0 .

It is denoted by: H ∼ PY(d, α, G 0 ).

Setting d = 0, the Pitman-Yor process reduces to the Dirichlet process with parameters α and G 0 , while the PY(d, 0, G 0 ) yields a measure whose random weights are based on a stable law with index 0 < d < 1. Definition 2 is a constructive definition of the PYP called stick-breaking representation. The PYP has also a characterization in terms of generalized Blackwell-MacQueen [BM73] urn scheme. Henceforth, K n will denote the random variable identifying the number of distinct values appearing in the sample (θ 1 . . . , θ n ) from a two-parameter Poisson-Dirichlet process, and k n a realization of K n .

Proposition 1 (Generalized Blackwell-MacQueen urn) Let d ∈ [0, 1), α > -d and G 0 a diffuse probability measure on Θ. Consider a sequence of (θ i ) i≥1 generated via the following predictive distributions

θ 1 ∼ G 0 (6) θ n+1 |θ 1 , . . . , θ n ∼ α + dk n α + n G 0 + kn j=1 n j -d α + n δ θ * j , (7)
where {θ * j , j = 1, . . . , k n } are the unique values among {θ i , i = 1, . . . , n} and n j the frequency of θ * j . The distribution of this sequence of exchangeable draws converges almost surely to a discrete distribution which is distributed according to a PY(d, α, G 0 ) when n goes to infinity: θ 1 , . . . , θ n |H ∼ H where H ∼ PY(d, α, G 0 ).

The predictive distributions (6)-(7) are the key components of marginal methods described in Section 3. They also make clear the role of the parameters in Pitman-Yor and Dirichlet process mixture models. In a DPM, the probability for θ n+1 to coincide with an already observed value, say θ * j , is determined by n j while the probability to sample a new value for θ n+1 from the base measure G 0 is proportional to the concentration parameter α. As a consequence, α is the only parameter which can be used to tune the number of clusters. [START_REF] Korwar | Contributions to the theory of Dirichlet processes[END_REF] showed that the asymptotic behaviour of the number of clusters K n that are induced by the Dirichlet process is

K n / log(n) → α.
For some applications, this feature can be too restrictive. In a PYM, together α and d control the formation of new clusters. The first controls the overall number of clusters whereas d determines the asymptotic growth of K n . The larger d, the more new values are generated. In [START_REF] Pitman | Combinatorial stochastic processes[END_REF], it is shown that

K n /n d → S d
where S d is a strictly positive random variable, with continuous density

P d,α (S ∈ ds) = g d,α (s) := Γ(α + 1) Γ( α d + 1) s α d g d (s)ds (s > 0)
where

g d = g d,0
is the density function of a positive stable random variable with parameter d. These effects can be visualized by looking at the induced prior distribution of the number of distinct clusters in a sample of size running from 1 to 10, 000, for both DP and PYP with respectively α = 10 and (d = 0.5, α = 10). Results are displayed in Figure 1. Another characteristic of the DP and the PYP is a reinforcement mechanism (named after "rich-get-richer") that tends to reinforce among the observed clusters those having higher frequencies. As noted above, in a DP the probability of joining an existing cluster j is proportional to the size n j of that cluster. In a PYP however, this probability is proportional to n j -d, with d acting as a discount parameter that reduces the probability of adding a new observation to an existing cluster. This yields a power-law behaviour for the PYP which greatly influences the clustering structure and makes this prior well-suited to natural language processing applications, unlike the DP [START_REF] Teh | A hierarchical Bayesian language model based on Pitmanyor processes[END_REF]. We refer to [DBFL + ] for more details on the role of d in the combined effect of the reinforcement mechanism and the increase in the rate at which new values are generated.

After having recalled the basis of Pitman-Yor processes, we present in the next section the MCMC samplers of both marginal and conditional type for the Dirichlet and Pitman-Yor mixture models.

Sampling from DPM and PYM

To be self-contained, in this section we recall basis of the two classes of MCMC algorithms for sampling under the posterior of a DPM and a PYM and briefly describe each of the algorithms that are used to compare our proposed method.

MCMC algorithms for PYM

Posterior distributions are intractable in DPM and PYM. Posterior inference is performed using approximation techniques such MCMC methods. There are many sampling MCMC algorithms which can roughly be divided into two categories: marginals and conditionals.

Marginal methods

These methods are called marginal since the infinite dimensional random component, namely the mixing measure H, is integrated out of the model and the predictive distributions are used within a Gibbs sampler to get posterior samples. Marginal methods can be sub-categorized into conjugate or non-conjugate models. By conjugacy, we mean that the mixture kernel p(•|θ) and the base distribution G 0 form a conjugate pair. In this case, calculations in conditional posterior distributions are simplified and can be performed analytically ( [START_REF] Neal | Bayesian mixture modeling[END_REF], [START_REF] Escobar | Estimating normal means with a Dirichlet process prior[END_REF], [START_REF] West | Hierarchical priors and mixture models, with application in regression and density estimation[END_REF], [START_REF] Escobar | Bayesian density estimation and ence using mixtures[END_REF] and [START_REF] Bush | A semiparametric Bayesian model for randomised block designs[END_REF]). In non-conjugate models however, posteriors can not be easily calculated. The sampling scheme is more difficult and requires elaborated techniques ([MM98], [START_REF] Walker | Sampling methods for Bayesian nonparametric inference involving stochastic processes[END_REF] and [START_REF] Green | Modelling heterogeneity with and without the Dirichlet process[END_REF]). The reader is referred to [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] for a more complete overview and discussions about these methods. Neal [Nea00] also proposes two novel sampling schemes for nonconjugate models: the first (referred to as "algorithm 7" in Neal's paper) uses a combination of Metropolis-Hastings steps with Gibbs updates. The second, named after "algorithm 8", is based on an augmentation scheme and extends the model to include auxiliary components which exist temporarily. We briefly detail this algorithm we will use to contrast our sampler since, to our knowledge, it achieves the best mixing properties in marginal methods. This algorithm has been developed for Dirichlet process mixture models in Neal's paper. Here, we slightly modify it by adding the second parameter in order to infer Pitman-Yor mixture models.

Algorithm 8 of [Nea00]:

The idea behind "algorithm 8" of [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] is to add auxiliary components (representing potential future components) when updating classification variables.

To be more precise, let c i such that

c i = k iff θ i = θ * k , where Θ * = (θ * 1 , θ * 2 , . . . ) denote the unique values among (θ 1 , . . . , θ n ).
Since data are exchangeable and labels of components completely arbitrary, each datum x i can be treated as the last. Using the predictive distribution (7), x i is thus assigned to an already represented component or to an auxiliary component. If we denote by k - n the number of active components disregarding observation i and n -i,k the number of c j for j = i that are equal to k, the prior probability to allocate x i to an active component is

n -i,k -d
α+n-1 and the probability to create a new component is α+dk - n α+n-1 , which will be equally distributed among the m auxiliary components.

The choice of m is left to the user. It is governed by a balance between computational considerations and mixing properties. Combining with the likelihood, updating classification variables is done via the conditional probabilities:

Pr(c i = k|x i , c -i , Θ * ) ∝        n -i,k -d n -1 + α p(x i |θ * k ) for k = 1, . . . , k - n , α + dk - n m(n -1 + α) p(x i |θ * k ) for k = k - n + 1, . . . , k - n + m.
After this step, parameters for non-empty components are updated according to their posterior law based on the prior G 0 and the likelihood of all data currently allocated to:

p(θ * k |X, c) ∝ G 0 (dθ * k ) {i:ci=k} p(x i |θ * k ).

Conditional methods

In contrast to marginal methods, conditional samplers retain the mixing measure H and explicitly represent it using, for example, the stick-breaking construction stated in Equation ( 5). The issue is to treat the infinite dimensionality of H. There are methods that approximate H by a deterministic truncation of the number of its components, and those that use a finite but random number of masses.

The truncated blocked Gibbs sampler [IJ01]

To remove the infinite dimensional aspect of the problem, Ishwaran and James ([IJ01], [START_REF] Ishwaran | Some further developments for stickbreaking priors: Finite and infinite clustering and classification[END_REF]) truncate H at a chosen integer value N . It is necessary to set v N = 1 in the stick-breaking construction of H to ensure that the truncated measure H N is a probability measure, with distribution:

H N (•) = N k=1 w k δ θ * k (•).
Under this truncated framework, and using classification variables c = (c 1 , . . . , c n ), the hierarchical model (4) can be rewritten as follows

x i |c i , Θ * ∼ p(x i |θ * ci ) c i |w ∼ N k=1 w k δ k (•) (8) w|d, α ∼ GEM(d, α) θ * k |G 0 ∼ G 0 ,
where X = (x 1 , . . . , x n ) represent the data, Θ * = (θ * 1 , . . . , θ * N ) the parameters of the mixture components and w = (w 1 , . . . , w N ) the components weights, with

w N = N -1 k=1 (1 -v k ) = 1 - N -1 k=1 w k .
Rewriting the model in this form makes direct posterior inference possible since one has to treat a finite number of components. If we denote by k n the number of currently non-empty components, the full conditionals involved in the Gibbs sampler are given by the following:

1. Conditional for c: c i |w, Θ * , X ∼ N k=1 w k,i δ k for i = 1, . . . , n where (w 1,i , . . . , w N,i ) ∝ (w 1 p(x i |θ * 1 ), . . . , w N p(x i |θ * N )) .
2. Conditional for w :

w 1 = v * 1 and w k = v * k k-1 l=1 (1 -v * l ) for k = 2, . . . , N - 1, with v * k ∼ Beta(1 -d + n k , α + kd + N l=k+1 n l ) for k = 1, . . . , N -1, where n k = #{i : c i = k}. 3. Conditional for Θ * : • θ * k ∼ G 0 for k = k n + 1, • • • , N . • θ * k |c, X has density proportional to G 0 (dθ * k ) {i:ci=k} p(x i |θ * k ) for k = 1, . . . , k n .
This sampler is easy to implement since the truncation allows it to be similar to standard Gibbs samplers in finite dimensional models. However, even if methods for controlling the truncation accuracy have been proposed ([IJ01], [START_REF] Ishwaran | Some further developments for stickbreaking priors: Finite and infinite clustering and classification[END_REF]), it would be better to avoid any hard approximation. To this purpose, algorithms that only use a finite number of elements at any iteration while allowing inference to the true infinite-dimensional prior have been proposed by Papaspiliopoulos and Roberts [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF] and Walker [Wal07]. This later uses an elegant strategy called slice sampling and is based on auxiliary variables. Walker's slice sampling was improved by Papaspiliopoulos [START_REF] Papaspiliopoulos | A note on posterior sampling from Dirichlet mixture models[END_REF] and Kalli et al. [START_REF] Kalli | Slice sampling mixture models[END_REF]. This later is named after "slice efficient".

The slice sampler ([Wal07], [KGW11])

The idea under the slice sampler proposed by [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF] for inference in DPM is to introduce uniform auxiliary variables which make the mixture model conditionally finite. To make it precise, let us consider the model (4) with kernel p and mixing measure H being constructed using the stick-breaking representation (5). The density of a single observation x i , given w and Θ * , is

f (x i ) = ∞ k=1 w k p (x i |θ * k ) . (9)
The infinite dimensional aspect of (9) is tackled by introducing u = (u 1 , u 2 , . . . , u n ) uniform auxiliary variables such that the joint density of any (x i , u i ) is

f (x i , u i ) = ∞ k=1 1(u i < w k ) p (x i |θ * k ) = ∞ k=1 w k p (x i |θ * k ) U (u i |0, w k ) ,
where

U (•|a, b) denotes the uniform density function over [a, b]. The conditional density of x i , given u i , is f (x i |u i ) = 1 ∞ k=1 1(u i < w k ) k∈{j:wj >ui} p (x i |θ * k ) .
So, given u i , the number of components of the mixture model f is finite. One can then complete the model by introducing a further assignment variable c i and considering the joint density function

f (x i , u i , c i ) = 1(u i < w ci ) p x i |θ * ci = w ci U (u i |0, w ci ) p x i |θ * ci .
The full conditionals required to implement this Gibbs sampler are those of: the slice variables u i , the indicators c i , the components parameters θ * k and the stick-breaking weights v k . In the slice sampler of Walker [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF], each of these variables is sampled independently and the sampling of the v k is quite hard. This was handled in the efficient version of the slice sampler ("slice efficient") proposed in [START_REF] Kalli | Slice sampling mixture models[END_REF]. In this version, the stick weights and the slice variables are blocked during iterations, which, by integrating out slice variables, dramatically simplifies generation of the weights and results in a more efficient sampler compared to Walker's algorithm. So, we will compare our methods with the "slice efficient" one. The required full conditionals of this later algorithm are now given.

Conditional for u

i : u i ∼ U(u i |0, w ci ).

Conditionals for θ *

k and w k : they are the same that in the truncated blocked Gibbs sampler of [IJ01] previously described.

Conditional for c

i : Pr(c i = k) ∝ 1(k : w k > u i ) p(x i |θ * k ).
To sample from this probability mass function, one needs to know the exact number of components that are required at each iteration of the sampler. It is given by the smallest K such that

K k=1 w k > 1 -u * ,
where u * = min{u 1 , . . . , u n }.

In this section, we have outlined a marginal and two conditional methods for inference in Pitman-Yor mixture models. In the next section, the merits and limitations of each class of algorithms are summarized. This motivates the development of a new sampling approach which falls within the class of conditional approaches.

A new sampling method

By integrating mixture components out of the model, marginal algorithms make the allocation step very sequential since they need to condition on all previously allocated data. These incremental updates make marginal samplers not easily parallelizable. This feature is prejudicial when working with large data sets since the allocation step forms the most time consuming part of the algorithm. However, marginal samplers have the advantage that they deal with exchangeable prediction rules and exhibit most of the time better mixing properties.

The stick-breaking representation of the Pitman-Yor process allows to explicitly represent the weights in terms of independent Beta random variables. As a consequence, updating indicator variables in the allocation step is done without conditioning on the other indicators. This makes conditional algorithms using this representation easy to implement in a parallel computer. However this simplicity in the representation comes at a cost of slower mixing, since the stick-breaking prior has a weak preference for components to be sorted by decreasing mass. Consequently, additional moves in the MCMC algorithm are necessary to improve mixing over clusters labels as suggested in [START_REF] Porteous | Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation[END_REF] and [START_REF] Papaspiliopoulos | A note on posterior sampling from Dirichlet mixture models[END_REF].

In this context, we propose a conditional algorithm which is rather different from the others discussed so far. Our sampler is an attempt to combine the main advantages of marginal and conditional algorithms. The underlying idea is to integrate out the explicit order of clusters labels like in marginal methods hence collapsing the model to a lower dimensional space while keeping components weights as done in conditional approaches. To these aims, we propose to replace the standard posterior updating of the mixing measure based on the stick-breaking representation, with a posterior updating of Pitman

w j δ θ * j + r kn H kn , (10)
where

(w 1 , . . . , w kn , r kn ) ∼ Dir(n 1 -d, . . . , n kn -d, α + dk n ) H kn ∼ PY(d, α + dk n , G 0 ),
and H kn independent of (w 1 , . . . , w kn , r kn ), with E(H kn ) = G 0 .

The posterior characterization (10) allows us to work on the space of equivalence classes of clusters θ * j and, due to exchangeability, to integrate out the order of cluster labels as in the marginal samplers. Indeed, Pitman showed in [START_REF] Pitman | Random discrete distributions invariant under size-biased permutation[END_REF] equivalence between exchangeability of the random partition generated by sampling from a discrete distribution and symmetry in the law characterizing the limiting frequencies of occupied components given the data. We can easily check that exchangeability is ensured in equation (10) since it sums to a symmetric Dirichlet distribution and an unconditional Pitman-Yor process (independent of the observed data). So, our sampler lives in the space of equivalence classes over clusters labels. These labels are then exchangeable and no mix over them is needed. This property has important consequences on the algorithm mixing. As opposed, in the conditional algorithms using the stick-breaking representation, exchangeability is lost when using the usual updating rule:

H(•)|θ 1 , . . . θ n = k∈c * w * k δ θ * k (•) + k / ∈c * w k δ Z k (•), (11)
where c * = (c * 1 , . . . , c * kn ) are the unique values of the classification variables c = (c 1 , . . . , c n ), the weights w * k follow a GEM distribution with updated parameters:

w * 1 = v * 1 , w * 2 = v * 2 (1 -v * 1 ), . . . , w * n = v * n n-1 i=1 (1 -v * i ) where v * l ∼ Beta(1 -d + n l , α + ld + ∞ m=l+1 n m ), and for all k / ∈ c * , Z k iid ∼ G 0 .
This clearly illustrates that the posterior distribution of a random probability measure constructed via the stick-breaking representation depends on which explicit atoms labels observations are allocated to. This property is not necessary and has the impact of bothering the Gibbs sampler.

Given the posterior characterization (10), we are now in position to set up a new Gibbs sampling scheme for simulating from the posterior of a Pitman-Yor mixture model.

Proposed variants

For simulating H kn , the continuous part of the posterior given in equation ( 10), we propose two variants. The first makes use of a thresholded version of the "slice efficient dependent" of [START_REF] Kalli | Slice sampling mixture models[END_REF]. The second is based on a truncation as originally suggested in [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF].

Exchangeable Thresholded Slice Sampler

We first propose a slice sampler inspired from [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF] and [START_REF] Kalli | Slice sampling mixture models[END_REF]. The main steps are now summarized.

We augment the state with additional slice variables u = (u 1 , u 2 , . . . , u n ) such that the joint density for any (x i , u i ), given a collection w of random masses and component parameters Θ * , is

f (x i , u i ) = ∞ k=1 w k p (x i |θ * k ) U (u i |0, ξ k ) , (12)
where ξ k is a dependent variable such that for all k,

ξ k = min (w k , ζ) , (13)
with ζ ∈ ]0, 1] and is independent of w k . Here, ζ is a threshold that we propose in order to improve mixing properties of the sampler compared to [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF] and [START_REF] Kalli | Slice sampling mixture models[END_REF]. The threshold ζ can be a random or deterministic variable. Here, the role of ζ is to ensure that on average at each iteration, all occupied clusters and at least a non-occupied one are proposed by the algorithm. For example, a deterministic typical value of ζ that gives rise to good trade-off between mixing properties and computational burden is the mean weight of the first atom (in size-biased order of H kn ) with no data allocated to, which can be expressed in the two-parameter case as

ζ = (α + d E α,d (K n ))(1 -d) (α + n)(α + 1) , where E α,d (K n ) is the expected value of K n : E α,d (K n ) = n i=1 (α + d) i-1↑ (α + 1) i-1↑ ,
where (x) a↑ = Γ(x + a)/Γ(x) is the Pochhammer symbol.

In the case of the Dirichlet process (d = 0),

E α,d (K n ) = n i=1 α α + i -1 = α log 1 + n α .
For d = 0, it can be easily checked that,

E α,d (K n ) = α d (α + d) n↑ (α) n↑ -1 .
For sufficiently large n, this expectation can be fairly approximated using Stirling's formula

E α,d (K n ) ≈ Γ(α + 1) d Γ(α + d) n d .
Coming back to the slice sampling formulation, using equation ( 13), we can rewrite equation (12) as follows:

f (x i , u i ) = 1 (ζ > u i ) ζ -1 w k >ζ w k p (x i |θ * k ) + w k ≤ζ 1 (w k > u i ) p (x i |θ * k ) ,
where both sums are finite since #{j : w j > ε} < ∞, for all ε > 0. The use of u allows to sample a finite number K * of weights and locations for H kn .

Let us here denote w = (w 1 , w 2 , . . . , w kn , w kn ) where w 1 , w 2 , . . . , w kn are the k n Dirichlet random weights in the posterior characterization (10), and w kn is a collection of random variables distributed according to a two-parameter GEM(d, α + dK n ) distribution; these are the stick-breaking random weights of H kn . The Gibbs sampler allows to generate variables from the joint posterior of (Θ * , c, w, u|X), by sampling iteratively from each full conditional. As in [START_REF] Kalli | Slice sampling mixture models[END_REF], we jointly sample w, u|c. The full conditional distributions involved in the steps of the sampler are then:

• p(c|θ * , w, u),

• p(θ * |c, w, u),

• p(w, u|c, θ * ) = p(u|w, c, θ * ) p(w|c, θ * ).

We now provide a way of simulating from each conditional.

Conditional for (w, u) :

We jointly sample w, u|c in three steps by first sampling w 1 , w 2 , . . . , w kn |c, then u|w 1 , w 2 , . . . , w kn , c, and finally w Kn |u. The mains steps are now given.

• Sample w k for k ≤ k n :

w 1 , . . . , w Kn ,r kn |c ∼ Dir (n 1 -d, . . . , n kn -d, α + k n d) .
• Sample u i |w 1 , w 2 , . . . , w kn , c:

u i |w 1 , w 2 , . . . , w kn , c ind. ∼ U (u i |0, min (w ci , ζ)) .
Set u * = min{u 1 , . . . , u n }.

• Sample w k for k > k n . While r k-1 > u * , v k ∼ Beta (1 -d, α + k d) , w k = v k r k-1 , r k = r k-1 (1 -v k ) .
Set K * = min ({k : r k < u * }).

Clearly, w k < u * for all k > K * , that is why we only have to sample a finite set of w K * .

Note that, at each iteration, non-empty clusters are re-labeled according to their order of appearance in the sampling. We operate in the space of equivalence classes over non-empty clusters labels which are thus exchangeable. The stick-breaking prior only concerns empty clusters for the given iteration of the Gibbs sampler. As pointed out, this encourages good mixing over clusters.

Conditional for c:

As underlined, sampling of classification variables requires the computation of a normalizing constant which becomes feasible using auxiliary variables since the choice of c i is from a finite set:

c i |w, u, Θ * , X ind ∼ K * k=1 w k,i δ k (•) , where w k,i ∝ 1 (w k > u i ) max (w k , ζ) p (x i |θ * k ) , and
K *
j=1 w k,i = 1. Note also that, in order to speed up computations, it is convenient to sort weights w k , k > k n in decreasing order. By this, we can avoid tests for all k > κ as soon as w κ < u i .

Conditional for Θ * :

• Updating parameters for non-empty components from the density proportional to:

G 0 (dθ * k ) i:ci=k p(x i |θ * k ) for all k ≤ k n .
• Sampling parameters for unallocated components from their priors:

θ * k iid ∼ G 0 , for k n < k ≤ K * .
The blocked Gibbs sampler structure allows easy implementation of the algorithm on a parallel computer.

Exchangeable Truncated Gibbs Sampler

The second variant of the algorithm we propose is an alternative of the first one.

It is still based on the posterior given in equation (10). But instead of using the slice sampling strategy to sample the continuous part H kn , we resort to an approximation by taking a fixed level L. This truncation eliminates the need of auxiliary variables. This scheme was suggested in [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF]. We approximate equation ( 10) by kn j=1

w j δ θ * j + r kn H * kn ,
where H * kn is an approximation of H kn , i.e a truncation of H kn at level L. The total number of represented components is then K * = k n + L. The main steps are now given.

• Sample classification variables:

(c i |w, u, Θ * , X) ind ∼ K k=1 w k,i δ k (•) , where w k,i ∝ w k p (x i |θ * k ) and K * k=1 w k,i = 1. • Sample w k for k ≤ k n : (w 1 , w 2 , . . . , w kn , r kn |c) ∼ Dir (n 1 -d, n 2 -d, . . . , n kn -d, α + k n d) . • Sample w k for k n < k ≤ K * : v k ∼ Beta (1 -d, α + k d) , w k = v k r k-1 , r k = r k-1 (1 -v k ) . Set w K * = r K * -1 such that v K * = 1.
• Sample components parameters using the density proportional to

G 0 (dθ * k ) i:ci=k p(x i |θ k ) for non-empty components (i.e. k ≤ k n ),
the priors for unallocated components:

θ * k iid ∼ G 0 , for k n < k ≤ K * .

Comparisons of algorithms

In this section, we carry out a comparative study that involves a variety of data sets, both real and simulated. We evaluate the performance of the samplers described in the previous sections and our new sampling method. We thus compare these following algorithms:

• algorithm 8 of [Nea00] ("Algo. 8"),

• the slice efficient of [START_REF] Kalli | Slice sampling mixture models[END_REF] ("Slice efficient"),

• the truncated blocked Gibbs sampler of [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF] ("Trunc."),

• the two variants of the proposed sampling scheme based on an exchangeable model ("Slice exch. thres." and "Trunc. exch.").

We also investigate the gain in the mixing performances of the algorithms due to the exchangeability property of the model on one hand, and to the proposed threshold on the other hand. For this reason, we implement in addition our slice sampler using the exchangeable model but without the threshold ("Slice exch. without thres.") and the slice efficient of [START_REF] Kalli | Slice sampling mixture models[END_REF] which uses an nonexchangeable model with the introduction of the threshold ("Slice eff. thres."). Note that the "Slice efficient" is referred to as "Slice efficient dependent" in [START_REF] Kalli | Slice sampling mixture models[END_REF], in contrast to their independent version which makes use of a deterministic slice function.

Data specification:

We tested the algorithms with p(•|θ) being a normal kernel with parameters θ * = (µ, σ 2 ) and G 0 a normal-inverse Gamma distribution i.e, G 0 (µ, σ -2 ) = N (µ|η, κ 2 ) ×G(σ -2 |γ, β) where G(•|γ, β) denotes the Gamma distribution with density proportional to x γ-1 e -x/β .

For comparison purposes, we considered the same real and simulated data sets as in [START_REF] Kalli | Slice sampling mixture models[END_REF].

1. The simulated data were generated from the following mixtures of Gaussians.

• A bimodal mixture (bimod): 0.5 N (-1, 0.5 2 ) + 0.5 N (1, 0.5 2 ).

• An unimodal lepto-kurtic mixture (lepto): 0.67 N (0, 1) + 0.33 N (0.3, 0.25 2 ).

These simulated densities are shown in Fig. 2.

In order to gauge algorithms performance for small and large data sets, we generated n = 100, n = 1, 000 and n = 10, 000 draws from each of these two mixtures.

The real data are

• Galaxy data, which are the velocities (in 10 3 km/s) of 82 distant galaxies diverging from our own. It is a popular data set in density estimation problems and is also used by [START_REF] Escobar | Bayesian density estimation and ence using mixtures[END_REF], [START_REF] Green | Modelling heterogeneity with and without the Dirichlet process[END_REF] for instance.

• S&P: this consists of 2023 daily index returns. This data set is unimodal, asymmetric, and heavy-tailed. We monitored the convergence of two quantities: the deviance of the estimated density and the number of occupied clusters. The deviance is a global function of all parameters of the model and is defined as

D = -2 n i=1 log   j n j n p(x i |θ * j )   ,
where n j is the size of cluster j.

The performance of competing samplers in their stationary regime was judged by looking at the integrated autocorrelation time (IAT) for each monitored quantity. IAT is defined in [START_REF] Sokal | Monte carlo methods in statistical mechanics: Foundations and new algorithms[END_REF] as,

τ = 1 + 2 ∞ j=1 ρ j ,
where ρ j is the sample autocorrelation at lag j. This quantity is an indicator of mixing behaviour of algorithms and measures effectiveness of MCMC samples. As such, it has also been used by other authors to compare MCMC methods (for example [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF], [START_REF] Green | Modelling heterogeneity with and without the Dirichlet process[END_REF], [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF], [START_REF] Kalli | Slice sampling mixture models[END_REF]). IAT controls the statistical error in Monte Carlo measurements. In fact, the correlated samples generated by a Markov chain at equilibrium cause a variance that is 2τ larger than in independent sampling [START_REF] Sokal | Monte carlo methods in statistical mechanics: Foundations and new algorithms[END_REF]. If we denote by τ j the integrated autocorrelation time produced by algorithm j for a given quantity, then τ 1 /τ 2 = k > 1 means that algorithm 1 requires k more iterations than algorithm 2 to produce the same Monte Carlo error [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF]. So, when comparing two alternative Monte Carlo algorithms for the same problem, the most efficient is the one that produces the smaller IAT since it provides better estimates.

However, the calculation of IAT is difficult in practice. Following [START_REF] Sokal | Monte carlo methods in statistical mechanics: Foundations and new algorithms[END_REF], an estimator of τ can be obtained by summing the estimated autocorrelations up to a fixed lag L:

τ = 1 + 2 L j=1 ρj .
The choice of the cut-off point L is left to the user.

One can also estimate the standard error of τ using this formula from [START_REF] Sokal | Monte carlo methods in statistical mechanics: Foundations and new algorithms[END_REF]:

std(τ ) ≈ 2(2L + 1) M τ 2
where M is the Monte-Carlo size.

Algorithms parametrization:

At first, we set the discount parameter d of the PYM to zero in order to reduce it to a DPM. The strength parameter of the PYM, that is now the precision parameter of the DPM, was respectively set to α = {1, 5}. Secondly, we investigated the behaviour of the competing algorithms in a power-law case (Pitman-Yor). The values d = 0.3 and α = 1 were chosen for the PYM. The expected number of a priori components are then much larger than in the DPM case. We report in Tables 1-3 the expected number of a priori components, for each data set length and each parametrization.

Data

E(K n ) Galaxy (n = 82) 4.4 Lepto/bimod (n = 100) 4.6 Lepto/bimod (n = 1, 000) 6.9 S&P 500 (n = 2023) 7.6 Lepto/bimod (n = 10, 000) 9.2 The hyperparameters have been fixed in a data-driven way according to [START_REF] Green | Modelling heterogeneity with and without the Dirichlet process[END_REF] and set as follows: if R is the range of the data we take η = R/2 (mid-range), κ 2 = 1/R 2 , γ = 2 and β = 0.02R 2 .

The blocked Gibbs sampler of [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF] was truncated at level N = 3αlog(n), where n is the data size. This induces a truncation error that stands for the L 1 distance between the marginal density of the data under the truncated model and the marginal density under the full model, (see [START_REF] Ishwaran | Gibbs sampling methods for stickbreaking priors[END_REF]). The corresponding truncation errors for the different data sets are reported in the following table. Data Galaxy (n = 82) 7.4139e-04 Lepto/bimod (n = 100) 9.0413e-04 Lepto/bimod (n = 1, 000) 8.2446e-06 S&P 500 (n = 2023) 2.2572e-06 Lepto/bimod (n = 10, 000) 7.5181e-08

We also truncated the second variant of our sampler at level L = 2αlog(n). The algorithm 8 of [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] was tested with m = 2 auxiliary components.

We followed the instructions of [START_REF] Sokal | Monte carlo methods in statistical mechanics: Foundations and new algorithms[END_REF] who recommends running the samplers for a sufficient number of iterations. For each of the data sets, we ran 2, 000, 000 iterations for each algorithm and discarded the first 200, 000 for the burn-in period. We believe that these numbers are sufficient to obtain reliable results.

Results and comments:

We report in Tables 4-9 the results of our comparisons for each set of data in the DPM case with α = 1. The other results are postponed in the appendix. Each table contains respectively, for each algorithm, the estimated IAT for the mean number of clusters and for the deviance, the estimated mean number of clusters and the estimated deviance. Estimated IAT are obtained by integrating autocorrelation values for each monitored quantity up to a fixed lag (L D for deviance and L C for the mean number of clusters). The estimates of standard errors are put inside parentheses.

For visual comparison purposes, the autocorrelation curves are displayed in Figures 4 and5 for Galaxy data. We also show in Figure 3 the histogram of the data and the density estimates when using each algorithm. By looking at the curves of the estimated densities, the values of the estimated deviances and the mean number of clusters, we made sure that all algorithms perform the estimation correctly and then they can be assessed through their mixing performance.

In the overall experiments, it turns out that:

• As expected, algorithm 8 performs better than all conditional algorithms since it works in an unidentifiable allocation structure. Furthermore, integrating out the mixture components speeds up the convergence since the dimensionality of the space is drastically reduced. One can refer to [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF] and [START_REF] Porteous | Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation[END_REF] for more details about why conditional approaches are outperformed by marginals.

• On the other hand, the two variants of our method are superior to all other competitors in conditional algorithms using the stick-breaking representation, thanks to exchangeability in the model and the introduction of the threshold we propose. The "Slice efficient" gives the worst performance.

We believe that the poor-mixing due to non-exchangeability in the posterior stick-breaking representation is emphasized by the lack of the weights in slice samplers. This could often hinder the Gibbs sampler in the allocation step, for changing an observation from a component associated with a few observations to a component associated with many. Introducing our threshold would facilitate this change. To validate this conjecture, we experimented the effect of the threshold in the "Slice efficient". This algorithm is referred as "Slice eff. thres." Furthermore, the threshold makes little difference between the thresholded slice efficient ("Slice eff. thres.) and the truncated blocked Gibbs sampler (Trunc.). This later considers the weights of the mixture components when updating classification variables. On the flip side, removing the threshold in our sampler ("Slice exch. without thres.") increases the IAT. Overall, it was observed on all data sets that the threshold causes a rapid decrease of autocorrelation curves in the first lags. However, it slightly increases the computation time per iteration. We underline that all algorithms have been implemented without any parallelization. All of them, excluding "Algo. 8", may be easily parallelized.

We now turn our attention to the benefits we reap thanks to the exchangeability property of the model. This is notable in differences between "Slice exch. thres." and "Slice eff. thres." and in differences between "Trunc. exch." and "Trunc.". We also notice on the curves that the autocorrelations obtained by "Slice exch. without thres." decrease and reach zero faster than in algorithms using non-exchangeable models ("Trunc.", "Slice eff. thres" and "Slice efficient"). This behaviour was observed on all data sets.

It is worth noting that the two variants of our algorithms and algorithm 8 of [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] were stable in all experiments: for various simulations, we always obtained the same results in each data set and in each size of data. On the contrary, the algorithms using non-exchangeable models ([IJ01] and [START_REF] Kalli | Slice sampling mixture models[END_REF]) did not always give the same results. We also observed erratic convergence behaviour of the Gibbs sampler in these two algorithms, particularly for large data sets (for example lepto with n = 10, 000).

Results for d = 0 and α = 1

In the following tables, n stands for the data set length, L D and L C are respectively the number of autocorrelation lags for deviance and for the clusters number. Values inside parentheses correspond to standard deviations of estimates.

IAT

Conclusion and discussion

When models become more and more complex, due to an increase in dimension, the poor mixing of a MCMC algorithm can be inhibiting. Therefore, it seems important to develop samplers that allow to improve mixing while experimenting strategies to reduce the computational cost. The present paper has aimed at providing a simple, efficient and easy to use Gibbs sampler for posterior simulation under Pitman-Yor and Dirichlet mixture models that satisfies the constraint of efficient parallelization ability while maintaining mixing properties closed to Pólya urn approaches. We have attempted to combine blocking properties of conditional approaches which retain the random distribution in the sampling, and exchangeability of the model which is maintained in Pólya urn based algorithms. The proposed approach combines an update but unused formula from [START_REF] Pitman | Some developments of the Blackwell-MacQueen urn scheme[END_REF] with the slice sampling sampling strategy and a tricky threshold. Our comparative study on both real and simulated data sets support our belief that the two novel variants of our conditional Gibbs sampler have the potential to be a useful addition to the menu of samplers for DP and PYP.

A difference between the two proposed variants is that for the truncated version ("Trunc. exch."), the fixed length of approximation has to be decided before effective sampling. This is most of the time not a crux for Dirichlet processes, but for the two-parameter case the fixed approximation may give rise to biased estimates for moderate truncation lengths. For large lengths, the computational burden is emphasized especially for large data sets. The exchangeable thresholded slice version ("Slice exch. thres.") achieves adaptive truncation at each iteration and maintain nice trade-off between IAT and time cost. This latter variant gives then rise to convenient trade-off between IAT and computation time while avoiding any hard truncation.

On one hand, as mentionned our samplers are applicable to Pitman-Yor and Dirichlet mixture models. On the other hand, since the "Slice efficient" of [START_REF] Kalli | Slice sampling mixture models[END_REF] has been developed for more general stick-breaking priors and that the introduction of the proposed threshold improves its mixing property, one can consider such a combination when working with mixtures based on general stick-breaking processes other than Dirichlet and Pitman-Yor processes. In this case, an interesting perspective could be to introduce also mixing moves over clusters labels as suggested in [START_REF] Porteous | Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation[END_REF] and [START_REF] Papaspiliopoulos | Retrospective Markov chain sampling Monte Carlo methods for Dirichlet process hierarchical models[END_REF]. As mentioned, the ordering of clusters labels matters in the stick-breaking representation. A step of labels permutation could result in a better mixing chain.

In our experimental study, it appeared that particularly for the two-parameter class, standard conditional algorithms may present unexpected biased results. This drawback is reinforced for large data sets. On the other hand, Pólya urn based algorithms and our proposed sampling schemes exhibit stable behaviour in all situations.

IAT

Figure 1 :

 1 Figure 1: Number of clusters as a function of the sample size for a DP(α = 10) (left) and a PYP(d = 0.5, α = 10) (right). The rate at which the number of clusters Kn increases is slower in the DP than in the PYP, being respectively log(n) and n d .

Figure 2 :

 2 Figure 2: Bimodal (bimod) and unimodal lepto-kurtic(lepto) mixtures

Figure 3 :

 3 Figure 3: Histogram of data and estimated densities (Galaxy data).

Figure 4 :

 4 Figure 4: Autocorrelation curves used to estimate the IAT for the number of clusters

Figure 5 :

 5 Figure 5: Autocorrelation curves used to estimate the IAT for the deviance (Galaxy data).

  -Yor processes under the class of Poisson-Kingman models introduced by [Pit03]. The next proposition summarizes this posterior characterization. Proposition 2 [Pit96b], Corollary 20 Let H ∼ PY(d, α, G 0 ) where G 0 is a diffuse probability measure s.t. E(H) = G 0 . Consider a sample θ 1 , . . . , θ n |H ∼ H. Let {θ * j } kn j=1 denote the set of unique values of {θ i } n i=1 and n j the number of occurrences of θ * j in the sample. Then the posterior of H can be expressed as follows H|θ 1 , . . . , θ n

d = kn j=1

Table 1 :

 1 E(Kn) for each data set length in a DP(α = 1)

	Data	E(K n )
	Galaxy (n = 82)	14.3
	Lepto/bimod (n = 100)	15.2
	Lepto/bimod (n = 1, 000)	26.5
	S&P 500 (n = 2023)	30
	Lepto/bimod (n = 10, 000)	38

Table 2 :

 2 E(Kn) for each data set length in a DP(α = 5)

	Data	E(K n )
	Galaxy (n = 82)	10.63
	Lepto/bimod (n = 100)	11.48
	Lepto/bimod (n = 1, 000)	25.5
	S&P 500 (n = 2023)	36.4
	Lepto/bimod (n = 10, 000)	58.9

Table 3 :

 3 E(Kn) for each data set length in a PYP(α = 1, d = 0.3)

Table 4 :

 4 Galaxy data n = 82, L D = 150, L M = 300.

		for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	14.48(0.37)	2.88(0.05)	3.986(0.93)	1561.14(21.61)
	Trunc. exch.	14.42(0.37)	2.94(0.05)	3.989(0.93)	1561.16(21.69)
	SE without thres.	35.52(0.92)	4.77(0.09)	3.989(0.93)	1561.15(21.61)
	Truncated	38.65(1.00)	3.63(0.07)	3.996(0.94)	1561.15(21.66)
	Slice efficient	60.65(1.57)	5.28(0.10)	3.991(0.93)	1561.15(21.62)
	Slice eff. thres.	37.82(0.98)	3.61(0.07)	3.986(0.93)	1561.08(22.17)
	Algo 8 (m = 2)	8.25(0.21)	2.57(0.05)	3.987(0.93)	1561.16(21.62)
		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	28.76(0.74)	5.85(0.11)	3.801(1.66)	287.59(8.46)
	Trunc. exch.	28.51(0.74)	6.00(0.11)	3.808(1.67)	287.58(8.48)
	SE without thres.	70.56(1.82)	9.54(0.17)	3.799(1.67)	287.58(8.43)
	Truncated	54.38(1.40)	5.89(0.11)	3.789(1.66)	287.58(8.42)
	Slice efficient	99.92(2.58)	8.76(0.16)	3.784(1.65)	287.58(8.42)
	Slice eff. thres.	55.41(1.43)	5.65(0.10)	3.794(1.66)	287.61(8.58)
	Algo 8 (m = 2)	15.59(0.40)	5.20(0.09)	3.794(1.66)	287.59(8.52)

Table 5 :

 5 Bimod data n = 100, L D = 150, L M = 300.

		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	93.28(3.93)	3.65(0.07)	3.806(1.73)	2735.14(8.66)
	Trunc. exch.	91.20(3.85)	3.69(0.07)	3.795(1.72)	2735.14(8.66)
	SE without thres.	228.64(9.64) 5.44(0.10)	3.809(1.73)	2735.15(8.67)
	Truncated	156.27(6.60) 3.71(0.07)	3.777(1.71)	2735.13(8.62)
	Slice efficient	257.25(10.85) 5.13(0.09)	3.766(1.68)	2735.12(8.61)
	Slice eff. thres.	150.13(6.33) 3.81(0.07)	3.798(1.71)	2735.15(8.72)
	Algo 8 (m = 2)	47.25(1.99)	3.06(0.06)	3.798(1.72)	2735.14(8.65)

Table 6 :

 6 Bimod data n = 1000, L D = 150, L M = 800.

		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	25.61(0.85)	13.53(0.29)	3.991(1.64)	239.74(11.38)
	Trunc. exch.	24.78(0.83)	13.46(0.28)	3.983(1.63)	239.75(11.31)
	SE without thres.	90.56(3.02)	42.74(0.90)	3.991(1.63)	239.73(11.26)
	Truncated	41.22(1.37)	17.03(0.36)	4.001(1.64)	239.72(11.31)
	Slice efficient	120.71(4.03) 46.28(0.98)	3.979(1.64)	239.77(11.29)
	Slice eff. thres.	44.73(1.49)	16.98(0.36)	3.989(1.64)	239.77(11.82)
	Algo 8 (m = 2)	14.79(0.49)	9.83(0.28)	3.994(1.63)	239.72(11.36)

Table 7 :

 7 Lepto data n = 100, L D = 200, L M = 500.

		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	235.49(9.93) 13.17(0.24)	4.006(2.05)	2400.51(18.68)
	Trunc. exch.	237.70(10.02) 13.66(0.25)	4.022(2.08)	2400.48(18.72)
	SE without thres.	462.09(19.49) 18.75(0.34)	3.973(1.99)	2400.53(18.73)
	Truncated	294.24(12.41) 12.67(0.23)	3.958(2.01)	2400.47(18.49)
	Slice efficient	472.95(19.95) 16.91(0.31)	3.864(1.92)	2400.45(18.26)
	Slice eff. thres.	302.50(12.76) 13.80(0.25)	3.978(2.07)	2400.53(18.93)
	Algo 8 (m = 2)	148.81(6.28) 11.55(0.21)	4.018(2.07)	2400.48(18.69)

Table 8 :

 8

Lepto data n = 1000, L D = 150, L M = 800.

(

  Galaxy data).

				Deviance: autocorrelation time
	1
		Slice exch. without thres.: 1561.152 (21.606) iat = 4.770 (0.087) T = 1.94e-03
			Slice exch. thres.: 1561.143 (21.611) iat = 2.885 (0.053) T = 2.74e-03
			Slice efficient: 1561.149 (21.619) iat = 5.283 (0.097) T = 2.27e-03
			Slice eff. thres.: 1561.084 (22.175) iat = 3.614 (0.066) T = 3.03e-03
	0.8			Trunc.: 1561.154 (21.661) iat = 3.627 (0.066) T = 4.33e-03 Algo. 8 (2): 1561.158 (21.616) iat = 2.568 (0.047) T = 1.58e-03
			Trunc. exch.: 1561.163 (21.692) iat = 2.941 (0.054) T = 1.57e-03
	0.6
	0.4
	0.2
	0
	-0.2
	0	20	40	60	80	100	120	140	160

Table 17 :

 17 S&P 500 data n = 2023, L D = 300, L M = 500.

		for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	17.09(0.57)	151.34(3.91) 7.476(2.65)	14989.81(53.22)
	Trunc. exch.	18.02(0.60)	143.45(3.71) 7.484(2.65)	14989.68(51.86)
	SE without thres.	56.57(1.89)	214.64(5.55) 7.473(2.64)	14989.66(52.30)
	Truncated	21.73(0.72)	150.94(3.90) 7.454(2.64)	14990.39(54.49)
	Slice efficient	66.67(2.22)	225.38(5.82) 7.481(2.65)	14990.43(54.78)
	Slice eff. thres.	23.14(0.77)	156.84(4.05) 7.475(2.65)	14990.18(56.20)
	Algo 8 (m = 2)	11.33(0.38)	97.28(2.51)	7.478(2.65)	14989.76(52.48)

Table 20 :

 20 Bimod data n = 1000, L D = 150, L M = 800.

		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	47.85(1.60)	27.00(0.57)	3.719(3.70)	223.92(8.55)
	Trunc. exch.	50.04(1.67)	26.91(0.57)	3.674(3.39)	223.86(8.61)
	SE without thres.	188.12(6.27) 70.11(1.48)	3.709(3.69)	223.94(8.50)
	Truncated	93.37(3.11)	35.96(0.76)	3.955(4.18)	223.78(8.75)
	Slice efficient	224.68(7.49) 74.51(1.57)	3.696(3.68)	223.94(8.49)
	Slice eff. thres.	85.39(2.85)	31.98(0.67)	3.732(3.73)	224.05(9.08)
	Algo 8 (m = 2)	27.28(0.91)	17.83(0.38)	3.720(3.71)	223.92(8.55)

Table 21 :

 21 Lepto data n = 100, L D = 200, L M = 500.

		IAT for # of	IAT for de-	Estimated #	Estimated
		clusters	viance	of clusters	deviance
	Slice exch. thres	156.21(7.37) 13.56(0.29)	4.255(3.22)	2371.35(16.11)
	Trunc. exch.	167.13(7.88) 13.08(0.28)	4.247(3.13)	2371.34(16.03)
	SE without thres.	341.73(16.11) 18.15(0.38)	4.252(3.18)	2371.34(15.97)
	Truncated	270.51(12.75) 17.04(0.36)	4.387(3.62)	2371.52(16.49)
	Slice efficient	422.09(19.90) 20.47(0.43)	4.291(3.30)	2371.50(16.53)
	Slice eff. thres.	217.35(10.25) 13.66(0.29)	4.226(3.19)	2371.54(16.60)
	Algo 8 (m = 2)	96.90(4.57)	10.90(0.23)	4.242(3.22)	2371.34(15.98)

Table 22 :

 22 Lepto data n = 1000, L D = 200, L M = 1000.

Appendix

In this appendix, we show the rest of the results from our comparative study, excepted results for lepto and bimod data with n = 10, 000.