Dependent Mixtures of Geometric Weights Priors
Spyridon J. Hatjispyros Christos Merkatas Theodoros Nicoleris Stephen G. Walker






Dependent Mixtures of Geometric Weights Priors




Spyridon J. Hatjispyros
     [Corresponding author. Tel.:+30 22730 82326

     E-mail address: schatz@aegean.gr ]^, *,
      Christos Merkatas^*,
      Theodoros Nicoleris^**,
      Stephen G. Walker^***

  




^* Department of Mathematics, University of the Aegean,
Karlovassi, Samos, GR-832 00, Greece. 
^** Department of Economics, National and Kapodistrian University of Athens,
Athens, GR-105 59, Greece.   
^***Department of Mathematics, University of Texas at Austin,
Austin, Texas 7812, USA. 



A new approach to the joint estimation of partially exchangeable observations is presented. This is achieved by constructing 
a model with pairwise dependence between random density functions, each of which is modeled as a mixture of geometric 
stick breaking processes. 
The claim is that mixture modeling with Pairwise Dependent Geometric Stick 
Breaking Process (PDGSBP) priors is sufficient for prediction and estimation purposes; that is, making the weights more exotic does not actually enlarge the support of the prior. Moreover, the corresponding Gibbs sampler for estimation
is faster and easier to implement than the Dirichlet Process counterpart. 


 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,…,X_n) 
coincides with the distribution of (X_π(1),…,X_π(n)), 
for all π∈ S(n), where S(n) is the set of permutations of  {1,…,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≤ i_j≤ n_j, 1≤ j≤ m, so for all π_j∈ S(n_j)


    ((X_1,i_1)_1≤ i_1≤ n_1,…,(X_m, i_m)_1≤ i_m≤ n_m)∼
    ((X_1,π_1(i_1))_1≤ i_1≤ n_1,…,(X_m, π_m(i_m))_1≤ i_m≤ n_m).


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 (_j)_1≤ j≤ m, where each _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) = ∫_ΘK(x|θ)(dθ); 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

    ( · ) = ∑_k=1^∞w_kδ_θ_k( · ),

where (θ_k)_k≥ 1 are independent and identically distributed from P_0 and (w_k)_k≥ 1 is
a stick breaking process; so if (v_k)_k≥ 1 are independent and identically distributed from 
Be(1, c), a beta distribution with mean (1+c)^-1, then w_1=v_1 and for k>1, w_k=v_k∏_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 
_z is constructed for each covariate z, 

    _z( · ) = ∑_k=1^∞ w_k(z)δ_θ_k(z)( · ),

where (w_k(z),θ_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≥ 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 (_j)_1≤ j≤ m, 
where each _j is a GSB random probability measure, such that there is a unique common component for 
each pair (_j,_j') with j≠ j'. 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≤ j≤ m. 
Hence we obtain a set of random densities (f_1,…,f_m), 
where marginally each f_j is a random density function

    f_j(x) = ∫_ΘK(x|θ) _j(dθ),

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≤ j≤ 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
_j = p_j_0 + (1-p_j)_j^*,
where _0 is the common component to all other distributions and {_j^*:j=1,…,m} are the  
idiosyncratic parts to each _j, and _0,_j^* DP(c,P_0). In Lijoi et al. (2014B)  
normalized random probability measures based on the σ–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≤ j≤ m, i≥ 1}
defined over a probability space (Ω, F,  P), 
that is partially exchangeable as in (<ref>). Let P denote
the set of probability measures over ; then de Finetti proved that there exists
a probability distribution Π over P^m, which satisfies  


    P{X_ji∈ A_ji:1≤ j≤ m,1≤ i≤ n_j}
        = ∫_ P^m P{X_ji∈ A_ji:1≤ j≤ m,1≤ i≤ n_j | _1,…,_m} Π(d_1,…,d_m)
        = ∫_ P^m∏_j=1^m P{X_ji∈ A_ji:1≤ i≤ n_j | _j} Π(d_1,…,d_m)
         = ∫_ P^m∏_j=1^m ∏_i=1^n_j_j(A_ji) Π(d_1,…,d_m) .


The de Finetti measure Π 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 Π, marginal measures _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≤ i≤ n_j:1≤ j≤ m},  

    x_ji|θ_ji    K( · |θ_ji)
    θ_ji|_j    _j( · )
    _j  =    ∑_l=1^mp_jl_jl,  ∑_l=1^mp_jl=1,  _jl=_lj
    _jl      DP(c, P_0),  1≤ j≤ l≤ m,

