A hybrid sampler for Poisson-Kingman mixture models
Mar\'ia Lomel\'i
Gatsby Unit
University College London
mlomeli@gatsby.ucl.ac.uk
This paper concerns the introduction of a new Markov Chain Monte Carlo scheme for posterior sampling in Bayesian nonparametric mixture models with priors that belong to the general Poisson-Kingman class. We present a novel compact way of representing the infinite dimensional component of the model such that while explicitly representing this infinite component it has less memory and storage requirements than previous MCMC schemes. We describe comparative simulation results demonstrating the efficacy of the proposed MCMC algorithm against existing marginal and conditional MCMC samplers.
NIPSPKpaper_resubmission_Arxivready_24sept2015.tex
\documentclass{article} %
\usepackage{nips15submit_e,times}
\usepackage{hyperref}
\usepackage{url}
\usepackage{amsopn}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{mathabx}
\usepackage{enumerate}
\usepackage{geometry}
 \usepackage{graphicx} 
 \usepackage{algpseudocode}
 \usepackage{algorithm}
\usepackage[numbers]{natbib} 


\title{A hybrid sampler for Poisson-Kingman mixture models}


\author{
Mar\'ia Lomel\'i \\
Gatsby Unit\\
University College London\\
\texttt{mlomeli@gatsby.ucl.ac.uk} \\
\And
Stefano Favaro \\
 Department of Economics and Statistics\\ University of Torino\\ and Collegio Carlo Alberto\\
\texttt{stefano.favaro@unito.it } \\
\And
Yee Whye Teh \\
Department of Statistics \\
University of Oxford \\
\texttt{y.w.teh@stats.ox.ac.uk } \\
}

\vspace{-5mm}

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}
\DeclareMathOperator*{\CRM}{CRM}
\DeclareMathOperator*{\NRM}{NRM}
\DeclareMathOperator*{\PK}{PK}
\DeclareMathOperator*{\NGG}{NGG}
\def\PP{\mathbb{P}}
\def\EE{\mathbb{E}}
 \def\simind{\stackrel{\mbox{\scriptsize{ind}}}{\sim}}
\def\simiid{\stackrel{\mbox{\scriptsize{iid}}}{\sim}}

\nipsfinalcopy %

\begin{document}

\graphicspath{ {./figures/} }

\maketitle

\begin{abstract}
This paper concerns the introduction of a new Markov Chain Monte Carlo scheme for posterior sampling in Bayesian nonparametric mixture models with priors that belong to the general Poisson-Kingman class. We present a novel compact way of representing the infinite dimensional component of the model such that while explicitly representing this infinite component it has less memory and storage requirements than previous MCMC schemes.  We describe  comparative  simulation  results  demonstrating  the  efficacy  of  the  proposed MCMC algorithm against existing marginal and conditional MCMC samplers.
\end{abstract}

\section{Introduction}

According to \citet{Zou15}, models that have a nonparametric component give us more flexiblity that could lead to better predictive performance. This is because their capacity to learn does not saturate hence their predictions should continue to improve as we get more and more data. Furthermore, we are able to fully consider our uncertainty about predictions thanks to the Bayesian paradigm. However, a major impediment to the widespread use of Bayesian nonparametric models is the problem of inference. Over the years, many MCMC methods have been proposed to perform inference which usually rely on a tailored representation of the underlying process \cite{EscWes1995a, Esc1994a, Nea1998b,PapRob2008a, Wal2007a, FavTeh2013a}.
This is an active research area since dealing with this infinite dimensional component forbids the direct use of standard simulation-based methods for posterior inference. These methods usually require a finite-dimensional representation and there are two main sampling approaches to facilitate simulation in the case of Bayesian nonparametric models: random truncation and marginalization. These two schemes are known in the literature as conditional and marginal samplers. 


In conditional samplers, the infinite-dimensional prior is replaced by a finite-dimensional representation chosen according to a truncation level. 
In marginal samplers, the need to represent the infinite-dimensional component can be bypassed by marginalising it out. Marginal samplers have less storage requirements than conditional samplers but could potentially have worst mixing properties. However, not integrating out the infinite dimensional compnent leads to a more comprehensive representation of the random probability measure, useful to compute expectations of interest with respect to the posterior.




In this paper, we propose a novel class of MCMC samplers for Poisson-Kingman mixture models, a very large class of Bayesian nonparametric mixture models that encompass all previously explored ones in the literature. 
 Our approach is based on a hybrid scheme that combines the main strengths of both conditional and marginal samplers.   In the flavour of probabilistic programming, we view our contribution as a step towards wider usage of flexible Bayesian nonparametric models, as it allows automated inference in probabilistic programs built out of a wide variety of Bayesian nonparametric building blocks.


\section{Poisson-Kingman processes}

Poisson-Kingman random probability measures (RPMs) have been introduced in \citet{Pit2003a} as a generalization of homogeneous Normalized Random Measures (NRMs) \cite{RegLijPru2003a, Jam2002a}. Let $\mathbb{X}$ be a complete and separable metric space endowed with the Borel $\sigma$-field $\mathcal{B}( \mathbb{X})$, let $\mu\sim\CRM(\rho,H_0)$ be a homogeneous  Completely Random Measure (CRM) with L\'evy measure$\rho$ and base distribution $H_0$ on this space, see \citet{Kin1967a} for a good overview about CRMs and references therein. Then, the corresponding total mass of $\mu$ is $T=\mu(\mathbb{X})$ and let it be finite, positive almost surely, and absolutely continuous with respect to Lebesgue measure. For any $t\in\mathbb{R}^+$, let us consider the conditional distribution of $\mu/t$ given that the total mass $T\in dt$. This distribution is denoted by $\PK(\rho,\delta_t,H_0)$, it is the distribution of a RPM, $\delta_t$ denotes the usual Dirac delta function. Poisson-Kingman RPMs form a class of RPMs whose distributions are obtained by mixing $\PK(\rho,\delta_t,H_0)$, over $t$, with respect to some distribution $\gamma$ on the positive real line. Specifically, a Poisson-Kingman RPM has following the hierarchical representation
\begin{align}
T &\sim \gamma \nonumber \label{PKP}\\
P|T=t &\sim \PK(\rho,\delta_t,H_0).
\end{align}
The RPM $P$ is referred to as the Poisson-Kingman RPM with L\'evy measure $\rho$, base distribution $H_0$ and mixing distribution $\gamma$. Throughout the paper we denote by $\PK(\rho,\gamma,H_0)$ the distribution of $P$ and, without loss of generality, we will assume that $\gamma(\text{d}t)\propto h(t) f_\rho(t)\text{d}t$ where $f_{\rho}$ is the density of the total mass $T$ under the CRM and $h$ is a non-negative function. Note that, when $\gamma(\text{d}t) = f_\rho(t)\text{d}t$ then the distribution $\PK(\rho,f_\rho,H_0)$ coincides with $\NRM(\rho,H_0)$. The resulting $P = \sum_{k\geq1}{p_k\delta_{\phi_{k}}}$ is almost surely discrete and since $\mu$ is homogeneous, the atoms $(\phi_k)_{k\ge 1}$ of $P$ are independent of their masses $(p_k)_{k\ge 1}$ and form a sequence of independent random variables identically distributed according to $H_{0}$. Finally, the masses of $P$ have distribution governed by the L\'evy measure $\rho$ and the distribution $\gamma$. 

