###### Abstract

The Bayesian approach to inference stands out for naturally allowing borrowing
information across heterogeneous populations, with different samples possibly
sharing the same distribution. A popular Bayesian nonparametric model for
clustering probability distributions is the nested Dirichlet process, which
however has the drawback of grouping distributions in a single cluster when
ties are observed across samples. With the goal of achieving a flexible and
effective clustering method for both samples and observations, we investigate
a nonparametric prior that arises as the composition of two different discrete
random structures and derive a closed-form expression for the induced
distribution of the random partition, the fundamental tool regulating the
clustering behavior of the model. On the one hand, this allows to gain a
deeper insight into the theoretical properties of the model and, on the other
hand, it yields an MCMC algorithm for evaluating Bayesian inferences of
interest. Moreover, we single out limitations of this algorithm when working
with more than two populations and, consequently, devise an alternative more
efficient sampling scheme, which as a by-product, allows testing homogeneity
between different populations. Finally, we perform a comparison with the
nested Dirichlet process and provide illustrative examples of both synthetic
and real data.

toltxlabel *Supp_au

Flexible clustering via hidden hierarchical Dirichlet priors

Antonio Lijoi 1, Igor Prünster 1, Giovanni Rebaudo2

1 Department of Decision Sciences and BIDSA, Bocconi University, via Röntgen
1, 20136 Milan, Italy

2 Department of Statistics and Data Sciences, University of Texas at Austin,
TX 78712-1823, USA

Key Words: Bayesian nonparametrics, clustering, dependent random partitions,
hierarchical Dirichlet process, mixture models, nested Dirichlet process,
vectors of random probabilities

##  1 Introduction

Dirichlet process (DP) mixtures are well-established and highly successful
Bayesian nonparametric models for density estimation and clustering, which
also enjoy appealing frequentist asymptotic properties (Lo,, 1984; Escobar,,
1994; Escobar and West,, 1995; Ghosal and van der Vaart,, 2017). However, they
are not suitable to model data
$\\{(X_{j,1},\ldots,X_{j,I_{j}}):j=1,\ldots,J\\}$ that are recorded under $J$
different, though related, experimental conditions. This is due to
exchangeability implying a common underlying distribution across populations,
a homogeneity assumption which is clearly too restrictive. To make things
concrete we consider the Collaborative Perinatal Project, which is a large
prospective epidemiologic study conducted from 1959 to 1974 (analyzed in
Section 5.3), where pregnant women were enrolled in 12 hospitals and followed
over time. Using a standard DP mixture on the patients enrolled across all 12
hospitals would correspond to ignoring the information on the specific center
$j$ where the data are collected and, thus, the heterogeneity across samples.
The opposite, also unrealistic, extreme case corresponds to modeling data from
each hospital independently, thus ignoring possible similarities among them.

A natural compromise between the aforementioned extreme cases is _partial
exchangeability_ (de Finetti,, 1938), which entails exchangeability within
each experimental condition (but not across) and dependent population–specific
distributions (thus allowing borrowing of information). See Kallenberg, (2005)
for a detailed account of the topic. In this framework the proposal of
dependent versions of the DP date back to the seminal papers of Cifarelli and
Regazzini, (1978) and MacEachern, (1999, 2000). Dependent DPs can be readily
used within mixtures leading to several success stories in topic modeling,
biostatistics, speaker diarization, genetics, fMRI analysis, and so forth. See
Dunson, (2010); Teh and Jordan, (2010); Foti and Williamson, (2015); Quintana
et al., (2022) and references therein.

Two hugely popular dependent nonparametric priors, which will also represent
the key ingredients of the present contribution, are the hierarchical
Dirichlet process (HDP) (Teh et al.,, 2006) and the nested Dirichlet process
(NDP) (Rodríguez et al.,, 2008). The HDP clusters observations within and
across populations. The NDP aims to cluster both population distributions and
observations, but as shown in Camerlenghi et al., 2019b , does not achieve
this goal. In fact, if there is a cluster of observations shared by different
samples, the model degenerates to exchangeability across samples. This issue
is successfully overcome in Camerlenghi et al., 2019b  by introducing latent
nested nonparametric priors. However, while this proposal has the merit of
being the first to solve the degeneracy problem, it suffers from other
limitations in terms of implementation and modeling: (a) with data from more
than two populations the analytical and computational burden implied by the
additive structure becomes overwhelming; (b) the model lacks the flexibility
needed to capture different weights that common clusters may feature across
different populations. More details can be found in the discussion to
Camerlenghi et al., 2019b .

The goal of this paper is thus to devise a principled Bayesian nonparametric
approach, which allows to cluster simultaneously distributions and
observations (within and across populations). We achieve this by blending
peculiar features of both the NDP and the HDP into a model, which we term
Hidden Hierarchical Dirichlet Process (HHDP). Importantly, the HHDP overcomes
the above-mentioned theoretical, modeling, and computational limitations since
it, respectively, does not suffer from the degeneracy flaw, is able to
effectively capture different weights of shared clusters and allows to handle
several populations as showcased in the real data application. Note that the
idea of the model was first hinted at in James, (2008) and, later, considered
in Agrawal et al., (2013) from a mere computational point of view without
providing results on distributional properties that are relevant for Bayesian
inference. Hence, as a by-product, our theoretical results shed also some
light on the topic modeling applications of Agrawal et al., (2013).
Additionally, the same model was independently applied in Balocchi et al.,
(2021) to successfully cluster urban areal units at different levels of
resolution simultaneously.

Section 2 concisely reviews the HDP and the NDP with a focus on the random
partitions they induce. In Section 3 we define the HHDP and investigate its
properties, foremost its clustering structure (induced by a partially
exchangeable array of observations). These findings lead to the development of
marginal and conditional Gibbs sampling schemes in Section 4. In Section 5 we
draw a comparison between HHDP and NDP on synthetic data and present a real
data application for our model. Finally, Section 6 is devoted to some
concluding remarks and possible future research.

##  2 Bayesian nonparametric priors for clustering

The assumption of exchangeability that characterizes widely used Bayesian
inferential procedures is equivalent to assuming data homogeneity. This is not
realistic in many applied contexts, for instance, for data recorded under $J$
different experimental conditions inducing heterogeneity. A natural assumption
that relaxes exchangeability and is suited for arrays of random variables
$\\{(X_{j,i})_{i\geq 1}:j=1,\ldots,J\\}$ is partial exchangeability, which
amounts to assuming homogeneity within each population, though not across
different populations. This is characterized by

$\\{(X_{j,i})_{i\geq
1}:j=1,\ldots,J\\}\overset{\text{d}}{=}\\{(X_{j,\sigma_{j}(i)})_{i\geq
1}:j=1,\ldots,J\\},$

for every finitary permutation $\\{\sigma_{j}:j=1,\ldots,J\\}$ with
$\overset{\text{d}}{=}$ henceforth denoting equality in distribution. Thanks
to de Finetti’s representation theorem for partially exchangeable arrays, the
dependence structure is effectively represented through the following
hierarchical formulation

$\begin{split}X_{j,i}\mid(G_{1},\ldots,G_{J})&\overset{\text{ind}}{\sim}G_{j},\quad\quad(i=1,\ldots,I_{j},j=1,\ldots,J)\\\
(G_{1},\ldots,G_{J})&\sim\mathcal{L}.\end{split}$ (1)

Here we focus on priors $\mathcal{L}$ defined as compositions of discrete
random structures and including, as special cases, both the HDP and the NDP.
More specifically, we consider $\mathcal{L}$ in (1) that is defined as follows

$G_{j}\,|\,Q\,\overset{\text{iid}}{\sim}\,\mathcal{L}(G_{j}|\,Q)\quad(j=1,\ldots,J);\qquad
Q\,|\,G_{0}\,\sim\,\mathcal{L}(Q|\,G_{0});\qquad
G_{0}\,\sim\,\mathcal{L}(G_{0}),$ (2)

with discrete random probability measures $G_{j}$ ($j=1,\ldots,J$), $Q$ and
$G_{0}$. The data are denoted by $\bm{X}=\\{\bm{X}_{1},\ldots,\bm{X}_{J}\\}$
with $\bm{X}_{j}=(X_{j,1},\ldots,X_{j,I_{j}})$ and $I_{j}$ the size of the
$j$th sample. Discreteness of these random structures entails that with
positive probability there are ties within each sample $\bm{X}_{j}$ and also
across samples $j=1,\ldots,J$, i.e. $\mathbb{P}(X_{j,i}=X_{j,\ell})>0$ for any
$i\neq\ell$, and $\mathbb{P}(X_{j,i}=X_{\kappa,\ell})>0$ for any
$j\neq\kappa$. Hence, $\bm{X}$ induces a random partition of the integers
$\\{1,2,\ldots,n\\}$ with $n=I_{1}+\,\cdots\,+I_{J}$, whose distribution
encapsulates the whole probabilistic clustering of the model and is,
therefore, the key quantity to study. Importantly, the random partition can be
characterized in terms of the partially exchangeable partition probability
function (pEPPF) as defined in Camerlenghi et al., 2019a . The pEPPF is the
natural generalization of the concept of exchangeable partition probability
function (EPPF) for the exchangeable case (see e.g. Pitman,, 2006). More
precisely, $D$ is the number of distinct values among the
$n=\sum_{j=1}^{J}I_{j}$ observations in the overall sample $\bm{X}$. The
vector of frequency counts is denoted by $\bm{n}_{j}=(n_{j,1},\ldots,n_{j,D})$
with $n_{j,d}$ indicating the number of elements in the $j$th sample that
coincide with the $d$th distinct value in order of arrival. Clearly,
$n_{j,d}\geq 0$ and $\sum_{i=1}^{J}n_{i,d}\geq 1$. One may well have
$n_{j,d}=0$, which implies that the $d$th distinct value is not recorded in
the $j$th sample, though by virtue of $\sum_{i=1}^{J}n_{i,d}\geq 1$ it must be
recorded at least in one of the samples. The $d$th distinct value is shared by
any two samples $j$ and $j^{\prime}$ if and only if
$n_{j,d}\,n_{j^{\prime},d}\geq 1$. The probability law of the random partition
is characterized by the pEPPF defined as

