Importance conditional sampling for Pitman-Yor mixtures

TL;DR

This paper introduces importance conditional sampling (ICS), a new efficient and stable method for Pitman-Yor mixture models, enhancing density estimation and clustering with better parameter robustness and extendability.

Contribution

The paper proposes ICS, a novel sampling strategy for Pitman-Yor mixtures that improves efficiency, interpretability, and stability over existing methods, and extends to other complex models.

Findings

ICS shows stable performance across various parameter settings.

The method is computationally efficient and parallelizable.

ICS can be extended to other nonparametric mixture models.

Abstract

Nonparametric mixture models based on the Pitman-Yor process represent a flexible tool for density estimation and clustering. Natural generalization of the popular class of Dirichlet process mixture models, they allow for more robust inference on the number of components characterizing the distribution of the data. We propose a new sampling strategy for such models, named importance conditional sampling (ICS), which combines appealing properties of existing methods, including easy interpretability and a within-iteration parallelizable structure. An extensive simulation study highlights the efficiency of the proposed method which, unlike other conditional samplers, shows stable performances for different specifications of the parameters characterizing the Pitman-Yor process. We further show that the ICS approach can be naturally extended to other classes of computationally demanding models, such as nonparametric mixture models for partially exchangeable data.

Tables

Table 1. Table B.1: Relative frequency of the of times that the bound 10 5 superscript 10 5 10^{5} on the number of jumps is reached when σ = 0.4 𝜎 0.4 \sigma=0.4 for the independent slice-efficient algorithm (I-SE) and the dependent slice-efficient algorithm (D-SE).

		I-SE	D-SE
$\vartheta =1$	$n=100$	0.00	0.00
	$n=250$	0.00	0.00
	$n=1000$	0.00	0.00
$\vartheta =10$	$n=100$	0.06	0.07
	$n=250$	0.8	0.14
	$n=1000$	0.16	0.31
$\vartheta =25$	$n=100$	0.30	0.39
	$n=250$	0.39	0.54
	$n=1000$	0.63	0.89

Equations

$$\tilde{f}(x)={\int }_{\Theta }\mathcal{K}(x;\theta )\mathrm{d}\tilde{p}(\theta ),$$

$$\begin{array}{rl}{X}_i\mid {\theta }_i& \stackrel{\text{mbox}ind}{\sim }\mathcal{K}(X_i;{\theta }_i),i=1,\dots ,n\\ {\theta }_i\mid \tilde{p}& \stackrel{\text{mbox}iid}{\sim }\tilde{p},\\ \tilde{p}& \sim PY(\sigma ,\vartheta ;P_0).\end{array}$$

$$\text{mbox}P({\theta }_{i+1}\in \mathrm{d}t\mid {\theta }_1,\dots ,{\theta }_i)=\frac{\vartheta +k_i\sigma }{\vartheta +i}{P}_0(\mathrm{d}t)+\sum _{j=1}^{k_i}\frac{n_j-\sigma }{\vartheta +i}{\delta }_{{\theta }_j^{\ast }}(\mathrm{d}t),$$

$$\tilde{p}=\sum _{j=1}^{\infty }{p}_j{\delta }_{{\tilde{\theta }}_j}.$$

$$p_1$$

$$p_j$$

$$V_j$$

$$p(\mathbf{X},\mathbf{\theta },\tilde{p})=\prod _{i=1}^n\mathcal{K}(X_i;{\theta }_i)\prod _{j=1}^{k_n}\tilde{p}(\mathrm{d}{\theta }_j^{\ast }{)}^{n_j}Q(\tilde{p}),$$

$$p_0\tilde{q}(\cdot )+\sum _{j=1}^{k_n}{p}_j{\delta }_{t_j^{\ast }}(\cdot ),$$

$$\text{mbox}P({\theta }_i\in \mathrm{d}t\mid X_i,\tilde{p})\propto p_0\mathcal{K}(X_i;t)\tilde{q}(\mathrm{d}t)+\sum _{j=1}^{k_n}{p}_j\mathcal{K}(X_i;t_j^{\ast }){\delta }_{t_j^{\ast }}(\mathrm{d}t),$$

$$w_{\ell }=\frac{\mathcal{K}(X_i,s_{\ell })\tilde{q}(s_{\ell })}{\tilde{q}(s_{\ell })}=\mathcal{K}(X_i,s_{\ell }),$$

$$\text{mbox}P({\theta }_i\in \mathrm{d}t\mid X_i,\tilde{p})\propto p_0\sum _{j=1}^{r_m}\frac{m_j}{m}\mathcal{K}(X_i;s_j^{\ast }){\delta }_{s_j^{\ast }}(\mathrm{d}t)+\sum _{j=1}^{k_n}{p}_j\mathcal{K}(X_i;t_j^{\ast }){\delta }_{t_j^{\ast }}(\mathrm{d}t).$$

$$\text{mbox}P(s_{\ell +1}\in \mathrm{d}s\mid s_1,\dots ,s_{\ell })=\frac{\vartheta +\sigma (k_n+r_{\ell })}{\vartheta +\sigma k_n+\ell }{P}_0(\mathrm{d}s)+\sum _{j=1}^{r_{\ell }}\frac{m_j-\sigma }{\vartheta +\sigma k_n+\ell }{\delta }_{s_j^{\ast }}(\mathrm{d}s),$$

$$\text{mbox}P(s_{\ell +1}^{(r)}\in \cdot \mid s_1^{(r)},\dots ,s_{\ell }^{(r)})=\frac{\vartheta +\sigma (k_n^{(r-1)}+r_{\ell }^{(r)})}{\vartheta +\sigma k_n^{(r-1)}+\ell }{P}_0(\cdot )+\sum _{j=1}^{r_{\ell }^{(r)}}\frac{m_j^{(r)}-\sigma }{\vartheta +\sigma k_n^{(r-1)}+\ell }{\delta }_{s_j^{\ast (r)}}(\cdot );$$

