Dependent Mixtures of Geometric Weights Priors

TL;DR

This paper introduces a novel pairwise dependent mixture model using geometric stick breaking processes, offering faster computation and effective estimation for partially exchangeable data.

Contribution

It develops a new dependent mixture model based on geometric weights, improving computational efficiency over traditional Dirichlet process methods.

Findings

Faster Gibbs sampler implementation.

Effective prediction and estimation with dependent mixtures.

Successful application on synthetic and real data.

Abstract

A new approach on the joint estimation of partially exchangeable observations is presented by constructing pairwise dependence between $m$ random density functions, each of which is modeled as a mixture of geometric stick breaking processes. This approach is based on a new random central masses version of the Pairwise Dependent Dirichlet Process prior mixture model (PDDP) first introduced in Hatjispyros et al. (2011). The idea is to create pairwise dependence through random measures that are location-preserving-expectations of Dirichlet random measures. Our contention is that mixture modeling with Pairwise Dependent Geometric Stick Breaking Process (PDGSBP) priors is sufficient for prediction and estimation purposes; moreover the associated Gibbs sampler is much faster and easier to implement than its Dirichlet Process based counterpart. To this respect, we provide a-priori-synchronized comparison studies under sparse $m$-scalable synthetic and real data examples.

Tables

Table 1. Table 1: Hellinger distances for the case m = 4 𝑚 4 m=4 .

$i$	${\mathscr{H}}_{𝒢}(f_i^{\left(4\right)},{\hat{f}}_i^{\left(4\right)})$	${\mathscr{H}}_{𝒟}(f_i^{\left(4\right)},{\hat{f}}_i^{\left(4\right)})$
$1$	$0.17$	$0.17$
$2$	$0.19$	$0.18$
$3$	$0.22$	$0.22$
$4$	$0.20$	$0.20$

Table 2. Table 2: Mean execution times in seconds per 10 3 superscript 10 3 10^{3} iterations.

$m$	Model	Sample size	MET
2	PDGSBP	$n_j^{\left(2\right)}=60$	0.57
	rPDDP		1.68
3	PDGSBP	$n_j^{\left(3\right)}=120$	2.16
	rPDDP		6.57
4	PDGSBP	$n_j^{\left(4\right)}=200$	5.30
	rPDDP		17.87

Table 3. Table 3: Hellinger distances between true and estimated densities for the case m = 10 𝑚 10 m=10 of the sparse scalable data example.

$i$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$
${\mathscr{H}}_{𝒢}(f_i^{\left(10\right)},{\hat{f}}_i^{\left(10\right)})$	$0.08$	$0.10$	$0.09$	$0.14$	$0.14$	$0.13$	$0.14$	$0.09$	$0.11$	$0.22$
${\mathscr{H}}_{𝒟}(f_i^{\left(10\right)},{\hat{f}}_i^{\left(10\right)})$	$0.09$	$0.11$	$0.10$	$0.15$	$0.12$	$0.10$	$0.14$	$0.09$	$0.09$	$0.22$

Table 4. Table 4: Hellinger distance between the true and the estimated densities.

$i$	${\mathscr{H}}_{𝒢}(f_i,{\hat{f}}_i)$	${\mathscr{H}}_{𝒟}(f_i,{\hat{f}}_i)$
$1$	$0.19$	$0.18$
$2$	$0.18$	$0.15$

Table 5. Table 5: Hellinger distances for the gamma mixture data model.

$i$	${\mathscr{H}}_{𝒢}(f_i,{\hat{f}}_i)$	${\mathscr{H}}_{𝒟}(f_i,{\hat{f}}_i)$
$1$	$0.13$	$0.11$
$2$	$0.19$	$0.18$

Table 6. Table 6: Hellinger distances between the true and the estimated densities for the three scenario example.

$s$	${\mathscr{H}}_{𝒢}(f_1^{\left(s\right)},{\hat{f}}_1^{\left(s\right)})$	${\mathscr{H}}_{𝒢}(f_2^{\left(s\right)},{\hat{f}}_2^{\left(s\right)})$	${\mathscr{H}}_{𝒢}(f_3^{\left(s\right)},{\hat{f}}_3^{\left(s\right)})$
$1$	$0.14$	$0.19$	$0.13$
$2$	$0.15$	$0.22$	$0.15$
$3$	$0.12$	$0.26$	$0.12$

Equations

$$((X_{1,i_1}{)}_{1\le i_1\le n_1},\dots ,(X_{m,i_m}{)}_{1\le i_m\le n_m})\sim ((X_{1,{\pi }_1(i_1)}{)}_{1\le i_1\le n_1},\dots ,(X_{m,{\pi }_m(i_m)}{)}_{1\le i_m\le n_m}).$$

$$\mathbb{P}(\cdot )=\sum _{k=1}^{\infty }{w}_k{\delta }_{{\theta }_k}(\cdot ),$$

$${\mathbb{P}}_z(\cdot )=\sum _{k=1}^{\infty }{w}_k(z){\delta }_{{\theta }_k(z)}(\cdot ),$$

$$f_j(x)={\int }_{\Theta }K(x|\theta ){\mathbb{G}}_j(d\theta ),$$

$$\mathrm{P}\{X_{ji}\in A_{ji}:1\le j\le m,1\le i\le n_j\}$$

$$={\int }_{{\mathcal{P}}^m}\mathrm{P}\{X_{ji}\in A_{ji}:1\le j\le m,1\le i\le n_j|{\mathbb{Q}}_1,\dots ,{\mathbb{Q}}_m\}\Pi (d{\mathbb{Q}}_1,\dots ,d{\mathbb{Q}}_m)$$

$$={\int }_{{\mathcal{P}}^m}\prod _{j=1}^m\mathrm{P}\{X_{ji}\in A_{ji}:1\le i\le n_j|{\mathbb{Q}}_j\}\Pi (d{\mathbb{Q}}_1,\dots ,d{\mathbb{Q}}_m)$$

$$={\int }_{{\mathcal{P}}^m}\prod _{j=1}^m\prod _{i=1}^{n_j}{\mathbb{Q}}_j(A_{ji})\Pi (d{\mathbb{Q}}_1,\dots ,d{\mathbb{Q}}_m).$$

$$x_{ji}|{\theta }_{ji}$$

$${\theta }_{ji}|{\mathbb{Q}}_j$$

$${\mathbb{Q}}_j=$$

$${\mathbb{P}}_{jl}\stackrel{\mathrm{iid}}{\sim }$$

$$f_j(x):=f_j(x|{\mathbb{Q}}_j)=\sum _{l=1}^m{p}_{jl}{g}_{jl}(x|{\mathbb{G}}_{jl}),{\mathbb{Q}}_j=\sum _{l=1}^m{p}_{jl}{\mathbb{G}}_{jl},1\le j\le m.$$

$${\mathbb{G}}_{jl}=\sum _{k=1}^{\infty }{q}_{jlk}{\delta }_{{\theta }_{jlk}}\mathrm{with}{q}_{jlk}={\lambda }_{jl}(1-{\lambda }_{jl}{)}^{k-1},{\lambda }_{jl}\sim h(\cdot |{\xi }_{jl}),{\theta }_{jlk}\stackrel{\mathrm{iid}}{\sim }{G}_0,$$

$$\mathbb{Q}=\left(p\otimes \mathbb{G}\right)1,$$

$$f_j(x_{ji},N_{ji}=r)=r^{-1}\sum _{l=1}^m{p}_{jl}{f}_N(r|{\lambda }_{jl})\sum _{k=1}^{r}K(x_{ji}|{\theta }_{jlk}),$$