One nice property is that $P$ is almost surely discrete: if we obtain a sample $\left\{Y_i\right\}_{i=1}^n$ from it, there is a positive probability of $Y_i=Y_j$ for each pair of indexes $i\neq j$.  This induces a random partition $\Pi$ on $\mathbb{N}$, where $i$ and $j$ are in the same block in $\Pi$ if and only if $Y_i=Y_j$.  \citet{Kin1978a} showed that $\Pi$ is exchangeable, this property will be one of the main tools for the derivation of our hybrid sampler.

\subsection{Size-biased sampling Poisson-Kingman processes}

A  second object induced by a Poisson-Kingman RPM is a size-biased permutation of its atoms.  Specifically, order the blocks in $\Pi$ by increasing order of the least element in each block, and for each $k\in\mathbb{N}$ let $Z_k$ be the least element of the $k$th block.  $Z_k$ is the index among $(Y_i)_{i\ge 1}$ of the first appearance of the $k$th unique value in the sequence.  Let $\tilde{J}_k=\mu(\{Y_{Z_k}\})$ be the mass of the corresponding atom in $\mu$.  Then $(\tilde{J}_k)_{k\ge 1}$ is a size-biased permutation of the masses of atoms in $\mu$, with larger masses tending to appear earlier in the sequence.  It is easy to see that $\sum_{k\ge 1}\tilde{J}_k=T$, and that the sequence can be understood as a stick-breaking construction: starting with a stick of length $T_0=T$; break off the first piece of length $\tilde{J}_1$; the surplus length of stick is $T_1=T_0-\tilde{J}_1$; then the second piece with length $\tilde{J}_2$ is broken off, etc. 

Theorem 2.1 of \citet{PerPitYor1992a} states that the sequence of surplus masses  $(T_k)_{k\ge 0}$ forms a Markov chain and gives the corresponding initial distribution and transition kernels. The corresponding generative process for the sequence $(Y_i)_{i\ge 1}$ is as follows:  

\begin{itemize}
\item[i)] Start with drawing the total mass from its distribution $\PP_{\rho,h,H_0}(T\in dt) \propto h(t) f_\rho(t)dt$.%
\item[ii)] The first draw $Y_1$ from $P$ is a size-biased pick from the masses of $\mu$.  The actual value of $Y_1$ is simply $Y^*_1\sim H_0$, while the mass of the corresponding atom in $\mu$ is $\tilde{J}_1$, with conditional distribution
$$\PP_{\rho,h,H_0}(\tilde{J}_1\in ds_1|T\in dt) = \frac{s_1}{t}\rho(ds_1) \frac{f_\rho(t-s_1)}{f_\rho(t)},\hspace{5mm}\text{with surplus mass }\hspace{5mm}T_1=T-\tilde{J}_1.$$
\item[iii)] For subsequent draws $i\geq2$:
\begin{itemize}
\item Let $K$ be the current number of distinct values among $Y_1,\ldots,Y_{i-1}$, and $Y^*_1,\ldots,Y^*_K$ the unique values, i.e., atoms in $\mu$.  The masses of these first $K$ atoms are denoted by $\tilde{J}_1,\ldots,\tilde{J}_K$ and the surplus mass is $T_K=T-\sum_{k=1}^K \tilde{J}_k$.
\item For each $k\le K$, with probability $\tilde{J}_k/T$, we set $Y_i=Y^*_k$.
\item With probability $T_K/T$, $Y_i$ takes on the value of an atom in $\mu$ besides the first $K$ atoms.  The actual value  $Y^*_{K+1}$ is drawn from $H_0$, while its mass is
drawn from
$$\PP_{\rho,h,H_0}(\tilde{J}_{K+1}\in ds_{K+1}|T_K\in dt_K) 
= \frac{s_{K+1}}{t_K}\rho(ds_{K+1}) \frac{f_\rho(t_K-s_{K+1})}{f_\rho(t_K)}, \hspace{5mm} T_{K+1}=T_K-\tilde{J}_{K+1}.$$
\end{itemize}
\end{itemize}
By multiplying the above infinitesimal probabilities one obtains the joint distribution of the random elements $T$, $\Pi$,  $(\tilde{J}_{i})_{i\geq1}$ and $(Y^{\ast}_{i})_{i\geq1}$:

\begin{align}\label{dist_prop}
&\PP_{\rho,h,H_0}( \Pi_n=(c_k)_{k\in[K]}, Y_k^*\in dy_k^*, \tilde{J}_k\in ds_k\text{ for }k\in[K], T\in \text{d}t)\\
&\notag\quad= t^{-n} f_\rho(t-\textstyle\sum_{k=1}^{K} s_k)h(t)\text{d}t
\displaystyle \prod_{k=1}^{K} s_k^{|c_k|} \rho(ds_k)H_0(dy_k^*),
\end{align}
where $(c_k)_{k\in[K]}$ denotes a particular partition of $[n]$ with $K$ blocks, $c_1,\ldots,c_K$, ordered by increasing least element  and $|c_k|$ is the cardinality of block $c_k$. The distribution \eqref{dist_prop} is invariant to the size-biased order. Such a joint distribution was first obtained in \citet{Pit2003a} , see also \citet{Pit2006a} for further details. 