$$\text{mbox}P({\theta }_i^{(r)}=t\mid \cdots )\propto \{\begin{array}{ll}{p}_0^{(r)}\frac{m_{\ell }^{(r)}}{m}\mathcal{K}(X_i;s_{\ell }^{\ast (r)})& \text{ if }t\in \{s_1^{\ast (r)},\dots ,s_{r_m^{(r)}}^{\ast (r)}\}\\ p_j^{(r)}\mathcal{K}(X_i;t_j^{\ast (r)})& \text{ if }t\in \{t_1^{\ast (r)},\dots ,t_{k_n^{(r-1)}}^{\ast (r)}\}\\ 0& \text{ otherwise}\end{array}$$

$$\text{mbox}P({\theta }_j^{\ast }\in \mathrm{d}t\mid \mathbf{X})\propto P_0(\mathrm{d}t)\prod _{i\in C_j}\mathcal{K}(X_i;t),$$

$${\tilde{f}}_m(x)=p_0\sum _{l=1}^{r_m}\frac{m_l}{m}\mathcal{K}(x;s_l^{\ast })+\sum _{j=1}^{k_n}{p}_j\mathcal{K}(x;t_j^{\ast }).$$

$$\hat{f}(x)=\frac{1}{R-R_b}\sum _{r=R_b+1}^R{\tilde{f}}_m^{(r)}(x),$$

$$\tilde{p}\mid \mathbf{\theta }\sim DP\left(\vartheta +n;\frac{\vartheta }{\vartheta +n}{P}_0+\sum _{j=1}^{k_n}\frac{n_j}{\vartheta +n}{\delta }_{{\theta }_j^{\ast }}\right),$$

$$\text{mbox}dev(\mathbf{X},{\mathbf{\theta }}^{(r)})=-2\sum _{i=1}^n\left(\sum _{j=1}^{k_n^{(r)}}\frac{n_j^{(r)}}{n}\mathcal{K}(X_i;{\theta }_j^{\ast (r)})\right),$$

$${\mathbf{n}}_1=(n_{1,1},n_{2,1},\dots ,n_{12,1})=(236,51,59,38,92,56,67,51,61,187,81,44)$$

$${\mathbf{n}}_2=(n_{1,2},n_{2,2},\dots ,n_{12,2})=(245,73,91,39,113,98,74,90,56,197,70,144)$$

$$X_{i,l}\mid {\mathbf{\theta }}_1,\dots ,{\mathbf{\theta }}_L$$

$${\theta }_{i_l,l}\mid ({\tilde{p}}_1,\dots ,{\tilde{p}}_L)$$

$$({\tilde{p}}_1,\dots ,{\tilde{p}}_L)$$

$${\tilde{p}}_l={\gamma }_l{w}_l+{\gamma }_0(1-w_l),$$

$$N_n=\min \left\{l\ge 1:\sum _{j\le l}{p}_j>1-\min (U_1{p}_{c_1},\dots ,U_n{p}_{c_n})\right\},$$

$$M_n$$

$$\le \min \left\{l\ge 1:\sum _{j\le l}{p}_j>1-\min (U_1{p}_{c_1},\dots ,U_n{p}_{c_n})\right\}=N_n.$$

$$M_n-1\stackrel{\text{mbox}a.s.}{\sim }{\left(\frac{B_n{T}_{\sigma ,\vartheta }}{\sigma }\right)}^{-\sigma /(1-\sigma )},$$

$$M(ϵ)-1\stackrel{\text{mbox}a.s.}{\sim }{\left(\frac{ϵT_{\sigma ,\vartheta }}{\sigma }\right)}^{-\sigma /(1-\sigma )},$$

$$A$$

$$B$$

Full text

Keywords: Bayesian nonparametrics, Dependent Dirichlet process, Importance conditional sampling, Nonparametric mixtures, Pitman-Yor process, Sampling-importance resampling, Slice sampler.

1 Introduction

Bayesian nonparametric mixtures are flexible models for density estimation and clustering, nowadays a well-established modelling option for applied statisticians (Frühwirth-Schnatter et al., , 2019). The first of such models to appear in the literature was the Dirichlet process (DP) (Ferguson, , 1973) mixture of Gaussian kernels by Lo, (1984), a contribution which paved the way to the definition of a wide variety of nonparametric mixture models. In recent years, increasing interest has been dedicated to the definition of mixture models based on nonparametric mixing random probability measures that go beyond the DP (e.g. Nieto-Barajas et al., , 2004; Lijoi et al., 2005b, ; Lijoi et al., 2005a, ; Lijoi et al., , 2007; Argiento et al., , 2016). Among these measures, the Pitman-Yor process (PY) (Perman et al., , 1992; Pitman, , 1995) stands out for conveniently combining mathematical tractability, interpretability, and modelling flexibility (see, e.g., De Blasi et al., , 2015).

Let ${\bm{X}}=(X_{1},\dots,X_{n})$ be an $n$-dimensional sample of observations defined on some probability space $(\Omega,\mathscr{A},\mbox{P})$ and taking values in $\mathds{X}$, and $\mathscr{F}$ denote the space of all probability distributions on $\mathds{X}$. A Bayesian nonparametric mixture model is a random distribution taking values in $\mathscr{F}$, defined as

[TABLE]

where $\mathcal{K}(x;\theta)$ is a kernel and $\tilde{p}$ is a discrete random probability measure. In this paper we focus on $\tilde{p}\sim PY(\sigma,\vartheta;P_{0})$, that is we assume that $\tilde{p}$ is distributed as a PY process with discount parameter $\sigma\in[0,1)$, strength parameter $\vartheta>-\sigma$, and diffuse base measure $P_{0}\in\mathscr{F}$. The DP is recovered as a special case when $\sigma=0$. Model (1.1) can alternatively be written in hierarchical form as

[TABLE]

The joint distribution of $\bm{\theta}=(\theta_{1},\dots,\theta_{n})$ is characterized by the predictive distribution of the PY, which, for any $i=1,2,\ldots$, is given by

[TABLE]

where $k_{i}$ is the number of distinct values $\theta_{j}^{*}$ observed in the first $i$ draws and $n_{j}$ is the number of observed $\theta_{l}$, for $l=1,\ldots,i$, coinciding with $\theta_{j}^{*}$, such that $\sum_{j=1}^{k_{i}}n_{j}=i$.

Markov chain Monte Carlo (MCMC) sampling methods represent the gold standard for carrying out posterior inference based on nonparametric mixture models. Resorting to the terminology adopted by Papaspiliopoulos and Roberts, (2008), most of the existing MCMC sampling methods for nonparametric mixtures can be classified into marginal and conditional, the two classes being characterized by different ways to deal with the infinite-dimensional random probability measure $\tilde{p}$. While marginal methods rely on the possibility of analytically marginalizing $\tilde{p}$ out, the conditional ones exploit suitable finite-dimensional summaries of $\tilde{p}$.

Marginal methods for nonparametric mixtures were first devised by Escobar, (1988) and Escobar and West, (1995), contributions which focused on DP mixtures of univariate Gaussian kernels. Extensions of such proposal include the works of Müller et al., (1996), MacEachern, (1994), MacEachern and Müller, (1998), Neal, (2000), Barrios et al., (2013), Favaro and Teh, (2013), and Lomelí et al., (2017). It is worth noting that, despite being the first class of MCMC methods for Bayesian nonparametric mixtures appeared in the literature, marginal methods are still routinely used in popular packages such as the DPpackage (Jara et al., , 2011), the de facto standard software for many Bayesian nonparametric models. Alternatively, conditional methods rely on the use of summaries—of finite and possibly random dimension—of realizations of $\tilde{p}$. To this end, the stick-breaking representation for the PY (Pitman and Yor, , 1997) turns out to be very convenient. The almost sure discreteness of the PY allows $\tilde{p}$ to be written as an infinite sum of random jumps $\{p_{j}\}_{j=1}^{\infty}$ occurring at random locations $\{\tilde{\theta}_{j}\}_{j=1}^{\infty}$, that is

[TABLE]

The distribution of the locations is independent of that of the jumps and, while $\tilde{\theta}_{j}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}P_{0}$, the distribution of the jumps is characterized by the following construction:

[TABLE]

A first example of conditional approach can be found in Ishwaran and James, (2001) and Ishwaran and Zarepour, (2002), contributions that consider a fixed truncation of the stick-breaking representation of a large class of random probability measures, and provide a bound for the introduced truncation error. Along similar lines, Muliere and Tardella, (1998) and Arbel et al., (2019) make the truncation level of the DP and the PY, respectively, random so to make sure that the resulting error is smaller than a given threshold. Exact solutions avoiding introducing truncation errors are the slice samplers of Walker, (2007) and Kalli et al., (2011), the improved slice sampler of Ge et al., (2015), and the retrospective sampler of Papaspiliopoulos and Roberts, (2008). It is worth noticing that, although originally introduced for the case of DP mixture models, the ideas behind slice and retrospective sampling algorithms are naturally extended to the more general class of mixture models for which the mixing random probability measure admits a stick-breaking representation (Ishwaran and James, , 2001), thus including the PY mixture model as a special case. In this context Favaro and Walker, (2013) propose a general framework for slice sampling the class of mixtures of $\sigma$-stable Poisson–Kingman model. Henceforth we will use the term slice sampling to refer to the proposals of Walker, (2007) and Kalli et al., (2011), and not to the general definition of slice sampling.