$\Pi_{D}^{(n)}(\bm{n}_{1},\ldots,\bm{n}_{J})=\mathbb{E}\int_{\mathbb{X}^{D}_{\ast}}\prod_{d=1}^{D}\\{G_{1}(\mathrm{d}x_{d})\\}^{n_{1,d}}\ldots\\{G_{J}(\mathrm{d}x_{d})\\}^{n_{J,d}},$
(3)

with the constraint $\sum_{d=1}^{D}n_{j,d}=I_{j}$, for each $j=1,\ldots,J$ and
where $\mathbb{X}$ is the space in which the $X_{j,i}$’s take values and
$\mathbb{X}^{D}_{\ast}$ is the collection of vectors in $\mathbb{X}^{D}$ whose
entries are all distinct. We stress that the expected value in (3) is computed
with respect to the joint law of the vector of random probabilities
$(G_{1},\ldots,G_{J})$, that is the de Finetti measure $\mathcal{L}$ in (1).
Hence, the pEPPF may also be interpreted as a marginal likelihood when
$(G_{1},\ldots,G_{J})$ directly model the observations according to (1).
Obviously, for a single population, that is $J=1$, the standard EPPF is
recovered and (3) is further interpretable as an extension of a product
partition model to a multiple samples framework. As such, it provides an
alternative approach to popular covariate–dependent product partition models.
See, e.g., Müller et al., (2011), Page and Quintana, (2016) and Page and
Quintana, (2018).

If we specify $\mathcal{L}(\,\cdot\,|Q)$ and $Q$ such that they give rise to
an NDP, then one may have ties also among the population probability
distributions $G_{1},\ldots,G_{J}$, i.e. $\mathbb{P}(G_{j}=G_{\kappa})>0$ for
any $j\neq\kappa$. Therefore, in the framework of (1) and (2), one may
investigate two types of clustering: (i) distributional clustering, which is
related to $G_{1},\ldots,G_{J}$ and (ii) observational clustering, which
refers to $\bm{X}$. The composition of these two clustering structures is the
main tool we rely on to devise a simple, yet effective, model that
considerably improves over existing alternatives.

###  2.1 Hierarchical Dirichlet process

Probably the most popular nonparametric prior for the partially exchangeable
case is the HDP of Teh et al., (2006), which can be nicely framed in the
composition scheme (2) as

$\mathcal{L}(G_{j}|Q)=\textsc{DP}(G_{j}|\beta,Q),\quad\mathcal{L}(Q|G_{0})=\delta_{G_{0}}(Q),\quad\mathcal{L}(G_{0})=\textsc{DP}(G_{0}|\beta_{0};H),$
(4)

where $\textsc{DP}(\,\cdot\,|\alpha,P)$ denotes the law of a DP with
concentration parameter $\alpha>0$ and baseline probability measure $P$. Here
we assume that $H$ is a non–atomic probability measure on $\mathbb{X}$ and we
refer to such prior as the $J$-dimensional HDP denoted by
$(G_{1},\ldots,G_{J})\sim\textsc{HDP}(\beta,\beta_{0};H)$. Hence, the
$G_{j}$’s share the atoms through $G_{0}$ and this leads to the creation of
shared clusters of observations (or latent features) across the $J$ groups.
The pEPPF induced by a partially exchangeable array in (1) with
$\mathcal{L}=\textsc{HDP}(\beta,\beta_{0};H)$ has been determined in
Camerlenghi et al., 2019a . It is important to stress that the model is not
suited for comparing populations’ distributions since
$\mathbb{P}(G_{j}=G_{\kappa})=0$ for any $j\neq\kappa$ (unless the $G_{j}$’s
are degenerate at $G_{0}$, in which case all distributions are equal). Similar
compositions have been considered in Camerlenghi et al., 2019a  and, later, in
Argiento et al., (2020) and Bassetti et al., (2020). Hierarchically dependent
mixture hazards have been introduced in Camerlenghi et al., (2021). Anyhow,
the HDP and its variations cannot be used to cluster both populations and
observations. To achieve this, one has to rely on priors induced by nested
structures, the most popular being the NDP.

###  2.2 Nested Dirichlet process

The NDP, introduced by Rodríguez et al., (2008), is the most widely used
nonparametric prior allowing to cluster both observations and populations.
However, as proved in Camerlenghi et al., 2019b , it suffers from a
_degeneracy issue_ , because even a single tie shared across samples is enough
to group the $J$ population distributions into a single cluster.

Like the HDP, also the NDP can be framed in the composition structure (2) as

$\mathcal{L}(G_{j}|Q)=Q(G_{j}),\quad\mathcal{L}(Q|G_{0})=\textsc{DP}(Q|\alpha;G_{0}),\quad\mathcal{L}(G_{0})=\delta_{\textsc{DP}(\beta;H)}(G_{0}),$
(5)

where $Q$ is a random probability measure on the space
$\mathscr{P}_{\mathbb{X}}$ of probability measures on $\mathbb{X}$ and $G_{0}$
is degenerate at the atom $\textsc{DP}(\beta;H)$, which is the law of a DP on
the sample space $\mathbb{X}$. As in (4), $H$ is assumed to be a non-atomic
probability measure on $\mathbb{X}$. Henceforth, we write
$(G_{1},\ldots,G_{J})\sim\textsc{NDP}(\alpha,\beta;H)$. By virtue of the
well–known stick–breaking representation of the DP (Sethuraman,, 1994) one has

$Q=\sum_{k\geq
1}\pi^{\ast}_{k}\delta_{G^{\ast}_{k}},\quad(\pi^{\ast}_{k})_{k\geq
1}\sim\textsc{GEM}(\alpha),\quad
G^{\ast}_{k}\overset{\text{iid}}{\sim}\textsc{DP}(\beta;H),$ (6)