\subsection{Relationship to the usual Stick-breaking construction}
In the generative process above, we mentioned that it is reminiscent of the well known stick breaking construction from \citet{IshJam2001a}, where you break a stick of length one, but it is not exactly the same. However, by starting with equation \eqref{dist_prop}, we can recover the usual construction due to two useful identities in distribution: $P_j\stackrel{d}{=}\frac{\tilde{J}_j}{T-\sum_{\ell<j}{\tilde{J}_{\ell}}}$ and $V_j\stackrel{d}{=}\frac{P_j}{1-\sum_{\ell<j}{P_{\ell}}}$ for $j=1,\hdots,K$. Indeed, we can reparameterize the model using these identities and then obtain the corresponding joint in terms of $K$ $(0,1)$-valued stick-breaking weights $\left\{V_j\right\}_{j=1}^K$ which correspond to the usual stick-breaking representation. Note that this joint distribution is for a general L\'evy measure $\rho$, density $f_{\rho}$ and it is conditioned on the valued of the random variable $T$. Even so, we can recover the standard Stick breaking representations for the Dirichlet and Pitman-Yor processes, for a specific choice of $\rho$ and if we integrate out $T$. However, in general, these stick-breaking random variables form a sequence of dependent random variables with a complicated distribution, except for the two previously mentioned processes, see \citet{Pit96} for details. 
\subsection{Poisson-Kingman mixture model}
We are mainly interested in using Poisson-Kingman RPMs as a building block for an infinite mixture model. Indeed, we can use Equation \eqref{PKP} as the top level of the following hierarchical specification:

\vspace{-7mm}

\begin{align}\label{infmixturemodel}
T&\sim \gamma \nonumber \\
P|T &\sim \PK(\rho_{\sigma}, \delta_T, H_0)\nonumber\\
Y_i \mid P &\stackrel{_{\textrm{iid}}} {\sim}P\nonumber\\
X_i \mid Y_i &\stackrel{\textrm{ind}}{\sim} F(\cdot \mid Y_i)
\end{align}

\vspace{-4mm}

where $F(\cdot \mid Y)$ is the likelihood term for each mixture component, and our dataset consists of  $n$ observations $(x_i)_{i\in[n]}$ of the corresponding variables $(X_i)_{i\in[n]}$.  We will assume that $F(\cdot \mid Y)$ is smooth. After specifying the model we would like to carry out inference for clustering and/or density estimation tasks. We can do it exactly and more efficiently than with known MCMC samplers with our novel approach. In the next section, we present our main contribution and in the following one we show how it outperforms other samplers.
\vspace{-4mm}
\section{Hybrid Sampler}
Equation's \eqref{dist_prop} joint distribution is written in terms of the first $K$ size-biased weights. In order to obtain a complete representation of the RPM, we would need to size-bias sample from it for a countably infinite number of times. Succesively, some way of representing exactly this object in a computer with finite memory and storage is needed. 

We introduce the following novel strategy: starting from equation \eqref{dist_prop}, we exploit the generative process of section 2.1 when reassigning observations to clusters. In addition to this, we reparameterize the model in terms of a surplus mass random variable $V= T-\sum_{k=1}^K \tilde{J}_k$ and end up with the following joint distribution:%
\begin{align}\label{dist_prop2}
&\PP_{\rho,h,H_0}( \Pi_n=(c_k)_{k\in[K]}, Y_k^*\in dy_k^*, \tilde{J}_k\in ds_k\text{ for }k\in[K], T-\sum_{k=1}^K \tilde{J}_k\in dv,X_i\in dx_i\text{ for } i\in[n])\\
&\notag\quad= (v+\sum_{k=1}^K{s_k})^{-n}h\left(v+\sum_{k=1}^K{s_k}\right) f_\rho(v)
\displaystyle \prod_{k=1}^{K} s_k^{|c_k|} \rho(ds_k)H_0(dy_k^*) \prod_{i\in c_k} F(dx_i|y^*_k).
\end{align}
For this reason, while having a complete representation of the infinite dimensional part of the model we only need to explicitly represent those size-biased weights associated to occupied clusters plus a surplus mass term which is associated to the rest of the empty clusters, as  Figure 1 shows. The cluster reassignment step can be seen as a lazy sampling scheme since we explicitly represent and update the weights associated to occupied clusters  and create a size-biased weight only when a new cluster appears. To make this possible we use the induced partition and we call Equation \eqref{dist_prop2} the varying table size Chinese restaurant representation because the size-biased weights can be thought as the sizes of the tables in our restaurant. In the next subsection, we compute the complete conditionals of each random variable of interest to implement an overall Gibbs sampling MCMC scheme.  
 \begin{figure}
\begin{center}
\fbox{\includegraphics[scale=0.7]{hybrid_sampler_nipspic.pdf}}
\end{center}
\caption{Varying table size Chinese restaurant representation for observations $\left\{X_i\right\}_{i=1}^9$ }%
\end{figure}


