Gibbs sampling for mixtures in order of appearance: the ordered allocation sampler

Mixture models represent one of the most successful applications of Bayesian methods. Bayesian inference proceeds by placing a prior on the mixing distribution, whose atoms and their sizes represent the mixture component parameters and the weights, respectively. An important issue is that the number of components is rarely known in advance. The nonparametric approach consists in modelling the mixing distribution with infinitely many support points, and infer the number of components through the number of groups observed in the data (e.g. Escobar and West, 1995). Another alternative is to assign a prior to this unknown quantity (cf. Richardson and Green, 1997; Miller and Harrison, 2018). Posterior inference for mixture models is customarily based on Gibbs sampling methods which have been broadly classified into two categories: marginal and conditional samplers. Marginal methods (Escobar and West, 1995; Neal, 2000) are termed this way because they partially integrate out the mixing distribution and exploit the generalized Pólya urn scheme representation (Blackwell and MacQueen, 1973; Pitman, 2006) of the prediction rule of a sample from the mixing distribution. By doing so marginal samplers avoid dealing with a potentially infinite or random model dimension. They are well suited for models such as the Dirichlet (Ferguson, 1973; Escobar and West, 1995), the Pitman-Yor (Pitman and Yor, 1997; Ishwaran and James, 2001) and mixtures of finite mixtures (Miller and Harrison, 2018). However, they are challenging to adapt to mixing priors without a tractable prediction rule. In alternative, one can use conditional methods which include the mixing distribution and update it as a component of the sampler. While being more widely applicable they bring some issues in their design and implementation. In infinite mixture models, some sort of truncation either deterministic or random is necessary to avoid dealing with infinitely many mixture components. Finite dimensional approximations of the mixing prior were proposed in (Ishwaran and James, 2001). Of random truncation type are the exact conditional samplers derived by Walker (2007); Kalli et al. (2011) and Papaspiliopoulos and Roberts (2008). It has been observed that they require Metropolis-Hastings steps that swap components’ labels so to speed up mixing (Porteous et al., 2006; Papaspiliopoulos and Roberts, 2008). As for mixtures with a random number of components, the main challenge of conditional samplers is the need to explore parameter spaces of different dimensions. The standard method is the reversible jump MCMC algorithm (Richardson and Green, 1997) but this can be difficult to implement.

In this paper we contribute both methodologically and computationally to mitigate these issues by developing a novel conditional Gibbs sampling method named the ordered allocation sampler. The sampler works with the mixture components in the random order in which they are discovered. To derive it we use in depth the theory of species sampling models set forth in Pitman (1995, 1996a, 1996b) where, in particular, it is established that the law of the weights in order of appearance corresponds to the distribution of the weights that is invariant under size-biased permutations. This one admits a simple stick-breaking representation for the Dirichlet and the Pitman-Yor processes. Working with this specific rearrangement of mixture components allows us to exploit a (conditional) prediction rule in the sampler. Thus, it bears similarities with marginal methods such as the fact that mixing takes place in the space of partitions and not in the space of cluster’s labels as it occurs in other conditional samplers (Porteous et al., 2006). Empirical studies confirm that this endows our sampler with nice mixing properties. A second major advantage of our proposal is that, since $n$ data points can not be generated from more than $n$ distinct components, at most the sampler needs to update the first $n$ components in order of appearance. This is especially relevant for infinite mixture models as it avoids truncation. In particular, for Pitman-Yor processes with slowly decaying weights, our sampler proves to be very convenient computationally wise. A third important consequence is that marginalization over the weights in order of appearance yields exactly the exchangeable partition probability function (EPPF Pitman, 1996b). For mixtures with random dimension, this translates to a simple way of updating the number of components without resorting to a reversible jump step. Finally, as other conditional methods, the ordered allocation sampler allows direct inference on the mixing distribution and can be adapted to a wide range of mixing priors.