$$f_j(x_{ji},N_{ji}=r,d_{ji}=k|{\delta }_{ji}=l)=\frac{1}{r}{f}_N(r|{\lambda }_{jl})\mathcal{I}(k\le r)K(x_{ji}|{\theta }_{jlk}).$$

$$f_j(x_{ji}|N_{ji}=r)=\sum _{l=1}^m\mathcal{W}(r|{\lambda }_{jl})\sum _{k=1}^r{r}^{-1}K(x_{ji}|{\theta }_{jlk}),$$

$$\mathcal{W}(r|{\lambda }_{jl})=\frac{p_{jl}{f}_N(r|{\lambda }_{jl})}{{\sum }_{l^{\prime }=1}^m{p}_{j{l}^{\prime }}{f}_N(r|{\lambda }_{j{l}^{\prime }})}.$$

$$f_j(x_{ji}|N_{ji}=r,{\delta }_{ji}=l)=\sum _{k=1}^r{r}^{-1}K(x_{ji}|{\theta }_{jlk}).$$

$$f_j(x_{ji}|{\delta }_{ji}=l)=\sum _{k=1}^{\infty }\left(\sum _{r=k}^{\infty }{r}^{-1}{f}_N(r|{\lambda }_{jl})\right)K(x_{ji}|{\theta }_{jlk}).$$

$$f_N(r|{\lambda }_{jl})=r{\lambda }_{jl}^2(1-{\lambda }_{jl}{)}^{r-1}\mathcal{I}(r\ge 1),$$

$$f_j(x_{ji}|{\delta }_{ji}=l)=\sum _{k=1}^{\infty }{q}_{jlk}K(x_{ji}|{\theta }_{jlk})\mathrm{with}{q}_{jlk}={\lambda }_{jl}(1-{\lambda }_{jl}{)}^{k-1},$$

$$f_j(x_{ji})=\sum _{l=1}^m{p}_{jl}\sum _{k=1}^{\infty }{q}_{jlk}K(x_{ji}|{\theta }_{jlk}).$$

$$f(\mathbf{x},\mathbf{N},\mathbf{d}|\mathbf{\delta })$$

$$x_{ji},N_{ji}|d_{ji},{\delta }_{ji},({\theta }_{jr{\delta }_{ji}}{)}_{1\le r\le m},{\lambda }_{j{\delta }_{ji}}\stackrel{\mathrm{ind}}{\sim }\prod _{r=1}^m{\left\{{\lambda }_{jr}^2(1-{\lambda }_{jr}{)}^{N_{ji}-1}K(x_{ji}|{\theta }_{jr{d}_{ji}})\right\}}^{{\delta }_{ji}^r}\mathcal{I}(N_{ji}\ge d_{ji})$$

$$d_{ji}|N_{ji}\stackrel{\mathrm{ind}}{\sim }\mathcal{DU}({\mathcal{S}}_{ji}),\mathrm{P}\{{\delta }_{ji}={\mathbf{e}}_l\}=p_{jl}$$

$$q_{jik}={\lambda }_{ji}(1-{\lambda }_{ji}{)}^{k-1},{\theta }_{jik}\stackrel{\mathrm{iid}}{\sim }{G}_0,k\in \mathbb{N}.$$

$$\mathbb{E}[g_{\mathbb{G}}(x{)}^2]=\left(\frac{1}{2-\lambda }\right)\left\{\lambda {\int }_{\Theta }K(x|\theta {)}^2{G}_0(d\theta )+2(1-\lambda ){\left({\int }_{\Theta }K(x|\theta )G_0(d\theta )\right)}^2\right\}.$$

$$\mathrm{Cov}(f_j(x),f_i(x))=p_{ji}{p}_{ij}\mathrm{Var}\left({\int }_{\Theta }K(x|\theta ){\mathbb{G}}_{ji}(d\theta )\right),$$

$$\mathrm{Var}\left({\int }_{\Theta }K(x|\theta ){\mathbb{G}}_{ji}(d\theta )\right)=\frac{{\lambda }_{ji}}{2-{\lambda }_{ji}}\mathrm{Var}(K(x|\theta )).$$

$$\mathrm{Corr}(f_j^{\mathcal{G}}(x),f_i^{\mathcal{G}}(x))=\frac{{\lambda }_{ji}{p}_{ji}{p}_{ij}}{2-{\lambda }_{ji}}{\left(\sum _{l=1}^m\sum _{r=1}^m\frac{p_{jl}^2{p}_{ir}^2{\lambda }_{jl}{\lambda }_{ir}}{(2-{\lambda }_{jl})(2-{\lambda }_{ir})}\right)}^{-1/2},$$

$$\mathrm{Corr}(f_j^{\mathcal{D}}(x),f_i^{\mathcal{D}}(x))=\frac{p_{ji}{p}_{ij}}{1+c_{ji}}{\left(\sum _{l=1}^m\sum _{r=1}^m\frac{p_{jl}^2{p}_{ir}^2}{(1+c_{jl})(1+c_{ir})}\right)}^{-1/2}.$$

$$\mathrm{Corr}(f_j^{\mathcal{G}}(x),f_i^{\mathcal{G}}(x))=\mathrm{Corr}(f_j^{\mathcal{D}}(x),f_i^{\mathcal{D}}(x))=p_{ji}{p}_{ij}{\left(\sum _{l=1}^m\sum _{r=1}^m{p}_{jl}^2{p}_{ir}^2\right)}^{-1/2}.$$

Full text

Keywords: Bayesian nonparametric inference; Mixture of Dirichlet process; Geometric stick breaking weights; Geometric Stick Breaking Mixtures; Dependent Process.

1. Introduction. In Bayesian nonparametric methods, the use of priors such as the Dirichlet process (Ferguson, 1973), is justified from the assumption that the observations are exchangeable, which means the distribution of $(X_{1},\ldots,X_{n})$ coincides with the distribution of $(X_{\pi(1)},\ldots,X_{\pi(n)}),$ for all $\pi\in S(n)$, where $S(n)$ is the set of permutations of $\{1,\ldots,n\}$. However, in real life applications, data are often partially exchangeable. For example, they may consist of observations sampled from $m$ populations, or may be sampled from an experiment conducted in $m$ different geographical areas. This means that the joint law is invariant under permutations within the $m$ subgroups of observations $(X_{j,i_{j}})_{1\leq i_{j}\leq n_{j}},\,1\leq j\leq m$, so for all $\pi_{j}\in S(n_{j})$

[TABLE]

When the exchangeability assumption fails one needs to use non–exchangeable priors. There has been substantial research interest following the seminal work of MacEachern (1999) in the construction of suitable dependent stochastic processes. Such then act as priors in Bayesian nonparametric models. These processes are distributions over a collection of measures indexed by values in some covariate space, such that the marginal distribution is described by a known nonparametric prior. The key idea is to induce dependence between a collection of random probability measures $(\mathbb{P}_{j})_{1\leq j\leq m}$, where each $\mathbb{P}_{j}$ comes from a Dirichlet process (DP) with concentration parameter $c>0$ and base measure $P_{0}$. Such random probability measures typically are used in mixture models to generate random density functions $f(x)=\int_{\Theta}K(x|\theta)\mathbb{P}(d\theta)$; see Lo (1984).

There is a variety of ways that a DP can be extended to dependent DP. Most of them use the stick-breaking representation (Sethuraman, 1994), that is

[TABLE]