\subsection{Complete conditionals}
Starting from equation \eqref{dist_prop2}, we obtain the following complete conditionals for the Gibbs sampler:
\begin{align}\label{compcond}
\mathbb{P}\left(	V\in \text{d}v\mid \text{Rest}\right)&\propto \left(v+\sum_{k=1}^K{s_k}\right)^{-n} f_\rho(v)h\left(v+\sum_{k=1}^K{s_k}\right)\text{d}v\\
\mathbb{P}\left(	\tilde{J}_i\in \text{d}s_i\mid \text{Rest}\right)&\propto  \left(v+s_i+\sum_{k\neq i}{s_k}\right)^{-n} h\left(v+s_i+\sum_{k\neq i}{s_k}\right) s_i^{|c_i|} \rho(ds_i)  \mathbb{I}_{\left(0,\text{Surpmass}_i\right)}(s_i)\text{d}s_{i}\nonumber
\end{align}
where Surpmass$_i = V+\sum_{j=1}^k\tilde{J}_j-\sum_{j<i}\tilde{J}_j $. \begin{align}
\mathbb{P}(c_i=c\mid \mathbf{c}_{-i},\textrm{Rest})&\propto \begin{cases} s_cF(d x_i\mid \left\{X_j\right\}_{j\in c} Y_c^*) &\text{if $i$ is assigned to existing cluster c}\\\frac{v}{M}F(dx_i\mid Y_c^*)&\text{if $i$ is assigned to a new cluster c}\end{cases}\nonumber
\end{align} 
According to the rule above, the $i$th observation will be either reassigned to an existing cluster or to one of the $M$ new clusters in the ReUse algorithm as in \citet{FavTeh2013a}. If it is assigned to a new cluster, then we need to sample a new size-biased weight from the following:
\begin{align}\label{newtablesizeden2}
\mathbb{P}\left(\tilde{J}_{k+1}\in\text{d}s_{k+1}\mid \text{Rest}\right)&\propto f_{\rho}(v-s_{k+1})\rho(s_{k+1})s_{k+1}\mathbb{I}_{(0,v)}(s_{k+1})\text{d}s_{k+1}.
\end{align}
Every time a new cluster is created we need to obtain its corresponding size-biased weight  which could happen $1\leq R \leq n$ times per iteration hence, it has a significant contribution to the overall computational cost. For this reason, an independent and identically distributed (i.i.d.) draw from its corresponding complete conditional \eqref{newtablesizeden2} is highly desirable. In the next subsection we present a way to achieve this.
  Finally, for updating cluster parameters $\left\{Y^*_k\right\}_{k\in[K]}$,  in the case where $H_0$ is non-conjugate to the likelihood, we use an extension of \citet{FavTeh2013a}'s ReUse algorithm, see Algorithm 2 in the supplementary material for details. 
  
  The complete conditionals in Equation \eqref{compcond} do not have a standard form but a generic MCMC method can be applied to sample from each within the Gibbs sampler. We use slice sampling from \citet{Nea2003a} to update the size-biased weights and the surplus mass. However, there is a class of priors where the total mass's density is intractable so an additional step needs to be introduced to sample the surplus mass. In the next subsection we present two alternative ways  to overcome this issue.



 %
\subsection{Example of classes of Poisson-Kingman priors}
\textbf{a)} $\mathbf{\sigma}$\textbf{-Stable Poisson-Kingman processes \citep{Pit2003a}.}  For any $\sigma\in(0,1)$, let  $f_{\sigma}(t) = \frac{1}{\pi}\sum_{j=0}^{\infty}{\frac{(-1)^{j+1}}{j!}}\textrm{sin}(\pi \sigma j)\frac{\Gamma(\sigma j +1)}{t^{\sigma j +1}}$ be the density function of a positive $\sigma$-Stable random variable and $\rho(dx)=\rho_{\sigma}(dx):=\frac{\sigma}{\Gamma(1-\sigma)}x^{-\sigma-1}dx$.
This class of RPMs is denoted by $\PK(\rho_{\sigma},h_T,H_0)$  where $h$ is a function that indexes each member of the class. For example, in the experimental section, we picked 3 choices of the $h$ function that index the following processes: Pitman-Yor, Normalized Stable and Normalized Generalized Gamma processes. This class includes all Gibbs type priors with parameter $\sigma\in(0,1)$, so other choices of $h$ are possible, see \citet{GnePit2006a} and \citet{DeBlasi2015} for a noteworthy account of this class of Bayesian nonparametric priors.
 In this case, the total mass's density is intractable and we propose two ways of dealing with this. Firstly, we used \citet{Kan1975a}'s integral representation for the $\sigma$-Stable density as in \citet{LomFavTeh2014a}, introduce an auxiliary variable $Z$ and slice sample each variable:
\begin{align}
 %
 \mathbb{P}\left(V\in \text{d}v\mid \text{Rest}\right)&\propto \left(v+\sum_{i=1}^{k}s_i\right)^{-n}v^{-\frac{\sigma}{1-\sigma}}\exp\left[-v^{\frac{-\sigma}{1-\sigma}}A(z)\right] h\left(v+\sum_{i=1}^{k}s_i\right)\text{d}v\nonumber\\