for some kernel density K( · | · ), concentration parameter c>0 and parametric 
central measure P_0 for which (_jl(dθ))=P_0(dθ).
So, we have assumed that the random densities f_j(x) are dependent mixtures 
of the dependent random measures _j via f_j(x|_j)=∫_Θ K(x|θ)_j(dθ),
or equivalently, dependent mixtures of the m independent mixtures 
g_jl(x| _jl)=∫_Θ K(x| ) _jl(d), l=1,…,m.
To introduce the rPDDP model, we randomize the PDDP model by sampling the 
_jl measures from the independent Dirichlet 
processes DP(c_jl,P_0) and then impose gamma priors on the concentration masses, i.e.
_jl  DP(c_jl, P_0),   c_jl G(a_jl,b_jl),  1≤ j≤ l≤ m.


 B.
To develop a pairwise dependent geometric stick breaking version, 
we let the random density functions f_j(x) generated via

    f_j(x):=f_j(x| _j)  =  ∑_l=1^m p_jl g_j l(x| _jl),  _j  = ∑_l=1^mp_jl_jl,   1≤ j≤ m.

The g_jl(x):=g_jl(x| _jl)=∫_Θ K(x| ) _jl(d) random densities are 
now independent mixtures of GSB processes, satisfying g_jl=g_lj, under the slightly altered definition

    _jl=∑_k=1^∞ q_jlkδ_θ_jlk   with   q_jlk=λ_jl(1-λ_jl)^k-1,  λ_jl∼ h( · |ξ_jl),   θ_jlk G_0,

where h is a parametric density supported over the interval (0,1) depending on some parameter ξ_jl∈Ξ, 
and G_0 is the associated parametric central measure. 

The independent GSB processes {_jl: 1≤ j, l≤ m} form a matrix 
 of random distributions with _jl=_lj.
In matrix notation

    ℚ = (p⊗) 1,

where p=(p_jl) is the matrix of random selection weights, and p⊗ is the Hadamard
product of the two matrices defined as (p⊗)_jl=p_jl_jl. By letting
1 to denote the m× 1 matrix of ones it is that the jth element of vector 
ℚ is given by equation (<ref>).


 C.