where $(\theta_{k})_{k\geq 1}$ are independent and identically distributed from $P_{0}$ and $(w_{k})_{k\geq 1}$ is a stick breaking process; so if $(v_{k})_{k\geq 1}$ are independent and identically distributed from ${\cal B}e(1,c)$, a beta distribution with mean $(1+c)^{-1}$, then $w_{1}=v_{1}$ and for $k>1$, $w_{k}=v_{k}\prod_{l<k}(1-v_{l})$. Dependence is introduced through the weights and/or the atoms. A classical example of the use of dependent DP’s is the Bayesian nonparametric regression problem where a random probability measure $\mathbb{P}_{z}$ is constructed for each covariate $z$,

[TABLE]

where $(w_{k}(z),\theta_{k}(z))$ is a collection of processes indexed in $z$–space. Extensions to dependent DP models can be found in De Iorio et al. (2004), Griffin and Steel (2006), and Dunson and Park (2008).

Recently there has been growing interest for the use of simpler random probability measures which while simpler are yet sufficient for Bayesian nonparametric density estimation. The geometric stick breaking (GSB) random probability measure (Fuentes–García, et al. 2010) has been used for density estimation and has been shown to provide an efficient alternative to DP mixture models. Some recent papers extend this nonparametric prior to a dependent nonparametric prior. In the direction of covariate dependent processes, GSB processes have been seen to provide an adequate model to the traditional dependent DP model. For example, for Bayesian regression, Fuentes–Garcia et al. (2009) propose a covariate dependent process based on random probability measures drawn from a GSB process. Mena et al. (2011) used GSB random probability measures in order to construct a purely atomic continuous time measure–valued process, useful for the analysis of time series data. In this case, the covariate $z\geq 0$ denotes the time that each observation is (discretely) recorded and conditionally on each observation is drawn from a time–dependent nonparametric mixture model based on GSB processes. However, to the best of our knowledge, random probability measures drawn from a GSB process, for modeling related density functions when samples from each density function are available, has not been developed in the literature.

In this paper we will construct pairwise dependent random probability measures based on GSB processes. That is, we are going to model a finite collection of $m$ random distribution functions $(\mathbb{G}_{j})_{1\leq j\leq m}$, where each $\mathbb{G}_{j}$ is a GSB random probability measure, such that there is a unique common component for each pair $(\mathbb{G}_{j},\mathbb{G}_{j^{\prime}})$ with $j\neq j^{\prime}$. We are going to use these measures in the context of GSB mixture models, generating a collection of $m$ GSB pairwise dependent random densities $(f_{j}(x))_{1\leq j\leq m}$. Hence we obtain a set of random densities $(f_{1},\ldots,f_{m})$, where marginally each $f_{j}$ is a random density function

[TABLE]

thus generalizing the GSB priors to a multivariate setting for partially exchangeable observations.

In the problem considered here, these random density functions $(f_{j})_{1\leq j\leq m}$ are thought to be related or similar, e.g. perturbations of each other, and so we aim to share information between groups to improve estimation of each density, especially for those densities $f_{j}$ for which the corresponding sample size $n_{j}$ is small. In this direction, the main references include the work of Müller et al. (2014), Bulla et al. (2009), Kolossiatis et al. (2013) and Griffin et al. (2013); more rigorous results can be found in Lijoi et al. (2014A, 2014B). All these models have been proposed for the modeling of an arbitrary but finite number of random distribution functions, via a common part and an index specific idiosyncratic part so that for $0<p_{j}<1$ we have $\mathbb{P}_{j}=p_{j}\mathbb{P}_{0}+(1-p_{j})\mathbb{P}_{j}^{*},$ where $\mathbb{P}_{0}$ is the common component to all other distributions and $\{\mathbb{P}_{j}^{*}:j=1,\ldots,m\}$ are the idiosyncratic parts to each $\mathbb{P}_{j}$, and $\mathbb{P}_{0},\mathbb{P}_{j}^{*}\stackrel{{\scriptstyle\rm iid}}{{\sim}}{\cal DP}(c,P_{0})$. In Lijoi et al. (2014B) normalized random probability measures based on the $\sigma$–stable process are used for modeling dependent mixtures. Although similar (all models coincide only for the $m=2$ case), these models are different from our model which is based on pairwise dependence of a sequence of random measures (Hatjispyros et al. 2011, 2016A).

We are going to provide evidence through numerical experiments that dependent GSB mixture models are an efficient alternative to pairwise dependent DP (PDDP) priors. First, we will randomize the existing PDDP model of Hatjispyros et al. (2011, 2016A), by imposing gamma priors on the concentration masses (leading to the more efficient rPDDP model). Then, for the objective comparison of execution times, we will conduct a-priori synchronized density estimation comparison studies between the randomized PDDP and the pairwise dependent GSB process (PDGSBP) models using synthetic and real data examples.

This paper is organized as follows. In Section 2 we will demonstrate the construction of pairwise dependent random densities, using a dependent model suggested by Hatjispyros et al. (2011). We also demonstrate how specific choices of latent random variables can recover the model of Hatjispyros et al. and the dependent GSB model introduced in this paper. These latent variables will form the basis of a Gibbs sampler for posterior inference, given in Section 3. In Section 4 we resort to simulation. We provide comparison studies between the randomized version of the PDDP model and our newly introduced dependent GSB model, involving five cases of synthetic data and a real data set. Finally, Section 5 concludes with a summary and future work.

2. Preliminaries. We consider an infinite real valued process $\{X_{ji}:1\leq j\leq m,\,i\geq 1\}$ defined over a probability space $(\Omega,{\cal F},{\rm P})$, that is partially exchangeable as in (1). Let ${\cal P}$ denote the set of probability measures over $\mathbb{R}$; then de Finetti proved that there exists a probability distribution $\Pi$ over ${\cal P}^{m}$, which satisfies

[TABLE]

The de Finetti measure $\Pi$ represents a prior distribution over partially exchangeable observations.

We start off by describing the PDDP model, with no auxiliary variables, using only the de Finetti measure $\Pi$, marginal measures $\mathbb{Q}_{j}$, then, we proceed to the definition of a randomized version of it, and to the specific details for the case of the GSB random measures.

A. In Hatjispyros et al. (2011), the following hierarchical model was introduced. For $m$ subgroups of observations $\{(x_{ji})_{1\leq i\leq n_{j}}:1\leq j\leq m\}$,

[TABLE]

for some kernel density $K(\,\cdot\,|\,\cdot\,)$, concentration parameter $c>0$ and parametric central measure $P_{0}$ for which $\mathbb{E}(\mathbb{P}_{jl}(d\theta))=P_{0}(d\theta)$. So, we have assumed that the random densities $f_{j}(x)$ are dependent mixtures of the dependent random measures $\mathbb{Q}_{j}$ via $f_{j}(x|\mathbb{Q}_{j})=\int_{\Theta}K(x|\theta)\mathbb{Q}_{j}(d\theta)$, or equivalently, dependent mixtures of the $m$ independent mixtures $g_{jl}(x|\,\mathbb{P}_{jl})=\int_{\Theta}K(x|\,\theta)\,\mathbb{P}_{jl}(d\theta)$, $l=1,\ldots,m$. To introduce the rPDDP model, we randomize the PDDP model by sampling the $\mathbb{P}_{jl}$ measures from the independent Dirichlet processes ${\cal DP}(c_{jl},P_{0})$ and then impose gamma priors on the concentration masses, i.e. $\mathbb{P}_{jl}\stackrel{{\scriptstyle\rm ind}}{{\sim}}\,{\cal DP}(c_{jl},P_{0}),\quad c_{jl}\stackrel{{\scriptstyle\rm ind}}{{\sim}}{\cal G}(a_{jl},b_{jl}),\,\,1\leq j\leq l\leq m.$