The rest of the paper is organized as follows. In Section 2 we provide background theory on species sampling priors in mixture models. It will set the stage for the ordered allocation sampler. In Section 3 we first derive the sampler for models with tractable size-biased permuted weights, and later we show how to adapt it when the law of this arrangement of the weights is not available in explicit form. In Section 4 we illustrate the performance of our sampler with well-known real and simulated datasets. Some concluding remarks and discussion points are brought in Section 5. Proofs, technical details and additional illustrations are in the Appendix.

In Bayesian mixture models we model exchangeable data, $(y_{i})=(y_{i})_{i=1}^{n}$, taking values in a Borel space, $(\mathbb{Y},\mathcal{B}(\mathbb{Y}))$, as conditionally independent and identically distributed (iid) from

where $g(\cdot\mid x)$ is a density for each $x$, and the mixing distribution, $P=\sum_{j=1}^{m}p_{j}\delta_{x_{j}}$, is an almost surely discrete random probability measure over the Borel parameter space $(\mathbb{X},\mathcal{B}(\mathbb{X}))$. Species sampling models, introduced and studied by Pitman (1996b), constitute a very general class of random probability measures that provide a convenient prior specification for the mixing distribution $P$. In species sampling models, $P=\sum_{j=1}^{m}p_{j}\delta_{x_{j}}$, the atoms, $(x_{j})=(x_{1},\ldots,x_{m})$, are iid from a diffuse distribution, $\nu$, over $(\mathbb{X},\mathcal{B}(\mathbb{X}))$, and are independent of $(m,(p_{j}))$. The weights $(p_{j})=(p_{1},\ldots,p_{m})$ are positive random variables with $\sum_{j=1}^{m}p_{j}=1$ almost surely, and the number of support points, $m$, can be finite, infinite or random. A sequence, $(\theta_{i})=(\theta_{i})_{i=1}^{\infty}$, is a species sampling sequence driven by $P$ if it is exchangeable and the almost sure limit of the empirical distributions, $P=\lim_{n\to\infty}n^{-1}\sum_{i=1}^{n}\delta_{\theta_{i}},$ is a species sampling model. By de Finetti’s theorem, the latter is equivalent to the existence of a species sampling model, $P$, such that given $P$, $(\theta_{i})$ are conditionally iid according to $P$ (cf. Theorem 1.1 and Proposition 1.4 of Kallenberg, 2005). The laws of $(\theta_{i})$ and $P$ determine each other, and both are fully determined by that of $(m,(p_{j}))$, and the diffuse distribution, $\nu$. A key aspect to note is that the law of $P$ is invariant under weights permutations. That is, $P=\sum_{j=1}^{m}p_{j}\delta_{x_{j}}$ is equal in distribution to $\sum_{j=1}^{m}p_{\rho(j)}\delta_{x_{j}}$ for every permutation $\rho$ of $\{1,\ldots,m\}$, which means that working with an ordering of the weights or another does not change the mixing prior. This is reflected through the so called exchangeable partition probability function (EPPF), given by

where the sum ranges over all $k$-tupples of distinct positive integers, and $p_{j}=0$ for $j>m$. In fact, $\pi(n_{1},\ldots,n_{k})$ describes the probability that a sample, $(\theta_{1},\ldots,\theta_{n})$, of size $n=\sum_{j=1}^{k}n_{j}$, exhibits exactly $k$ distinct values, with corresponding frequencies, $n_{1},\ldots,n_{k}$ (Pitman, 1996b, 2006). Whenever $\pi(n_{1},\ldots,n_{k})$ can be computed in closed form, the prediction rule, $\mathbb{P}[\theta_{i+1}\in\cdot\mid\theta_{1},\ldots,\theta_{i}]$, of $(\theta_{i})$ is available and can be described in terms of a generalized Pólya urn scheme (cf Blackwell and MacQueen, 1973).

The invariance under permutations of $\sum_{j=1}^{m}p_{\rho(j)}\delta_{x_{j}}$, as well as the complexity of computing the unordered sum in (2), have motivated the study of weights permutations that simplify the analysis. An ordering of the weights of paramount importance is the size-biased permutation, $(\tilde{p}_{j})=(\tilde{p}_{1},\ldots,\tilde{p}_{m})$, given by $\tilde{p}_{j}=p_{\alpha_{j}}$, and $(\alpha_{j})=(\alpha_{1},\ldots,\alpha_{m})$ defined by