Recent contributions have proposed hybrid strategies for posterior sampling nonparametric mixture models, that combine steps of marginal and conditional algorithms and therefore cannot be classified as either type of algorithm. Notable examples are the hybrid sampler of Lomelí et al., (2015) for the general class of Poisson–Kingman mixture models, and the hybrid approach proposed by Dubey et al., (2020) for a wide range of Bayesian nonparametric models based on completely random measures.

Marginal methods are appealing for their simplicity and for the fact that the number of random elements that must be drawn at each iteration of the sampler, i.e. the components of $\bm{\theta}$, is deterministic and thus bounded. At the same time, quantifying the posterior uncertainty, e.g. via posterior credible sets, by using the output of marginal methods is in general not straightforward since marginal methods do not generate realizations of the posterior distribution of $\tilde{f}$, but only of its conditional expectation $\mathds{E}[\tilde{f}\mid\bm{\theta},{\bm{X}}]$, where the expectation is taken with respect to $\tilde{p}$. To this end, convenient strategies have been proposed, which typically exploit the possibility of sampling approximate realizations of $\tilde{p}$ conditionally on the values of $\bm{\theta}$ generated by the marginal algorithm (see discussions in Gelfand and Kottas, , 2002; Taddy and Kottas, , 2012; Arbel et al., , 2016). Conditional methods, instead, produce approximate trajectories from the posterior distribution of $\tilde{f}$, which can be readily used to quantify posterior uncertainty. Moreover, by exploiting the conditional independence of the parameters $\theta_{i}$’s, given $\tilde{p}$ or a finite summary of it, conditional methods conveniently avoid sequentially updating the components of $\bm{\theta}$ at each iteration of the MCMC, thus leading to a fully parallelizable updating step within each iteration. On the other hand, the random truncation at the core of conditional methods such as slice and retrospective samplers makes the number of atoms and jumps that must be drawn at each iteration of the algorithm, random and unbounded. By confining our attention to the slice sampler of Walker, (2007) and, equivalently, its dependent slice-efficient version (Kalli et al., , 2011), we observe that, while its sampling routines are efficient and reliable when the DP case is considered, the same does not hold for the more general class of PY mixtures, specially when large values of $\sigma$ are considered. In practice, we noticed that, even for small sample sizes, the number of random elements that must be drawn at each iteration of the algorithm can be extremely large, often so large to make an actual implementation of the slice sampler for PY mixture models unfeasible. It is clear-cut that this limitation represents a major problem as the discount parameter $\sigma$ greatly impacts the robustness of the prior with respect to model-based clustering (see Lijoi et al., , 2007; Canale and Prünster, , 2017). In order to shed some light on this aberrant behaviour, we investigate the distribution of the random number $N_{n}$ of jumps that must be drawn at each iteration of a slice sampler, implemented to carry out posterior inference based on a sample of size $n$. We can define—see Appendix A for details—a data-free lower bound for $N_{n}$, that is a random variable $M_{n}$ such that $N_{n}(\omega)\geq M_{n}(\omega)$ for every $\omega\in\Omega$ and for every sample of size $n$. $M_{n}$ is distributed as $\min\left\{l\geq 1\,:\,\prod_{j\leq l}(1-V_{j})<B_{n}\right\}$, where the $V_{j}$’s are defined as in (1.5) and $B_{n}\sim\text{Beta}(1,n)$: studying the distribution of the lower bound $M_{n}$ will provide useful insight on $N_{n}$. Note that, in addition, $M_{n}$ coincides with the number of jumps to be drawn in order to generate a sample of size $n$ by adapting to the PY case the retrospective sampling idea introduced for the DP by Papaspiliopoulos and Roberts, (2008).

Figure 1 shows the empirical distribution of $M_{n}$, with $n=100$, for various combinations of $\vartheta$ and $\sigma$. The estimated median of the distribution of $M_{n}$ grows with $\sigma$ and, for any given value of $\sigma$, with $\vartheta$. It can be appreciated that the size of the values taken by $M_{n}$, and thus by $N_{n}$, explodes when $\sigma$ grows beyond $0.5$, fact which leads to the aforementioned computational bottlenecks in routine implementations of the slice sampler. For example, when $\sigma=0.8$, the estimated probability of $M_{n}$ exceeding $10^{9}$ is equal to $0.35$, $0.42$ and $0.63$, for $\vartheta$ equal to $0.1$, $1$ and $10$, respectively. From an analytic point of view, following Muliere and Tardella, (1998), it is easy to show that in the DP case (i.e. $\sigma=0$), $(M_{n}-1)\sim\mbox{Poisson}(\vartheta\log(1/B_{n}))$. Beyond the DP case (i.e. $\sigma\in(0,1)$), an application of Arbel et al., (2019) allows us to derive an analogous asymptotic result, which corroborates our empirical findings on the practical impossibility of using the slice sampler for PY mixtures with $\sigma\geq 0.5$. See Proposition A.1 and related discussion in the Appendix.

Herein, we propose a new sampling strategy, named importance conditional sampling (ICS), for PY mixture models, which combines the appealing features of both conditional and marginal methods, while avoiding their weaknesses, including the computational bottleneck depicted in Figure 1. Like marginal methods, the ICS has a simple and interpretable sampling scheme, reminiscent of Blackwell-MacQueen’s Pólya urn (Blackwell and MacQueen, , 1973), and allows to work with the update of a bounded number of random elements per iteration; at the same time, being a conditional method, it allows for fully parallelizable parameters update and it accounts for straightforward approximate posterior quantification. Our proposal exploits the posterior representation of the PY process, derived by Pitman, (1996) in combination with an efficient sampling-importance resampling idea. The structure of Pitman, (1996)’s representation makes it suitable for numerical implementations of PY based models, as indicated in Ishwaran and James, (2001), and nicely implemented by Fall and Barat, (2014).

The rest of the paper is organized as follows. The ICS is described in Section 2. Section 3 is dedicated to an extensive simulation study, comparing the performance of the ICS with state-of-the-art marginal and conditional sampling methods. Section 4 proposes (reports) an illustrative application, where the proposed algorithm is used to analyse a data set from the Collaborative Perinatal Project (Klebanoff, , 2009). In this context, Section 4.2 is dedicated to illustrate how the ICS approach can be extended to the case of nonparametric mixture models for partially exchangeable data. Section 5 concludes the paper with a discussion. Additional results are presented in the Appendix.

2 Importance conditional sampling

The random elements involved in a PY mixture model defined as in (1.2) are observations ${\bm{X}}$, latent parameters $\bm{\theta}$ and the PY random probability measure $\tilde{p}$. The joint distribution of $({\bm{X}},\bm{\theta},\tilde{p})$ can be written as

[TABLE]

where $\bm{\theta}^{*}=(\theta_{1}^{*},\dots,\theta_{k_{n}}^{*})$ is the vector of unique values in $\bm{\theta}$, with frequencies $(n_{1},\dots,n_{k_{n}})$ such that $\sum_{j=1}^{k_{n}}n_{j}=n$, and $Q$ is the distribution of $\tilde{p}\sim PY(\sigma,\vartheta;P_{0})$. In line of principle, the full conditional distributions of all random elements can be derived from (2.1) and used to devise a Gibbs sampler. Given that the vector ${\bm{X}}$, conditionally on $\bm{\theta}$, is independent of $\tilde{p}$, the update of $\bm{\theta}$ is the only step of the Gibbs sampler which works conditionally on a realization of the infinite-dimensional $\tilde{p}$. The conditional distribution $p(\bm{\theta}\mid{\bm{X}},\tilde{p})$ therefore will be the main focus of our attention: its study will allow us to identify a finite-dimensional summary of $\tilde{p}$, sufficient for the purpose of updating $\bm{\theta}$ from its full conditional distribution. As a result, as far as $\tilde{p}$ is concerned, only the update of its finite-dimensional summary will need to be included in the Gibbs sampler. Our proposal exploits a convenient representation of the posterior distribution of a PY process (Pitman, , 1996), reported in the next proposition.