B. To develop a pairwise dependent geometric stick breaking version, we let the random density functions $f_{j}(x)$ generated via

[TABLE]

The $g_{jl}(x):=g_{jl}(x|\,\mathbb{G}_{jl})=\int_{\Theta}K(x|\,\theta)\,\mathbb{G}_{jl}(d\theta)$ random densities are now independent mixtures of GSB processes, satisfying $g_{jl}=g_{lj}$, under the slightly altered definition

[TABLE]

where $h$ is a parametric density supported over the interval $(0,1)$ depending on some parameter $\xi_{jl}\in\Xi$, and $G_{0}$ is the associated parametric central measure.

The independent GSB processes $\{\mathbb{G}_{jl}:\,1\leq j,\,l\leq m\}$ form a matrix $\mathbb{G}$ of random distributions with $\mathbb{G}_{jl}=\mathbb{G}_{lj}$. In matrix notation

[TABLE]

where $p=(p_{jl})$ is the matrix of random selection weights, and $p\otimes\mathbb{G}$ is the Hadamard product of the two matrices defined as $(p\otimes\mathbb{G})_{jl}=p_{jl}\mathbb{G}_{jl}$. By letting $\bf 1$ to denote the $m\times 1$ matrix of ones it is that the $j$th element of vector $\mathbb{Q}$ is given by equation (2).

C. Following a univariate construction of geometric slice sets (Fuentes–García et al. 2010), we define the stochastic variables ${\mathbf{N}}=(N_{ji})$ for $1\leq i\leq n_{j}$ and $1\leq j\leq m$, where $N_{ji}$ is an almost surely finite random variable of mass $f_{N}$ possibly depending on parameters, associated with the sequential slice set ${\cal S}_{ji}=\{1,\ldots,N_{ji}\}$. Following Hatjispyros et al. (2011, 2016a) we introduce:

The GSB mixture selection variables ${\bm{\delta}}=(\delta_{ji})$; for an observation $x_{ji}$ that comes from $f_{j}(x)$, $\delta_{ji}$ selects one of the mixtures $\{g_{jl}(x):l=1,\ldots,m\}$. Then the observation $x_{ji}$ came from mixture $g_{j\delta_{ji}}(x)$. 2. 2.

The GSB clustering variables ${\bm{d}}=(d_{ji})$; for an observation $x_{ji}$ that comes from $f_{j}(x)$, given $\delta_{ji}$, $d_{ji}$ allocates the component of the GSB mixture $g_{j\delta_{ji}}(x)$ that $x_{ji}$ came from. Then the observation $x_{ji}$ came from component $K(x|\theta_{j\delta_{ji}d_{ji}})$.

In what follows, unless otherwise specified, the random densities $f_{j}(x)$ are mixtures of independent GSB mixtures.

**Proposition 1. **Suppose that the clustering variables $(d_{ji})$ conditionally on the slice variables $(N_{ji})$ are having the discrete uniform distribution over the sets $({\cal S}_{ji})$ that is $d_{ji}|N_{ji}\sim{\cal DU}({\cal S}_{ji})$, and ${\rm P}\{N_{ji}=r|\delta_{ji}=l\}=f_{N}(r|\lambda_{jl})$, then

[TABLE]

and

[TABLE]

The proof is given in Appendix A.

The following proposition gives a multivariate analogue of equation ($2$) in Fuentes–García, et al. (2010):

**Proposition 2. **Given the random set ${\cal S}_{ji}$, the random functions in (2) become finite mixtures of a.s. finite equally weighted mixtures of the $K(\,\cdot\,|\,\cdot\,)$ probability kernels, that is

[TABLE]

where the probability weights $\{{\cal W}(r|\lambda_{jl}):1\leq l\leq m\}$ are given by

[TABLE]

The proof is given in Appendix A.

Note that, the one–dimensional model introduced in Fuentes–García et al. (2010), under our notation attains the representation

[TABLE]

2.1 The model. Marginalizing (6) with respect to the variable $(N_{ji},d_{ji})$, we obtain

[TABLE]

The quantity inside the parentheses on the right-hand side of the previous equation is $f_{j}(d_{ji}|\delta_{ji}=l)$. Following Fuentes–García, et al. (2010), we substitute $f_{N}(r|\lambda_{jl})$ with the negative binomial distribution ${\cal NB}(r|2,\lambda_{jl})$, i.e.

[TABLE]

so equation (8) becomes

[TABLE]

and the $f_{j}$ random densities take the form of a finite mixture of GSB mixtures

[TABLE]

We denote the set of observations along the $m$ groups as ${\bm{x}}=(x_{ji})$ and with ${\bm{x}}_{j}$ the set of observations in the $j$th group. The three sets of latent variables in the $j$th group will be denoted as ${\bm{N}}_{j}$ for the slice variables, ${\bm{d}}_{j}$ for the clustering variables, and finally ${\bm{\delta}}_{j}$ for the set of GSB mixture allocation variables. From now on, we are going to leave the auxiliary variables unspecified; especially for $\delta_{ji}$ we use the notation $\delta_{ji}=(\delta_{ji}^{1},\ldots,\delta_{ji}^{m})\in\left\{\mathbf{e}_{1},\ldots,\mathbf{e}_{m}\right\}\,\,\,{\rm with}\,\,\,{\rm P}\{\delta_{ji}=\mathbf{e}_{l}\}=p_{jl},$ where $\mathbf{e}_{l}$ denotes the usual basis vector having its only nonzero component equal to $1$ at position $l$. Hence, for a sample of size $n_{1}$ from $f_{1}$, a sample of size $n_{2}$ from $f_{2}$, etc., a sample of size $n_{m}$ from $f_{m}$ we can write the full likelihood as a multiple product:

[TABLE]

In a hierarchical fashion, using the auxiliary variables, we have for $j=1,\ldots,m\text{ and }i=1,\ldots,n_{j},$

[TABLE]

2.2 The PDGSBP covariance and correlation. In this sub–section we find the covariance and the correlation between $f_{j}(x)$ and $f_{i}(x)$. First we provide the following lemma.

**Lemma 1. **Let $g_{\mathbb{G}}(x)=\int_{\Theta}K(x|\theta)\mathbb{G}(d\theta)$ be a random density, with $\mathbb{G}=\lambda\sum_{j=1}^{\infty}(1-\lambda)^{j-1}\delta_{\theta_{j}}$ and $\theta_{j}\stackrel{{\scriptstyle\rm iid}}{{\sim}}G_{0}$, then

[TABLE]

The proof is given in Appendix A.

**Proposition 3. **It is that

[TABLE]

with

[TABLE]

The proof is given in Appendix A.

Suppose now that $(f_{j}^{\cal D}(x))_{1\leq j\leq m}$ and $(f_{j}^{\cal G}(x))_{1\leq j\leq m}$ are two collections of $m$ DP and $m$ GSB pairwise dependent random densities respectively, i.e. $f_{j}^{\cal D}(x)=\sum_{l=1}^{m}p_{jl}g_{jl}^{\cal D}(x)$ with $g_{jl}^{\cal D}(x)=g_{jl}(x|\mathbb{P}_{jl})$, and $f_{j}^{\cal G}(x)=\sum_{l=1}^{m}p_{jl}g_{jl}^{\cal G}(x)$ with $g_{jl}^{\cal G}(x)=g_{jl}(x|\mathbb{G}_{jl})$. Then we have the following proposition:

**Proposition 4. **For given parameters $(\lambda_{ji})$, $(c_{ji})$, and matrix of selection probabilities $(p_{ji})$ it is that

The PDGSBP and rPDDP correlations are given by

[TABLE]

and

[TABLE] 2. 2.

When $\lambda_{ji}=\lambda$ and $c_{ji}=c$ for all $1\leq j\leq i\leq m$, the expressions for the rPDDP and PDGSBP correlations simplify to

[TABLE]

The proof is given in Appendix A.

It is clear that, irrespective of the model, the random densities $f_{j}(x)$ and $f_{i}(x)$ are positively correlated whenever $p_{ji}=p_{ij}=1$. Similarly, the random densities $f_{j}(x)$ and $f_{i}(x)$ are independent (have no common part) whenever $p_{ji}=p_{ij}=0$. Another, less obvious feature, upon synchronization, is the ability of controlling the correlation among the models. For example, suppose that for $m=2$, the random densities $f_{1}(x)$ and $f_{2}(x)$ are dependent, and that $\lambda_{ji}=(1+c_{ji})^{-1}$; then consider the expression

[TABLE]

Since correlations are positive, $D_{12}\geq 0$ whenever ${\rm Corr}(f_{1}^{\cal G}(x),f_{2}^{\cal G}(x))\leq{\rm Corr}(f_{1}^{\cal D}(x),f_{2}^{\cal D}(x)),$ and that $D_{12}<0$ whenever ${\rm Corr}(f_{1}^{\cal G}(x),f_{2}^{\cal G}(x))>{\rm Corr}(f_{1}^{\cal D}(x),f_{2}^{\cal D}(x))$. Then, it is not difficult to see that

[TABLE]

with $r_{k}=(2-\lambda_{12})/(2-\lambda_{kk})$, $k=1,2$. We have the following cases:

$\lambda_{12}>\max\{\lambda_{11},\lambda_{22}\}\,\Leftrightarrow\,r_{1}<1,r_{2}<1\,\Leftrightarrow\,{\rm Corr}(f_{1}^{\cal G}(x),f_{2}^{\cal G}(x))>{\rm Corr}(f_{1}^{\cal D}(x),f_{2}^{\cal D}(x))$. 2. 2.

$\lambda_{12}<\min\{\lambda_{11},\lambda_{22}\}\,\Leftrightarrow\,r_{1}>1,r_{2}>1\,\Leftrightarrow\,{\rm Corr}(f_{1}^{\cal G}(x),f_{2}^{\cal G}(x))<{\rm Corr}(f_{1}^{\cal D}(x),f_{2}^{\cal D}(x))$. 3. 3.

$\lambda_{12}=\lambda_{11}=\lambda_{22}\,\Leftrightarrow\,r_{1}=r_{2}=1\,\Leftrightarrow\,{\rm Corr}(f_{1}^{\cal G}(x),f_{2}^{\cal G}(x))={\rm Corr}(f_{1}^{\cal D}(x),f_{2}^{\cal D}(x))$.

3. The PDGSBP Gibbs sampler. In this section we will describe the PDGSBP Gibbs sampler for estimating the model. The details for the sampling algorithm of the PDDP model can be found in Hatjispyros et al. (2011, 2016A). At each iteration we will sample the variables,

[TABLE]

with $N^{*}=\max_{j,i}N_{ji}$ being almost surely finite.

1. For the locations of the random measures for $k=1,\ldots,d^{*}$ where $d^{*}=\max_{j,i}d_{ji}$, it is that

[TABLE]

If $N^{*}>d^{*}$ we sample additional locations $\theta_{jl,d^{*}+1},\ldots,\theta_{jl,N^{*}}$ independently from the prior.

2. Here we sample the allocation variables $d_{ji}$ and the mixture component indicator variables $\delta_{ji}$ as a block. For $j=1,\ldots,m$ and $i=1,\ldots,n_{j}$, we have

[TABLE]

3. The slice variables $N_{ji}$ have full conditional distributions given by

[TABLE]

which are truncated geometric distributions over the set $\{l,l+1,\ldots\}$.

4. The full conditional for $j=1,\ldots,m$ for the selection probabilities ${\bm{p}}_{j}=(p_{j1},\ldots,p_{jm})$, under a Dirichlet prior $f({\bm{p}}_{j}\,|\,{\bm{a}}_{j})\propto\prod_{l=1}^{m}p_{jl}^{a_{jl}-1}$, with hyperparameter ${\bm{a}}_{j}=(a_{j1},\dots,a_{jm})$, is Dirichlet

[TABLE]

5. Here we update the geometric probabilities $(\lambda_{jl})$ of the GSB measures. For $1\leq j\leq l\leq m$, it is that

[TABLE]

To complete the model, we assign priors to the geometric probabilities. For a fair comparison of the execution time between the two models, we apply $\lambda_{jl}=(1+c_{jl})^{-1}$ transformed priors. So, by placing gamma priors $c_{jl}\sim{\cal G}(a_{jl},b_{jl})$ over the concentration masses $c_{jl}$ of the PDDP model, we have

[TABLE]

In the Appendix, we give the full conditionals for $\lambda_{jl}$’s, their corresponding embedded Gibbs sampling schemes, and the sampling algorithm for the concentration masses.

3.1 The complexity of the rPDDP and PDGSBP samplers. The main difference between the two samplers in terms of execution time, comes from the blocked sampling of the clustering and the mixture indicator variables $d_{ji}$ and $\delta_{ji}$.

The rPDDP model: The state space of the variable $(d_{ji},\delta_{ji})$ conditionally on the slice variable $u_{ji}$ is $(d_{ji},\delta_{ji})(\Omega)=\cup_{l=1}^{m}\left(A_{w_{jl}}(u_{ji})\times\{\mathbf{e}_{l}\}\right),$ where $A_{w_{jl}}(u_{ji})=\{r\in{\mathbb{N}}:u_{ji}<w_{jlr}\}$ is the a.s. finite slice set corresponding to the observation $x_{ji}$ (Walker, 2007). At each iteration of the Gibbs sampler, we have $m(m+1)/2$ vectors of stick-breaking weights $\mathbf{w}_{jl}$, each of length $N_{jl}^{*}$; where $N_{jl}^{*}\sim 1+{\rm Poisson}(-c_{jl}\log u_{jl}^{*})$ with $c_{jl}$ being the concentration parameter of the Dirichlet process $\mathbb{P}_{jl}$ and $u_{jl}^{*}$ being the minimum of the slice variables in densities $f_{j}$ and $f_{l}.$ Algorithm $1$ gives the blocked sampling procedure of the clustering and mixture indicator variables. An illustration of the effect of the slice variable $u_{ji}$ is given in Figure 1(a).

Since the weights forming the stick-breaking representation are not in an ordered form, the construction of the slice sets in step 5 of Algorithm 1 requires a complete search in the array where the weights are stored. This operation is done in ${\cal O}(N_{jl}^{*})$ time. For the sampling of the $d_{ji}$ and $\delta_{ji}$ variables in step 6, the choice of their value is an element from the union $\cup_{l=1}^{m}\left(A_{w_{jl}}(u_{ji})\times\{\mathbf{e}_{l}\}\right).$ This means that the rPDDP algorithm for each $j$, must create $m$ slice sets which require $N_{jl}^{*}$ comparisons each. The worst case scenario is that the sampled $(d_{ji},\delta_{ji})$ is the last element of $\cup_{l=1}^{m}\left(A_{w_{jl}}(u_{ji})\times\{\mathbf{e}_{l}\}\right)$. Thus, the DP based procedure of sampling $(d_{ji},\delta_{ji})$ is of order