Following a univariate construction of geometric slice sets (Fuentes–García et al. 2010), we define
the stochastic variables 𝐍=(N_ji) for 1≤ i≤ n_j and 1≤ j≤ 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 S_ji={1,…,N_ji}. 
Following Hatjispyros et al. (2011, 2016a) we introduce: 


  
  * The GSB mixture selection variables δ=(_̣ji); 
        for an observation x_ji that comes from f_j(x), 
        _̣ji selects one of the mixtures {g_jl(x):l=1,…,m}. 
        Then the observation x_ji came from mixture g_j δ_ji(x).
        
  
  * The GSB clustering variables d=(d_ji); for an observation 
        x_ji that comes from f_j(x), given δ_ji, d_ji allocates the 
        component of the GSB mixture g_j δ_ji(x) that x_ji came from. 
        Then the observation x_ji came from component K(x|_j_̣jid_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 ( S_ji) that is d_ji|N_ji∼ DU( S_ji), and P{N_ji=r|_̣ji=l}=f_N(r|λ_jl), then

    f_j(x_ji,N_ji=r)=r^-1∑_l=1^mp_jlf_N(r|λ_jl)∑_k=1^r K(x_ji|θ_jlk),

and 

    f_j(x_ji,N_ji=r,d_ji=k|δ_ji=l) = 1 rf_N(r|λ_jl)  I(k≤ r) K(x_ji|θ_jlk).


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 S_ji, the random functions in (<ref>) become finite mixtures of a.s. finite equally weighted mixtures 
of the K( · | · ) probability kernels, that is

    f_j(x_ji|N_ji=r)=∑_l=1^m W(r|λ_jl)∑_k=1^rr^-1K(x_ji|θ_jlk),

where the probability weights { W(r|λ_jl):1≤ l≤ m} are given by

    W(r|λ_jl)=p_jl f_N(r|λ_jl)∑_l'=1^m p_jl' f_N(r|λ_jl').


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

    f_j(x_ji|N_ji=r,_̣ji=l)=∑_k=1^rr^-1K(x_ji|θ_jlk).





 2.1 The model. Marginalizing (<ref>) with respect to the variable (N_ji, d_ji), we obtain

    f_j(x_ji|_̣ji=l)=∑_k=1^∞(∑_r=k^∞
        r^-1f_N(r|λ_jl))K(x_ji|θ_jlk).

The quantity inside the parentheses on the right-hand side of the previous equation is f_j(d_ji|_̣ji=l).
Following Fuentes–García, et al. (2010), we substitute f_N(r|ł_jl) with the negative binomial distribution 
NB(r|2,ł_jl), i.e. 

    f_N(r|ł_jl)  = r ł_jl^2(1-ł_jl)^r-1 I(r≥ 1),

so equation (<ref>) becomes

    f_j(x_ji|_̣ji=l)=∑_k=1^∞ q_jlk K(x_ji|θ_jlk)    with   q_jlk=ł_jl(1-ł_jl)^k-1,
 
and the f_j random densities take the form of a finite mixture of GSB mixtures

    f_j(x_ji)=∑_l=1^m p_jl∑_k=1^∞ q_jlk K(x_ji|θ_jlk).

We denote the set of observations along the m groups as x=(x_ji) and with 
x_j the set of observations in the jth group.
The three sets of latent variables in the jth group will be denoted as N_j for the slice variables, 
d_j for the clustering variables, and finally _j for the set of GSB mixture allocation variables. 
From now on, we are going to leave the auxiliary variables unspecified; especially for _̣ji we use the notation
_̣ji=(_̣ji^1,…,_̣ji^m)∈{e⃗_1,…,e⃗_m}    with    P{_̣ji=e⃗_l}=p_jl,
where 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:

    f( x, N, d | ) 
               =    ∏_j=1^m  f( x_j, N_j, d_j | _j)
        =    ∏_j=1^m∏_i=1^n_j I(d_ji≤ N_ji)∏_l=1^m
                 {ł_jl^2(1-ł_jl)^N_ji-1K(x_ji| θ_j l d_ji)}^δ_ji^l.

In a hierarchical fashion, using the auxiliary variables, we have for j=1,…,m  and  i = 1,…, n_j,

    x_ji, N_ji |  d_ji, _̣ji, (_jr_̣ji)_1≤ r≤ m,λ_j_̣ji  ∏_r=1^m{ł_jr^2(1-ł_jr)^N_ji-1K(x_ji|θ_jr d_ji)}^δ_ji^r I(N_ji≥ d_ji)
        d_ji | N_ji DU( S_ji),    P{δ_ji=e⃗_l} = p_jl
        q_jik=λ_ji(1-λ_ji)^k-1,   _jik G_0,   k∈.





 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_(x)=∫_Θ K(x|)(d) be a random density, with 
=ł∑_j=1^∞ (1-ł)^j-1_̣_j and _j G_0, then

    [g_(x)^2] = (1 2-ł){ł∫_Θ K(x|θ)^2G_0(dθ) 
     + 2(1-ł)(∫_Θ K(x|θ)G_0(dθ))^2}.


The proof is given in Appendix A.



 Proposition 3. It is that

    Cov(f_j(x),f_i(x))  =  p_ji p_ij Var(∫_Θ K(x|)_ji(d)),

with

    Var(∫_Θ K(x|)_ji(d))=ł_ji 2-ł_ji Var(K(x|)).


The proof is given in Appendix A.

Suppose now that (f_j^ D(x))_1≤ j≤ m and 
(f_j^ G(x))_1≤ j≤ m are two collections
of m DP and m GSB pairwise dependent random densities respectively, i.e.
f_j^ D(x)=∑_l=1^m p_jlg_jl^ D(x) with g_jl^ D(x)=g_jl(x|_jl), and
f_j^ G(x)=∑_l=1^m p_jlg_jl^ G(x) with g_jl^ G(x)=g_jl(x|_jl).
Then we have the following proposition:


 Proposition 4. For given parameters (ł_ji), (c_ji), and matrix of selection
probabilities (p_ji) it is that


  * The PDGSBP and rPDDP correlations are given by  

    Corr(f_j^ G(x),f_i^ G(x))  = λ_jip_jip_ij 2-λ_ji(∑_l=1^m∑_r=1^m  p_jl^2 p_ir^2λ_jlλ_ir (2-λ_jl)(2-λ_ir))^-1/2,

and 

    Corr(f_j^ D(x),f_i^ D(x)) = p_jip_ij 1+c_ji(∑_l=1^m∑_r=1^m  p_jl^2 p_ir^2 (1+c_jl)(1+c_ir))^-1/2.



  * When λ_ji=λ and c_ji=c for all 1≤ j≤ i≤ m, 
the expressions for the rPDDP and PDGSBP correlations simplify to

    Corr(f_j^ G(x),f_i^ G(x))= Corr(f_j^ D(x),f_i^ D(x))
     = p_jip_ij( ∑_l=1^m ∑_r=1^m  p_jl^2 p_ir^2 )^-1/2.



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 ł_ji=(1+c_ji)^-1; then consider the expression

    D_12:=ł_12^2  p_12^2  p_21^2 { Corr(f_1^ G(x),f_2^ G(x))^-2
    - Corr(f_1^ D(x),f_2^ D(x))^-2}.

Since correlations are positive,
D_12≥ 0 whenever
Corr(f_1^ G(x),f_2^ G(x))≤ Corr(f_1^ D(x),f_2^ D(x)), and that
D_12<0 whenever Corr(f_1^ G(x),f_2^ G(x))>
 Corr(f_1^ D(x),f_2^ D(x)). Then, it is not difficult to see that

    D_12=(p_12^2ł_12+r_1p_11^2ł_11)
           (p_21^2ł_12+r_2 p_22^2ł_22)
          -(p_12^2ł_12+p_11^2ł_11)
           (p_21^2ł_12+p_22^2ł_22)

with r_k=(2-ł_12)/(2-ł_kk), k=1,2. We have the following cases:


  * ł_12>max{ł_11,ł_22} ⇔ r_1<1, r_2<1 ⇔  Corr(f_1^ G(x),f_2^ G(x))> Corr(f_1^ D(x),f_2^ D(x)).

  * ł_12<min{ł_11,ł_22} ⇔ r_1>1, r_2>1 ⇔  Corr(f_1^ G(x),f_2^ G(x))< Corr(f_1^ D(x),f_2^ D(x)).

  * ł_12=ł_11=ł_22 ⇔ r_1=r_2=1 ⇔  Corr(f_1^ G(x),f_2^ G(x))= Corr(f_1^ D(x),f_2^ 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,

    θ_jlk, 1≤ j ≤ l ≤ m, 1≤ k ≤ N^*,
        d_ji,N_ji,δ_ji, 1≤ j ≤ m, 1≤ i ≤ n_j,
        p_jl, 1≤ j ≤ m, 1≤ l ≤ m,

with N^*=max_j,iN_ji being almost surely finite. 

1. 
For the locations of the random measures for k=1,…,d^* where d^*=max_j,id_ji, it is that

    f(θ_jlk|⋯) ∝ f(θ_jlk)
    ∏_i=1^n_jK(x_ji|θ_jlk)^ I(_̣ji=e⃗_l, d_ji=k)∏_i=1^n_lK(x_li|θ_jlk)^ I(_̣li=e⃗_j, d_li=k)      l>j ,
    ∏_i=1^n_jK(x_ji|θ_jjk)^ I(_̣ji=e⃗_j, d_ji=k)      l=j .

If N^*>d^* we sample additional locations _jl,d^*+1,…,_jl,N^* independently from the prior. 

2. 
Here we sample the allocation variables d_ji and the mixture 
component indicator variables δ_ji as a block. For j=1,…,m and i=1,…,n_j, we have

    P(d_ji=k,_̣ji=e⃗_l |N_ji=r,⋯) ∝  p_jl K(x_ji|θ_jlk)  I(l≤ m)  I(k≤ r).



3. 
The slice variables N_ji have full conditional distributions given by

    P(N_ji = r | _̣ji=e⃗_l,d_ji=l,⋯)∝(1-λ_jl)^r  I(r≥ l),

which are truncated geometric distributions over the set {l, l+1,…}.


4. 
The full conditional for j=1,…,m for the selection probabilities p_j=(p_j1,…,p_jm), 
under a Dirichlet prior f( p_j |  a_j)∝∏_l=1^m p_jl^a_jl-1, with hyperparameter
a_j=(a_j1,…,a_jm), is Dirichlet

    f( p_j |⋯) ∝ ∏_l=1^m p_jl^a_jl+∑_i=1^n_l I(δ_ji = e⃗_l) - 1.



5. 
Here we update the geometric probabilities (ł_jl) of the GSB measures. For 1≤ j≤ l≤ m, it is that

    f(ł_jl|⋯) ∝ f(ł_jl)
    ∏_i=1^n_j{ł_jl^2(1-ł_jl)^N_ji-1}^ I(_̣ji=e⃗_l)∏_i=1^n_l{ł_jl^2(1-ł_jl)^N_li-1}^ I(_̣li=e⃗_j)      l>j
    ∏_i=1^n_j{ł_jj^2(1-ł_jj)^N_ji-1}^ I(_̣ji=e⃗_j)      l=j .

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 ł_jl=(1+c_jl)^-1
transformed priors. So, by placing gamma priors c_jl∼ G(a_jl,b_jl) 
over the concentration masses c_jl of the PDDP model, we have

    f(ł_jl)= TG(ł_jl | a_jl,b_jl)∝ł_jl^-(a_jl+1)e^-b_jl/ł_jl(1-ł_jl)^a_jl-1  I(0<ł_jl<1).

In the Appendix, we give the full conditionals for ł_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 _̣ji.

The rPDDP model:
The state space 
of the variable (d_ji,δ_ji) conditionally on the slice variable u_ji is
(d_ji,δ_ji)(Ω)=
∪_l=1^m(A_w_jl(u_ji) ×{e⃗_l}),
where A_w_jl(u_ji)={r∈ℕ: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 w⃗_jl, each of length N_jl^*; 
where N_jl^*∼ 1 +  Poisson(-c_jllog u_jl^*) with c_jl being the concentration parameter of the Dirichlet 
process ℙ_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 <ref> requires a complete search in the array where the weights are stored. This operation 
is done in O(N_jl^*) time. For the sampling of the d_ji and _̣ji variables in step 6, the choice of 
their value is an element from the union ∪_l=1^m(A_w_jl(u_ji)×{e⃗_l}). 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,_̣ji) is the last element of ∪_l=1^m(A_w_jl(u_ji)×{e⃗_l}). 
Thus, the DP based procedure of sampling (d_ji,_̣ji) is of order 

    O( m^2 n_j N_jl^* ∑_l=1^m|A_w_jl(u_ji)| ) 
    =  O(N_jl^*∑_l=1^m|A_w_jl(u_ji)|).





The PDGSBP model:
The state space 
of the variable (d_ji,δ_ji) conditionally on the slice variable N_ji is
(d_ji, _̣ji)(Ω)=∪_l=1^m( S_ji×{e⃗_l}).
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 <ref> the worst case scenario is that the sampled (d_ji,_̣ji) will be the last element 
of ∪_l=1^m( S_ji×{e⃗_l}). 
Thus, the GSB based procedure of sampling (d_ji,_̣ji) is of order  
𝒪(m^2n_jN_jl) = 𝒪(N_jl).










 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|θ) = 𝒩(x|θ) where θ=(μ,τ^-1) and τ=σ^-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

    P_0(dμ,dτ)=G_0(dμ,dτ)=  N(μ | μ_0,τ_0^-1)
      G(τ | ϵ_1,ϵ_2) dμ dτ.

Attempting a noninformative prior specification (unless otherwise specified),
we took μ_0=0 and τ_0=ϵ_1=ϵ_2=10^-3.
For the concentration masses of the rPDDP model, a-priori, we set c_jl∼ 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 ł_jl∼ TG(a_jl,b_jl) – the transformed gamma density given in equation (<ref>). 
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 (_jl) of the Dirichlet 
priors over the matrix of the selection probabilities p=(p_jl) has been set to _jl=1.

In all cases, we measure the similarity between probability distributions with the Hellinger distance.
So for example, H_ G(f,f̂) and H_ D(f,f̂), will denote the Hellinger 
distance between the true density f and the predictive density f̂ of the PDGSBP and rPDDP algorithms,
respectively. The Gibbs samplers run for 11× 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 
M_m={f_j^(m):1≤ j≤ m} of dimension m∈{2,3,4}, based on a 
4× 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  

    M=
      [ 1 1 1 1 0 0 0 0 0 0; 0 0 1 0 1 0 0 1 1 0; 0 1 0 0 0 1 0 1 0 1; 1 0 0 0 0 0 1 0 1 1;   ],

and for m∈{2,3,4}, we define 

    M_m: f_j^(m)(x)∝∑_k=5-m^2(m+1)M_jk  N(x|10(k-6),1),  1≤ j≤ m,

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

    f_j^(m)(x) ∝  N(x|(j-1)ξ,1)  I(1≤ j<m)+
    ∑_k=1^m-1 N(x|(k-1) ξ,1)  I(j=m),

with sample sizes 
n_j^(m)=n{ I(1≤ j<m)+(m-1)  I(j=m)}.
We have chosen ξ=10 and an average sample size per mode of n=20, for m∈{2,…,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
H_G(f_10^(10),f̂_10^(10))≈ H_D(f_10^(10),f̂_10^(10))≈ 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 

    (f_1,f_2)=(1 2  g_11 + 1 2  g_12,  4 7  g_21 + 3 7 g_22),

with

    g_11   = 2/7 N(-8,0.25^2) + 3/7 N(1,0.5^2) + 2/7 N(10,1) 
    
    g_12   = 1/7 N(-10,0.5^2) + 3/7 N(-3,0.75^2) + 1/7 N(3,0.25^2) + 2/7 N(7, 0.25^2)
    
    g_21   = 2/8 N(-10,0.5^2) + 3/8 N(-3,0.75^2) + 2/8 N(3,0.25^2) + 1/8 N(7, 0.25^2)
    
    g_22   = 1/3 N(-6,0.5^2) + 1/3 N(-1,0.25^2) + 1/3 N(5,0.5^2).

For this case, a-priori we took (μ_0,τ_0,ϵ_1,ϵ_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 

    (f_1,f_2)=(2  5  g_11 + 3  5 g_12,  7 10  g_12 + 3 10 g_22),

with

    g_11   =  2/ 3 G(  2,1.1) + 1/ 3 G( 80,2)  
    
    g_12   =  8/14 G( 10,0.9) + 6/14 G(200,8.1)
    
    g_22    = 2/ 3 G(105,3)   + 1/ 3 G(500,10),

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|θ) = ℒ𝒩(x|θ) with θ=(μ,σ^2), 
is the log-normal density with mean exp(μ+σ^2/2). For this case, a-priori we set 

    (μ_0,τ_0,ϵ_1,ϵ_2)=(S̅,0.5,2,0.01),  S̅=1 n_1+n_2(∑_j=1^n_1log x_1j+∑_j=1^n_2log x_2j).

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

    _ G(p | (x_ji)) = [ 0.42 0.58; 0.64 0.36 ],  _ D(p | (x_ji)) = [ 0.42 0.58; 0.69 0.31 ],  
                     p_ true=[ 0.4 0.6; 0.7 0.3 ].





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∈{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

    (f_1^(s),f_2^(s),f_3^(s))=
    ((1-q^(s))f+q^(s)g_1,  f,  (1-q^(s))f+q^(s)g_2),

where 

    f      =3/10 N(-10,1) + 2/10 N(-6,1) 
          +2/10 N(6,1) + 3/10 N(10,1)
    
    g_1     = 1/2 N(-4,1) + 1/2 N(4,1)
    
    g_2     = 1/2 N(-12,1) + 1/2 N(12,1).


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.


  * 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.


  * 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 f̂_2^(s) and the true density f_2^(s) increases as the borrowing of strength 
weakens, it is that
H_ G(f_2^(1),f̂_2^(1))< H_ G(f_2^(2),f̂_2^(2))
< H_ G(f_2^(3),f̂_2^(3)).








4.4 Real data example.  The data set is to be found at  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.

  * The individual had a transplant.

  * 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 

    _jl=  10,     j=l=1  j=l=3
     1,




  * 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.
      

  * The estimated a-posteriori selection probabilities are given below 

    _ G(p | (x_ji)) = [ 0.61 0.23 0.16; 0.34 0.10 0.56; 0.08 0.12 0.80 ],  _ D(p | (x_ji))  = [ 0.67 0.16 0.17; 0.29 0.15 0.56; 0.10 0.12 0.78 ].


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 
S_ji={1,…,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 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
ł_ji=(1+c_ji)^-1, the correlation between the models can be controlled by
imposing further restrictions among the ł_ji parameters.

Finally, an interesting research path would be the generalization of the pairwise dependent
_j measures to include all possible interactions, in the sense that

    _j( · )=p_j _j( · )+∑_l=2^m∑_η ∈  C_j,l,mp_j,η _η_(j)( · )   with  
    p_j+∑_l=2^m∑_η ∈  C_j,l,mp_j,η=1,

where the _j and the _η_(j)'s are independent GSB processes,
C_j,l,m={(k_1,…,k_l-1):1≤ k_1<⋯<k_l-1≤ m, k_r≠ j, 1≤ r≤ m-1}
and η_(j) is the ordered vector of the elements of the vector η 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,…,f_m)
will be 2^m-1.









 Appendix A

Proof of Proposition 1.
Starting from the N_ji-augmented random densities we have

    f_j(x_ji,N_ji=r)     =    ∑_l=1^mf_j(x_ji,N_ji=r,δ_ji=l) = ∑_l=1^mp_jl f_j(x_ji,N_ji=r|δ_ji = l) 
        =    ∑_l=1^mp_jl∑_k=1^∞f_j(x_ji,N_ji=r,d_ji=k|δ_ji=l) 
        =    ∑_l=1^mp_jlf_j(N_ji=r|δ_ji=l)∑_k=1^∞f_j(d_ji=k|N_ji=r)f_j(x_ji|d_ji=k,δ_ji=l).


Because f_j(N_ji=r|_̣ji=l)=f_N(r|λ_jl) and
f_j(x_ji|d_ji=k,δ_ji=l)=K(x_ji|_jlk), the last equation gives

    f_j(x_ji,N_ji=r)  =  ∑_l=1^m p_jl f_N(r|λ_jl)∑_k=1^∞1 r I(k≤ r)K(x_ji|_jlk)
                            =  1/r∑_l=1^mp_jlf_N(r|λ_jl)∑_k=1^r K(x_ji|θ_jlk).

Augmenting further with the variables d_ji and _̣ji yields

    f_j(x_ji,N_ji=r,d_ji=k,_̣ji=l)=1 r p_jl f_N(r|λ_jl)  I(k≤ r) K(x_ji|_jlk).

Because P(_̣ji=l)=p_jl, the last equation leads to equation (<ref>) and the proposition follows. □




Proof of Proposition 2.
Marginalizing the joint of x_ji and N_ji with respect to x_ji we obtain

    f_j(N_ji=r)=∑_l=1^m p_jlf_N(r|λ_jl).

Then dividing equation (<ref>) with the probability that N_ji equals r we obtain equation (<ref>). □




Proof of Lemma 1.
Because g_(x)=ł∑_j=1^∞(1-ł)^j-1K(x|_j), we have

    {g_(x)^2}=ł^2 {(∑_j=1^∞(1-ł)^j-1K(x|θ_j))^2}
        =ł^2{∑_j=1^∞(1-ł)^2j-2 [K(x|_j)^2]+2∑_k=2^∞∑_j=1^k-1(1-ł)^j+k-2 [K(x|_j)K(x|_k)]}
        =ł^2{∑_j=1^∞(1-ł)^2j-2[K(x|)^2]+2 ∑_k=2^∞∑_j=1^k-1(1-ł)^j+k-2[K(x|)]^2}
        =ł^2{1ł(2-ł)[K(x|)^2]+2 1-łł^2(2-ł)[K(x|)]^2},

which gives the desired result.□



Proof of Proposition 3.
The random densities f_i(x)=∑_l=1^m p_il g_il(x) and f_j(x)=∑_l=1^m p_jl g_jl(x) 
depend to each other through the random measure _ji, therefore 

    [f_i(x)f_j(x)]= [ (f_i(x)f_j(x)|_ji) ]={ [f_i(x)|_ji] [f_j(x)|_ji] },

and

    [f_j(x)|_ji]=∑_l≠ ip_jl [g_jl(x)]+p_jig_ji(x)
                             =(1-p_ji) [K(x|)]+p_jig_ji(x)
       [f_i(x)|_ji]=∑_l≠ jp_il [g_il(x)]+p_ijg_ji(x)
                             =(1-p_ij) [K(x|)]+p_ijg_ji(x) .

Substituting back to equation (<ref>) one obtains

    [f_i(x)f_j(x)]=(1-p_ijp_ji) [K(x|)]^2+p_ijp_ji [g_ji(x)^2].

Using lemma 1, the last equation becomes

    [f_i(x)f_j(x)]=ł_jip_jip_ij 2-ł_ji{[K(x|)^2]-[K(x|)]^2}+[K(x|)]^2,

or that

    Cov(f_j(x),f_i(x))  = ł_jip_ji p_ij 2-ł_ji Var(K(x|)).

The desired result, comes from the fact that

    Var(∫_Θ K(x|)_ji(d)) 
               =
    {ł_ji 2-ł_ji[K(x|)^2]+2(1-ł_ji) 2-ł_ji[K(x|)]^2}-[K(x|)]^2
        = ł_ji 2-ł_ji([K(x|)^2]-[K(x|)]^2).

□


Proof of Proposition 4.

(1.) From equation (<ref>) and proposition 3, we have that

    Var(f_j^ G(x))
    =   Var(∑_l=1^m p_jlg_jl^ G(x))
    =∑_l=1^m p_ji^2λ_ji 2-λ_ji Var(K(x|)).

Normalizing the covariance in equation (<ref>) with the associated
standard deviations, yields

    Corr(f_j^ G(x),f_i^ G(x))  = λ_jip_jip_ij 2-λ_ji(∑_l=1^m∑_r=1^m  p_jl^2 p_ir^2λ_jlλ_ir (2-λ_jl)(2-λ_ir))^-1/2.

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

    Var(f_j^ D(x))
    =∑_l=1^m p_ji^2 1+c_ji Var(K(x|)),

and

    Corr(f_j^ D(x),f_i^ D(x))  = p_jip_ij 1+c_ji(∑_l=1^m∑_r=1^m  p_jl^2 p_ir^2λ_jlλ_ir (1+c_jl)(1+c_ir))^-1/2.



(2.)
When λ_ji=λ and c_ji=c for all 1≤ j≤ i≤ m, from equations (<ref>)
and (<ref>), it is clear that 

    Corr(f_j^ G(x),f_i^ G(x))= Corr(f_j^ D(x),f_i^ D(x))
     = p_jip_ij( ∑_l=1^m ∑_r=1^m  p_jl^2 p_ir^2 )^-1/2.








 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|_jl), where _jl∼ DP(c_jl,P_0). 
We randomize the concentrations by letting c_jl∼ 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 x_j and the number of unique clusters for which 
_̣ji= e_j. Letting 

    ρ_jj=#{d_jj:_̣ji= e_j,1≤ i≤ n_j},
 
we have

    ∼̱ Be(c_jj+1, n_j)
        c_jj | ,̱ρ_jj ∼ π_ G(a_jj+ρ_jj, b_jj-log)̱ +
             (1-π_)̱  G(a_jj+ρ_jj-1, b_jj-log)̱

with the weights π_$̱ satisfyingπ_/1-π_=a_jj+ρ_jj-1/n_j(b_jj-log)̱.

 B. Forj≠l, the posteriorc_jl's will be affected by the size 
of the data setsx_jandx_land the cumulative number 
of unique clustersd_jifor which_̣ji=e_land the
unique clustersd_lifor which_̣li=e_j. Lettingρ_jl=#{d_ji:_̣ji= e_l,1≤ i≤ n_j}+
          #{d_li:_̣li= e_j,1≤ i≤ n_l},it is that

    ∼̱ Be(c_jl+1, n_j+n_l)
        c_jl | ,̱ρ_jl ∼ π_𝒢(a_jl+ρ_jl, b_jl-log)̱ +
             (1-π_)̱ 𝒢(a_jl+ρ_jl-1, b_jl-log)̱,

with the weightsπ_$̱ satisfying 
π_/1-π_=a_jl+ρ_jl-1/(n_j+n_l)(b_jl-log)̱.

Bear in mind that ρ_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 λ_jl 
under beta conjugate and transformed gamma nonconjugate priors.
We let 

    S_jl =∑_i=1^n_j I(_̣ji= e_l)   and  
    S_jl'=∑_i=1^n_j I(_̣ji= e_l)(N_ji-1).

 
A. For the choice of prior λ_jl∼ Be(a_jl,b_jl), 
for l=j it is that

    f(ł_jj|⋯)= Be(ł_jl| a_jj + 2 S_jj, b_jj + S_jj'),

also, for l≠ j we have

    f(ł_jl|⋯) =  Be(ł_jl| a_jl + 2(S_jl+S_lj), b_jl + S_jl' + S_lj').


B. For the choice of prior  ł_jl∼ TG(a_jl,b_jl), for l=j it is that

    f(ł_jj|…) ∝ł_jj^2S_jj - a_jj-1(1-ł_jj)^S_jj'+ a_jj-1e^-b_jj/ł_jj  I(0<ł_jj<1).

To sample from this density, we include the positive auxiliary random variables ν_1 and ν_2 such that

    f(ł_jj,ν_1,ν_2|⋯) ∝ł_jj^2S_jj - a_jj-1 I(ν_1<(1-ł_jj)^S_jj'+ a_jj-1) 
     I(ν_2<e^-b_jj/ł_jj) I(0<ł_jj<1).

The full conditionals for ν_1,ν_2 are uniforms

    f(ν_1|⋯) = U(ν_1|0, (1-ł_jj)^S_jj'+ a_jj-1)   and  
    f(ν_2|⋯) = U(ν_2|0, e^-b_jj/ł_jj),

and the new full conditional for λ_jj becomes

    f(ł_jj|ν_1,ν_2,…) ∝ł_jj^2S_jj - a_jj-1 I(-b_jjlogν_2<ł_jj<1-ν_1^1/L_jj)        L_jj≥ 0 
     I(max{-b_jjlogν_2,1-ν_1^1/L_jj}<ł_jj<1)     L_jj<0,

where we have set L_jj=S_jj'+ a_jj-1.
We can sample from this density using the inverse cumulative distribution function technique.
Also, for l≠ j we apply the same embedded Gibbs sampling technique to the full conditional 
density

    f(ł_jl|⋯)∝ł_jl^2(S_jl+S_lj)-a_jl-1(1-ł_jl)^S_jl'+S_lj'+ 
    a_jl-1e^-b_jl/ł_jl  I(0<ł_jl<1).






 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.