In other words, $(\alpha_{j})$ and $(\tilde{p}_{j})$ are sampled without replacement from $\{1,\ldots,m\}$ and $(p_{j})$, respectively, with probabilities $(p_{j})$. By construction the distribution of $(\tilde{p}_{j})$ is invariant under size-biased permutations. As shown by Pitman (1995, 1996a), the EPPF (2) can be computed through

hence, if the distribution of $(\tilde{p}_{j})$ is available, it becomes easier to compute $\pi(n_{1},\ldots,n_{k})$. Another advantage of working with size-biased permutations is that $\tilde{p}_{j}$ coincides with the long-run proportion of indexes $i$ such that $\theta_{i}=\tilde{x}_{j}$, where $\tilde{x}_{j}$ is the $j$th distinct value to appear in $(\theta_{i})$. Furthermore, the conditional law of $(\theta_{i})$ given $(\tilde{x}_{j})=(\tilde{x}_{1},\ldots,\tilde{x}_{m})$ and $(\tilde{p}_{j})$ admits a simple prediction rule as detailed next.

Theorem 1 is based on theory laid down in Pitman (1995, 1996b), nonetheless, we provide a self-contained proof in Appendix A, due to the crucial role it plays in the derivation of the new sampler.

The canonical example of species sampling models in Bayesian nonparametric statistics is the Dirichlet process (Ferguson, 1973). It has $m=\infty$ support points and its size-biased permuted weights, $(\tilde{p}_{j})$, admit the stick-breaking representation

where $(v_{j})=(v_{j})_{j=1}^{\infty}$ are iid from the Beta distribution $\mathsf{Be}(1,\theta)$ (Sethuraman, 1994). The Dirichlet model can be generalized to the two-parameter $(\sigma,\theta)$-model (Pitman, 2006) which features size-biased permuted weights as in (5) with independent $v_{j}\sim\mathsf{Be}(1-\sigma,\theta+j\sigma)$ according to one of the following two regimes:

1. a)

   $\sigma\in[0,1)$ and $\theta>-\sigma$. In this case $m=\infty$ and the species sampling model $P$ has been named the Pitman-Yor process by Ishwaran and James (2001) after Pitman and Yor (1997). Evidently the choice $\sigma=0$ reduces to a Dirichlet process.

2. b)

   Given $m\in\mathbb{N}$, $\sigma=-\gamma<0$, and $\theta=m\gamma$. In agreement with the notation we have established $m$ stands for the number of support points of $P$. It turns out that the law of $(\tilde{p}_{j})$ corresponds to that of the size-biased permutation of symmetric Dirichlet weights, $(p_{1},\ldots,p_{m})\sim\mathsf{Dir}(\gamma,\ldots,\gamma)$ (Pitman, 1996a). When $\gamma$ is fixed and $m$ is random $P$ belongs to the class of Gibbs-type priors (see De Blasi et al., 2015, for a recent review), while the allied mixture model corresponds to the mixture of finite mixtures of Miller and Harrison (2018).

Another type of species sampling models for which a stick-breaking characterization of size-biased weights is available are homogeneous normalized random measures with independent increments (cf. Regazzini et al., 2003; Favaro et al., 2016). Unfortunately, such characterization remains elusive for most species sampling models used in mixture modelling, examples are finite dimensional approximations of the Pitman-Yor process (Ishwaran and James, 2001), the Geometric process (Fuentes-García et al., 2010), the probit stick-breaking process (Rodríguez and Dunson, 2011) and exchangeable stick-breaking processes studied by Gil–Leyva and Mena (2021). For all these species sampling priors the weights can be defined in terms of a stick-breaking decomposition, $p_{j}=v_{j}\prod_{l=1}^{j-1}(1-v_{l})$, for some sequence of random variables $(v_{j})$ with values in $[0,1]$, yet $(p_{j})$ is not invariant under size-biased permutations.