where the weights $(\pi^{\ast}_{k})_{k\geq 1}$ and the random distributions
$(G^{\ast}_{k})_{k\geq 1}$ are independent. Recall that GEM stands for the
distribution of probability weights after Griffiths, Engen, and McCloskey,
according to the well-established terminology of Ewens, (1990). Given a
sequence $(V_{i})_{i\geq 1}$ such that
$V_{i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mbox{Beta}(1,\alpha)$,
this means that $\pi_{1}^{*}=V_{1}$ and
$\pi_{k}^{*}=V_{k}\,\prod_{i=1}^{k-1}(1-V_{i})$, for any $k\geq 2$. Since
$\mathbb{P}(G_{j}=G_{\kappa})=1/(\alpha+1)$ for any $j\neq\kappa$, $Q$
generates ties among the random distributions $G_{j}$’s with positive
probability and, thus, clusters populations. Furthermore, a structure similar
to the one displayed in (6) holds for each $G^{\ast}_{k}$, i.e.

$G^{\ast}_{k}=\sum_{l\geq
1}\omega_{k,l}\delta_{X^{\ast}_{k,l}},\quad(\omega_{k,l})_{l\geq
1}\overset{\text{iid}}{\sim}\textsc{GEM}(\beta),\quad
X^{\ast}_{k,l}\overset{\text{iid}}{\sim}H,$

and, due to the non–atomicity of $H$, the $X_{k,l}^{\ast}$ are all distinct
values.

The discrete structure of the $G_{k}^{\ast}$’s generates ties across the
samples $\\{\bm{X}_{j}:j=1,\ldots,J\\}$ with positive probability. For
example,
$\mathbb{P}(X_{j,i}=X_{j^{\prime},i^{\prime}})=1/\\{(\alpha+1)(\beta+1)\\}$
for any $j\neq j^{\prime}$. Hence, the $G_{k}^{\ast}$’s induce the clustering
of the observations $\bm{X}$.

If the data $\bm{X}$ are modeled as in (1), with
$(G_{1},\ldots,G_{J})\sim\textsc{NDP}(\alpha,\beta;H)$, conditional on a
partition of the $G_{j}$’s the observations from populations allocated to the
same cluster are exchangeable and those from populations allocated to distinct
clusters are independent. This potentially appealing feature of the NDP is
however the one responsible for the above-mentioned degeneracy issue. For
exposition clarity, consider the case of $J=2$ populations. If the two
populations belong to different clusters, i.e. $G_{1}\neq G_{2}$, they cannot
share even a single atom $X_{k,l}^{\ast}$ due to the non–atomicity of $H$.
Hence, $\mathbb{P}(X_{1,l}=X_{2,l^{\prime}}|G_{1}\neq G_{2})=0$ for any $l$
and $l^{\prime}$. Therefore there is neither clustering of observations nor
borrowing of information across different populations. On the contrary,
$\mathbb{P}(X_{1,i}=X_{2,i^{\prime}}|G_{1}=G_{2})=1/(\beta+1)>0$. These two
findings are quite intuitive. Indeed, $G_{1}\neq G_{2}$ means they are
independent realizations of a DP with atoms iid from the same non-atomic
probability distribution $H$ and, thus, they are almost surely different.
Instead, $G_{1}=G_{2}$ corresponds to all observations coming from the same
population distribution, more precisely from the same DP, and ties occur with
positive probability. A less intuitive fact is that when a single atom, say
$X_{k,l}^{\ast}$, is shared between $G_{1}$ and $G_{2}$ the model degenerates
to the exchangeable case, namely
$\mathbb{P}(G_{1}=G_{2}|X_{1,i}=X_{2,i^{\prime}})=1$ and the two populations
have (almost surely) equal distributions. Hence, the NDP is not an appropriate
specification when aiming at clustering both populations and observations
across different populations. This was shown in Camerlenghi et al., 2019b
where, spurred by this anomaly of the NDP, a novel class of priors named
latent nested processes (LNP) designed to ensure that $\mathbb{P}(G_{1}\neq
G_{2}|X_{1,i}=X_{2,i^{\prime}})>0$ is proposed. However, while this formally
solves the problem, it has computational and modeling limitations. On the one
hand, the implementation of LNPs with more than two samples is not feasible
due to severe computational hurdles. On the other hand, LNPs have limited
flexibility since the weights of the common clusters of observations across
different populations are the same. This feature is not suited to several
applications and the discussion to Camerlenghi et al., 2019b  provides
interesting examples. See also Soriano and Ma, (2019); Christensen and Ma,
(2020); Denti et al., (2021); Beraha et al., (2021) for further stimulating
contributions to this literature.

Hence, within the composition structure framework (2), our goal is to obtain a
prior distribution able to infer the clustering structure of both populations
and observations, which is highly flexible and implementable for a large
number of populations and associated samples.

##  3 Hidden hierarchical Dirichlet process

Our proposal consists in blending the HDP and the NDP in a way to leverage on
their strengths, namely clustering data across multiple heterogeneous samples
for the HDP and clustering different populations (or probability
distributions) for the NDP. More precisely we combine these two models in a
structure (2) as

$\mathcal{L}(G_{j}|Q)=Q(G_{j}),\quad\mathcal{L}(Q|G_{0})=\textsc{DP}(Q|\alpha;\textsc{DP}(\beta;G_{0})),\quad\mathcal{L}(G_{0})=\textsc{DP}(G_{0}|\beta_{0};H).$

This leads to the following definition.

######  Definition 1.

The vector of random probability measures $(G_{1},\ldots,G_{J})$ is a _hidden
hierarchical Dirichlet process_ (HHDP) if

$G_{j}\mid Q\overset{\text{iid}}{\sim}Q,\quad Q=\sum_{k\geq
1}\pi^{\ast}_{k}\delta_{G^{\ast}_{k}},\quad(\pi^{\ast}_{k})_{k\geq
1}\sim\textsc{GEM}(\alpha),\quad(G^{\ast}_{k})_{k\geq
1}\sim\textsc{HDP}(\beta,\beta_{0};H),$

with $(\pi^{\ast}_{k})_{k\geq 1}$ and $(G^{\ast}_{k})_{k\geq 1}$ independent.
In the sequel we write
$(G_{1},\ldots,G_{J})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$.

In terms of a graphical model, the HHDP can be represented as in Figure 1.

$G_{j}$$z_{j}$$\alpha$$G^{\ast}_{k}$$G_{0}$$\beta$$H$$\beta_{0}$$J$$\infty$
Figure 1: Graphical model representing the dependencies for a
$\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$. Here the $z_{j}$’s are auxiliary
integer–valued random variables that assign each $G_{j}$ to a specific atom
$G_{k}^{\ast}$ of $Q$.

The sequence $(G^{\ast}_{k})_{k\geq 1}$ acts as a hidden, or latent, component
that is crucial to avoid the bug of the NDP, namely clustering of populations
when they share some observations. Moreover, by extending (4) to $J=\infty$,
it can be more conveniently represented as

$\begin{split}G_{k}^{\ast}=\sum_{l\geq 1}\omega_{k,l}\,\delta_{Z_{k,l}},\quad
Z_{k,l}|G_{0}&\overset{\text{iid}}{\sim}G_{0},\quad G_{0}=\sum_{l\geq
1}\omega_{0,l}\,\delta_{X_{l}^{\ast}},\quad
X_{\ell}^{\ast}\overset{\text{iid}}{\sim}H,\\\ (\omega_{k,l})_{l\geq
1}\overset{\text{iid}}{\sim}\textsc{GEM}(\beta),&\quad(\omega_{0,l})_{l\geq
1}\sim\textsc{GEM}(\beta_{0}),\end{split}$ (7)

where independence holds true between the sequences $(\omega_{k,l})_{l\geq 1}$
and $(Z_{k,l})_{l\geq 1}$ and between $(\omega_{0,l})_{l\geq 1}$ and
$(X_{l}^{\ast})_{l\geq 1}$. Combining the stick-breaking representation and a
closure property of the DP with respect to grouping, one further has

$G^{\ast}_{k}=\sum_{l\geq
1}\omega^{\ast}_{k,l}\delta_{X_{l}^{\ast}},G_{0}=\sum_{l\geq
1}\omega_{0,l}\delta_{X_{l}^{\ast}},$

where $((\omega^{\ast}_{k,l})_{l\geq
1}\mid\bm{\omega}_{0})\overset{\text{iid}}{\sim}\textsc{DP}(\beta;\bm{\omega}_{0})$,
$\bm{\omega}_{0}=(\omega_{0,l})_{l\geq 1}\sim\textsc{GEM}(\beta_{0})$ and
$X_{l}^{\ast}\overset{\text{iid}}{\sim}H$, for $l\geq 1$.

In this scheme, the clustering of populations is governed, a priori, by the
NDP layer $Q$ through $(\pi^{\ast}_{k})_{k\geq 1}\sim\textsc{GEM}(\alpha)$.
However, the aforementioned degeneracy issue of the NDP, a posteriori, is
successfully avoided. The intuition is quite straightforward: unlike for the
NDP, the distinct distributions $G^{\ast}_{k}$ in the HHDP are dependent and
have a common random discrete base measure $G_{0}$, which leads to shared
atoms across the $G_{k}^{\ast}$’s and thus borrowing of information, similarly
to the HDP case.

###  3.1 Some distributional properties

Given the discreteness of
$(G_{1},\ldots,G_{J})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$, the key
quantity to derive is the induced random partition, which controls the
clustering mechanism of the model. However, it is useful to start with a
description of pairwise dependence of the elements of the vector
$(G_{1},\ldots,G_{J})$, which allows a better understanding of the model and
intuitive parameter elicitation. To this end, as customary, we evaluate the
correlation between $G_{j}(A)$ and $G_{j^{\prime}}(A)$: whenever it does not
depend on the specific measurable set $A\subset\mathbb{X}$, it is used as a
measure of overall dependence between $G_{j}$ and $G_{j^{\prime}}$.

######  Proposition 1.

If $(G_{1},\ldots,G_{J})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$ and $A$
is a measurable subset of $\mathbb{X}$, then

$\displaystyle\mathrm{Var}[G_{j}(A)]$
$\displaystyle=\dfrac{H(A)[1-H(A)](\beta_{0}+\beta+1)}{(\beta+1)(\beta_{0}+1)}$
$\displaystyle(j=1,\ldots,J),$
$\displaystyle\mathrm{Corr}[G_{j}(A),G_{j^{\prime}}(A)]$
$\displaystyle=1-\dfrac{\alpha\beta_{0}}{(\alpha+1)(\beta+\beta_{0}+1)}$
$\displaystyle(j\neq j^{\prime}).$

Arguments similar to those in the proof of Proposition 1 lead to determine the
correlation between observations, either from the same or from different
samples.

######  Proposition 2.

If $\\{\bm{X}_{j}:j=1,\ldots,J\\}$ are from
$(G_{1},\ldots,G_{J})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$ according to
(1), then

$\displaystyle\mathrm{Corr}(X_{j,i},X_{j^{\prime},i^{\prime}})=\mathbb{P}(X_{j,i}=X_{j,i^{\prime}})=\begin{cases}\dfrac{1}{\beta_{0}+1}+\dfrac{\beta_{0}}{(1+\alpha)(1+\beta)(1+\beta_{0})}&(j\neq
j^{\prime})\\\\[7.0pt]
\dfrac{\beta+\beta_{0}+1}{(\beta+1)(\beta_{0}+1)}&(j=j^{\prime}).\end{cases}$

Figure 2: Correlations as functions of the hyperparameters $\beta$ and
$\beta_{0}$ with $\alpha=1$. The left plot represents the correlation between
random probabilities $G_{j}(A)$, the middle one between observations collected
in the same population and the right one between observations from different
populations.

The correlation between observations of the same sample depends only on the
parameters of the underlying HDP$(\beta,\beta_{0};H)$ that governs the atoms
$G_{k}^{\ast}$: this is not surprising since, whatever the value of the
parameter $\alpha$ at the NDP layer, observations from the same sample are
exchangeable. Moreover, an appealing feature is that such a correlation is
higher than for the case of observations from different samples, i.e. $j\neq
j^{\prime}$. As for the dependence on the hyperparameters
$(\alpha,\beta_{0},\beta)$, when $\alpha\to\infty$ the $G_{j}$’s ar forced to
equal different unique distributions $G_{k}^{\ast}$, similarly to the NDP
case. However, unlike the NDP, this does not imply that the distributions are
independent, and the correlation is controlled by the hyperparameters $\beta$
and $\beta_{0}$ (increasing in $\beta$ and decreasing in $\beta_{0}$). In Fig.
2 we report the aforementioned correlations as functions of $\beta$ and
$\beta_{0}$ with $\alpha$ set equal $1$. Finally, if $\alpha\rightarrow 0$ the
a priori probability to degenerate to the exchangeable case, i.e. all
$G_{j}$’s coincide a.s., tends to $1$ and so does also
$\mathrm{Cor}[G_{j}(A),G_{j^{\prime}}(A)]$.

We now investigate the random partition structure associated with a HHDP,
namely the partition of $\\{1,\ldots,n\\}$, with $n=\sum_{j=1}^{J}I_{j}$,
induced by a partially exchangeable sample $\bm{X}$ modeled as in (1). Since a
$\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$ arises from the composition of two
discrete random structures, it is clear that the partition induced by $\bm{X}$
will depend on the partition, say $\Psi^{(J)}$, of the random probability
measures $G_{1},\ldots,G_{J}$. As for the latter, the $G_{i}$’s are drawn from
a discrete random probability measure on $\mathscr{P}_{\mathbb{X}}$ whose
weights have a $\textsc{GEM}(\alpha)$ distribution and whose atoms are almost
surely different since they are sampled from an
$\textsc{HDP}(\beta,\beta_{0};H)$. Then the probability distribution of
$\Psi^{(J)}$ is the well–known Ewens sampling formula, namely

$\mathbb{P}[\Psi^{(J)}=\\{B_{1},\ldots,B_{R}\\}]=\phi_{R}^{(J)}(m_{1},\ldots,m_{R})=\frac{\alpha^{R}}{\alpha_{(J)}}\prod_{r=1}^{R}(m_{r}-1)!,$

where $\\{B_{1},\ldots,B_{R}\\}$ is a partition of $\\{1,\ldots,J\\}$, with
$1\leq R\leq J$, the frequencies $m_{r}=\mbox{card}(B_{r})$ are such that
$\sum_{r=1}^{R}m_{r}=J$ and $\alpha_{(J)}=\Gamma(\alpha+J)/\Gamma(\alpha)$.
This structure a priori implies, as in the NDP case, that
$\mathbb{P}(G_{j}=G_{\kappa})\in(0,1)$ for any $j\neq\kappa$. However, unlike
the NDP, a posteriori the HHDP yields
$\mathbb{P}(G_{j}=G_{\kappa}\mid\bm{X})<1$, regardless of the shared clusters
across the samples $\bm{X}$. Moreover, let
$\Phi_{D,R}^{(n)}(\,\cdots\,;\beta,\beta_{0})$ denote the pEPPF of a
$\textsc{HDP}(\beta,\beta_{0};H)$, namely

$\Phi_{D,R}^{(n)}(\bm{n}_{1}^{\ast},\ldots,\bm{n}_{R}^{\ast};\beta,\beta_{0})=\mathbb{E}\int_{\mathbb{X}_{*}^{D}}\prod_{d=1}^{D}\hat{G}_{1}(\,\textrm{d}{x}_{d})^{n_{1,d}^{\ast}}\>\cdots\>\hat{G}_{R}(\,\textrm{d}{x}_{d})^{n_{R,d}^{\ast}},$

where $(\hat{G}_{1},\ldots,\hat{G}_{R})\sim\textsc{HDP}(\beta,\beta_{0};H)$,
$D\in\\{1,\ldots,n\\}$ and $\sum_{r=1}^{R}\sum_{d=1}^{D}n_{r,d}^{*}=n$. An
explicit expression of $\Phi_{D,R}^{(n)}$ has been established in Camerlenghi
et al., 2019a , even beyond the DP case. Now we can state the pEPPF induced by
$\\{\bm{X}_{j}:j=1,\ldots,J\\}$ in (1), where $\mathcal{L}$ is the law of a
$\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$.

######  Theorem 1.

The random partition induced by the partially exchangeable array
$\\{\bm{X}_{j}:j=1,\ldots,J\\}$ drawn from
$(G_{1},\ldots,G_{J})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$, according
to (1), is characterized by the following pEPPF

$\Pi_{D}^{(n)}(\bm{n}_{1},\ldots,\bm{n}_{J})=\sum\phi_{R}^{(J)}(m_{1},\ldots,m_{R};\alpha)\Phi_{D,R}^{(n)}(\bm{n}_{1}^{\ast},\ldots,\bm{n}^{\ast}_{R};\beta,\beta_{0}),$
(8)

where the sum runs over all partitions $\\{B_{1},\ldots,B_{R}\\}$ of
$\\{1,\ldots,J\\}$ and $n^{\ast}_{r,d}=\sum_{j\in B_{r}}n_{j,d}$ for each
$r\in\\{1,\ldots,R\\}$, $d\in\\{1,\ldots,D\\}$.

Given the composition structure underlying the
$\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$, the pEPPF (8) unsurprisingly is a
mixture of pEPPF’s induced by different HDPs. For ease of interpretation
consider the case of $J=2$ populations and note that the pEPPF boils down to

$\Pi_{D}^{(n)}(\bm{n}_{1},\bm{n}_{2})=\dfrac{1}{\alpha+1}\Phi_{D,1}(\bm{n}_{1}+\bm{n}_{2})+\dfrac{\alpha}{\alpha+1}\Phi_{D,2}(\bm{n}_{1},\bm{n}_{2}),$
(9)

where $\Phi^{(n)}_{D,1}$ is the EPPF of a single
$\textsc{HDP}(\beta,\beta_{0};H)$, namely $J=1$, while $\Phi^{(n)}_{D,2}$ is
the pEPPF of a $\textsc{HDP}(\beta,\beta_{0};H)$ with two samples, namely
$J=2$. Clearly (9) arises from mixing with respect to partitions of
$\\{G_{1},G_{2}\\}$ in either $R=1$ and $R=2$ groups, where the former
corresponds to exchangeability across the two populations. Still for the case
$J=2$, a straightforward application of the pEPPF leads to the posterior
probability of gathering the two probability curves, $G_{1}$ and $G_{2}$, in
the same cluster thus making the two samples exchangeable, or homogeneous.

######  Proposition 3.

If the sample $\\{\bm{X}_{j}:j=1,2\\}$ is from
$(G_{1},G_{2})\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$, according to (1),
the posterior probability of degeneracy is

$\mathbb{P}(G_{1}=G_{2}\mid\bm{X})=\dfrac{\Phi^{(n)}_{D,1}(\bm{n}_{1}+\bm{n}_{2})}{\Phi^{(n)}_{D,1}(\bm{n}_{1}+\bm{n}_{2})+\alpha\,\Phi^{(n)}_{D,2}(\bm{n}_{1},\bm{n}_{2})},$
(10)

where $\Phi^{(n)}_{D,1}$ and $\Phi^{(n)}_{D,2}$ are the EPPF and the pEPPF
induced by the $\textsc{HDP}(\beta,\beta_{0};H)$ for a single exchangeable
sample and for two partially exchangeable samples, respectively.

The pEPPF is a fundamental tool in Bayesian calculus and it plays, in the
partially exchangeable framework, the same role of the EPPF in the
exchangeable case. Indeed, the pEPPF governs the learning mechanism, e.g. the
strength of the borrowing information, clustering, and, in view of Proposition
3, it allows to perform hypothesis testing for distributional homogeneity
between populations. Finally, one can obtain a Pólya urn scheme that is
essential for inference and prediction, See LABEL:subsec:MargGibbs in the
Supplementary Material. In the next section, we provide a characterization of
the $\textsc{HHDP}(\alpha,\beta,\beta_{0};H)$ that is reminiscent of the
popular Chinese restaurant franchise metaphor for the HDP and allows us to
devise a suitable sampling algorithm and further understand the model
behavior.

###  3.2 The hidden Chinese restaurant franchise

The marginalization of the underlying random probability measures, as
displayed in Theorem 1, can be characterized in terms of a hidden Chinese
restaurant franchise (HCRF) metaphor. This representation sheds further light
on the HHDP and clarifies the sense in which it generalizes the well-known
Chinese restaurant (CRP) and franchise (CRF) processes induced by the DP and
the HDP, respectively. For simplicity we consider the case $J=2$.

As with simpler sampling schemes, all restaurants of the franchise share the
same menu, which has an infinite number of dishes generated by the non–atomic
base measure $H$. However, unlike the standard CRF, the restaurants of the
franchise are merged into a single one if $G_{1}=G_{2}$, while they differ if
$G_{1}\neq G_{2}$. Moreover, each $X_{j,i}$ identifies the label of the dish
that customer $i$ from the $j$–th population chooses from the shared menu
$(X_{d}^{\ast})_{d\geq 1}$, with the unique dishes
$X^{\ast}_{d}\overset{\text{iid}}{\sim}H$. If $G_{1}\neq G_{2}$, customers may
be assigned to different restaurants and when $G_{1}=G_{2}$, they are all
seated in the same restaurant. Given such a grouping of the restaurants, the
customers are, then, seated according to the CRF applied either to a single
restaurant or to two distinct restaurants (Teh et al.,, 2006; Camerlenghi et
al.,, 2018). Furthermore, each restaurant has infinitely many tables. The
first customer $i$ who arrives at a previously unoccupied table chooses a dish
that is shared by all the customers who will join the table afterward. It is
to be noted that distinct tables within each restaurant and across restaurants
may share the same dish. An additional distinctive feature, compared to the
CRF, is that tables can be shared across populations when they are assigned to
the same restaurant, i.e. when $G_{1}=G_{2}$. Accordingly, the allocation of
each customer $X_{j,i}$ to a specific restaurant clearly depends on having
either $G_{1}=G_{2}$ or $G_{1}\neq G_{2}$.

The sampling scheme simplifies if latent variables $T_{j,i}$’s, denoting the
tables’ labels for customer $i$ from population $j$, are introduced. We stress
that, if $G_{1}\neq G_{2}$, the number of shared tables across the two
populations is zero, given the populations $j=1,2$ are assigned to different
restaurants, labeled $r=1,2$, respectively. Conversely, if $G_{1}=G_{2}$, one
may have shared tables across populations, since they are assigned to the same
restaurant $r=1$.

Now define $q_{r,t,d}$ as the frequencies of observations sitting at table $t$
eating the $d$th dish, for a table specific to restaurant $r$. Moreover,
$D_{t}$ is the dish label corresponding to table $t$ and $\ell_{r,d}$ the
frequency of tables serving dish $d$ in restaurant $r$. Marginal frequencies
are represented with dots, e.g. $\ell_{r,\cdot}$ is the number of tables in
restaurant $r$. Throughout the symbol $\bm{x}^{-i}$ identifies either a set or
a frequency obtained upon removing the element $i$ from $\bm{x}$. Finally,
$\Delta$ stands for an indicator function such that $\Delta=1$ if
$G_{1}=G_{2}$, while $\Delta=0$ if $G_{1}\neq G_{2}$.

The stepwise structure of the sampling procedure reflects the composition of
the three layers $\mathcal{L}(G_{j}|Q)$, $\mathcal{L}(Q|G_{0})$ and
$\mathcal{L}(G_{0})$ in (7) relying on a conditional CRF. First, one sample
the populations’ clustering $\Delta$ and, given the allocations of the
populations to the restaurants, one has a CRF. Hence, the algorithm becomes

  1. (1)

Sample the population assignments to the restaurants from
$\mathbb{P}(\Delta=1)=1/(\alpha+1)$.

  2. (2)

Sequentially sample the table assignments $T_{j,i}$ and corresponding dishes
$D_{T_{j,i}}$ from

$\displaystyle
p(T_{j,i},D_{T_{j,i}}\mid\bm{T}^{-(ji+)},\bm{X}^{-(ji+)},\Delta)\propto\begin{cases}T_{j,i}=t&\frac{q_{r,t,\cdot}^{-(ji+)}}{q_{r,\cdot,\cdot}^{-(ji+)}+\beta}\\\
T_{j,i}=t^{\text{new}},D_{t^{\text{new}}}=d&\frac{\beta}{q_{r,\cdot,\cdot}^{-(ji+)}+\beta}\frac{\ell_{\cdot,d}^{-(ji+)}}{\ell_{\cdot,\cdot}^{-(ji+)}+\beta_{0}}\\\
T_{j,i}=t^{\text{new}},D_{t^{\text{new}}}=d^{\text{new}}&\frac{\beta}{q_{r,\cdot,\cdot}^{-(ji+)}+\beta}\frac{\beta_{0}}{\ell_{\cdot,\cdot}^{-(ji+)}+\beta_{0}},\end{cases}$

where $(ji+)=\\{(ji^{\prime}):i^{\prime}\geq
i\\}\cup\\{(j^{\prime}i^{\prime}):j^{\prime}\geq j\\}$ is the index set
associated to the future random variables not yet sampled.

##  4 Posterior Inference for HHDP mixture models

Thanks to the results of Section 3, we now devise MCMC algorithms for drawing
posterior inferences with mixture models driven by a HHDP. Though the samplers
are tailored to mixture models, they are easily adapted to other inferential
problems such as e.g. survival analysis and species sampling. Henceforth,
$\mathcal{K}$ is a density kernel and we consider

$\displaystyle X_{j,i}\mid\theta_{j,i}$
$\displaystyle\overset{\text{ind}}{\sim}\mathcal{K}(\cdot|\theta_{j,i}),\quad\quad$
$\displaystyle(i=1,\ldots,I_{j}\quad j=1,\ldots,J),$ (11)
$\displaystyle\theta_{j,i}\mid G_{j}$
$\displaystyle\overset{\text{ind}}{\sim}G_{j},\quad\quad\qquad\>\>$
$\displaystyle(i=1,\ldots,I_{j},\quad j=1,\ldots,J),$
$\displaystyle(G_{1},\ldots,G_{J})$
$\displaystyle\sim\textsc{HHDP}(\alpha,\beta,\beta_{0};H).$

We develop two samplers: (i) a marginal algorithm that relies on the posterior
degeneracy probability (Proposition 3) in LABEL:subsec:MargGibbs of the
Supplementary Material; (ii) a conditional blocked Gibbs sampler, in the same
spirit of the sampler proposed for the NDP by Rodríguez et al., (2008), in
Section 4.1. As for (i), the underlying random probability measures $G_{0}$
and $G_{k}^{\ast}$’s are integrated out leading to urn schemes that extend the
class of Blackwell-MacQueen Pólya urn processes. In such a way we generalize
the a posteriori sampling scheme of the Chinese restaurant process for the DP
mixture Neal, (2000) and the one of the Chinese restaurant franchise for the
HDP mixture (Teh et al.,, 2006). In the Supplementary Material, we describe
the marginal sampler for the case of $J=2$ populations. Even if in principle
it can be generalized in a straightforward way, it is computationally
intractable for a larger number of populations. Similarly to the hidden
Chinese restaurant franchise situation, one has to evaluate the posterior
probability of all possible groupings of $G_{1},\ldots,G_{J}$, which boils
down to $\mathbb{P}(G_{1}=G_{2}|\bm{X})$ when $J=2$ but becomes involved for
$J>2$.

This shortcoming is overcome by the conditional algorithm we derive in Section
4.1, which relies on finite–dimensional approximations of the trajectories of
the underlying random probability measure. Its effectiveness in dealing with
$J>2$ populations is further illustrated in the synthetic data example 5.2 and
in the application of Section 5.3.

###  4.1 A conditional blocked Gibbs sampler

A more effective algorithm is based on a simple blocked conditional procedure.
To this end, we use a finite approximation of the DP in the spirit of Muliere
and Tardella, (1998) and Ishwaran and James, (2001). However, instead of
truncating the stick–breaking representation of the DP, we use a finite
Dirichlet approximation. See Ishwaran and Zarepour, (2002). Therefore, we
approximate $\bm{\pi}^{\ast},\bm{\omega}_{0}^{\ast}$, with a $K$– and an
$L$–dimensional Dirichlet distribution, respectively. More precisely, we
consider the following approximation

$\bm{\pi}^{\ast}\sim\textsc{Dir}(\alpha/K,\ldots,\alpha/K),\qquad\bm{\omega}_{0}^{\ast}\sim\textsc{Dir}\big{(}\beta_{0}/L,\ldots,\beta_{0}/L\big{)}\\\
$ (12)

implying that
$(\bm{\omega}^{\ast}_{k}\mid\bm{\omega}^{\ast}_{0})\overset{\text{iid}}{\sim}\textsc{Dir}(\beta\,\bm{\omega}_{0}^{\ast})$,
for $k\geq 1$.

Introduce the auxiliary variables $z_{j}$ and $\zeta_{j,i}$ which represent
the distributional and observational cluster memberships, respectively, such
that $z_{j}=k$ and $\zeta_{j,i}=l$ if and only if $G_{j}=G^{\ast}_{k}$ and
$\theta_{j,i}=\theta^{\ast}_{l}$. Henceforth,
$\bm{S}=\\{(\theta^{\ast}_{l})_{l=1}^{L},\bm{\pi}^{\ast},\bm{\omega}^{\ast}_{0},(\bm{\omega}^{\ast}_{k})_{k=1}^{K},(z_{j})_{j=1}^{J},(\zeta_{j,i})_{j,i},(X_{j,i})_{j,i}\\}$
and, in order to identify the full conditionals of the Gibbs sampler, we note
that under the finite Dirichlet approximation (12)

$p(\bm{S})=p(\bm{\pi}^{\ast})p(\bm{\omega}^{\ast}_{0})\bigg{[}\prod_{l=1}^{L}p(\theta^{\ast}_{l})\bigg{]}\bigg{[}\prod_{k=1}^{K}p(\bm{\omega}^{\ast}_{k}\mid\bm{\omega}^{\ast}_{0})\bigg{]}\bigg{\\{}\prod_{j=1}^{J}p(z_{j}\mid\bm{\pi}^{\ast})\bigg{[}\prod_{i=1}^{I_{j}}p(X_{j,i}\mid\theta^{\ast}_{\zeta_{j,i}})p(\zeta_{j,i}\mid\bm{\omega}^{\ast}_{z_{j}})\bigg{]}\bigg{\\}}.$

This leads to the following

  1. (1)

Sample the unique $\theta^{\ast}_{l}$ from

$p(\theta^{\ast}_{l}\mid\bm{S}^{-\theta^{\ast}_{l}})\propto
H(\theta^{\ast}_{l})\prod_{\\{j,i:\zeta_{j,i}=l\\}}\mathcal{K}(X_{j,i}\mid\theta^{\ast}_{l}).$

  2. (2)

Sample distributional cluster probabilities from

$p(\bm{\pi}^{\ast}\mid\bm{S}^{-\bm{\pi}^{\ast}})=\textsc{Dir}(\bm{\pi}^{\ast}\mid\alpha/K+m_{1},\ldots,\alpha/K+m_{K}),$

with $m_{k}=\sum_{j=1}^{J}\mathbb{1}\\{z_{j}=k\\}$.

  3. (3)

Sample probability weights of the base DP from

$\displaystyle\begin{split}p(\bm{\omega}^{\ast}_{0}\mid\bm{S}^{-\bm{\omega}^{\ast}_{0}})\propto\prod_{l=1}^{L}\bigg{[}\dfrac{(\omega_{0,l}^{\ast})^{\beta_{0}/L-1}\xi_{l}^{\beta\omega_{0,l}^{\ast}}}{\Gamma(\beta_{0}\omega_{0,l}^{\ast})^{K}}\bigg{]},\end{split}$
(13)

with $\xi_{l}=\prod_{k=1}^{K}\omega_{k,l}^{\ast}$.

  4. (4)

Sample the observational cluster probabilities independently from

$p(\bm{\omega}^{\ast}_{k}\mid\bm{S}^{-\bm{\omega}^{\ast}_{k}})=\textsc{Dir}(\bm{\omega}^{\ast}_{k}\mid\beta\bm{\omega}^{\ast}_{0}+\bm{n}_{k}),$

with
$n_{k,l}=\sum_{\\{j:z_{j}=k\\}}\sum_{i=1}^{I_{j}}\mathbb{1}\\{\zeta_{j,i}=l\\}$.

  5. (5)

Sample distributional and observational cluster membership from

$\displaystyle p(z_{j}=k\mid\bm{S}^{-\\{z_{j},\bm{\zeta}_{j}\\}})$
$\displaystyle\propto\pi^{\ast}_{k}\prod_{i=1}^{I_{j}}\sum_{l=1}^{L}\omega^{\ast}_{k,l}\mathcal{K}(X_{j,i}\mid\theta^{\ast}_{l})$
$\displaystyle(k=1,\ldots,K),$ $\displaystyle
p(\zeta_{j,i}=l\mid\bm{S}^{-\zeta_{j,i}})$
$\displaystyle\propto\omega^{\ast}_{z_{j}l}\,\mathcal{K}(X_{j,i}\mid\theta^{\ast}_{l})$
$\displaystyle(l=1,\ldots,L).$

Importantly, all the full conditional distributions are available in simple
closed forms, with the exception of the distributions of
$\bm{\omega}_{0}^{\ast}$ and, possibly, of $\theta^{\ast}_{l}$. To update
$\bm{\omega}_{0}^{\ast}$ we perform a Metropolis-Hastings step, where we work
on the unconstrained space $\mathbb{R}^{L-1}$ after the transformation
$[\log(\omega_{0,1}/\omega_{0,L}),\ldots,\log(\omega_{0,L-1}/\omega_{0,L})]$
and we adopt a component–wise adaptive random walk proposal following Roberts
and Rosenthal, (2009). The update of the unique atoms $\theta^{\ast}_{l}$ is
standard, as with the DP mixture model in the exchangeable case.

In Section 5 we assume a Gaussian kernel
$\mathcal{K}(\cdot|\theta)=\textsc{N}(\cdot|\mu,\sigma^{2})$ and a conjugate
Normal-inverse-Gamma base measure
$H(\cdot)=\textsc{NIG}(\cdot\mid\mu_{0},\lambda_{0},s_{0},S_{0})$ and obtain

$p(\theta^{\ast}_{l}\mid\bm{S}^{-\theta^{\ast}_{l}})=\textsc{NIG}(\theta^{\ast}_{l}\mid\mu_{l},\lambda_{l},s_{l},S_{l}),$

with $\mu_{l}=\dfrac{n_{l}\bar{y}_{l}+\lambda_{0}\mu_{0}}{\lambda_{0}+n_{l}}$,
$S_{l}=S_{0}+\dfrac{1}{2}\bigg{(}e^{2}_{l}+\dfrac{n_{l}\lambda_{0}(\bar{y}_{l}-\mu_{0})^{2}}{\lambda_{0}+n_{l}}\bigg{)}$,
$\lambda_{l}=\lambda_{0}+n_{l}$, and $s_{l}=n_{l}/2+s_{0}$, where
$n_{l}=\sum_{j=1}^{J}\sum_{i=1}^{I_{j}}\mathbb{1}\\{\zeta_{j,i}=l\\}$,
$\bar{y}_{l}=\sum_{\\{j,i:\zeta_{j,i}=l\\}}X_{j,i}/n_{l}$, and
$e^{2}_{l}=\sum_{\\{j,i:\zeta_{j,i}=l\\}}(X_{j,i}-\bar{y}_{l})^{2}$ are the
observational cluster sizes, means and deviances, respectively.

##  5 Illustration

In this section, we compare the performance of our proposal (11) with the same
model where the HHDP is replaced by a NDP as in (5), on synthetic data
involving $J=2$ and $J=4$ populations. Note that for the latter, the
implementation of the latent nested prior process mixture of Camerlenghi et
al., 2019b  is not feasible, while the proposed HHDP mixture model can easily
handle that level of complexity. The inferential results that we display are
obtained by relying on the blocked Gibbs sampler of Section 4.

###  5.1 Inference with two populations

The data are simulated from the same scenarios considered in Camerlenghi et
al., 2019b . More precisely, we consider two populations and the data in each
population are iid from a mixture of two normals:

  * Scen 1.

We simulate the data from the two populations independently from the same
density

$X_{1,i}\overset{\text{d}}{=}X_{2,i^{\prime}}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.5\textsc{N}(0,1)+0.5\textsc{N}(0,1).$

  * Scen 2.

We simulate the data in the two populations independently from mixtures of two
normals with one shared component

$X_{1,i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.9\textsc{N}(5,0.6)+0.1\textsc{N}(10,0.6)\quad\quad
X_{2,i^{\prime}}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.1\textsc{N}(5,0.6)+0.9\textsc{N}(0,0.6).$

  * Scen 3.

We simulate the data in the two populations independently from mixtures of two
normals having the same components, though with different weights

$X_{1,i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.8\textsc{N}(5,1)+0.2\textsc{N}(0,1)\quad\quad
X_{2,i^{\prime}}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.2\textsc{N}(5,1)+0.8\textsc{N}(0,1).$

In all these scenarios we consider balanced sample sizes $I_{1}=I_{2}=100$ and
an HHDP mixture model (11), with $\alpha=1$, $\beta=1$, $\beta_{0}=1$ and
$H(\cdot)=\textsc{NIG}(\cdot\mid\mu_{0},\lambda_{0},s_{0},S_{0}).$ We set
standard values of the hyperparameters in terms of the mean $\bar{y}$ and
variance $\mathrm{Var}(y)$ of the data, i.e. $\mu_{0}=\bar{y}$,
$\lambda_{0}=1/(3\,\mathrm{Var}(y))$, $s_{0}=1$ and $S_{0}=4$. In drawing the
comparison between (11) and the $\textsc{NDP}(\alpha,\beta;H)$, we further set
$\alpha=\beta=1$. Furthermore, we set the concentration parameters all equal
to 1. In LABEL:app:hyperpar we perform a sensitivity analysis with respect to
hyperparameters’ specifications as done, for instance, by Zuanetti et al.,
(2018) for the NDP. The mean measure of the marginal underlying random
distributions $\mathbb{E}[G_{j}(A)]=H(A)$ is the same for all populations.
Also variances are comparable (see Proposition 1) since
$\mathrm{Var}[G_{j}(A)]$ equals $H(A)[1-H(A)]/2$ for the NDP and
$3H(A)[1-H(A)]/4$ for the HHDP. The sensitivity analysis leads, for all the
considered settings, to the same conclusions in terms of comparison of the two
models. Moreover, we fix the dimensions of the finite approximations $L=K=50$
in (12) and we do the same for the truncation levels in the algorithm of
Rodríguez et al., (2008). In the Supplementary Material, we perform an
empirical analysis trying different levels of $L$ and $K$ which corroborates
the fact that the approximation error is negligible in terms of inferential
results.

Inference is based on $10\,000$ iterations with the first half discarded as
burn-in. As for the output, besides obtaining density estimates for the two
populations we also determine the point estimate of the clustering of
observations that minimizes the variation of information (VI) loss function.
See Meilǎ, (2007) and Wade and Ghahramani, (2018) for detailed discussions on
VI and point summaries of probabilistic clustering. Additionally, we estimate
the probability that observations co-cluster, namely
$\mathbb{P}(\zeta_{j,i}=\zeta_{j^{\prime},i^{\prime}}\mid\bm{X})$ through the
average over MCMC draws

$\dfrac{\sum_{b=1}^{B}\mathbb{1}\\{\zeta_{j,i}^{b}=\zeta_{j^{\prime},i^{\prime}}^{b}\\}}{B},$

where $B$ is the number of MCMC iterations. These are visualized through
heatmaps as in Fig. 4, with colors ranging from white, if the probability is
$0$, to dark red, if the probability is $1$. Our analysis is completed by
reporting the estimated distributions of the numbers of mixture components in
each scenario.

As expected, both models yield accurate estimates of the true densities in all
scenarios. In Fig. 3 we report the true and estimated models under the third
scenario. In terms of clustering, in the first scenario both models correctly
cluster together the two populations, thus degenerating to the exchangeable
case as they should. However, in the second and third scenarios the NDP makes
the two samples $\bm{X}_{1}$ and $\bm{X}_{2}$ independent, therefore
preventing borrowing of information across the two populations. As the
distributions have a shared component, the only way for the NDP to recover
correctly the true densities is by missing such a component. Had it been
detected, the density estimates of the two populations would have been equal
and, thus, far from the truth. The point estimate of the observations’
clustering in Table 2, the heatmaps of the posterior co-clustering
probabilities in Fig. 4 and the posterior distributions of the overall number
of occupied components in Table 1 showcase the theoretical findings, namely
that the NDP in the second and third scenarios cannot learn the shared
components. Hence, it overestimates the total number of occupied components
and does not cluster observations across populations. In contrast, the HHDP
model is able to cluster observations across populations, learns the shared
components and borrows information also when the model does not degenerate to
the exchangeable case.

NDP HHDP  Pop 1 Pop 2 Pop 1 Pop 2  Scen III

Figure 3: True (dashed lines), posterior mean (solid lines) densities and 95% point-wise posterior credible intervals (shaded gray) estimated under the third scenario.  |  | Overall number of components
---|---|---
Scen | Model | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |  $\geq$10
I | NDP | 0 | 0.4090 | 0.3615 | 0.1647 | 0.0492 | 0.0136 | 0.0020 | 0 | 0 | 0
HHDP | 0 | 0.5374 | 0.3743 | 0.0788 | 0.0080 | 0.0016 | 0 | 0 | 0 | 0
II | NDP | 0 | 0 | 0 | 0.2959 | 0.3906 | 0.2151 | 0.0700 | 0.0256 | 0.0024 | 0.0004
HHDP | 0 | 0 | 0.5742 | 0.3339 | 0.0796 | 0.0116 | 0.0008 | 0 | 0 | 0
III | NDP | 0 | 0 | 0 | 0.1331 | 0.3055 | 0.2947 | 0.1743 | 0.0608 | 0.0232 | 0.0084
HHDP | 0 | 0.5010 | 0.3966 | 0.0856 | 0.0164 | 0.0004 | 0 | 0 | 0 | 0
Table 1: Posterior distributions of the number of overall occupied components estimated with the two models under different scenarios. | Scenario I | Scenario II | Scenario III
---|---|---|---
| NDP | HHDP | NDP | HHDP | NDP | HHDP
Population | 1 | 2 | 1 | 2 | 1 | 2 | 3 | 4 | 1 | 2 | 3 | 1 | 2 | 3 | 4 | 1 | 2
1 | 56 | 44 | 56 | 44 | 87 | 13 | 0 | 0 | 87 | 13 | 0 | 85 | 15 | 0 | 0 | 85 | 15
2 | 48 | 52 | 48 | 52 | 0 | 0 | 88 | 12 | 12 | 0 | 88 | 0 | 0 | 80 | 20 | 21 | 79
Table 2: Frequencies of observations in the two populations allocated to the
point estimate of the clustering that minimizes the VI loss with the two
models under different scenarios.

True HHDP NDP  Pop 1 Pop 2 Pop 1 Pop 2 Pop 1 Pop 2  Scenario I Pop 1 Pop 2
Scenario II Pop 1 Pop 2
Scenario III Pop 1 Pop 2

Figure 4: Heatmaps of the true and estimated posterior probability of co-
clustering of observations, ordered by population memberships, under the HHDP
and the NDP models, for the three different scenarios in Section 5.1.

###  5.2 Inference with more than two populations

Here we consider $J=4$ populations and deal with the same scenario discussed
in Beraha et al., (2021). More precisely, we simulate independently across
populations $I_{j}=100$ (for $j=1,\ldots,4$) observations as follows

$\displaystyle
X_{1,i}\overset{\text{d}}{=}X_{2,i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.5\textsc{N}(0,1)+0.5\textsc{N}(5,1)\quad
X_{3,i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.5\textsc{N}(0,1)+0.5\textsc{N}(-5,1)\quad
X_{4,i}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}0.5\textsc{N}(-5,1)+0.5\textsc{N}(5,1)$

Our prior corresponds to a Gaussian mixture model with the same specification
for the HHDP used in the previous Section with $J=2$ population. Fig. 5 shows
that the HHDP mixture model is able to recover the data generating densities
also in this scenario.

Pop 1 Pop 2 Pop 3 Pop 4

Figure 5: True (dashed lines), posterior mean (solid lines) densities and 95%
point-wise posterior credible intervals (shaded gray) estimated under the
fourth scenario.

In terms of clustering of populations the point estimate that minimizes the VI
loss coincides with the data generating truth. Fig. 6 reports the heatmaps of
the posterior co-clustering probabilities of the four populations that show
little uncertainty around the correct point estimate, e.g. the estimated
probability that populations $1$ and $2$ are correctly clustered together is
$0.9858$.

Figure 6: Heatmap of the estimated posterior probabilities of co-clustering of
the population estimated with the HHDP mixture model under the fourth scenario
in Section 5.2.

Finally, the point estimate of the observations’ clustering in Table 3 shows
the HHDP model is able to cluster observations across populations, learns the
shared components and borrows information also when there are more than two
populations.

observational cluster | 1 | 2 | 3
---|---|---|---
Pop 1 | 53 | 47 | 0
Pop 2 | 56 | 44 | 0
Pop 3 | 48 | 0 | 52
Pop 4 | 0 | 52 | 48
Table 3: Frequencies of observations in the four populations allocated to the
point estimate of the clustering that minimizes the VI loss with HHDP under
the fourth scenario.

###  5.3 Collaborative perinatal project data

A multi-center application is the focus of this section. We consider a data
set from the Collaborative Perinatal Project (CPP), a large prospective
epidemiologic study conducted from 1959 to 1974. Pregnant women were enrolled
in 12 hospitals between 1959 and 1966 and were followed over time. Among
several pre–pregnancy measurements, we focus on the birth weight $X_{j,i}$ for
non-smoking woman $i$ in center $j$. We assume the following Gaussian mixture
model:

$\displaystyle X_{j,i}\mid\mu_{j,i},\sigma_{j,i}$
$\displaystyle\overset{\text{ind}}{\sim}\textsc{N}(\mu_{j,i},\sigma_{j,i})$
$\displaystyle(i=1,\ldots,I_{j},\quad j=1;\ldots,12),$
$\displaystyle\mu_{j,i},\sigma_{j,i}\mid G_{j}$
$\displaystyle\overset{\text{ind}}{\sim}G_{j}$
$\displaystyle(i=1,\ldots,I_{j},\quad j=1;\ldots,12).$

The same HHDP prior used for the previous synthetic data is placed the vector
of random distributions. This model specification is coherent with what is
suggested by Dunson, (2010) for the CPP data. Indeed, it is known that the
pregnancy outcome can vary substantially for women from different ethnicity
and socioeconomic groups. Therefore, we specify a model allowing to capture
differences between the centers since different groups of hospitals can serve
different women. Canale et al., (2019) provide further analysis of the CPP
data.

The heatmap of the co-clustering posterior probability for the 12 hospitals is
shown in Fig. 7. Such probabilities imply that the clustering point estimate
of the hospitals that minimizes the VI loss has two blocks and, in the same
figure, the mean posterior densities associated with the two clusters are
reported. Given the partition of the hospitals, the posterior mean densities
are evaluated based on all patients belonging to hospitals in each of the two
partition groups. The heatmap shows the posterior distribution of the
clustering of the hospitals and can be used to perform uncertainty
quantification. As expected, the lack of well-separated data generating
mixtures of Gaussians entails more uncertainty around the point estimate of
the clustering of the populations with respect to the numerical experiments.
However, the heatmap shows that the point estimate of the clustering of
distributions is a reliable summary. More precisely, the point estimate that
minimizes the VI loss entails that the first cluster of hospitals includes the
hospitals with (reordered) labels $1,2,3$: these are well-separated from the
remaining hospitals according to the posterior probabilities of co-clustering
in the heatmap. The heatmap shows also that another meaningful point estimate
of the clustering of the hospitals is the finer partition
$\\{\\{1,2,3\\},\\{4,5,6,7\\},\\{8,9,10,11,12\\}\\}$. However, the VI loss
suggests a more parsimonious clustering of the hospitals in two blocks, that
is $\\{\\{1,2,3\\},\\{4,5,6,7,8,9,10,11,12\\}\\}$. Note that in the second
cluster of hospitals (red dashed density in Fig. 7) the distribution of the
birth weights is slightly shifted on lower values and the two mean densities
are similar in the two clusters of populations. Coherently the proposed model
allows to borrow information across clusters of hospitals for estimating the
posterior mean densities of the birth weights. Furthermore the model can be
used to identify clusters of women shared in the two different clusters of
hospitals. Indeed, Table 4 shows that some clusters of observations are shared
across different clusters of hospitals, thus allowing the borrowing of
information for estimating the densities of the birth weights in the two
groups.

Figure 7: Heatmap of the estimated posterior probability of co-clustering of hospitals and estimated population cluster-specific posterior densities for the CPP data. number of observational clusters | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
only in the second cluster of hospitals | 0.3530 | 0.3670 | 0.2040 | 0.0640 | 0.0100 | 0.0020
only in the first cluster of hospitals | 0.7750 | 0.1850 | 0.0340 | 0.0060 | 0 | 0
shared across clusters of hospitals | 0 | 0.1680 | 0.4800 | 0.2660 | 0.0780 | 0.0080
Table 4:  Posterior distributions of the number of clusters shared and not
shared across the two clusters of hospitals.

##  6 Discussion

As highlighted in the recent literature, NDP mixture models are often not an
appropriate tool for clustering simultaneously population distributions and
observations. In contrast, the HHDP, overcomes the issues plaguing the NDP,
while preserving tractability and clustering flexibility even when the number
of populations $J$ is larger than $2$. We have further devised sampling
schemes allowing for efficient inference and prediction. This work paves the
way for future intriguing research directions that we plan to address in
forthcoming work. First, it is natural to move beyond DPs and consider models
based on alternative discrete nonparametric priors, such as the Pitman-Yor
process and normalized completely random measures, while studying the induced
clustering. The characterization of the HHDP in terms of the induced random
partition suggests a nice connection of our work with recent and exciting
advances on time-dependent random partition models such as those proposed,
e.g., in Page et al., (2022) and Zanini et al., (2019). Indeed, these papers
define a general framework that can be tailored to HHDP priors for generating
time-dependent models suited for analyzing, e.g., longitudinal data thus
allowing for the investigation of the joint evolution of observational and
distributional clustering through time. The theory we have developed in
Sections 3 and 4 provides the necessary tools for successfully carrying out
such a program. Moreover, the general composition scheme, where we have
embedded the HHDP, seems a promising and effective approach for addressing
other interesting inferential problems, beyond density estimation and
clustering. Finally, the general scheme that we have introduced in (2) seems
an appropriate specification for capturing the inherent complexity and
heterogeneity of data that arise when drawing predictions with multivariate
species sampling models and when performing inferences in survival and
functional data analysis. These will be the object of forthcoming work.

## REFERENCES

  * Agrawal et al., (2013) Agrawal, P., Tekumalla, L. S., and Bhattacharya, I. (2013).  Nested Hierarchical Dirichlet process for nonparametric entity-topic analysis.  In Jt. Eur. Conf. Mach. Learn. Knowl. Discov. Databases, volume 8189 LNAI, pages 564–579.
  * Argiento et al., (2020) Argiento, R., Cremaschi, A., and Vannucci, M. (2020).  Hierarchical normalized completely random measures to cluster grouped data.  J. Amer. Statist. Assoc., 115(529):318–333.
  * Balocchi et al., (2021) Balocchi, C., George, E. I., and Jensen, S. T. (2021).  Clustering areal units at multiple levels of resolution to model crime incidence in Philadelphia.  Preprint arXiv: 2112.02059.
  * Bassetti et al., (2020) Bassetti, F., Casarin, R., and Rossini, L. (2020).  Hierarchical species sampling models.  Bayesian Anal., 15(3):809–838.
  * Beraha et al., (2021) Beraha, M., Guglielmi, A., and Quintana, F. A. (2021).  The semi-hierarchical dirichlet process and its application to clustering homogeneous distributions.  Bayesian Anal., 16(4):1187–1219.
  * (6) Camerlenghi, F., Dunson, D. B., Lijoi, A., Prünster, I., and Rodríguez, A. (2019b).  Latent nested nonparametric priors.  Bayesian Anal., 14:1303–1356.  (With discussion).
  * (7) Camerlenghi, F., Lijoi, A., Orbanz, P., and Prünster, I. (2019a).  Distribution theory for hierarchical processes.  Ann. Stat., 47(1):67–92.
  * Camerlenghi et al., (2018) Camerlenghi, F., Lijoi, A., and Prünster, I. (2018).  Bayesian nonparametric inference beyond the Gibbs-type framework.  Scand. J. Stat., 45(4):1062–1091.
  * Camerlenghi et al., (2021) Camerlenghi, F., Lijoi, A., and Prünster, I. (2021).  Survival analysis via hierarchically dependent mixture hazards.  Ann. Stat., 49(2):863–884.
  * Canale et al., (2019) Canale, A., Corradin, R., and Nipoti, B. (2019).  Importance conditional sampling for Bayesian nonparametric mixtures.  Preprint at arXiv: 1906.08147.
  * Christensen and Ma, (2020) Christensen, J. and Ma, L. (2020).  A Bayesian hierarchical model for related densities using Pólya trees.  J. R. Stat. Soc. Ser. B, 82(1):127–153.
  * Cifarelli and Regazzini, (1978) Cifarelli, D. M. and Regazzini, E. (1978).  Problemi statistici non parametrici in condizioni di scambiabilita parziale e impiego di medie associative.  Quaderni Istituto Matematica Finanziaria dell’Universita di Torino.
  * Denti et al., (2021) Denti, F., Camerlenghi, F., Guindani, M., and Mira, A. (2021).  A common atom model for the Bayesian nonparametric analysis of nested data.  J. Am. Stat. Assoc., (in press).
  * Dunson, (2010) Dunson, D. B. (2010).  Nonparametric Bayes applications to biostatistics.  In Bayesian Nonparametrics, pages 223–273. Cambridge University Press.
  * Escobar, (1994) Escobar, M. D. (1994).  Estimating normal means with a Dirichlet process prior.  J. Am. Stat. Assoc., 89(425):268–277.
  * Escobar and West, (1995) Escobar, M. D. and West, M. (1995).  Bayesian density estimation and inference using mixtures.  J. Am. Stat. Assoc., 90(430):577–588.
  * Ewens, (1990) Ewens, W. J. (1990).  Population genetics theory - the past and the future.  In Math. Stat. Dev. Evol. Theory, pages 177–227. Springer, Dordrecht.
  * de Finetti, (1938) de Finetti, B. (1938).  Sur la condition d’equivalence partielle.  Actual. Sci. Ind., 739:5–18.
  * Foti and Williamson, (2015) Foti, N. J. and Williamson, S. A. (2015).  A survey of non-exchangeable priors for Bayesian nonparametric models.  IEEE Trans. Pattern Anal. Mach. Intell., 37(2):359–371.
  * Ghosal and van der Vaart, (2017) Ghosal, S. and van der Vaart, A. (2017).  Fundamentals of nonparametric Bayesian inference.  Cambridge University Press.
  * Ishwaran and James, (2001) Ishwaran, H. and James, L. F. (2001).  Gibbs sampling methods for stick-breaking priors.  J. Am. Stat. Assoc., 96(453):161–173.
  * Ishwaran and Zarepour, (2002) Ishwaran, H. and Zarepour, M. (2002).  Exact and approximate sum representations for the Dirichlet process.  Can. J. Stat., 30(2):269–283.
  * James, (2008) James, L. (2008).  Discussion of Nested Dirichlet Process paper by Rodríguez, Dunson and Gelfand.  J. Am. Stat. Assoc., 483:1131.
  * Kallenberg, (2005) Kallenberg, O. (2005).  Probabilistic symmetries and invariance principles.  Springer.
  * Lo, (1984) Lo, A. Y. (1984).  On a class of Bayesian nonparametric estimates: I. density estimates.  Ann. Stat., 12(1):351–357.
  * MacEachern, (1999) MacEachern, S. N. (1999).  Dependent nonparametric processes.  In ASA Proc. Sect. Bayesian Stat. Sci., pages 50–55.
  * MacEachern, (2000) MacEachern, S. N. (2000).  Dependent Dirichlet processes.  Technical report, The Ohio State University.
  * Meilǎ, (2007) Meilǎ, M. (2007).  Comparing clusterings-an information based distance.  J. Multivar. Anal., 98(5):873–895.
  * Muliere and Tardella, (1998) Muliere, P. and Tardella, L. (1998).  Approximating distributions of random functionals of Ferguson-Dirichlet priors.  Can. J. Stat., 26(2):283–297.
  * Müller et al., (2011) Müller, P., Quintana, F., and Rosner, G. L. (2011).  A product partition model with regression on covariates.  J. Comput. Graph. Statist., 20(1):260–278.
  * Neal, (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  J. Comput. Graph. Stat., 9(2):249–265.
  * Page and Quintana, (2016) Page, G. L. and Quintana, F. A. (2016).  Spatial product partition models.  Bayesian Anal., 11(1):265–298.
  * Page and Quintana, (2018) Page, G. L. and Quintana, F. A. (2018).  Calibrating covariate informed product partition models.  Stat. Comput., 28(5):1009–1031.
  * Page et al., (2022) Page, G. L., Quintana, F. A., and Dahl, D. B. (2022).  Dependent modeling of temporal sequences of random partitions.  J. Comput. Graph. Stat., (in press).
  * Pitman, (2006) Pitman, J. (2006).  Combinatorial stochastic processes.  Springer.
  * Quintana et al., (2022) Quintana, F. A., Müller, P., Jara, A., and MacEachern, S. N. (2022).  The dependent dirichlet process and related models.  Stat. Sci., (in press).
  * Roberts and Rosenthal, (2009) Roberts, G. O. and Rosenthal, J. S. (2009).  Examples of adaptive MCMC.  J. Comput. Graph. Stat., 18(2):349–367.
  * Rodríguez et al., (2008) Rodríguez, A., Dunson, D. B., and Gelfand, A. E. (2008).  The nested Dirichlet process.  J. Am. Stat. Assoc., 103(483):483–1131.
  * Sethuraman, (1994) Sethuraman, J. (1994).  A constructive definition of Dirichlet priors.  Stat. Sin., 4(2):639–650.
  * Soriano and Ma, (2019) Soriano, J. and Ma, L. (2019).  Mixture modeling on related samples by $\psi$-stick breaking and kernel perturbation.  Bayesian Anal., 14(1):161–180.
  * Teh and Jordan, (2010) Teh, Y. W. and Jordan, M. I. (2010).  Hierarchical Bayesian nonparametric models with applications.  In Bayesian Nonparametrics, pages 158–207. Cambridge University Press.
  * Teh et al., (2006) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006).  Hierarchical Dirichlet processes.  J. Am. Stat. Assoc., 101(476):1566–1581.
  * Wade and Ghahramani, (2018) Wade, S. and Ghahramani, Z. (2018).  Bayesian cluster analysis: point estimation and credible balls.  Bayesian Anal., 13(2):559–626.
  * Zanini et al., (2019) Zanini, C. T. P., Müller, P., Ji, Y., and Quintana, F. A. (2019).  A Bayesian random partition model for sequential refinement and coagulation.  Biometrics, 75(3):988–999.
  * Zuanetti et al., (2018) Zuanetti, D. A., Müller, P., Zhu, Y., Yang, S., and Ji, Y. (2018).  Clustering distributions with the marginalized nested Dirichlet process.  Biometrics, 74(2):584–594.