Proposition 1**.**

*(Corollary 20 in *Pitman, , 1996)**. Let $t_{1},\ldots,t_{n}\mid\tilde{p}\sim\tilde{p}$ and $\tilde{p}\sim PY(\sigma,\vartheta;P_{0})$, and denote by $(t_{1}^{*},\ldots,t_{k_{n}}^{*})$ and $(n_{1},\ldots,n_{k_{n}})$ the set of $k_{n}$ distinct values and corresponding frequencies in $(t_{1},\ldots,t_{n})$. The conditional distribution of $\tilde{p}$, given $(t_{1},\ldots,t_{n})$, coincides with the distribution of

[TABLE]

where $(p_{0},p_{1},\ldots,p_{k_{n}})\sim\text{Dirichlet}(\vartheta+k_{n}\sigma,n_{1}-\sigma,\ldots,\allowbreak n_{k_{n}}-\sigma)$ and $\tilde{q}\sim PY(\sigma,\vartheta+k_{n}\sigma;P_{0})$ is independent of $(p_{0},p_{1},\ldots,p_{k_{n}})$.

In the context of mixture models, Pitman’s result implies that the full conditional distribution of $\tilde{p}$ coincides with the distribution of a mixture composed by a PY process $\tilde{q}$ with updated parameters, and a discrete random probability measure with $k_{n}$ fixed jump points at ${\bm{t}}=(t_{1}^{*},\ldots,t_{k_{n}}^{*})$. This means that, in the context of a Gibbs sampler, while, by conditional independence, the update of each parameter $\theta_{i}$ is done independently of the other parameters $(\theta_{1},\ldots,\theta_{i-1},\theta_{i+1},\ldots,\theta_{n})$, the distinct values $\bm{\theta}^{*}$ taken by the parameters at a given iteration, are carried on to the next iteration of the algorithm through $\tilde{p}$, in the form of fixed jump points ${\bm{t}}$. Specifically, if $\Theta^{*}=\Theta\setminus\{t^{*}_{1},\dots,t^{*}_{k_{n}}\}$, then, for every $i=1,\dots,n$, the full conditional distribution of the $i$-th parameter $\theta_{i}$ can be written as

[TABLE]

where $\tilde{q}$ is the restriction of $\tilde{p}$ to $\Theta^{*}$, $p_{0}=\tilde{p}(\Theta^{*})$ and $p_{j}=\tilde{p}(t_{j}^{*})$, for every $j=1,\ldots,k_{n}$. The full conditional in (2.2) is reminiscent of the Blackwell-MacQueen urn scheme characterizing the update of the parameters in marginal methods: the parameter $\theta_{i}$ can either coincide with one of the $k_{n}$ fixed jump points of $\tilde{p}$ or take a new value from a distribution proportional to $\mathcal{K}(X_{i};t)\tilde{q}(\mathrm{d}t)$. The key observation at the basis of the ICS is that, for the purpose of updating the parameters $\bm{\theta}$, there is no need to know the whole realization of $\tilde{p}$ but it suffices to know the vector ${\bm{t}}$ of fixed jump points of $\tilde{p}$, the value ${\bm{p}}=(p_{0},p_{1},\ldots,p_{k_{n}})$ taken by $\tilde{p}$ at the partition $(\Theta^{*},t^{*}_{1},\ldots,t^{*}_{k_{n}})$ of $\Theta$, and to be able to sample from a distribution proportional to $\mathcal{K}(X_{i},t)\tilde{q}(\mathrm{d}t)$. For the latter task, we adopt a sampling-importance resampling approach (see, e.g., Smith and Gelfand, , 1992) with proposal distribution $\tilde{q}$. It is remarkable that such solution allows us to approximately sample from the target distribution while avoiding the daunting task of simulating a realization of $\tilde{q}$ itself. Indeed, for any $m\geq 1$, a vector ${\bm{s}}=(s_{1},\ldots,s_{m})$ such that $s_{i}\mid\tilde{q}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\tilde{q}$ can be generated by means of an urn scheme exploiting (1.3). Given the almost sure discreteness of $\tilde{q}$, the generated vector will show ties with positive probability and thus will feature $r_{m}\leq m$ distinct values $(s_{1}^{*},\ldots,s_{r_{m}}^{*})$, with frequencies $(m_{1},\ldots,m_{r_{m}})$ such that $\sum_{j=1}^{r_{m}}m_{j}=m$. In turn, importance weights for the resampling step are computed, for any $\ell=1,\ldots,m$, as

[TABLE]

thus without requiring the evaluation of $\tilde{q}$. As a result, the full conditional (2.2) can be rewritten as

[TABLE]

Once more we highlight an interesting analogy between the conditional approach we propose and marginal methods: the introduction of the auxiliary random variables $s_{1}^{*},\ldots,s_{r_{m}}^{*}$ reminds of the augmentation introduced in Algorithm 8 of Neal, (2000), marginal algorithm proposed to deal with a non-conjugate specification of the mixture model. From (2.3) it is straightforward to identify $({\bm{s}},{\bm{t}},{\bm{p}})$ as a finite-dimensional summary of $\tilde{p}$, sufficient for the purpose of updating the parameters $\theta_{i}$ from their full conditionals. This means that, as far as $\tilde{p}$ is concerned, only its summary $({\bm{s}},{\bm{t}},{\bm{p}})$ must be included in the updating steps of the Gibbs sampler. To this end, Proposition 1 provides the basis for the update of $({\bm{s}},{\bm{t}},{\bm{p}})$. Indeed, conditionally on $\bm{\theta}$, the fixed jump points ${\bm{t}}$ coincide with the $k_{n}$ distinct values appearing in $\bm{\theta}$, while the random vectors ${\bm{p}}$ and ${\bm{s}}$ are independent with ${\bm{p}}\sim\text{Dirichlet}(\vartheta+\sigma k_{n},n_{1}-\sigma,\ldots,n_{k_{n}}-\sigma)$ and the joint distribution of ${\bm{s}}$ characterized by the predictive distribution of a PY$(\sigma,\vartheta+\sigma k_{n};P_{0})$, that is, for any $\ell=0,1,\ldots,m-1$,

[TABLE]

where $(s_{1}^{*},\ldots,s_{r_{\ell}}^{*})$ is the vector of $r_{\ell}$ distinct values appearing in $(s_{1},\ldots,s_{\ell})$, with corresponding frequencies $(m_{1},\ldots,m_{r_{\ell}})$ such that $\sum_{j=1}^{r_{\ell}}m_{j}=\ell$.

By combining the steps just described, as summarized in Algorithm 1, we can then devise a Gibbs sampler which we name ICS. In Algorithm 1 and henceforth, the superscript $(r)$ is used to denote the value taken by a random variable at the $r$-th iteration. In order to improve mixing, the ICS includes an acceleration step which consists in updating, at the end of each iteration, the distinct values $\bm{\theta}^{*}$ from their full conditional distributions. Namely, for every $j=1,\ldots,k_{n}$,

[TABLE]

where $C_{j}=\{i\in\{1,\ldots,n\}\,:\,\theta_{i}=\theta_{j}^{*}\}$.