[TABLE]

The PDGSBP model: The state space of the variable $(d_{ji},\delta_{ji})$ conditionally on the slice variable $N_{ji}$ is $(d_{ji},\delta_{ji})(\Omega)=\cup_{l=1}^{m}\left({\cal S}_{ji}\times\{\mathbf{e}_{l}\}\right).$ In the GSB case, the slice variable has a different rôlee. It indicates at which random point the search for the appropriate $d_{ji}$ will stop. In Figure 1(b) we illustrate this argument. In Algorithm 2 the worst case scenario is that the sampled $(d_{ji},\delta_{ji})$ will be the last element of $\cup_{l=1}^{m}\left({\cal S}_{ji}\times\{\mathbf{e}_{l}\}\right)$. Thus, the GSB based procedure of sampling $(d_{ji},\delta_{ji})$ is of order $\mathcal{O}\left(m^{2}n_{j}N_{jl}\right)=\mathcal{O}\left(N_{jl}\right).$

4. Illustrations. In this section we illustrate the efficiency of the PDGSBP model. For the choice of a normal kernel (unless otherwise specified) $K(x|\theta)=\mathcal{N}(x|\theta)$ where $\theta=(\mu,\tau^{-1})$ and $\tau=\sigma^{-2}$ is the precision. The prior over the means and precisions of the PDGSBP ($G_{0}$) and the rPDDP model ($P_{0}$) is the independent normal-gamma measure, given by

[TABLE]

Attempting a noninformative prior specification (unless otherwise specified), we took $\mu_{0}=0$ and $\tau_{0}=\epsilon_{1}=\epsilon_{2}=10^{-3}$. For the concentration masses of the rPDDP model, a-priori, we set $c_{jl}\sim{\cal G}(a_{jl},b_{jl})$. For an objective evaluation of the execution time, of the two algorithms under different scenarios, we choose a synchronized prior specification, namely, for the geometric probabilities, we set $\lambda_{jl}\sim{\cal TG}(a_{jl},b_{jl})$ – the transformed gamma density given in equation (14). In the appendix B, we show that such prior specifications are valid for $a_{jl}>1$. In all our numerical examples, we took $a_{jl}=b_{jl}=1.1$. For our numerical experiments (unless otherwise specified), the hyperparameters $(\alpha_{jl})$ of the Dirichlet priors over the matrix of the selection probabilities $p=(p_{jl})$ has been set to $\alpha_{jl}=1$.

In all cases, we measure the similarity between probability distributions with the Hellinger distance. So for example, ${\cal H}_{\cal G}(f,\hat{f})$ and ${\cal H}_{\cal D}(f,\hat{f})$, will denote the Hellinger distance between the true density $f$ and the predictive density $\hat{f}$ of the PDGSBP and rPDDP algorithms, respectively. The Gibbs samplers run for $11\times 10^{4}$ iterations leaving the first $10^{4}$ samples as a burn-in period.

4.1 Time execution efficiency of the PDGSBP model.

Nested normal mixtures with a unimodal common and idiosyncratic part: Here, we choose to include all pairwise and idiosyncratic dependences in the form of unimodal equally weighted normal mixture components. The mixture components are well separated with unit variance. We define each data model ${\cal M}_{m}=\{f_{j}^{(m)}:1\leq j\leq m\}$ of dimension $m\in\{2,3,4\}$, based on a $4\times 10$ matrix $M=(M_{jk})$, with entries in the set $\{0,1\}$, having at most two ones in each column and exactly four ones in each row. When there is exactly one entry of one, the column defines an idiosyncratic part. The appearance of exactly two ones in a column defines a common component. We let the matrix $M$ given by

[TABLE]

and for $m\in\{2,3,4\}$, we define

[TABLE]

We are taking independently samples of sizes $n_{j}^{(2)}=60$ from the $f_{j}^{(2)}$’s, $n_{j}^{(3)}=120$ from the $f_{j}^{(3)}$’s, and, $n_{j}^{(4)}=200$ from the $f_{j}^{(4)}$’s. In all cases, the PDGSBP and the rPDDP density estimations are of the same quality.

In Figures 2(a)–(d) we give the histograms of the data sets for the specific case $m=4$, which are overladed with the kernel density estimations (KDE’s) based on the predictive samples of the $f_{j}^{(4)}$’s coming from the PDGSBP (solid line) and the rPDDP (dashed line) models. The differences between the two models are nearly indistinguishable. The Hellinger distances between the true and the estimated densities for the case $m=4$ are given in table 1.

In Table 2 we summarize the mean execution times (MET’s) per $10^{3}$ iterations in seconds. The PDGSBP sampler is about three times faster than the rPDDP sampler. The corresponding MET ratios for $m=2,3$ and 4 are $2.96,3.04$ and 3.37 respectively. We can see that the PDGSBP Gibbs sampler gives slightly faster execution times with increasing $m$. This will become more clear in our next simulated data example, where the average sample size per mode is being kept constant.

Sparse $m$–scalable data set models: In this example, we attempt to create $m$-scalable normal mixture data sets of the lowest possible sample size. To this respect, we sample independently $m$ groups of data sets from the densities

[TABLE]

with sample sizes $n_{j}^{(m)}=n\{{\cal I}(1\leq j<m)+(m-1)\,{\cal I}(j=m)\}.$ We have chosen $\xi=10$ and an average sample size per mode of $n=20$, for $m\in\{2,\ldots,10\}$.

In Figure 3 we depict the average execution times as functions of the dimension $m$. We can see how fast the two MET-curves diverge with increasing $m$. In Figure 4(a)–(j), for the specific case $m=10$, we give the histograms of the data sets, overladed with the KDE’s based on the predictive samples of the $f_{j}^{(10)}$’s coming from the PDGSBP (solid line) and the rPDDP (dashed line) models. We can see that the PDGSBP and the rPDDP density estimations are of the same quality.

The Hellinger distances between the true and the estimated densities for the specific case $m=10$ are given in Table 3. The large values of the Hellinger distances ${\cal H}_{\cal{G}}(f_{10}^{(10)},\hat{f}_{10}^{(10)})\approx{\cal H}_{\cal{D}}(f_{10}^{(10)},\hat{f}_{10}^{(10)})\approx 0.22$, are caused by the enlargement of the variances of the underrepresented modes due to the small sample size.

4.2 Normal and gamma mixture models that are not well separated.

The normal mixture example: We will first consider a normal model for $m=2$, first appeared in Lijoi et. al (2014B). The data models for $f_{1}$ and $f_{2}$ are 7-mixtures. Their common part is a 4-mixture that is weighted differently between the two mixtures. More specifically, we sample two data sets of sample size $n_{1}=n_{2}=200$, independently from

[TABLE]

with

[TABLE]

For this case, a-priori we took $(\mu_{0},\tau_{0},\epsilon_{1},\epsilon_{2})=(0,10^{-3},1,10^{-2})$.

In Figure 5(a)–(b) we give the histograms of the data sets, with the predictive densities of the PDGSBP and rPDDP models superimposed in black solid and black dashed curves, respectively. We can see that the PDGSBP and the rPDDP density estimations are of the same quality. In Table 4, we give the Hellinger distance between the true and the estimated densities