\mathbb{P}\left(	Z\in \text{d}z\mid \text{Rest}\right)&\propto A(z)\exp{\left[-v^{\left(-\frac{\sigma}{1-\sigma}\right)}A(z)\right]}\text{d}z\nonumber,
\end{align}
see Algorithm 1 in the supplementary material for details. Alternatively, we can completely bypass the evaluation of the total mass's density by updating the surplus mass with a Metropolis-Hastings step with an independent proposal from a Stable or from an Exponentially Tilted Stable($\lambda$) . It is straight forward to obtain i.i.d draws from this proposals, see \citet{devroye} and \citet{Hof11}  for an improved rejection sampling method for the Exponentially tilted case. This leads to the following acceptance ratio:%
\begin{align}%
\frac{\mathbb{P}\left(	V'\in \text{d}v'\mid \text{Rest}\right)f_{\sigma}(v)\exp{(-\lambda v)}}{\mathbb{P}\left(V\in \text{d}v\mid \text{Rest}\right)f_{\sigma}(v')\exp{(-\lambda v')}}&=\frac{\left(v'+\sum_{i=1}^{k}{s_i}\right)^{-n}
h\left(v'+\sum_{i=1}^{k}{s_i}\right)\text{d}v'
\exp{(-v)}}{\left(v+\sum_{i=1}^{k}{s_i}\right)^{-n}
h\left(v+\sum_{i=1}^{k}{s_i}\right)\text{d}v\exp{(- v')}}\nonumber,
\end{align}
see Algorithm 4 in the supplementary material for details.
Finally, to sample a new size-biased weight:
\begin{align}\label{newtablesizeden3}
\mathbb{P}\left(\tilde{J}_{k+1}\in\text{d}s_{k+1}\mid \text{Rest}\right)&\propto f_{\sigma}(v-s_{k+1})s_{k+1}^{-\sigma}\mathbb{I}_{(0,v)}(s_{k+1})\text{d}s_{k+1}.
\end{align}
Fortunately, we can get an i.i.d. draw from the above due to an identity in distribution given by \citet{FavLomNip2014a} for the usual stick breaking weights for any prior in this class such that $\sigma=\frac{u}{v}$ where $u<v$ are coprime integers. Then we just reparameterize it back to obtain the new size-biased weight, see Algorithm 3 in the supplementary material for details.

\medskip

\textbf{b)}  \textbf{$\mathbf{-\log}$Beta-Poisson-Kingman processes \citep{RegLijPru2003a,Ren08}.} %
  Let $f_{\rho}(t)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\exp{(-at)}\left(1-\exp(-t)\right)^{b-1}$ be the density of a positive random variable $X\stackrel{\text{d}}{=}-\log Y$, where $Y\sim \text{Beta}(a,b)$ and $\rho(x)=\frac{\exp(-ax)\left(1-\exp(-bx)\right)}{x\left(1-\exp(-x)\right)}$.
  This class of RPMs generalises the Gamma process but has similar properties. Indeed, if we take $b=1$ and the density function for $T$ is $\gamma(t)=f_{\rho}(t)$ we recover the L\'evy measure and total mass's density function of a Gamma process.
 Finally, to sample a new size-biased weight:
\begin{align}\label{newtablesizeden4}
\mathbb{P}\left(\tilde{J}_{k+1}\in\text{d}s_{k+1}\mid \text{Rest}\right)
&\propto\frac{\left(1-\exp(s_{k+1}-v)\right)^{b-1}\left(1-\exp(-b s_{k+1})\right)}{1-\exp(-s_{k+1})}\textrm{d}s_{k+1}\mathbb{I}_{(0,v)}(s_{k+1})\nonumber
\end{align}
If $b>1$, this complete conditional is a monotone decreasing  unnormalised density with maximum at $b$. We can easily get an i.i.d. draw with a simple rejection sampler \citep{Dev1986a}  where the rejection constant is $bv$ and the proposal is $U(0,v)$. There is no other known sampler for this process.
\subsection{Relationship to marginal and conditional MCMC samplers}
Starting from equation \eqref{dist_prop}, another strategy would be to reparameterize the model in terms of the usual stick breaking weights. Next, we could choose a random truncation level and represent finitely many sticks as in \citet{favaroslice12}. Alternatively, we could integrate out the random probability measure and sample only the partition induced by it as in \citet{LomFavTeh2014a}. Conditional samplers have large memory requirements as often, the number of sticks needed can be very large. Furthermore, the conditional distributions of the stick lengths are quite involved so they tend to have slow running times. Marginal samplers have less storage requirements than conditional samplers but could potentially have worst mixing properties. For example, \citet{LomFavTeh2014a} had to introduce a number of auxiliary variables which worsen the mixing.

 Our novel hybrid sampler exploits marginal and conditional samplers advantages. It has less memory requirements since it just represents the size-biased weights of occupied as opposed to conditional samplers which represent both empty and occupied clusters. Also, it does not integrate out the size-biased weights thus, we obtain a more comprehensive representation of the RPM.

\section{Performance assesssment}

We illustrate the performance of our hybrid sampler on a range of Bayesian nonparametric mixture models, obtained by different
specifications of $\rho$ and $\gamma$, as in Equation \eqref{infmixturemodel}. At the top level of this hierarchical specification, different Bayesian nonparametric priors were chosen from both classes presented in the examples section. We chose the base distribution $H_0$ and the likelihood term $F$ for the $k$th cluster to be

$H_0(\textrm{d}\mu_k) = \mathcal{N}\left(\textrm{d}\mu_k\mid\mu_0,\sigma_0^2\right)$ and
$F(\textrm{d}x_1,\hdots, \textrm{d}x_{n_k}\mid \mu_k,\tau_1) = \prod_{i=1}^{n_k}\mathcal{N}\left(x_i\mid \mu_k, \sigma_1^2\right),$


  %
  

  
    \begin{table}
     \centering
  \scalebox{0.8}{
 \begin{tabular}{|cccc|}
 \hline
\textbf{Algorithm}&$\mathbf{\sigma}$ &\textbf{Running time}&\textbf{ESS($\pm \textrm{std}$)}\\
 \hline
 Pitman-Yor process ($ \theta = 10$)&&&\\
 \hline
Hybrid&0.3&7135.1(28.316)&\textbf{ 2635.488(187.335)}\\
Hybrid-MH ($\lambda=0$)&0.3&5469.4(186.066)&2015.625(152.030)\\
Conditional&0.3&NA&NA\\
Marginal&0.3&\textbf{4685.7(84.104)}&2382.799(169.359)\\
\hline
Hybrid&0.5&\textbf{3246.9(24.894)}&\textbf{3595.508(174.075)}\\ 
Hybrid-MH ($\lambda=50$)&0.5&4902.3(6.936)&3579.686(135.726)\\
Conditional&0.5&10141.6(237.735)& 905.444(41.475)\\
Marginal&0.5&4757.2(37.077)&2944.065(195.011)\\
\hline
Normalized Stable process&&&\\
 \hline
Hybrid&0.3&\textbf{5054.7(70.675)}&\textbf{5324.146(167.843)}\\
Hybrid-MH ($\lambda=0$)&0.3&7866.4(803.228)&5074.909(100.300)\\
Conditional&0.3&NA&NA\\
Marginal&0.3&7658.3(193.773)&2630.264(429.877)\\
\hline
Hybrid&0.5&5382.9(57.561)&\textbf{4877.378(469.794)}\\
Hybrid-MH ($\lambda=50$)&0.5&\textbf{4537.2(37.292)}&4454.999(348.356)\\
Conditional&0.5&10033.1(22.647)&912.382(167.089)\\
Marginal&0.5&8203.1(106.798)&3139.412(351.788)\\
\hline
Normalized Generalized Gamma process ($ \tau = 1$)&&&\\
 \hline
Hybrid&0.3&\textbf{4157.8(92.863)}&\textbf{5104.713(200.949)}\\
Hybrid-MH ($\lambda=0$)&0.3&4745.5(187.506)&4848.560(312.820)\\
Conditional&0.3&NA&NA\\
Marginal&0.3&7685.8(208.98)&3587.733(569.984)\\
\hline
Hybrid&0.5&6299.2(102.853)&\textbf{4646.987(370.955)}\\
Hybrid-MH ($\lambda=50$)&0.5&\textbf{4686.4(35.661)}&4343.555(173.113)\\
Conditional&0.5&10046.9(206.538)&1000.214(70.148)\\
Marginal&0.5&8055.6(93.164)&4443.905(367.297)\\
\hline
-logBeta ($a=1,b=2$)&&&\\
 \hline
 Hybrid& -&\textbf{2520.6(121.044)}&\textbf{3068.174(540.111)}\\
 %
 Conditional&-&\text{NA}&\text{NA}\\
Marginal&-&\text{NA}&\text{NA}\\
 \hline
 \end{tabular}}
  \caption{\footnotesize{Running times  in seconds and ESS averaged over 10 chains, 30,000 iterations, 10,000 burn in. \label{assess}}}
\end{table} 


where $\left\{X_j\right\}_{j=1}^{n_k}$ are the $n_k$ observations assigned to the $k$th cluster at some iteration. $\mathcal{N}$ denotes a Normal distribution with mean $\mu_k$ and variance $\sigma_1^2$, a common parameter among all clusters. The mean's prior distribution is Normal, centered at $\mu_0$ and with variance $\sigma_0^2$. Although the base distribution is conjugate to the likelihood  we treated it as non-conjugate case and sampled the parameters at each iteration rather than integrating them out.


We used the dataset from  \citet{galaxy} to test the algorithmic performance in terms of running time and effective sample size (ESS), as Table~\ref{assess} shows. The dataset consists of measurements of velocities in km/sec of $n=82$ galaxies from a survey of the Corona Borealis region. For the $\sigma$-Stable Poisson-Kingman class, we compared it against a variation of \citet{favaroslice12}'s conditional sampler and against the marginal sampler of \citet{LomFavTeh2014a}.We chose to compare our hybrid sampler against this existing approaches which follow the same general purpose paradigm.


Table \ref{assess} shows that different choices of $\sigma$ result in differences in the algorithm's running times and ESS. The reason for this is that in the $\sigma=0.5$ case there are readily available random number generators which do not increase the computational cost. In contrast, in the $\sigma=0.3$ case, a rejection sampler method is needed every time a new size-biased weight is sampled which increases the computational cost, see \citet{FavLomNip2014a} for details.  Even so, in most cases, we outperform both marginal and conditional MCMC schemes in terms of running times and in all cases, in terms of ESS. In the Hybrid-MH case, even thought the ESS and running times are competitive, we found that the acceptance rate is not optimal, we are currently exploring other choices of proposals. Finally, in Example b), our approach is the only one available and it has good running times and ESS. This qualitative comparison confirms our previous statements about our novel approach.




\section{Discussion}

Our main contribution is our Hybrid MCMC sampler as a general purpose tool for inference with a very large class of infinite mixture models. We argue in favour of an approach in which a generic algorithm can be applied to a very large class of models, so that the modeller has a lot of flexibility in choosing specific models suitable for his/her problem of interest. Our method is a hybrid approach since it combines the perks of the conditional and marginal schemes. Indeed, our experiments confirm that our hybrid sampler is more efficient since it outperforms both marginal and conditional samplers in running times in most cases and in ESS in all cases.


 We introduced a new compact way of representing the infinite dimensional component of the model such that it is feasible to perform inference and dealing with the corresponding intractabilities. However, there are still various intractabilities and challenges  that remain when dealing with this type of models.  For example, we would like to stress that there are some values for $\sigma$ where we are unable to perform inference with our novel sampler. Furthermore, there could be other ways to improve the mixing when a Metropolis-Hastings step is chosen in terms of better proposals. We consider these points to be an interesting avenue of future research.



\subsubsection*{Acknowledgments}
We thank Konstantina Palla for her insightful comments. Mar\'ia Lomel\'i is funded by the Gatsby Charitable Foundation, Stefano Favaro is supported by the European Research Council through StG N-BNP 306406 and Yee Whye Teh is supported by the European Research Council under the European Union�s Seventh Framework Programme (FP7/2007-2013) ERC grant agreement no. 617411. 



\small{
\bibliographystyle{authordate1}
\bibliography{refs_nipspkpaper}
}
\newpage
\appendix
\section{Pseudocode}


\begin{algorithm}
\caption{HybridSampler$\sigma$-PK$\left(K,V,\mathbf{c},\{X_i\}_{i\in[n]},\{Y^*_c\}_{c\in\Pi_n},H_0, M\right)$ }
\begin{algorithmic}
\For{$t = 2 \to iter$}
\State Update $v^{(t)}$: Slice sample $\tilde{\mathbb{P}}\left(V\in \textrm{d}v\mid  \textrm{rest}\right)$
\State Update $s_i^{(t)}$ for $i=1,\hdots,k$ : Slice sample $\tilde{\mathbb{P}}\left( \tilde{J}_i\in \textrm{d}s_i \mid  \textrm{rest}\right)$
\State Update \small{$\pi^{(t)}, \left\{y_c^*\right\}^{(t)} _{c\in \pi}, \left\{s_c\right\}^{(t)} _{c\in \pi}$: \textbf{AddTable\&ReUse}$\left(V,\Pi_n,M,\{X_i\}_{i\in[n]},\{Y^*_c\}_{c\in\Pi_n},\{\tilde{J}_c\}_{c\in\Pi_n},H_0\mid \textrm{rest}\right)$}
\EndFor 
\end{algorithmic}
\end{algorithm}


\begin{algorithm}
\caption{AddTable\&ReUse$\left(V,\Pi_n,M,\{X_i\}_{i\in[n]},\{Y^*_c\}_{c\in\Pi_n},\{S_c\}_{c\in\Pi_n},H_0\mid \textrm{rest}\right)$}
\begin{algorithmic}
\State Let $c\in \Pi_n$ be such that $i\in c$
\State $c\leftarrow c\setminus\{i\}$
\If {$c=\emptyset$}
\State $k\sim\textrm{UniformDiscrete}(\frac{1}{M})$
\State $Y_k^e \leftarrow Y^*_c$
\State $\Pi_n\leftarrow \Pi_n\setminus \{c\}$
 \State $V \leftarrow  V+ \tilde{J}_c$\Comment{Add back the discarded table size to the surplus}
\EndIf
\State Set  $c'$according to$\mathbb{P}(c_i=c\mid \mathbf{c}_{-i},\textrm{Rest})\propto\begin{cases} \tilde{J}_cF(x_i\mid \left\{X_i\right\}_{i\in c} Y_c^*) &\text{if existing}\\\frac{V}{M}F(x_i\mid Y_c^*)&\text{if new}\end{cases}$ 
\If{$c'\in[M]$ }
\State $ \tilde{J}_{\text{new}}\leftarrow  $\textbf{ExactSampleNewTableSize}$(\{ \tilde{J}_c\}_{c\in\pi},V,\Pi_n=\pi, \textrm{Rest})$%
\State $V \leftarrow  V- \tilde{J}_{\text{new}}$\Comment{Remove it from the old surplus}
\State $\Pi_n\leftarrow \Pi_n\cup \{\{i\}\}$
\State $Y^*_{\{i\}}\leftarrow Y_{c'}^e$ 
\State $Y_{c'}^e\sim H_0$ 
\Else
\State $c' \leftarrow c' \cup \{i\}$ 
\EndIf
\State Draw $\{Y_j^e\}_{j=1}^M \stackrel{\textrm{i.i.d.}}{\sim} H_0$
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{ExactSampleNewTableSize$(V,\sigma, \text{Rest})$ }
\begin{algorithmic}
\State
\If{$\sigma =0.5$}
\State $G\sim\text{Gamma}\left(\frac{3}{4},1\right)$
\State $IG\sim\text{Inverse Gamma}\left(\frac{1}{4},\frac{1}{4^3}V^{-2}\right)$
\State $V_{stick} = \frac{\sqrt{G}}{\sqrt{G}+\sqrt{IG}}$
\State $ \tilde{J}_{new}=V_{stick}V$
\Else
\If{$\sigma<0.5 \hspace{2mm}\& \& \hspace{2mm}\sigma=\frac{u_{\sigma}}{v_{\sigma}} ,u_{\sigma},v_{\sigma}\in \mathbb{N}$}
\State $\lambda=u_{\sigma}^2/ v_{\sigma}^{\frac{v_{\sigma}}{u_{\sigma}}}$
\State $IG\sim\text{Inverse Gamma}\left(1-\frac{u_{\sigma}}{v_{\sigma}},\lambda\right)$
\State $\frac{1}{G}\sim\mathcal{E_T}\left(\lambda,L^{-1/u}_{\frac{u_{\sigma}}{v_{\sigma}}}\right)$
\State $V_{stick} = \frac{G}{G+IG}$
\State $ \tilde{J}_{new}=V_{stick}V$
\EndIf
\EndIf
\end{algorithmic}
\end{algorithm}


\begin{algorithm}
\caption{HybridSampler-MH-$\sigma$PK$\left(K,\mathbf{S},V,\mathbf{c},\mathbf{\theta},\mathbf{Y},\mathbf{X},M\right)$ }
\begin{algorithmic}
\For{$t = 2 \to iter$}
\State Update $s_i^{(t)}$ for $i=1,\hdots,k$ : Slice sample $\tilde{\mathbb{P}}\left( \tilde{J}_i\in \textrm{d}s_i \mid  \textrm{rest}\right)$
\State Update $v^{(t)}$: MH step for $\tilde{\mathbb{P}}\left(V\in \textrm{d}v\mid  \textrm{rest}\right)$ with independent proposal \texttt{Stablernd}($\sigma$) or  \texttt{ExpTiltStablernd}($\lambda,\sigma$) . %
\State  Update \small{$\pi^{(t)}, \left\{y_c^*\right\}^{(t)} _{c\in \pi}, \left\{s_c\right\}^{(t)} _{c\in \pi}$: \textbf{AddTable\&ReUse}$\left(V,\Pi_n,M,\{X_i\}_{i\in[n]},\{Y^*_c\}_{c\in\Pi_n},\{\tilde{J}_c\}_{c\in\Pi_n},H_0\mid \textrm{rest}\right)$}
\EndFor 
\end{algorithmic}
\end{algorithm}



\end{document}

NIPSPKpaper_resubmission_Arxivready_24sept2015.bbl
\begin{thebibliography}{}

\bibitem[\protect\citename{De~Blasi {\em et~al.\ }\relax, }2015]{DeBlasi2015}
De~Blasi, P., Favaro, S., Lijoi, A., Mena, R.~H., Pr\"uenster, I., \& Ruggiero,
  M. 2015.
\newblock Are Gibbs-type priors the most natural generalization of the
  {D}irichlet process?
\newblock {\em Pages  212--229 of:} {\em IEEE Transactions on Pattern Analysis
  \& Machine Intelligence},  vol. 37.

\bibitem[\protect\citename{Devroye, }1986]{Dev1986a}
Devroye, L. 1986.
\newblock {\em Non-Uniform Random Variate Generation}.
\newblock Springer-Verlag.

\bibitem[\protect\citename{Devroye, }2009]{devroye}
Devroye, L. 2009.
\newblock Random variate generation for exponentially and polynomially tilted
  {S}table distributions.
\newblock {\em ACM Transactions on Modelling and Computer Simulation}, {\bf
  19}, 1--20.

\bibitem[\protect\citename{Escobar, }1994]{Esc1994a}
Escobar, M.~D. 1994.
\newblock Estimating normal means with a {D}irichlet process prior.
\newblock {\em Journal of the American Statistical Association}, {\bf 89},
  268--277.

\bibitem[\protect\citename{Escobar \& West, }1995]{EscWes1995a}
Escobar, M.~D., \& West, M. 1995.
\newblock {B}ayesian density estimation and inference using mixtures.
\newblock {\em Journal of the American Statistical Association}, {\bf 90},
  577--588.

\bibitem[\protect\citename{Favaro \& Teh, }2013]{FavTeh2013a}
Favaro, S., \& Teh, Y.~W. 2013.
\newblock {MCMC} for Normalized Random Measure Mixture Models.
\newblock {\em Statistical Science}, {\bf 28}(3), 335--359.

\bibitem[\protect\citename{Favaro \& Walker, }2012]{favaroslice12}
Favaro, S., \& Walker, S.~G. 2012.
\newblock Slice sampling {$\sigma$-S}table {Poisson-Kingman} mixture models.
\newblock {\em Journal of Computational and Graphical Statistics}, {\bf 22},
  830--847.

\bibitem[\protect\citename{Favaro {\em et~al.\ }\relax, }2014]{FavLomNip2014a}
Favaro, S., Lomeli, M., Nipoti, B., \& Teh, Y.~W. 2014.
\newblock On the Stick-Breaking representation of $\sigma$-{S}table
  {P}oisson-{K}ingman models.
\newblock {\em Electronic Journal of Statistics}, {\bf 8}, 1063--1085.

\bibitem[\protect\citename{Ghahramani, }2015]{Zou15}
Ghahramani, Z. 2015.
\newblock Probabilistic Machine Learning and Artificial Inteligence.
\newblock {\em Nature}, {\bf 521}, 452�459.

\bibitem[\protect\citename{Gnedin \& Pitman, }2006]{GnePit2006a}
Gnedin, A., \& Pitman, J. 2006.
\newblock Exchangeable {G}ibbs partitions and {S}tirling triangles.
\newblock {\em Journal of Mathematical Sciences}, {\bf 138}, 5674--5684.

\bibitem[\protect\citename{Hofert, }2011]{Hof11}
Hofert, M. 2011.
\newblock Efficiently sampling nested Archimedean copulas.
\newblock {\em Comput. Statist. Data Anal.}, {\bf 55}, 57�70.

\bibitem[\protect\citename{Ishwaran \& James, }2001]{IshJam2001a}
Ishwaran, H., \& James, L.~F. 2001.
\newblock {G}ibbs Sampling Methods for Stick-Breaking Priors.
\newblock {\em Journal of the American Statistical Association}, {\bf 96}(453),
  161--173.

\bibitem[\protect\citename{James, }2002]{Jam2002a}
James, L.~F. 2002.
\newblock {\em Poisson process partition calculus with applications to
  exchangeable models and {B}ayesian nonparametrics}.
\newblock ArXiv:math/0205093.

\bibitem[\protect\citename{Kanter, }1975]{Kan1975a}
Kanter, M. 1975.
\newblock Stable densities under change of scale and total variation
  inequalities.
\newblock {\em Annals of Probability}, {\bf 3}, 697--707.

\bibitem[\protect\citename{Kingman, }1967]{Kin1967a}
Kingman, J. F.~C. 1967.
\newblock Completely Random Measures.
\newblock {\em Pacific Journal of Mathematics}, {\bf 21}, 59--78.

\bibitem[\protect\citename{Kingman, }1978]{Kin1978a}
Kingman, J. F.~C. 1978.
\newblock The representation of partition structures.
\newblock {\em Journal of the London Mathematical Society}, {\bf 18}, 374--380.

\bibitem[\protect\citename{Lomeli {\em et~al.\ }\relax, }2014]{LomFavTeh2014a}
Lomeli, M., Favaro, S., \& Teh, Y.~W. 2014.
\newblock A marginal sampler for {$\sigma$-Stable P}oisson-{K}ingman mixture
  models.
\newblock {\em Arxiv 1407.4211}.

\bibitem[\protect\citename{Neal, }1998]{Nea1998b}
Neal, R.~M. 1998.
\newblock {\em {M}arkov Chain Sampling Methods for {D}irichlet Process Mixture
  Models}.
\newblock Tech. rept. 9815. Department of Statistics, University of Toronto.

\bibitem[\protect\citename{Neal, }2003]{Nea2003a}
Neal, R.~M. 2003.
\newblock Slice sampling.
\newblock {\em Annals of Statistics}, {\bf 31}, 705--767.

\bibitem[\protect\citename{Papaspiliopoulos \& Roberts, }2008]{PapRob2008a}
Papaspiliopoulos, O., \& Roberts, G.~O. 2008.
\newblock {Retrospective Markov chain Monte Carlo methods for Dirichlet process
  hierarchical models}.
\newblock {\em Biometrika}, {\bf 95}, 169--186.

\bibitem[\protect\citename{Perman {\em et~al.\ }\relax, }1992]{PerPitYor1992a}
Perman, M., Pitman, J., \& Yor, M. 1992.
\newblock Size-biased sampling of {P}oisson point processes and excursions.
\newblock {\em Probability Theory and Related Fields}, {\bf 92}, 21--39.

\bibitem[\protect\citename{Pitman, }1996]{Pit96}
Pitman, J. 1996.
\newblock Random discrete distributions invariant under size-biased
  permutation.
\newblock {\em Advances in Applied Probability}, {\bf 28}, 525--539.

\bibitem[\protect\citename{Pitman, }2003]{Pit2003a}
Pitman, J. 2003.
\newblock {P}oisson-{K}ingman Partitions.
\newblock {\em Pages  1--34 of:} Goldstein, D.~R. (ed), {\em Statistics and
  Science: a Festschrift for Terry Speed}.
\newblock Institute of Mathematical Statistics.

\bibitem[\protect\citename{Pitman, }2006]{Pit2006a}
Pitman, J. 2006.
\newblock {\em Combinatorial Stochastic Processes}.
\newblock Lecture Notes in Mathematics.
\newblock Springer-Verlag, Berlin.

\bibitem[\protect\citename{Regazzini {\em et~al.\ }\relax,
  }2003]{RegLijPru2003a}
Regazzini, E., Lijoi, A., \& Pr\"uenster, I. 2003.
\newblock Distributional results for means of normalized random measures with
  independent increments.
\newblock {\em Annals of Statistics}, {\bf 31}, 560--585.

\bibitem[\protect\citename{Roeder, }1990]{galaxy}
Roeder, K. 1990.
\newblock Density estimation with confidence sets exemplified by super-clusters
  and voids in the galaxies.
\newblock {\em Journal of the American Statistical Association}, {\bf 85},
  617--624.

\bibitem[\protect\citename{von Renesse {\em et~al.\ }\relax, }2008]{Ren08}
von Renesse, M., Yor, M., \& Zambotti, L. 2008.
\newblock Quasi-invariance properties of a class of subordinators.
\newblock {\em Stochastic Processes and their Applications}, {\bf 118},
  2038--2057.

\bibitem[\protect\citename{Walker, }2007]{Wal2007a}
Walker, Stephen~G. 2007.
\newblock Sampling the {D}irichlet Mixture Model with Slices.
\newblock {\em Communications in Statistics - Simulation and Computation}, {\bf
  36}, 45.

\end{thebibliography}