Finally, a realization from the posterior distribution of $({\bm{s}},{\bm{t}},{\bm{p}})$ defines an approximate realization $f$ of the posterior distribution of the random density defined in (1.1), that is

[TABLE]

If the algorithm is run for a total of $R$ iterations, the first $R_{b}$ of which discarded as burn-in, then the posterior mean is estimated by

[TABLE]

where $\tilde{f}_{m}^{(r)}$ denotes the approximate density sampled from the posterior at the $r$-th iteration. The set of densities $\tilde{f}_{m}^{(r)}$ can be also used to quantify posterior uncertainty. It is worth remarking though that any such quantification is based on realizations of a finite dimensional summary of the infinite-dimensional $\tilde{p}$ and thus is, by its nature, approximated. For a quantification of the approximating error one could resort to Arbel et al., (2019).

It is instructive to consider how the ICS works for the special case of DP mixture models, that is when $\sigma=0$. In such case, the steps described in Algorithm 1 can be nicely interpreted by resorting to three fundamental properties characterizing the DP, namely conjugacy, self-similarity, and availability of finite-dimensional distributions. More specifically, when $\sigma=0$, step 4 of Algorithm 1 consists in generating the random weights ${\bm{p}}$ from a Dirichlet distribution of parameters $(\vartheta,n_{1},\ldots,\allowbreak n_{k_{n}})$. This follows by combining the conjugacy of the DP (Ferguson, , 1973), for which

[TABLE]

with the availability of finite-dimensional distributions of DP (Ferguson, , 1973), which provides the distribution of ${\bm{p}}$, defined as the evaluation of the conditional distribution of $\tilde{p}$ on the partition of $\Theta$ induced by $\bm{\theta}$. Moreover, when $\sigma=0$, according to the predictive distribution displayed in step 6 of Algorithm 1, the auxiliary random variables ${\bm{s}}$ are exchangeable from $\tilde{q}\sim DP(\vartheta;P_{0})$, with $\tilde{q}$ independent of ${\bm{p}}$. This is nicely implied by the self-similarity of the DP (see, e.g., Ghosal, , 2010), according to which $\tilde{q}=\tilde{p}|_{\Theta^{*}}$ is independent of $\tilde{p}|_{\Theta\setminus\Theta^{*}}$, and therefore of ${\bm{p}}$, and is distributed as a $DP(\vartheta P_{0}(\Theta^{*});P_{0}|_{\Theta^{*}})$, and by the diffuseness of $P_{0}$. As a result, in the DP case, the auxiliary random variables ${\bm{s}}$ are generated from the prior model.

3 Simulation study

We performed a simulation study to analyze the performance of the ICS algorithm and to compare it with marginal and slice samplers. For the latter, two versions proposed by Kalli et al., (2011) were considered, namely the dependent and the independent slice-efficient algorithms. The independent version of the algorithm requires the specification of a deterministic sequence $\xi_{1},\xi_{2},\ldots$, which in our implementation was set equal to $\mathds{E}[p_{1}],\mathds{E}[p_{2}],\ldots$, with the $p_{j}$’s defined in (1.5), in analogy with what was proposed by Kalli et al., (2011) for the DP (see Algorithm 5 in the Appendix for more details). All algorithms were written in C++ and are implemented in the BNPmix package (Corradin et al., , 2021), available on CRAN. Aware that different implementations can lead to a biased comparison (see Kriegel et al., , 2017, for an insightful discussion), we aimed at reducing such bias to a minimum by letting the four algorithms considered here share the same code for most sub-routines.

Throughout this section we consider synthetic data generated from a simple two-component mixture of Gaussians, namely $f_{0}(x)=0.75\phi(x;-2.5,1)+0.25\phi(x;2.5,1)$, with $\phi(\cdot;\mu,\sigma^{2})$ denoting the density of a Gaussian random variable with mean $\mu$ and variance $\sigma^{2}$. All data were analyzed by means of the nonparametric mixture model defined in (1.1) and specified by considering a univariate Gaussian kernel $\mathcal{K}(x,\theta)=\phi(x;\mu,\sigma^{2})$, with $\theta=(\mu,\sigma^{2})$, and by assuming a normal-inverse gamma base measure $P_{0}$ such that $\sigma^{2}\sim IG(2,1)$ and $\mu\mid\sigma^{2}\sim N(0,5\sigma^{2})$. Different combinations of values for the parameters $\sigma$ and $\vartheta$, and for the sample size $n$ were considered. The results of this section are then obtained as averages over a specified number of replicates. All algorithms were run for $1\,500$ iterations, of which the first $500$ discarded as burn-in. Convergence of the chains was checked by visual inspection of the trace plots of randomly selected runs, which did not provide any evidence against it. The analysis was carried out by running BNPmix on R 4.0.3 on a 64-bit Windows machine with a 3.4-GHz Intel quad-core i7-3770 processor and 16 GB of RAM.

The first part of our investigation is dedicated to the role of $m$, the size of the auxiliary sample generated for the sampling-importance resampling step within the ICS. To this end, we considered two sample sizes, namely $n=100$ and $n=1\,000$, and generated $10$ data sets per size. Such data were then analyzed by considering a combination of values for the PY parameters, namely $\sigma\in\{0,0.2,0.4,0.6,0.8\}$ and $\vartheta\in\{1,10\}$, and by running the ICS with $m\in\{1,10,100\}$. Estimated posterior densities, not displayed here, did not show any noticeable effect of $m$. More interesting findings were obtained when the analysis focused on the quality of the generated posterior sample: larger values for $m$ appear to lead to a better mixing of the Markov chain at the price of additional computational cost. These effects were measured by considering the effective sample size (ESS), computed by resorting to the CODA package (Plummer et al., , 2006), and the ratio between runtime, in seconds, and ESS (time/ESS), both averaged over 100 replicates. Following the algorithmic performance analyses of Neal, (2000), Papaspiliopoulos and Roberts, (2008) and Kalli et al., (2011), the ESS was computed on the number of clusters—$k_{n}$ as far as the ICS is concerned— and on the deviance of the estimated density, with the latter defined as

[TABLE]

for the $r$-th MCMC draw. The ratio time/ESS takes into account both quality of the generated sample and computational cost, and can be interpreted as the average time needed to sample one independent draw from the posterior.

The results show that larger values of $m$ lead, on average, to a larger ESS, that is to better quality posterior samples. This is displayed in the top row of Figure 2, which shows the estimated ESS for $k_{n}$. We observe that, when averaging over all the considered scenarios, the ESS obtained by setting $m=100$ is $1.82$ and $1.09$ times larger than the average ESS obtained by setting $m=1$ and $m=10$, respectively. At the same time, larger values of $m$ require drawing more random objects per iteration and thus, as expected, lead to longer runtimes. In this sense, the bottom row of Figure 2 clearly indicates that, as far as $k_{n}$ is concerned, the ratio time/ESS tends to be larger for larger values of $m$. This is particularly evident, for example, when $\sigma=0.8$ as the ratio time/ESS corresponding to $m=100$ is, on average, $1.81$ and $1.06$ times larger than the same ratio corresponding to $m=1$ and $m=10$, respectively. Similar conclusions can be drawn by looking at Figure B.1, presented in Appendix B and displaying time/ESS for the deviance of the estimated densities.

When implementing the ICS, the value of $m$ can be tuned based on the desired algorithm performance in terms of quality of mixing and runtime. As for the rest of the paper, redand for the ease of illustration, we will work with $m=10$, chosen as a sensible compromise between good mixing and controlled computational cost.