The gamma mixture example: In this example we took $m=2$. The data models for $f_{1}$ and $f_{2}$ are gamma 4-mixtures. The common part is a gamma 2-mixture, weighted identically among the two mixtures. More specifically, we sample two data sets of sample size $n_{1}=n_{2}=160$, independently from

[TABLE]

with

[TABLE]

Because we want to estimate the density of non negative observations, we find it more appropriate to take the kernel to be a log-normal distribution (Hatjispyros et al. 2016B). That is $K(x|\theta)=\mathcal{LN}(x|\theta)$ with $\theta=(\mu,\sigma^{2})$, is the log-normal density with mean $\exp(\mu+\sigma^{2}/2)$. For this case, a-priori we set

[TABLE]

In Figure 6(a)-(b), we display the KDE’s based on the predictive samples of the two models. We can see that the PDGSBP and the rPDDP density estimations are of the same quality. In Table 5, we give the Hellinger distances.

Because the common part is equally weighted among $f_{1}$ and $f_{2}$, it makes sense to display the estimations of the selection probability matrices under the two models

[TABLE]

4.3 Borrowing of strength of the PDGSBP model. In this example we consider three populations $\{D_{j}^{(s)}:j=1,2,3\}$, under three different scenarios $s\in\{1,2,3\}$. The sample sizes are always the same, namely, $n_{1}=200$, $n_{2}=50$ and $n_{3}=200$ – the second population is sampled only once. The three data sets $D_{1}^{(s)}$, $D_{2}^{(s)}$ and $D_{3}^{(s)}$, are sampled independently from the normal mixtures

[TABLE]

where

[TABLE]

More specifically, the three scenarios are:

For $s=1$, we set, $q^{(1)}=0$. This is the case where the three populations are coming from the same 4–mixture $f$. We depict the density estimations under the first scenario in Figures 7(a)–(c). This is the case where the small data set, benefits the most in terms of borrowing of strength. 2. 2.

For $s=2$, we set, $q^{(2)}=1/2$. The 2-mixtures $g_{1}$ and $g_{2}$ are the the idiosyncratic parts of the 6-mixtures $f_{1}^{(2)}$ and $f_{3}^{(2)}$, respectively. The density estimations under the second scenario are given in Figures 7(d)–(f). In this case, the strength of borrowing between the small data set and the two large data sets weakens. 3. 3.

For $s=3$ we set $q^{(3)}=1$. In this case the three populations have no common parts. The density estimations are given in Figures 7(g)–(i). This is the worst case scenario, where there is no borrowing of strength between the small and the two large data sets.

The Hellinger distances between the true and the estimated densities, for the three scenarios, are given in table 6. In the second column of the Table we can see how the Hellinger distance of the estimation $\hat{f}_{2}^{(s)}$ and the true density $f_{2}^{(s)}$ increases as the borrowing of strength weakens, it is that ${\cal H}_{\cal G}(f_{2}^{(1)},\hat{f}_{2}^{(1)})<{\cal H}_{\cal G}(f_{2}^{(2)},\hat{f}_{2}^{(2)})<{\cal H}_{\cal G}(f_{2}^{(3)},\hat{f}_{2}^{(3)}).$

4.4 Real data example. The data set is to be found at http://lib.stat.cmu.edu/datasets/pbcseq and involves data from 310 individuals. We take the observation as SGOT (serum glutamic-oxaloacetic transaminase) level, just prior to liver transplant or death or the last observation recorded, under three conditions on the individual

The individual is dead without transplantation. 2. 2.

The individual had a transplant. 3. 3.

The individual is alive without transplantation.

We normalize the means of all three data sets to zero. Since it is reasonable to assume the densities for the observations are similar for the three categories (especially for the last two), we adopt the models proposed in this paper with $m=3$. The number of transplanted individuals is small (sample size of 28) so it is reasonable to borrow strength for this density from the other two. In this example, we set the hyperparameters of the Dirichlet priors for the selection probabilities to

[TABLE]

In Figure 8(a)–(c) we provide histograms of the real data sets and superimpose the KDE’s based on the predictive samples of the PDDP and PDGSBP samplers. The two models give nearly identical density estimations. 2. 2.

The estimated a-posteriori selection probabilities are given below

[TABLE]

By comparing the second rows of the selection matrices, we conclude that the strength of borrowing is slightly larger in the case of PDGSBP model .

5. Discussion. In this paper we have generalized the GSB process to a multidimensional dependent stochastic process which can be used as a Bayesian nonparametric prior for density estimation in the case of partially exchangeable data sets. The resulting Gibbs sampler is as accurate as its DP based counterpart, yet faster and far less complicated. The main reason for this is that the GSB sampled value of the allocation variable $d_{ji}$ will be an element of the sequential slice set ${\cal S}_{ji}=\{1,\ldots,N_{ji}\}$. Thus, there is no need to search the arrays of the weights; we know the state space of the clustering variables in advance. On the other hand, the sampling of $d_{ji}$ in the DP based algorithm will always have one more step; the creation of the slice sets.

For an objective comparison of the execution times of the two models, we have run the two samplers in an a-priori synchronized mode. This, involves the placing of ${\cal G}(a_{jl},b_{jl})$ priors over the DP $c_{jl}$ concentration masses, leading to a more efficient version of the PDDP model introduced in Hatjispyros et al. (2011, 2016A).

We have show that when the PDGSBP and PDDP models are synchronized, i.e. their parameters satisfy $\lambda_{ji}=(1+c_{ji})^{-1}$, the correlation between the models can be controlled by imposing further restrictions among the $\lambda_{ji}$ parameters.

Finally, an interesting research path would be the generalization of the pairwise dependent $\mathbb{Q}_{j}$ measures to include all possible interactions, in the sense that

[TABLE]

where the $\mathbb{G}_{j}$ and the $\mathbb{G}_{\eta_{(j)}}$’s are independent GSB processes, ${\cal C}_{j,l,m}=\{(k_{1},\ldots,k_{l-1}):1\leq k_{1}<\cdots<k_{l-1}\leq m,k_{r}\neq j,1\leq r\leq m-1\}$ and $\eta_{(j)}$ is the ordered vector of the elements of the vector $\eta$ and $\{j\}$. Now the $f_{j}$ densities will be a mixture of $2^{m-1}$ GSB mixtures, and the total number of the independent GSB processes needed to model $(f_{1},\ldots,f_{m})$ will be $2^{m}-1$.

Appendix A

Proof of Proposition 1. Starting from the $N_{ji}$-augmented random densities we have

[TABLE]

Because $f_{j}(N_{ji}=r|\delta_{ji}=l)=f_{N}(r|\lambda_{jl})$ and $f_{j}(x_{ji}|d_{ji}=k,\delta_{ji}=l)=K(x_{ji}|\theta_{jlk})$, the last equation gives

[TABLE]

Augmenting further with the variables $d_{ji}$ and $\delta_{ji}$ yields

[TABLE]

Because P$(\delta_{ji}=l)=p_{jl}$, the last equation leads to equation (6) and the proposition follows. $\square$

Proof of Proposition 2. Marginalizing the joint of $x_{ji}$ and $N_{ji}$ with respect to $x_{ji}$ we obtain

[TABLE]

Then dividing equation (5) with the probability that $N_{ji}$ equals $r$ we obtain equation (7). $\square$

Proof of Lemma 1. Because $g_{\mathbb{G}}(x)=\lambda\sum_{j=1}^{\infty}(1-\lambda)^{j-1}K(x|\theta_{j})$, we have