As mentioned in Section 2, in mixture models data points $(y_{i})=(y_{i})_{i=1}^{n}$ are treated as conditionally iid from a random density $Q$ as in (1). Whenever the mixing distribution, $P$, is a species sampling model we can equivalently assume $y_{i}\mid\theta_{i}\sim g(\cdot\mid\theta_{i})$, independently for $i\leq n$, where $(\theta_{i})$ is a species sampling sequence driven by $P$. In this setting, marginal samplers integrate out $P$ and exploit the exchangeability of $(\theta_{i})$ as well as the prediction rule, $\mathbb{P}[\theta_{i+1}\in\cdot\mid\theta_{1},\ldots,\theta_{i}]$, to derive an algorithm for posterior inference (cf. Neal, 2000; Favaro and Teh, 2013; Miller and Harrison, 2018). Instead, conditional samplers (e.g. Ishwaran and James, 2001; Papaspiliopoulos and Roberts, 2008; Kalli et al., 2011) include the mixing distribution, $P$, and update its atoms, $(x_{j})$, and weights, $(p_{j})$, as components of the sampler. The ordered allocation sampler is a conditional sampler as it includes the mixing distribution, $P$, however similarly to marginal samplers it relies on a prediction rule for species sampling sequences. Explicitly, motivated by Theorem 1 we work with the atoms, $(\tilde{x}_{j})$, and weights, $(\tilde{p}_{j})$, of $P$ in the order in which they were discovered by $(\theta_{i})$. As commonly done in other samplers, we augment the model with latent allocation variables that identify each observation, $y_{i}$, with the mixture component it was sampled from. Here, in accordance with the order of appearance we introduce what we call ordered allocation variables, $(d_{i})=(d_{i})_{i=1}^{n}$, given by $d_{i}=j$ if and only if $y_{i}$ was sampled from $g(\cdot\mid\tilde{x}_{j})$, i.e. $\theta_{i}=\tilde{x}_{j}$. Thus, $\theta_{i}=\tilde{x}_{d_{i}}$, and $y_{i}\mid((\tilde{x}_{j}),d_{i})\sim g(\cdot\mid\tilde{x}_{d_{i}})$, independently for $i\leq n$. If $k_{i}$ denotes the number of distinct values in $(\theta_{1},\ldots,\theta_{i})$ then $k_{i}$ coincides with $\max\{d_{1},\ldots,d_{i}\}$, and $d_{i+1}$ necessarily takes a value in $\{1,\ldots,k_{i}+1\}$. More precisely, iii of Theorem 1 allows us to compute

for $i\leq n$, independently of elements in $(\tilde{x}_{j})$ (see also iv of Theorem 1). This yields the augmented likelihood

where $k_{n}=\max\{d_{1},\ldots,d_{n}\}$, $D_{j}=\{i\leq n:d_{i}=j\}$, $n_{j}=|D_{j}|$, and $\mathcal{D}$ is the event that $\{D_{1},\ldots,D_{k_{n}}\}$ is a partition of $\{1,\ldots,n\}$ with blocks in the least element order, in particular $D_{j}\neq\emptyset$, for $j\leq k_{n}$, and $\min\left(D_{1}\right)<\min\left(D_{2}\right)<\cdots<\min\big{(}D_{{k_{n}}}\big{)}$. The full conditional distributions required at each iteration of the sampler are proportional to the product of (7) times the prior distributions, $\mathbbm{p}[(\tilde{x}_{j})]$ and $\mathbbm{p}[(\tilde{p}_{j})]$, of the atoms and weights of $P$ in order of appearance. We first derive the ordered allocation sampler for those species sampling mixing distributions where the prior of $(\tilde{p}_{j})$ can be modelled directly as is the case of the $(\sigma,\theta)$-model. Latter we explain how to adapt the sampler for the more general case where the law of $(\tilde{p}_{j})$ is not available.