The second part of the simulation study compares the performance of ICS, marginal sampler, dependent and independent slice-efficient samplers. For the sake of clarity, pseudo-code of the implemented algorithms is provided in Appendix D. We considered the sample sizes $n=100$, $n=250$ and $n=1\,000$, and generated 10 data sets per size from $f_{0}$. These data were then analyzed by considering a combination of values for the PY parameters, namely $\sigma\in\{0,0.2,0.4,0.6,0.8\}$ and $\vartheta\in\{1,10,25\}$. The results we report are obtained, for each scenario, by averaging over the 10 replicates. As for the two slice samplers, due to the aforementioned explosion of the number of drawings per iteration when $\sigma$ takes large values, our analysis was forcefully confined to the case $\sigma\leq 0.4$. Moreover, the results referring to the case $\sigma=0.4$ are approximate as they were obtained by constraining the slice sampler to draw at most $10^{5}$ components at each iteration: such limitation of our study could not be avoided, given the otherwise unmanageable computational burden associated with this specific setting. Table B.1 in Appendix B shows that such bound was reached more often when large data sets were analyzed. For example, while for $n=100$ the bound was reached on average 12% and 15% of the iterations, for independent and dependent slice-efficient samplers respectively, the same happened on average 26% and 40% of the iterations when $n=1\,000$. For this reason, these specific results must be considered approximated and, as far as the runtime is concerned, conservative. The four algorithms were compared by using the same measures adopted in the first part of the simulation study, namely the ESS for the number of clusters, the ESS for the deviance of the estimated density, and the corresponding ratios time/ESS.

A clear trend can be appreciated in Figure 3 where the focus is on the ESS for the number of clusters: the marginal sampler displays, on average, a larger ESS than ICS, whose ESS appears, in turn, uniformly larger than the ones characterizing the two slice samplers. As for the latter two, while the displayed trend is similar, it can be appreciated that the independent algorithm is uniformly characterized by a better mixing. Results referring to the ratio time/ESS, for the variable number of clusters, are displayed in Figure 4. ICS and marginal sampler show in general similar performances. It is interesting to notice though that, while the ratio time/ESS for the marginal algorithm is rather stable over the values of $\sigma$ considered in the study, the same quantity for ICS indicates a slightly better performance when $\sigma$ takes large values. On the counterpart, the efficiency of the two slice samplers is heavily affected by the value of $\sigma$, with time/ESS exploding when $\sigma$ moves from 0 to 0.4 and when $\vartheta$ is increasing. On the basis of this study, the slice samplers appear competitive options when $\sigma\in\{0,0.2\}$ and a small $\vartheta$ are considered. On the contrary, it is apparent that larger values of $\sigma$ make the two slice samplers less efficient than ICS and marginal sampler. Similar considerations can be drawn when analyzing the performance in terms of deviance of the estimated densities, with the plots for ESS and the ratio time/ESS displayed in Figures B.2 and B.3.

4 Illustrations

We consider a data set from the Collaborative Perinatal Project (CPP), a large prospective study of the cause of neurological disorders and other pathologies in children in the United States. Pregnant women were enrolled between 1959 and 1966 when they showed up for prenatal care at one of 12 hospitals. While several measurements per pregnancy are available, our attention focuses on two main quantities: the gestational age (in weeks) and the logarithm of the concentration level of DDE, a persistent metabolite of the pesticide DDT, known to have adverse impact on the gestational age (Longnecker et al., , 2001). Our analysis has a two-fold goal. First, we focus on estimating and comparing the joint density of gestational age and DDE for two groups of women, namely smokers and non-smokers. This will also allow us to assess how the probability of premature birth varies conditionally on the level of DDE. Adopting a nonparametric mixture model will allow us to investigate the presence of clusters within the data. Second, we consider the data set partitioned in the 12 hospitals of the study and focus on the estimation of the hospital-specific distribution of the gestational age, by accounting for possible association across subsamples collected at different hospitals. For this analysis we adopt a nonparametric mixture model for partially exchangeable data and propose an extension of the ICS approach presented in Section 2.

4.1 Cross-hospital analysis

Smokers and non-smokers groups have sample size of $n_{1}=1023$ and $n_{2}=1290$, respectively. For the two groups we independently model the joint distribution of gestational age and DDE by means of a PY mixture model (1.2) with bivariate Gaussian kernel function $\mathcal{K}(x,\bm{\theta})=\phi(x,\bm{\theta})$, with $\bm{\theta}=(\bm{\mu},\bm{\Sigma})$, and with conjugate normal-inverse Wishart base measure $P_{0}=N\text{-}IW({\bm{m}}_{0},k_{0},\nu_{0},{\bm{S}}_{0})$. In absence of precise prior information on the density to be estimated, we specify a vague base measure following an empirical Bayes approach. Specifically we let ${\bm{m}}_{0}$ be equal to the sample average, ${\bm{S}}_{0}$ be equal to three times the empirical covariance, $k_{0}=1/10$, and $\nu_{0}=5$. These settings are equivalent to assuming that the scale parameter of the generic mixture component coincides with 1.5 times the empirical covariance, while the location parameter is centered on the sample mean with prior variance equal to 10 times the scale parameter. Next, we set the parameters $\vartheta$ and $\sigma$ on the basis of the prior distribution they imply on the number of clusters $k_{n}$, within each group. Specifically, we set the prior expectation and prior standard deviation for $k_{n}$ equal to 10 and 20, respectively. Our choice implies that a small probability ($\approx 0.05$) is assigned to the event $k_{n}\geq 50$. This argument leads to set $(\sigma,\vartheta)$ equal to $(0.548,-0.485)$ and $(0.5295,-0.4660)$ for the groups of smokers and non-smokers, respectively. The values specified for $\sigma$ are thus larger than $0.5$, situation that is conveniently tackled by the ICS, as displayed by the simulation study of Section 3. An alternative modelling strategy is achieved by introducing a hyperprior distribution for both $\sigma$ and $\theta$. While not explored in this illustration, it is worth stressing that this strategy might be conveniently implemented by adopting the ICS: if the prior on $\sigma$ is defined on $(0,1)$, an implementation of the model requires a sampler whose efficiency is not compromised by the specific values of $\sigma$ explored by the chain.

The analysis of both samples was carried out by running the ICS for $12\,000$ iterations, with the first $7\,000$ discarded as burn-in. Convergence of the chain was assessed as satisfactory by visually investigating the trace plots and by means of the Geweke’s diagnostics (Geweke, , 1992). Running the analysis of the two samples took less than two minutes in total. It is important to stress that, given the model specification, the same analysis could not be carried out by implementing the slice samplers described in Algorithms 4 and 5, as the value of $\sigma$ would make computation time endless. We could instead implement the marginal sampler described in Algorithm 3, which, as expected, took considerably longer than the ICS (about 11 minutes), due to the moderately large sample sizes.

The contour curves of the estimated joint densities of gestational age and DDE for the two groups are displayed in the left panel of Figure 5 and suggest different distributions between smokers and non-smokers, specially when large values for DDE are considered. Differences between the two groups are further highlighted by the right panel of Figure 5, which shows the estimated probability—along with corresponding pointwise 90% posterior credible bands—of premature birth (i.e. gestational age smaller than 37 weeks), conditionally on the value taken by DDE, for the two groups. Once again, a difference between smokers and non-smokers can be appreciated for large values of DDE, although a sizeable uncertainty is associated with posterior estimates, as displayed by the large credible bands.

4.2 Multi-hospital analysis

The same data set as in the previous section is considered here, with observations classified according to both smoking habits of women and the hospitals where they were enrolled. This leads to two samples stratified into $L=12$ strata, with cardinalities summarized by the vectors

[TABLE]

and

[TABLE]

for smokers and non-smokers, respectively. The focus of the analysis is modelling the distribution of gestational age.