[TABLE]

which gives the desired result. $\square$

Proof of Proposition 3. The random densities $f_{i}(x)=\sum_{l=1}^{m}p_{il}\,g_{il}(x)$ and $f_{j}(x)=\sum_{l=1}^{m}p_{jl}\,g_{jl}(x)$ depend to each other through the random measure $\mathbb{G}_{ji}$, therefore

[TABLE]

and

[TABLE]

Substituting back to equation (15) one obtains

[TABLE]

Using lemma $1$, the last equation becomes

[TABLE]

or that

[TABLE]

The desired result, comes from the fact that

[TABLE]

$\square$

**Proof of Proposition 4.

**(1.) From equation (11) and proposition 3, we have that

[TABLE]

Normalizing the covariance in equation (10) with the associated standard deviations, yields

[TABLE]

Similarly, from proposition 1 in Hatjispyros et al. (2011), it is that

[TABLE]

and

[TABLE]

(2.) When $\lambda_{ji}=\lambda$ and $c_{ji}=c$ for all $1\leq j\leq i\leq m$, from equations (16) and (17), it is clear that

[TABLE]

Appendix B

1. Sampling of the concentrations masses for the rPDDP model.

In this case, the random densities $(f_{j})$ are represented as finite mixtures of the DP mixtures $g_{jl}(x|\mathbb{P}_{jl})$, where $\mathbb{P}_{jl}\sim{\cal DP}(c_{jl},P_{0})$. We randomize the concentrations by letting $c_{jl}\sim{\cal G}(a_{jl},b_{jl})$. Following West (1992) we have the following two specific cases:

A. For $j=l$, the posterior $c_{jj}$’s will be affected only by the size of the data set ${\bm{x}}_{j}$ and the number of unique clusters for which $\delta_{ji}={\bf e}_{j}$. Letting

[TABLE]

we have

[TABLE]

with the weights $\pi_{\beta}$ satisfying $\frac{\pi_{\beta}}{1-\pi_{\beta}}=\frac{a_{jj}+\rho_{jj}-1}{n_{j}(b_{jj}-\log\beta)}$.

B. For $j\neq l$, the posterior $c_{jl}$’s will be affected by the size of the data sets ${\bm{x}}_{j}$ and ${\bm{x}}_{l}$ and the cumulative number of unique clusters $d_{ji}$ for which $\delta_{ji}={\bf e}_{l}$ and the unique clusters $d_{li}$ for which $\delta_{li}={\bf e}_{j}$. Letting

[TABLE]

it is that

[TABLE]

with the weights $\pi_{\beta}$ satisfying $\frac{\pi_{\beta}}{1-\pi_{\beta}}=\frac{a_{jl}+\rho_{jl}-1}{(n_{j}+n_{l})(b_{jl}-\log\beta)}$.

Bear in mind that $\rho_{jl}=0$ is always a possibility, so that we impose $a_{jl}>1$.

2. Sampling of the geometric probabilities for the PDGSBP model.

In this section we provide the full conditionals for the geometric probabilities $\lambda_{jl}$ under beta conjugate and transformed gamma nonconjugate priors. We let

[TABLE]

A. For the choice of prior $\lambda_{jl}\sim{\cal B}e(a_{jl},b_{jl})$, for $l=j$ it is that

[TABLE]

also, for $l\neq j$ we have

[TABLE]

B. For the choice of prior $\lambda_{jl}\sim{\cal TG}(a_{jl},b_{jl})$, for $l=j$ it is that

[TABLE]

To sample from this density, we include the positive auxiliary random variables $\nu_{1}$ and $\nu_{2}$ such that

[TABLE]

The full conditionals for $\nu_{1},\nu_{2}$ are uniforms

[TABLE]

and the new full conditional for $\lambda_{jj}$ becomes

[TABLE]

where we have set $L_{jj}=S_{jj}^{\prime}+a_{jj}-1$. We can sample from this density using the inverse cumulative distribution function technique. Also, for $l\neq j$ we apply the same embedded Gibbs sampling technique to the full conditional density

[TABLE]

References.

Bulla, P., Muliere, P. and Walker, S.G. (2009). A Bayesian nonparametric estimator of a multivariate survival function. Journal of Statistical Planning and Inference 139, 3639–3648.

De Iorio, M., Müller, P., Rosner, G.L. and MacEachern, S.N. (2004). An ANOVA model for dependent random measures. Journal of the American Statistical Association 99, 205–215.

Dunson, D.B. and Park, J.H. (2008). Kernel stick–breaking processes. Biometrika 95, 307–323.

Ferguson, T.S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics 1, 209–230.

Fuentes–Garcia, R., Mena, R.H., Walker, S.G. (2009). A nonparametric dependent process for Bayesian regression Statistics and Probability Letters 79, 1112–1119.

Fuentes–Garcia, R., Mena, R.H., Walker, S.G. (2010). A new Bayesian nonparametric mixture model. Comm.Statist.Simul.Comput 39, 669–682.

Griffin, J.E. and Steel, M.F.J. (2006). Order–based dependent Dirichlet processes. Journal of the American Statistical Association 101, 179–194.

Griffin, J.E., Kolossiatis, M. and Steel, M.F.J. (2013). Comparing distributions by using dependent normalized ranom–measure mixtures. Journal of the Royal Statistical Society, Series B 75, 499–529.

Hatjispyros, S.J., Nicoleris, T. and Walker, S.G. (2011). Dependent mixtures of Dirichlet processes. Computational Statistics and Data Analysis 55, 2011–2025.

Hatjispyros, S.J., Nicoleris, T. and Walker, S.G. (2016a). Dependent random density functions with common atoms and pairwise dependence. Computational Statistics and Data Analysis 101, 236–249.

Hatjispyros, S.J., Nicoleris, T. and Walker, S.G. (2016b). Bayesian nonparametric density estimation under length bias. Communications in Statistics

DOI: 10.1080/03610918.2016.1263735

Lijoi, A., Nipoti, B. and Prüenster, I. (2014a). Bayesian inference with dependent normalized completely random measures. Bernoulli, 20, 1260–1291.

Lijoi, A., Nipoti, B. and Prüenster, I. (2014b). Dependent mixture models: clustering and borrowing information. Computational Statistics and Data Analysis 71, 17–433.

Kolossiatis, M., Griffin, J.E. and Steel, M.F.J. (2013). On Bayesian nonparametric modelling of two correlated distributions. Statistics and Computing 23, 1–15.

Lo, A.Y. (1984). On a class of Bayesian nonparametric estimates I. Density estimates. Annals of Statistics 12, 351–357.

MacEachern, S.N. (1999). Dependent nonparametric processes. In “Proceedings of the Section on Bayesian Statistical Science” pp. 50-55. American Statistical Association.

Müller, P., Quintana, F., and Rosner, G., (2004). A method for combining inference across related nonparametric Bayesian models. Journal of the Royal Statistical Society, Series B 66, 735–749.

Mena, R.H., Ruggiero, M. and Walker, S.G. (2011). Geometric stick–breaking processes for continuous–time Bayesian nonparametric modeling. Journal of Statistical Planning and Inference 141 (9), 3217–3230.

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica 4 639–650.

Walker, S.G. (2007). Sampling the Dirichlet mixture model with slices Communications in Statistics 36 45–54.

West, M. (1992). Hyperparameter estimation in Dirichlet process mixture models. Technical report 92-A03, Duke University, ISDS.