4.2.1 A mixture model for partially exchangeable data

Smokers and non-smokers data are analyzed independently. For each group, heterogeneity across hospitals suggests to assume that data are partially exchangeable in the sense of de Finetti, (1938). To account for this assumption, we consider a mixture model for partially exchangeable data, where the stratum-specific mixing random probability measures form the components of a dependent Dirichlet process. Within this flexible class of processes (see Foti and Williamson, , 2015, and references therein), we consider the Griffiths-Milne dependent Dirichlet processes (GM-DDP), as defined and studied in Lijoi et al., 2014a ; Lijoi et al., 2014b . For an allied approach see Griffin et al., (2013). Let $X_{i,l}$ be the gestational age of the $i$-th woman in the $l$-th hospital, and $\bm{\theta}_{l}$ be the vector of latent variables $\theta_{i,l}$ referring to the $l$-th hospital. The mixture model can be represented in its hierarchical form as

[TABLE]

with $l=1,\ldots,L,i=1,\ldots,n_{l}$, $\vartheta>0$, $z\in(0,1)$, $P_{0}$ is a probability distribution on $\mathds{R}\times\mathds{R}^{+}$, and the GM-DDP distribution of the vector $(\tilde{p}_{1},\ldots,\tilde{p}_{L})$ coincides with the distribution of the vector of random probability measures whose components are defined, for every $l=1,\ldots,L$, as

[TABLE]

where $\gamma_{1},\ldots,\gamma_{L}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}DP(\vartheta z;P_{0})$ and $\gamma_{0}\sim DP(\vartheta(1-z);P_{0})$ is independent of $\gamma_{l}$, for any $l=1,\ldots,L$. Moreover, the vector of random weights ${\bm{w}}=(w_{1},\dots,w_{L})$, taking values in $[0,1]^{L}$, is distributed as a multivariate beta of parameters $(\vartheta z,\ldots,\vartheta z,\vartheta(1-z))$, as defined in Olkin and Liu, (2003), and its components are independent of the random probability measures $\gamma_{0},\gamma_{1},\ldots,\gamma_{L}$. As a result, the random probabilities $\tilde{p}_{l}$ are, marginally, identically distributed with $\tilde{p}_{l}\sim DP(\vartheta;P_{0})$ (see Lijoi et al., 2014a, , for details).

4.2.2 ICS for GM-DDP mixture model and its application

The ICS can be easily adapted to a variety of models. For example, it naturally fits the partially exchangeable framework of model (4.2.1). The ICS algorithm for GM-DDP mixture models is described in Algorithm 2, and consists of three main steps. First, conditionally on the allocation of observations to clusters referring to either the idiosyncratic process $\gamma_{l}$, with $l=1,\ldots,L$, or the common process $\gamma_{0}$, summaries of all the processes, that is $({\bm{s}}_{l},{\bm{t}}_{l},{\bm{p}}_{l})$, for $l=0,\dots,L$, are updated as done in Section 2 for a single process, with the proviso that $\sigma=0$. Second, the latent variables $\theta_{i,l}$ are updated for every $l=1,\dots,L$ and $1\leq i\leq n_{l}$; and, third, the components of ${\bm{w}}$ are sampled. The full conditional distributions for $\theta_{i,l}$ and ${\bm{w}}$ are provided in Appendix C.

Model (4.2.1) is specified by assuming a univariate Gaussian kernel and normal-inverse gamma base measure $P_{0}=N\text{-}IG(0,5,4,1)$. Moreover, the specification $\vartheta=1$ and $z=0.5$ is adopted, with the latter choice corresponding to equal prior weights assigned to idiosyncratic and common components $\gamma_{l}$ and $\gamma_{0}$. The ICS algorithm for the GM-DDP mixture model was run for 10 000 iterations, the first 5 000 of which were discarded as burn-in. Estimating posterior densities for smokers and non-smokers, required a total runtime of less than two and a half minutes. Convergence of the chains was assessed by visually investigating the trace plots, which did not provide any evidence against it.

Figure 6 shows the estimated densities of the gestational age, for each stratum, with a comparison between smokers and non-smokers. The distribution for smokers is globally more skewed and shifted to the left than the one for non-smokers, indicating an expected more adverse effect of smoke on gestational age.

5 Discussion

We proposed a new sampling strategy for PY mixture models, named ICS, which combines desirable properties of existing marginal and conditional methods: the ICS shares easy interpretability with marginal methods, while allowing, likewise conditional samplers, for a parallelizable update of the latent parameters $\bm{\theta}$, and for a straightforward quantification of posterior uncertainty. The simulation study of Section 3 showed that the ICS overtakes some of the computational bottlenecks characterizing the conditional methods considered in the comparison. Specifically, the ICS can be implemented for any value of the discount parameter $\sigma$, with its efficiency being stable to the specification of $\sigma$. This is appealing as the discount parameter plays a crucial modelling role when PY mixture models are used for model-based clustering: the ICS allows for an efficient implementation of such models, without the need of setting artificial constraints on the value of $\sigma$. As far as the comparison of the performances of ICS and other algorithms is concerned, it is important to remark that the independent slice-efficient algorithm proposed by Kalli et al., (2011) is more general than the one considered in Section 3 as other specifications of the deterministic sequence $\xi_{1},\xi_{2},\ldots$ are possible. As nicely discussed by Kalli et al., (2011), the choice of such sequence “is a delicate issue and any choice has to balance efficiency and computational time”. Alternative specifications of the sequence may be explored on a case-by-case basis but, in our experience, the computational time can be reduced only at the cost of worsening the mixing of the algorithm. It is also worth remarking that the ICS does not rely on any assumption of conjugacy between base measure and kernel, and thus it can be considered by all means a valid alternative to the celebrated Algorithm 8 of Neal, (2000) when non-conjugate mixture models are to be implemented. Finally, while originally introduced to overtake computational problems arising in the implementation of algorithms for PY mixture models, the idea behind the ICS approach can be naturally extended to other classes of computationally demanding models. As an example, we implemented the same idea to deal with posterior inference based on a flexible class of mixture models for partially exchangeable data. Other extensions are also possible and are currently subject of ongoing research.

Acknowledgments. The first author is supported by the University of Padova under the STARS Grant. The second and the third authors are grateful to the DEMS Data Science Lab for supporting this work by providing computational resources.

Appendix A On the number of jumps to be drawn with the slice sampler

Let $N_{n}$ be the random number of jumps which need to be drawn at each iteration of a slice sampler (Walker, , 2007) or, equivalently, its dependent slice-efficient version (Kalli et al., , 2011), implemented to carry out posterior inference based on a sample of size $n$. Conditionally on the cluster assignment variables $c_{1},\ldots,c_{n}$ and on the weights $p_{c_{1}},\ldots,p_{c_{n}}$ of the non-empty components of the mixture, $N_{n}$ is given by

[TABLE]

where the random weights $p_{j}$’s are defined as in (1.5) and $U_{1},\ldots,U_{n}$ are independent uniform random variables, independent of the weights $p_{j}$’s. We next define a second random variable $M_{n}$, function of the same uniform random variables $U_{1},\ldots,U_{n}$, as

[TABLE]

The random number $M_{n}$ is thus a data-free lower bound for $N_{n}$, where the inequality $M_{n}(\omega)\leq N_{n}(\omega)$ holds for every $\omega\in\Omega$. Studying the distribution of $M_{n}$ will shed light on the distribution of its upper bound $N_{n}$. Interestingly, $M_{n}$ represents also the random number of jumps to be drawn in order to generate a sample of size $n$ from a PY by adapting the retrospective sampling idea of Papaspiliopoulos and Roberts, (2008), described in their Section 2 for the DP case. The distribution of $M_{n}$ coincides with the distribution of $\min\left\{l\geq 1\,:\,\prod_{j\leq l}(1-V_{j})<B_{n}\right\}$, where the stick-breaking variables $(V_{j})_{j\geq 1}$ are defined as in (1.5) and $B_{n}$ is a beta random variable with parameters $1$ and $n$. Following Muliere and Tardella, (1998), it is easy to show that, when $\sigma=0$, then $M_{n}-1$ is distributed as a mixture of Poisson distributions, specifically $(M_{n}-1)\sim\text{Poisson}(\vartheta\log(1/B_{n}))$. This leads to $\mathds{E}[M_{n}]=\vartheta H_{n}+1$, where $H_{n}=\sum_{l=1}^{n}l^{-1}$ is the $n$-th harmonic number. It is worth noting that, for $n\rightarrow\infty$, $\mathds{E}[M_{n}]\approx\vartheta\log(n)$, that is the growth is logarithmic in $n$, while the contribution of $\vartheta$ is linear. As for the PY process, we resort to Arbel et al., (2019), where the asymptotic distribution of the minimum number of jumps of a PY, needed to guarantee that the truncation error is smaller than a deterministic threshold, is studied. We introduce the notation $a_{n}\stackrel{{\scriptstyle\mbox{\scriptsize{a.s.}}}}{{\sim}}b_{n}$ to indicate that $\mbox{P}(\lim_{n\rightarrow\infty}a_{n}/b_{n}=1)=1$ and, by exploiting Theorem 2 in Arbel et al., (2019), we prove the following proposition.

Proposition A.1**.**

Let $M_{n}=\min\left\{l\geq 1\,:\,\prod_{j\leq l}(1-V_{j})<B_{n}\right\}$ where the sequence $(V_{j})_{j\geq 1}$ is defined as in (1.5) and $B_{n}$ is a beta random variable with parameters $1$ and $n$. Then, for $n\rightarrow\infty$,

[TABLE]

where $T_{\sigma,\vartheta}$, independent of $B_{n}$, is a polynomially tilted stable random variable (Devroye, , 2009), with probability density function proportional to $t^{-\vartheta}f_{\sigma}(x)$, where $f_{\sigma}$ is the density function of a unilateral stable random variable with Laplace transform equal to $\exp\{-\lambda^{\sigma}\}$.

Proof.

Define $M(\epsilon)=\min\left\{l\geq 1\,:\,\prod_{j\leq l}(1-V_{j})<\epsilon\right\}$. Following Arbel et al., (2019),

[TABLE]

as $\epsilon\rightarrow 0$. Observe that $M_{n}=M(B_{n})$ and that $B_{n}\stackrel{{\scriptstyle\mbox{\scriptsize{a.s.}}}}{{\sim}}0$ as $n\rightarrow\infty$. We then define the events

[TABLE]

and observe that $C\subset A\cup B$. Which implies that $\mbox{P}(C)\leq\mbox{P}(A\cup B)\leq\mbox{P}(A)+\mbox{P}(B)=0$. ∎

If we define $L_{n}=\left(B_{n}T_{\sigma,\vartheta}/\sigma\right)^{-\sigma/(1-\sigma)}$, for any positive integer $n$, the statement of Proposition A.1 is tantamount to $M_{n}-1\stackrel{{\scriptstyle\mbox{\scriptsize{a.s.}}}}{{\sim}}L_{n}$ as $n\rightarrow\infty$. The random variable $L_{n}$ has finite mean if and only if $\sigma\in(0,1/2)$, case in which $\mathds{E}[L_{n}]=c_{\sigma,\vartheta}\Gamma(n+1)/\Gamma(n+2-1/(1-\sigma))$, where

[TABLE]

which implies that $\mathds{E}[L_{n}]\approx c_{\sigma,\theta}n^{\sigma/(1-\sigma)}$, when $n\rightarrow\infty$. A simple simulation experiment was run to empirically investigate the quality of the asymptotic approximation of $M_{n}$ provided by $L_{n}$. The random variable $T_{\sigma,\vartheta}$ appearing in the defintion of $L_{n}$ was sampled by resorting to Hofert, (2011). Figure A.1 displays the estimated probability of the events $M_{n}>10^{6}$ and $L_{n}>10^{6}$, as a function of $\sigma\in(0,1)$, for $\vartheta\in\{0.1,1,10\}$ and for different sample sizes $n\in\{100,1\,000,10\,000\}$.

Appendix B Additional details on the simulation study

This section provides additional results of the simulation study presented in Section 3. Table B.1 reports on the number of times the upper bound for the number of jumps drawn at each iteration of dependent and independent slice-efficient samplers was reached. Figures B.1, B.2 and B.3 focus on the functional deviance and display results analogous to those presented in Section 3 for the random variable number of clusters.

Appendix C ICS for GM-DDP

In order to describe the full conditional distributions of $\theta_{i,l}$ and ${\bm{w}}$, and to provide the pseudo-code of the ICS for the GM-DDP mixture model, some notation needs to be introduced. Let $r_{m,0}$ and $r_{m,l}$, for $l=1,\ldots,L$, represent the number of distinct values $s_{j,0}^{*}$ and $s_{j,l}^{*}$ appearing in the vectors ${\bm{s}}_{0}$ and ${\bm{s}}_{l}$, respectively. The corresponding frequencies are given by $m_{j,0}$ and $m_{j,l}$, and are such that $\sum_{j=1}^{k_{m,l}}m_{j,l}=m$, for every $l=0,1,\ldots,L$. Let $\bm{\theta}_{0}^{*}$ be the vector of distinct values appearing in $(\bm{\theta}_{1},\ldots,\bm{\theta}_{L})$ coinciding with either the $k_{\mathbf{n},0}$ fixed jump points ${\bm{t}}_{0}$ of the common process $\gamma_{0}$ or with any of the $r_{m,0}$ values appearing in ${\bm{s}}_{0}$. Similarly, for any $l=1,\ldots,L$, $\bm{\theta}_{l}^{*}$ denotes the vector of distinct values appearing in $(\bm{\theta}_{1},\ldots,\bm{\theta}_{L})$ coinciding with either the $k_{\mathbf{n},l}$ fixed jump points ${\bm{t}}_{l}$ of the idiosyncratic process $\gamma_{l}$ or with any of the $r_{m,l}$ values appearing in ${\bm{s}}_{l}$. Finally, we let $\mathcal{C}_{j,0}=\{(i,l)\;:\;\theta_{i,l}=\theta^{*}_{j,0}\}$ and, for $l=1,\ldots,L$, $\mathcal{C}_{j,l}=\{i\;:\;\theta_{i,l}=\theta^{*}_{j,l}\}$.

The full conditional distribution of $\theta_{i,l}$, for every $l=1,\dots,L$ and $1\leq i\leq n_{l}$, is given, up to a proportionality constant, by

[TABLE]

The full conditional for ${\bm{w}}$ is given, up to a proportionality constant, by

[TABLE]

where

[TABLE]

[TABLE]

The pseudo-code of the ICS for the GM-DDP mixture model is presented in Algorithm 2.

Appendix D Implementation of the algorithms in Section 3

This section reports the pseudo-code of marginal, dependent slice-efficient and independennt slice-efficient samplers, algorithms which were implemented for the performance comparison described in Section 3. For the sake of simplicity, all the algorithms are described without specifying prior distributions for the hyperparameters. Algorithm 3 is based on Escobar and West, (1995). Algorithms 4 and 5 are implemented by following the dependent and independent slice-efficient versions of the slice sampler described in Kalli et al., (2011).