# MCMC Sampling Estimation of Poisson-Dirichlet Process Mixture Models

Authors: Xiang Qiu; Linlin Yuan; Xueqin Zhou
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this article, we aim to estimate the parameters of Poisson-Dirichlet mixture model with multigroup data structure by empirical Bayes. The number of mixture components with Bayesian nonparametric process priors is not fixed in advance and it can grow with the increase of data. Empirical Bayes is the useful method to estimate the mixture components without information on them in advance. We give the procedure to construct smooth estimates of base distributionG0 and estimates of the two parameters α,θ. The performances of estimations for parameters under multigroup data are better than those of the single-group data with the same total size of individuals in the perspectives of bias, standard deviations, and mean squared errors by numerical simulation. Also, we applied Poisson-Dirichlet mixture models to well-known real datasets.


---

## Body

## 1. Introduction

Mixture models are typically used to model data in which each observation belongs to one of a set of distributions, being finite or infinite, with certain probabilities. In many situations, some elements or all of the mixturing distributions are unknown. However, accuracy of prior information plays an important role in statistical inferences, and the deviation of prior information inevitably results in biased inference which may cause serious problems in subsequent decisions. A possible solution to the inference on parameters of the prior is the empirical Bayes (EB) methods with the price of using more data sources than pure Bayesian inference. EB is a kind of modern statistical analysis methodologies that blend frequentist inference into Bayesian ideas with the parameters in prior distributions estimated from the observed data (or Bayes decisions estimated from observed data). EB has been widely applied in many fields in its development process of more than half a century. Relevant contributions are, among others, Seidel [1] on estimating partial prior information by EB, George and Foster [2] on EB variable selection in regression analysis, Clyde and George [3] on EB estimation for wavelets, Yang and Wu [4] on EB estimates of the precision parameter and baseline distribution with monotone missing data subject to Dirichlet process priors, Bi and Davuluri [5] developing a novel nonparametric empirical Bayesian-based approach to model the RNA-seq data, and Zhou et al. [6] presenting careful comparisons among several empirical Bayes estimates to the precision parameter of Dirichlet process prior. As of EB inference for nonparametric problems when priors are modeled by mixture models, one can refer to, for example, the works of Liu [7] and McAuliffe et al. [8].A Dirichlet process (DP) [9] is a distribution of distributions and, thanks to its characteristics, can act as a prior distribution in nonparametric problems with possibly tied observations. Therefore, a DP mixture model allows for common distributions being shared by subsets of observations and thus forms generalizations of the finite mixture models (FMM). In FMM, the number of unknown distributions is supposed to be known and fixed and every individual is sampled from one of those distributions in different probability, whereas, in DP mixture models, the number of clusters is not needed to be given in advance and inferred based on the data. The number of clusters in a DP mixture model is random and grows in log(sample size). The research efforts on DP mixture model in Bayesian nonparametric are quite large, which can be found in, for example, the works of Ferguson [9], MacEachern and Müller [10], Neal [11], Ishwaran and James [12], Quintana et al. [13], and so on.Being proposed by Pitman and Yor [14], Poisson-Dirichlet process (PDP) is also known as Pitman-Yor process, which can be seen as the extension of DPs. In PDP model, a new observation is either from the previous ones in some cluster with a probability which is an increasing function of the cluster size or a new draw from the base distribution.As an effective and widely applied density estimation method, mixture modeling covers a great deal of phenomena in real-world practice. But it is usually hard to determine the number of mixture components. There are many methods treating the number of components as an unknown constant which is based on the data beforehand, for example, model selection methods. However, the assumption could be too rigid or wrong for the new individual may be from the previous cluster or new cluster as mentioned above. There are a lot of researches on the DP mixture model, such as the works of Berry and Christensen [15], McAuliffe et al. [8], Mostofi and Kharrati-Kopaei [16], and Hu et al. [17]. But there are few literatures about the PDP mixture model and its parameters’ estimation.In this article, we are interested in the estimation of base distribution under the multigroup data, where the density of base distribution is continuous. There are two parameters and base distribution in PDP mixture model, different from DP mixture model. One difficulty is the initial value ofα,θ and its impact on the estimation of base distribution, where we suppose that α∈0,1,θ>−α in the context. Another one is that the base distribution could not be directly observed, which is the parameter of infinite dimensional distribution. The individuals of every group follow the random probability distribution in view of the base. We handle this by the nonparametric empirical Bayes; and the estimations of α,θ are referred to in the works, among others, of Carlton [18] and Favaro et al. [19]. The performance of estimations for α,θ under multigroup data is better than that of the single-group data with the same total size of individuals in the perspectives of bias, standard deviations, and mean squared errors.The remainder of this paper is structured as follows. In Section2, we review PDP mixture models. In Section 3, under the framework of multigroup data, the procedure of estimation on the base density is described and the algorithm with the nonconjugate Gibbs sampler is given. In Section 4, some simulations with different base distributions and pairs of α,θ are shown and we compare the results between single-group and multigroup data. In Section 5, the PDP mixture model is applied to the well-known real-world datasets.


## 2. Poisson-Dirichlet Process Mixture Model

This section formulates the Poisson-Dirichlet processes mixture models we are going to discuss.The fundamental setting is a set of independent random couplesX1,η1,X2,η2…,Xn,ηn such that, for every i=1,2,…,n, Xi|ηi∼fx|ηi for certain known density function f and η1,η2,…,ηn are identically distributed as an unknown probability measure G, as discussed in, for example, in the works of Neal [11] and McAuliffe et al. [8]. There are quite a few real-world applications that can be set in this framework. Examples include such general insurance practice where the observations Xi are the claims of an insured at n periods distributed identically as fx|ηi with ηi reflecting the overall profile of the economic, political, and other situations at period i, such as presence of extreme weather, earthquake, war, and epidemics. In contrast to the most extensively discussed model under which X1,X2,…,Xn|η∼iidfx|η,η∼G, the parameters ηi are of evolutionary feature from one period to another [20]. In many real scenarios, the values of evolutionary ηi,i=1,2,…,n exhibit clustering feature and the number and the profiles of the clusters are unknown, so that it is appropriate to use stick-breaking distributions to model ηi,i=1,2,…,n.Pitman-Yor processes [14, 21] are a class of stick-breaking distributions covering Dirichlet processes as special elements. Let α∈0,1 (referred to as the discount parameter), let θ>−α (referred to as the concentration parameter), and let Vnn≥1 be a sequence of mutually independent random variables with Vn∼Beta1−α,θ+nα. Set V˜1=V1andV˜n=Vn∏i=1n−11−Vi,n≥2,sothat∑n=1∞V˜n=1. Let ξnn≥1 be a sequence of i.i.d. random variables following a common probability measure G0, taking values in a measurable space, X,ℬ, for example, and independent also of Vnn≥1. Let δξ,⋅ stand for the Dirac measure degenerated at point ξ. Then the stick-breaking distribution G⋅=∑n=1∞V˜nδξn,⋅ on X,ℬ is referred to as a Poisson-Dirichlet process (PDP) with parameters α,θ and base distribution G0, writing G∼PDG0;α,θ. Note that parameters α and θ determine the characteristic of power law of a PDP and the case α=0 corresponds to a Dirichlet process.Letη=η1,η2,…,ηn be a sequence of i.i.d. observations drawn from G so as to reflect the clustering feature of distributions. Write K for the number of clusters of η, η1∗,η2∗,…,ηK∗ for the distinct observations among η1,η2,…,ηn and mk is the number of elements in cluster k. Then, by Pitman [21], given η1,η2,…,ηn, the probability distribution of a new future ηn+1 is(1)Prηn+1∈A|η1,η2,…,ηn=Kα+θn+θG0A+∑k=1Kmk−αn+θδηk∗,A,that is, taking values ηk∗ with probability mk−α/n+θ and a new one from G0 with probability Kα+θ/n+θ. This is the generalized Pólya urn scheme to produce ηi,i=1,2,….A PDP mixture model is defined by(2)G∼PDG0;α,θ,η1,…,ηn|G∼i.i.d.G,i=1,2,…,n,Xi|ηi∼f⋅|ηi,i=1,2,…,n,where f⋅|η is a known density function characterized by parameter η, η=η1,η2,…,ηn consists of the latent parameters, and G0 is absolutely continuous with respect to certain σ-finite measure with an entirely unknown density function g0, as depicted in Figure 1. Thanks to the clustering feature of elements of η from PD processes, under this model, the observations X1,X2,…,Xn are randomly divided into a random number of clusters so that all the elements of every cluster share a common parameter picked from η1∗,η2∗,…,ηK∗ and a new observation Xn+1 either is generated from a distribution with a new parameter, meaning that Xn+1 comes from a new cluster, or falls into one of the previously existing clusters with parameter from η1∗,η2∗,…,ηK∗. Note that α=0 corresponds to the special case of DP mixture models, on which some recent efforts can be found in, for example, the work of Neal [11] who discussed the methods of sampling from the posterior distribution of a DP mixture model and the work of McAuliffe et al. [8] who studied the estimation of the base distribution G0 in DP mixture model with a single group of data X1,X2,…,Xn.Figure 1
PDP mixture models with a single-group data.In this paper, however, we consider the estimation of the parameters in the distributionG when data is organized in independent and identically distributed groups. To be specific, we consider the following PDP mixture model: the data are divided into I groups and the parameters ηij,j=1,2,…,ni, of each group form an i.i.d. sample from G. A graphical representation of a PDP mixture model is shown in Figure 2, and its hierarchical structure in terms of conditional distributions is clearly expressed by(3)G1,…,GI∼i.i.d.PDG0;α,θ,ηi1,…,ηini|Gi∼i.i.d.Gi,i=1,2,…,I,Xij|ηij∼f⋅|ηij,j=1,2,…,ni,i=1,2,…,I.Figure 2
PDP mixture models with multigroup data.Under this setting, the individuals from different groups are mutually independent and the random probabilitiesGi,i=1,2,…,I, are independent and identically distributed as a PDG0;α,θ with unknown α,θ, and G0 to be estimated.On the one hand, the structure defined above arises quite frequently in real-world applications. In the context of general insurance, it is a typical structure of data extensively investigated by the so-called credibility theory [20]. On the other hand, it has been commonly recognized that, for the majority of EB problems, it is impossible to make efficient estimate of the unknown elements in the prior distributions if one has only a single group of conditionally independent and identically distributed observations given the random parameters following the priors. The only exception arises when the prior can be characterized by a so-called stick-breaking distribution such as Dirichlet processes (see [4, 22]). Even for stick-breaking priors, existing research efforts have exhibited the possibly low efficiency in estimating the parameters in prior distributions with a single group of data; see, for example, the works of Zehnwirth [22], Yang and Wu [4], and Zhang et al. [23].Fori=1,2,…,I, denote by Ki the number of distinct values among ηi=ηi1,ηi2,…,ηini, ηi1∗,ηi2∗,…,ηiKi∗ those distinct values, and mij, j=1,2…,Ki, the number of elements in cluster j of group i. Then, similar to Theorem 2.5 of Korwar and Hollander [24], given Ki, the distinct values ηi1∗,ηi2∗,…,ηiKi∗∼i.i.d.G0 and are also independent of mi1,mi2,…,miKi, so that we can retrieve information as of the base distribution G0 from the information delivered by the i.i.d. ηi1∗,ηi2∗,…,ηiKi∗,i=1,2,…,I and of the parameters θ and α from Ki,mi1…,miKi,i=1,2,…,I.Generally, denotemi=mi1,…,miKi,i=1,2,…,I and rearrange Xi1,…,Xini into Ki groups X11i,…,Xmi11i,…,X1Kii,…,XmiKiKii so that, for every k=1,2,…,Ki, the observations in Xjki,j=1,2,…,mki share the common random parameter ηik∗ of cluster k in group i. Then, by means of marginalization, the likelihood under framework (3) with data X=Xij,j=1,2,…,ni;i=1,2,…,I can be expressed as(4)Lg0,α,θ=∏i=1IEηi∏j=1nifXij|ηij=∏i=1IEηi∏j=1mi1fXj1i|ηi1∗…∏j=1miKifXjKii|ηiKi∗=∏i=1IEKi,mi∫∏j=1mi1fXj1i|ug0udu⋯∫∏j=1miKifXjKii|ug0udu,where the expectation Eηi is taken with respect to marginal distribution of ηi with Xij taken fixed, EKi,mi is understood in the same way, and g0 is the unknown density of G0. Imaginably, the likelihood function (4) of g0,α,θ is generally expressed as a product of I summations, of which every summand is, in turn, a product of the form expressed inside the parentheses in (4). So, it is too intricate to be maximized; even numerically and certain indirect alternatives to the maximum likelihood estimates of g0,α,θ need to be developed, of which the following is constructed with the help of the ideas from empirical Bayes.


## 3. Estimation of Poisson-Dirichlet Mixture Models

This section concentrates on developing reasonable estimates of the parametersg0,α, and θ and elucidates in detail algorithms deriving those estimates. We proceed in two subsections with one for the idea and the other for algorithms to realize the ideas.


### 3.1. Constructing Estimates

Should one have observedη=ηij,j=1,2,…,ni,i=1,2,…,I or, equivalently, K1,K2,…,KI together with ηij∗,mij,j=1,2,…,Ki,i=1,2,…,I. By the fact that ηij∗∼i.i.d.g0 for all j=1,2,…,Ki,i=1,2,…,I, a natural choice is to estimate g0 by certain maturely developed density estimators. We use kernel method to result in(5)g˜0⋅;η=1∑Ki∑i=1I∑j=1Kiωh⋅,ηij∗,where ωhx,ηij∗=1/hωx−ηij∗/h is defined by a certain kernel function ω⋅ together with a bandwidth h>0. For α and θ, similar to Carlton [18] who studied the estimation of PDP priors and Favaro et al. [19] who established the PDP model so as to predict species for the problem of species searching, we can use the maximum likelihood estimate derived in terms of the distribution of Ki,mi1,…,miKi,i=1,2,…,I; that is,(6)α˜η,θ˜η=Argmaxα,θ∏i=1I1θ+1ni−1⋅∏l=1Ki−1θ+lα⋅∏j=2ni1−αj−1aji,where 1+xk=1+x2+x…k−1+x and aji is the cardinality of mik:k=1,2,…,Kisuch thatmik=j, that is, each of the clusters of which has j items.Now that the latent variablesηij,j=1,2,…,ni,i=1,2,…,I are unobservable, we construct unbiased pseudoestimates:(7)g˜˜0⋅;g0,α,θ=Eg˜0⋅;η|X=E1∑i=1IKi∑i=1I∑j=1Kiωh⋅,ηij∗|X,α˜˜,θ˜˜g0,α,θ=Eα˜η,θ˜η|X,where the arguments g0,α,θ indicate that the “estimates” nonetheless depend on the unknown parameters g0, α, and θ. Note that these pseudoestimates can be computed by means of the posterior distribution for η given the observations X.By these pseudoestimates, an algorithm can be motivated, which, starting from certain appropriately selected initial values ofg0,α, and θ, iteratively computes for r=0,1,….(8)g^0r+1⋅=Eg^0r,α^r,θ^rg˜0⋅;η|X,α^r+1,θ^r+1=Eg^0r,α^r,θ^rα˜η,θ˜η|X,with the symbol Eg^0r,α^r,θ^r indicating that the posterior expectations of g˜0⋅;η,α˜η,θ˜η given X are computed with the unknown parameters g0,α,θ replaced by g^0r,α^r,θ^r. Consequently, estimate g0,α, and θ by the limit(9)g^0,α^,θ^=limr⟶∞g^0r,α^r,θ^r.In every iteration to computeg^0r+1⋅,α^r+1,θ^r+1 from g^0r,α^r,θ^r, numerical solutions are necessary because it is not an easy task to work out closed forms of the posterior expectations in the equations in (8). A natural way is to employ an MCMC technique that will be described in detail later on.To sum up, the estimates are approximated with the following workable procedure.


### 3.2. Algorithms

This subsection details Algorithm1 in the selection of initial values, the steps to perform the MCMC sampling, and the computation of the estimation proposed.Algorithm 1: Estimation procedure.
(1)
The procedure begins with certain appropriately set initial valuesg^00⋅,α^0,θ^0at round 0.(2)
Having computedg^0r⋅,α^r,θ^rat roundr, ther+1-th round proceeds in the following steps:(a) Start from certain initial valuesη0 of the latent variable η and produce a sequence of ηb,b=1,2,…,B by certain MCMC method for a sufficiently large positive integer B.(b) Computeg˜0⋅;ηb by (5) and α˜ηb,θ˜ηb by (6) for every b=1,2,…,B.(c) End by assigning a new set of valuesg^0r+1⋅,α^r+1,θ^r+1≈1/B∑b=1Bg˜0⋅;ηb,α˜ηb,θ˜ηbto the parametersg0⋅,α,θ.(3)
Repeat the algorithm until the sequenceg^0r⋅,α^r,θ^ris stable.


#### 3.2.1. Initial Values

To fix an initial estimateg^00, compute first the maximum likelihood estimate η^ij by fXij|η^ij=maxufXij|u and then compute a kernel estimate of g0 by(10)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,η^ij.Similar to the formula in (5), in many important cases, f is a location translation of some standard density with mode zero (e.g., a standard normal distribution) such that η^ij=Xij. In these cases, the initial values of g0 are given by(11)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,Xij.For parametersα,θ, the initial value of α could be 0.1 or 0.5 and θ could be 1 or 5 in case there is no further information on them. The simulation, which will be reported later on, showed that the initial values of α,θ had only quite minor effects on the convergence of the algorithm.


#### 3.2.2. The Gibbs Sampling

We next state the full conditional distribution used in the MCMC procedure, that is, the conditional distribution of a single latent variableηij given the data Xi=Xi1,Xi2,…,Xini,i=1,2,…,I, and the latent variables ηi∖ηij∪k≠iηk. Denote by ηi−j=ηi∖ηij the subvector of ηi by taking away ηij and Ki−j correspondingly the numbers of the distinct values in ηi−j, j=1,2,…,ni,i=1,2,…,I. With an abuse of notation, let ηi1∗,ηi2∗,…,ηiKi−j∗ be the Ki−j distinct values in ηi−j and mil=#r:r∈1,2,…,ni∖jandηir=ηil∗ the number of such elements in ηi−j which take the value ηil∗.Thanks to the mutual independence among the random vectorsGi,Xi1,ηi1…,Xini,ηini,i=1,2,…,I, it suffices to discuss the single-group situation, as stated in Theorem 3.1.Theorem 1.
For eachj=1,2,…,ni, givenηi−jandXi, the conditional distribution ofηijis(12)ηij|ηi−j,Xi;g0,α,θ∼q0Gη⋅|Xij+∑l=1Ki−jqlδηil∗,⋅,where Gη⋅|Xij is the posterior distribution of ηij derived from a single observation Xij with density gηu|Xij=fXij|ug0u/∫fXij|ug0du, qll=0Ki−j are such that q0=bθ+αKi−j∫fXij|ug0udu, ql=bmil−αfXij|ηil∗,l=1,…,Ki−j, and the normalization quantity b is determined by the equality ∑l=0Ki−jql=1.Proof.
By the generalized Pólya urn scheme (1), the posterior distribution under PDg0;α,θ is(13)Prηij∈⋅|ηi−j=θ+αKi−jθ+ni−1g0⋅+1θ+ni−1∑l=1Ki−jmil−αδηil∗,⋅.
Therefore, the conditional distribution ofηij given ηi−j,Xi can be computed by(14)dGηij|ηi−j,Xi==fXij|ηijdFηij|ηi−j∫fXij|ηijdFηij|ηi−j=θ+αKi−jfXij|ηijg0ηijdηij+∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηijθ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗.
Set(15)ΛXij=1θ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗,and rewrite the equation above as(16)dGηij|ηi−j,Xi=ΛXijθ+αKi−jfXij|ηijg0ηijdηij⋅fXij|ηijg0ηij∫fXij|ηijg0dηij+ΛXij⋅∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηij.
That is, the conditional distribution ofηij is(17)ηij|ηi−j,Xi;g0,α,θ∼q0Gηηij|Xij+∑l=1Ki−jqlδηil∗,ηij,
whereq0=ΛXijθ+αKi−j∫fXij|ηijg0ηijdηij and ql=ΛXijmil−αfXij|ηil∗. Obviously, q0 and qls satisfy ∑l=0Ki−jql=1. The proof is completed.


#### 3.2.3. The Iteration Procedure

With the theorem above and the initial valuesg^00,α^0,θ^0, we now describe the iteration to compute g^0r+1,α^r+1,θ^r+1 from g^0r,α^r,θ^r for all r=0,1,2…, which consists of two alternatively running phases: one for sampling from the posterior distribution of the latent variables η with the unknown parameters replaced by the estimates g^0r,α^r,θ^r and the other for estimation of the parameters using the sampled latent variables η, as depicted in what follows:(1) Sampling. The sampling phase is indeed a two-step procedure applied to every group i=1,2,…,I separately so as to produce the posterior samples ηi, starting with an arbitrarily selected ηi0=ηi10,ηi20,…,ηini0. The first step is the standard Gibbs sampling derived from the full conditionals presented in Theorem 1 and the second performs an acceleration.As frequently done in the literature (see, e.g., [11]), introduce an auxiliary vector cηi=Δc=c1,c2,…,cni to stand for the cluster indices of ηijs in ηi1∗,ηi2∗,…,ηiKi∗; that is, cj=k⇔ηij=ηik∗, so that ηi=ηic1∗,ηic2∗,…,ηicni∗. In this notation, ηi is equivalent to the combination of c and ηi1∗,ηi2∗,…,ηiKi∗. With the cluster indicators in c, in this two-step procedure, the first step updates the values of cj and in passing the values of ηij respecting the full conditionals in Theorem 1. For the case ∑l=1Ki−jql is much larger than q0, the Markovian chain directly derived from Theorem 1 may converge slowly because the chance to introduce new values for ηij is quite small; see equation (12). Hence, an acceleration method (e.g., [11]) may be helpful, which samples another round of ηi for acceleration. This is realized by sampling ηi conditional on Xi and the cluster indicators c just obtained in Step 1 and then sampling c given Xi,ηi. In other words, Step 2 is driven by the grouped conditional distribution of ηi given Xi,c and c given Xi,ηi. Because of the particular feature that c is a deterministic function of ηi, Step 2 only needs to draw values from the conditional distribution of ηi given Xi,c.Moreover, for everyj=1,…,ni, the integration ∫fXij|ηijg^0rηijdηij involved in computing Gη⋅|Xij and ql,l=0,1,…,Ki−j is approximated by an average(18)hrj=1d∑l=1dfXij|ul,where d is a positive integer and u1,…,ud∼iidg^0r. Because u1,…,ud∼iidg^0r depend only on the iteration indicator r, this approximation can be prepared at the beginning of each iteration.As a final point, we would like to note that, due to the presence of the acceleration Step 2, the update of the values ofηij in Step 1 essentially takes effect only on the computation of the probabilities ql in equation (12) for every j=1,2,…,ni.The details of this two-step procedure of generatingηib are presented in Algorithm 2. For each group i=1,2,…,I, this procedure is repeated to generate ηib=ηi1b,…,ηinib,b=1,2,…,B, for a sufficiently large B.Algorithm 2: Gibbs sampling under single group.
(1)
Drawdiid valuesu1,u2,…,udfromg^0rand computehrj,j=1,…,niby (13).(2)
Initiate the algorithm with an arbitrarily selectedηi0. Correspondingly, setc=cηi0. Letη∗=ηi1∗,ηi2∗,…,ηiKi∗,0,…,0be the vector of dimensionniin which the firstKielements store the distinct value ofηi0andm=mi1,mi2,…,miKi,0,…,0be anni-vector to represent the numbers ofηij∗occurring inηi0.(3)
Now assume that we have gotηib−1and the corresponding auxiliary vectorsc,mandη∗. The sampled valuesηibare generated by the following two-step procedure:Step 1. Forj=1,2,…,ni,(1) setmicj=micj−1, Ki−j=Ki−1ifmicj=0Kiotherwise;(2) draw a new value forcj from 1,2,…,ni according to the distribution lawqc=Prcj=c|Xi,c∖cj,ηi1∗,…,ηiKi−j∗=ℬmic−α^rfXij|ηic∗for1≤c≤ni,mic≠0b/ni−Ki−jθ^r+α^rKi−jhrjfor1≤c≤ni,mic=0and then, in casemicj=0, draw a value for ηicj from Gη⋅|Xij with g0 replaced by g^0r, where ℬ is the normalization constant subject to ∑c=0niqc=1; and then(3) setmicj=micj+1.Step 2. For allc∈c1,c2,…,cni which are obtained by above Step 1, draw ηb from the conditional densityfηic|c,Xi∝g^0rηic∏j:cj=cfXij|ηic.(2) Estimation. In phase estimation, with ηb,b=1,2,…,B, in hand, the estimates of α,θ, and g0 are obtained according to (10).By these newly obtained estimates, repeat the above two phases for the next iteration until the sequence of estimates is stable. The above computation process on the estimation ofα,θ,g0 under the framework in (3) can be briefly arranged as in Algorithm 3.Algorithm 3: Gibbs sampling under multigroup.
Step 1. Initializeα^0,θ^0andg^00;Step 2. At roundr+1,r=1,2,…, for alli=1,2,…,I,(i) arbitrarily selection initial values forηi0 and(ii) updateηi with Algorithm 2 to produce ηb,b=1,2,…,B;Step 3. Obtainα^r+1,θ^r+1andg^0r+1according to (10).We now conclude this section by the following remark on one of Neal’s [11] algorithms.Remark 1.
Algorithm2 is similar to, but with two minor differences, to Algorithm 8 of Neal [11] at Step 1. First, we only draw a single set of d i.i.d. values from g^0r in each iteration, while Neal [11] suggested drawing d i.i.d. values for each j=1,…,ni, hence needing ni−1 times more than the sampling computation here. The algorithm here hence may cause a significant reduction in the total computation load. Second, we draw a value for ηicj from Gη⋅|Xij if a new cluster for cj is produced, while Neal [11] simply drew a value from g0r. According to Theorem 1, we think Algorithm 2 draws a correct value.


## 3.1. Constructing Estimates

Should one have observedη=ηij,j=1,2,…,ni,i=1,2,…,I or, equivalently, K1,K2,…,KI together with ηij∗,mij,j=1,2,…,Ki,i=1,2,…,I. By the fact that ηij∗∼i.i.d.g0 for all j=1,2,…,Ki,i=1,2,…,I, a natural choice is to estimate g0 by certain maturely developed density estimators. We use kernel method to result in(5)g˜0⋅;η=1∑Ki∑i=1I∑j=1Kiωh⋅,ηij∗,where ωhx,ηij∗=1/hωx−ηij∗/h is defined by a certain kernel function ω⋅ together with a bandwidth h>0. For α and θ, similar to Carlton [18] who studied the estimation of PDP priors and Favaro et al. [19] who established the PDP model so as to predict species for the problem of species searching, we can use the maximum likelihood estimate derived in terms of the distribution of Ki,mi1,…,miKi,i=1,2,…,I; that is,(6)α˜η,θ˜η=Argmaxα,θ∏i=1I1θ+1ni−1⋅∏l=1Ki−1θ+lα⋅∏j=2ni1−αj−1aji,where 1+xk=1+x2+x…k−1+x and aji is the cardinality of mik:k=1,2,…,Kisuch thatmik=j, that is, each of the clusters of which has j items.Now that the latent variablesηij,j=1,2,…,ni,i=1,2,…,I are unobservable, we construct unbiased pseudoestimates:(7)g˜˜0⋅;g0,α,θ=Eg˜0⋅;η|X=E1∑i=1IKi∑i=1I∑j=1Kiωh⋅,ηij∗|X,α˜˜,θ˜˜g0,α,θ=Eα˜η,θ˜η|X,where the arguments g0,α,θ indicate that the “estimates” nonetheless depend on the unknown parameters g0, α, and θ. Note that these pseudoestimates can be computed by means of the posterior distribution for η given the observations X.By these pseudoestimates, an algorithm can be motivated, which, starting from certain appropriately selected initial values ofg0,α, and θ, iteratively computes for r=0,1,….(8)g^0r+1⋅=Eg^0r,α^r,θ^rg˜0⋅;η|X,α^r+1,θ^r+1=Eg^0r,α^r,θ^rα˜η,θ˜η|X,with the symbol Eg^0r,α^r,θ^r indicating that the posterior expectations of g˜0⋅;η,α˜η,θ˜η given X are computed with the unknown parameters g0,α,θ replaced by g^0r,α^r,θ^r. Consequently, estimate g0,α, and θ by the limit(9)g^0,α^,θ^=limr⟶∞g^0r,α^r,θ^r.In every iteration to computeg^0r+1⋅,α^r+1,θ^r+1 from g^0r,α^r,θ^r, numerical solutions are necessary because it is not an easy task to work out closed forms of the posterior expectations in the equations in (8). A natural way is to employ an MCMC technique that will be described in detail later on.To sum up, the estimates are approximated with the following workable procedure.


## 3.2. Algorithms

This subsection details Algorithm1 in the selection of initial values, the steps to perform the MCMC sampling, and the computation of the estimation proposed.Algorithm 1: Estimation procedure.
(1)
The procedure begins with certain appropriately set initial valuesg^00⋅,α^0,θ^0at round 0.(2)
Having computedg^0r⋅,α^r,θ^rat roundr, ther+1-th round proceeds in the following steps:(a) Start from certain initial valuesη0 of the latent variable η and produce a sequence of ηb,b=1,2,…,B by certain MCMC method for a sufficiently large positive integer B.(b) Computeg˜0⋅;ηb by (5) and α˜ηb,θ˜ηb by (6) for every b=1,2,…,B.(c) End by assigning a new set of valuesg^0r+1⋅,α^r+1,θ^r+1≈1/B∑b=1Bg˜0⋅;ηb,α˜ηb,θ˜ηbto the parametersg0⋅,α,θ.(3)
Repeat the algorithm until the sequenceg^0r⋅,α^r,θ^ris stable.


### 3.2.1. Initial Values

To fix an initial estimateg^00, compute first the maximum likelihood estimate η^ij by fXij|η^ij=maxufXij|u and then compute a kernel estimate of g0 by(10)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,η^ij.Similar to the formula in (5), in many important cases, f is a location translation of some standard density with mode zero (e.g., a standard normal distribution) such that η^ij=Xij. In these cases, the initial values of g0 are given by(11)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,Xij.For parametersα,θ, the initial value of α could be 0.1 or 0.5 and θ could be 1 or 5 in case there is no further information on them. The simulation, which will be reported later on, showed that the initial values of α,θ had only quite minor effects on the convergence of the algorithm.


### 3.2.2. The Gibbs Sampling

We next state the full conditional distribution used in the MCMC procedure, that is, the conditional distribution of a single latent variableηij given the data Xi=Xi1,Xi2,…,Xini,i=1,2,…,I, and the latent variables ηi∖ηij∪k≠iηk. Denote by ηi−j=ηi∖ηij the subvector of ηi by taking away ηij and Ki−j correspondingly the numbers of the distinct values in ηi−j, j=1,2,…,ni,i=1,2,…,I. With an abuse of notation, let ηi1∗,ηi2∗,…,ηiKi−j∗ be the Ki−j distinct values in ηi−j and mil=#r:r∈1,2,…,ni∖jandηir=ηil∗ the number of such elements in ηi−j which take the value ηil∗.Thanks to the mutual independence among the random vectorsGi,Xi1,ηi1…,Xini,ηini,i=1,2,…,I, it suffices to discuss the single-group situation, as stated in Theorem 3.1.Theorem 1.
For eachj=1,2,…,ni, givenηi−jandXi, the conditional distribution ofηijis(12)ηij|ηi−j,Xi;g0,α,θ∼q0Gη⋅|Xij+∑l=1Ki−jqlδηil∗,⋅,where Gη⋅|Xij is the posterior distribution of ηij derived from a single observation Xij with density gηu|Xij=fXij|ug0u/∫fXij|ug0du, qll=0Ki−j are such that q0=bθ+αKi−j∫fXij|ug0udu, ql=bmil−αfXij|ηil∗,l=1,…,Ki−j, and the normalization quantity b is determined by the equality ∑l=0Ki−jql=1.Proof.
By the generalized Pólya urn scheme (1), the posterior distribution under PDg0;α,θ is(13)Prηij∈⋅|ηi−j=θ+αKi−jθ+ni−1g0⋅+1θ+ni−1∑l=1Ki−jmil−αδηil∗,⋅.
Therefore, the conditional distribution ofηij given ηi−j,Xi can be computed by(14)dGηij|ηi−j,Xi==fXij|ηijdFηij|ηi−j∫fXij|ηijdFηij|ηi−j=θ+αKi−jfXij|ηijg0ηijdηij+∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηijθ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗.
Set(15)ΛXij=1θ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗,and rewrite the equation above as(16)dGηij|ηi−j,Xi=ΛXijθ+αKi−jfXij|ηijg0ηijdηij⋅fXij|ηijg0ηij∫fXij|ηijg0dηij+ΛXij⋅∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηij.
That is, the conditional distribution ofηij is(17)ηij|ηi−j,Xi;g0,α,θ∼q0Gηηij|Xij+∑l=1Ki−jqlδηil∗,ηij,
whereq0=ΛXijθ+αKi−j∫fXij|ηijg0ηijdηij and ql=ΛXijmil−αfXij|ηil∗. Obviously, q0 and qls satisfy ∑l=0Ki−jql=1. The proof is completed.


### 3.2.3. The Iteration Procedure

With the theorem above and the initial valuesg^00,α^0,θ^0, we now describe the iteration to compute g^0r+1,α^r+1,θ^r+1 from g^0r,α^r,θ^r for all r=0,1,2…, which consists of two alternatively running phases: one for sampling from the posterior distribution of the latent variables η with the unknown parameters replaced by the estimates g^0r,α^r,θ^r and the other for estimation of the parameters using the sampled latent variables η, as depicted in what follows:(1) Sampling. The sampling phase is indeed a two-step procedure applied to every group i=1,2,…,I separately so as to produce the posterior samples ηi, starting with an arbitrarily selected ηi0=ηi10,ηi20,…,ηini0. The first step is the standard Gibbs sampling derived from the full conditionals presented in Theorem 1 and the second performs an acceleration.As frequently done in the literature (see, e.g., [11]), introduce an auxiliary vector cηi=Δc=c1,c2,…,cni to stand for the cluster indices of ηijs in ηi1∗,ηi2∗,…,ηiKi∗; that is, cj=k⇔ηij=ηik∗, so that ηi=ηic1∗,ηic2∗,…,ηicni∗. In this notation, ηi is equivalent to the combination of c and ηi1∗,ηi2∗,…,ηiKi∗. With the cluster indicators in c, in this two-step procedure, the first step updates the values of cj and in passing the values of ηij respecting the full conditionals in Theorem 1. For the case ∑l=1Ki−jql is much larger than q0, the Markovian chain directly derived from Theorem 1 may converge slowly because the chance to introduce new values for ηij is quite small; see equation (12). Hence, an acceleration method (e.g., [11]) may be helpful, which samples another round of ηi for acceleration. This is realized by sampling ηi conditional on Xi and the cluster indicators c just obtained in Step 1 and then sampling c given Xi,ηi. In other words, Step 2 is driven by the grouped conditional distribution of ηi given Xi,c and c given Xi,ηi. Because of the particular feature that c is a deterministic function of ηi, Step 2 only needs to draw values from the conditional distribution of ηi given Xi,c.Moreover, for everyj=1,…,ni, the integration ∫fXij|ηijg^0rηijdηij involved in computing Gη⋅|Xij and ql,l=0,1,…,Ki−j is approximated by an average(18)hrj=1d∑l=1dfXij|ul,where d is a positive integer and u1,…,ud∼iidg^0r. Because u1,…,ud∼iidg^0r depend only on the iteration indicator r, this approximation can be prepared at the beginning of each iteration.As a final point, we would like to note that, due to the presence of the acceleration Step 2, the update of the values ofηij in Step 1 essentially takes effect only on the computation of the probabilities ql in equation (12) for every j=1,2,…,ni.The details of this two-step procedure of generatingηib are presented in Algorithm 2. For each group i=1,2,…,I, this procedure is repeated to generate ηib=ηi1b,…,ηinib,b=1,2,…,B, for a sufficiently large B.Algorithm 2: Gibbs sampling under single group.
(1)
Drawdiid valuesu1,u2,…,udfromg^0rand computehrj,j=1,…,niby (13).(2)
Initiate the algorithm with an arbitrarily selectedηi0. Correspondingly, setc=cηi0. Letη∗=ηi1∗,ηi2∗,…,ηiKi∗,0,…,0be the vector of dimensionniin which the firstKielements store the distinct value ofηi0andm=mi1,mi2,…,miKi,0,…,0be anni-vector to represent the numbers ofηij∗occurring inηi0.(3)
Now assume that we have gotηib−1and the corresponding auxiliary vectorsc,mandη∗. The sampled valuesηibare generated by the following two-step procedure:Step 1. Forj=1,2,…,ni,(1) setmicj=micj−1, Ki−j=Ki−1ifmicj=0Kiotherwise;(2) draw a new value forcj from 1,2,…,ni according to the distribution lawqc=Prcj=c|Xi,c∖cj,ηi1∗,…,ηiKi−j∗=ℬmic−α^rfXij|ηic∗for1≤c≤ni,mic≠0b/ni−Ki−jθ^r+α^rKi−jhrjfor1≤c≤ni,mic=0and then, in casemicj=0, draw a value for ηicj from Gη⋅|Xij with g0 replaced by g^0r, where ℬ is the normalization constant subject to ∑c=0niqc=1; and then(3) setmicj=micj+1.Step 2. For allc∈c1,c2,…,cni which are obtained by above Step 1, draw ηb from the conditional densityfηic|c,Xi∝g^0rηic∏j:cj=cfXij|ηic.(2) Estimation. In phase estimation, with ηb,b=1,2,…,B, in hand, the estimates of α,θ, and g0 are obtained according to (10).By these newly obtained estimates, repeat the above two phases for the next iteration until the sequence of estimates is stable. The above computation process on the estimation ofα,θ,g0 under the framework in (3) can be briefly arranged as in Algorithm 3.Algorithm 3: Gibbs sampling under multigroup.
Step 1. Initializeα^0,θ^0andg^00;Step 2. At roundr+1,r=1,2,…, for alli=1,2,…,I,(i) arbitrarily selection initial values forηi0 and(ii) updateηi with Algorithm 2 to produce ηb,b=1,2,…,B;Step 3. Obtainα^r+1,θ^r+1andg^0r+1according to (10).We now conclude this section by the following remark on one of Neal’s [11] algorithms.Remark 1.
Algorithm2 is similar to, but with two minor differences, to Algorithm 8 of Neal [11] at Step 1. First, we only draw a single set of d i.i.d. values from g^0r in each iteration, while Neal [11] suggested drawing d i.i.d. values for each j=1,…,ni, hence needing ni−1 times more than the sampling computation here. The algorithm here hence may cause a significant reduction in the total computation load. Second, we draw a value for ηicj from Gη⋅|Xij if a new cluster for cj is produced, while Neal [11] simply drew a value from g0r. According to Theorem 1, we think Algorithm 2 draws a correct value.


## 3.2.1. Initial Values

To fix an initial estimateg^00, compute first the maximum likelihood estimate η^ij by fXij|η^ij=maxufXij|u and then compute a kernel estimate of g0 by(10)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,η^ij.Similar to the formula in (5), in many important cases, f is a location translation of some standard density with mode zero (e.g., a standard normal distribution) such that η^ij=Xij. In these cases, the initial values of g0 are given by(11)g^00⋅=1∑i=1Ini∑i=1I∑j=1niωh⋅,Xij.For parametersα,θ, the initial value of α could be 0.1 or 0.5 and θ could be 1 or 5 in case there is no further information on them. The simulation, which will be reported later on, showed that the initial values of α,θ had only quite minor effects on the convergence of the algorithm.


## 3.2.2. The Gibbs Sampling

We next state the full conditional distribution used in the MCMC procedure, that is, the conditional distribution of a single latent variableηij given the data Xi=Xi1,Xi2,…,Xini,i=1,2,…,I, and the latent variables ηi∖ηij∪k≠iηk. Denote by ηi−j=ηi∖ηij the subvector of ηi by taking away ηij and Ki−j correspondingly the numbers of the distinct values in ηi−j, j=1,2,…,ni,i=1,2,…,I. With an abuse of notation, let ηi1∗,ηi2∗,…,ηiKi−j∗ be the Ki−j distinct values in ηi−j and mil=#r:r∈1,2,…,ni∖jandηir=ηil∗ the number of such elements in ηi−j which take the value ηil∗.Thanks to the mutual independence among the random vectorsGi,Xi1,ηi1…,Xini,ηini,i=1,2,…,I, it suffices to discuss the single-group situation, as stated in Theorem 3.1.Theorem 1.
For eachj=1,2,…,ni, givenηi−jandXi, the conditional distribution ofηijis(12)ηij|ηi−j,Xi;g0,α,θ∼q0Gη⋅|Xij+∑l=1Ki−jqlδηil∗,⋅,where Gη⋅|Xij is the posterior distribution of ηij derived from a single observation Xij with density gηu|Xij=fXij|ug0u/∫fXij|ug0du, qll=0Ki−j are such that q0=bθ+αKi−j∫fXij|ug0udu, ql=bmil−αfXij|ηil∗,l=1,…,Ki−j, and the normalization quantity b is determined by the equality ∑l=0Ki−jql=1.Proof.
By the generalized Pólya urn scheme (1), the posterior distribution under PDg0;α,θ is(13)Prηij∈⋅|ηi−j=θ+αKi−jθ+ni−1g0⋅+1θ+ni−1∑l=1Ki−jmil−αδηil∗,⋅.
Therefore, the conditional distribution ofηij given ηi−j,Xi can be computed by(14)dGηij|ηi−j,Xi==fXij|ηijdFηij|ηi−j∫fXij|ηijdFηij|ηi−j=θ+αKi−jfXij|ηijg0ηijdηij+∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηijθ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗.
Set(15)ΛXij=1θ+αKi−j∫fXij|ηijg0dηij+∑l=1Ki−jmil−αfXij|ηil∗,and rewrite the equation above as(16)dGηij|ηi−j,Xi=ΛXijθ+αKi−jfXij|ηijg0ηijdηij⋅fXij|ηijg0ηij∫fXij|ηijg0dηij+ΛXij⋅∑l=1Ki−jmil−αfXij|ηil∗δηil∗,ηij.
That is, the conditional distribution ofηij is(17)ηij|ηi−j,Xi;g0,α,θ∼q0Gηηij|Xij+∑l=1Ki−jqlδηil∗,ηij,
whereq0=ΛXijθ+αKi−j∫fXij|ηijg0ηijdηij and ql=ΛXijmil−αfXij|ηil∗. Obviously, q0 and qls satisfy ∑l=0Ki−jql=1. The proof is completed.


## 3.2.3. The Iteration Procedure

With the theorem above and the initial valuesg^00,α^0,θ^0, we now describe the iteration to compute g^0r+1,α^r+1,θ^r+1 from g^0r,α^r,θ^r for all r=0,1,2…, which consists of two alternatively running phases: one for sampling from the posterior distribution of the latent variables η with the unknown parameters replaced by the estimates g^0r,α^r,θ^r and the other for estimation of the parameters using the sampled latent variables η, as depicted in what follows:(1) Sampling. The sampling phase is indeed a two-step procedure applied to every group i=1,2,…,I separately so as to produce the posterior samples ηi, starting with an arbitrarily selected ηi0=ηi10,ηi20,…,ηini0. The first step is the standard Gibbs sampling derived from the full conditionals presented in Theorem 1 and the second performs an acceleration.As frequently done in the literature (see, e.g., [11]), introduce an auxiliary vector cηi=Δc=c1,c2,…,cni to stand for the cluster indices of ηijs in ηi1∗,ηi2∗,…,ηiKi∗; that is, cj=k⇔ηij=ηik∗, so that ηi=ηic1∗,ηic2∗,…,ηicni∗. In this notation, ηi is equivalent to the combination of c and ηi1∗,ηi2∗,…,ηiKi∗. With the cluster indicators in c, in this two-step procedure, the first step updates the values of cj and in passing the values of ηij respecting the full conditionals in Theorem 1. For the case ∑l=1Ki−jql is much larger than q0, the Markovian chain directly derived from Theorem 1 may converge slowly because the chance to introduce new values for ηij is quite small; see equation (12). Hence, an acceleration method (e.g., [11]) may be helpful, which samples another round of ηi for acceleration. This is realized by sampling ηi conditional on Xi and the cluster indicators c just obtained in Step 1 and then sampling c given Xi,ηi. In other words, Step 2 is driven by the grouped conditional distribution of ηi given Xi,c and c given Xi,ηi. Because of the particular feature that c is a deterministic function of ηi, Step 2 only needs to draw values from the conditional distribution of ηi given Xi,c.Moreover, for everyj=1,…,ni, the integration ∫fXij|ηijg^0rηijdηij involved in computing Gη⋅|Xij and ql,l=0,1,…,Ki−j is approximated by an average(18)hrj=1d∑l=1dfXij|ul,where d is a positive integer and u1,…,ud∼iidg^0r. Because u1,…,ud∼iidg^0r depend only on the iteration indicator r, this approximation can be prepared at the beginning of each iteration.As a final point, we would like to note that, due to the presence of the acceleration Step 2, the update of the values ofηij in Step 1 essentially takes effect only on the computation of the probabilities ql in equation (12) for every j=1,2,…,ni.The details of this two-step procedure of generatingηib are presented in Algorithm 2. For each group i=1,2,…,I, this procedure is repeated to generate ηib=ηi1b,…,ηinib,b=1,2,…,B, for a sufficiently large B.Algorithm 2: Gibbs sampling under single group.
(1)
Drawdiid valuesu1,u2,…,udfromg^0rand computehrj,j=1,…,niby (13).(2)
Initiate the algorithm with an arbitrarily selectedηi0. Correspondingly, setc=cηi0. Letη∗=ηi1∗,ηi2∗,…,ηiKi∗,0,…,0be the vector of dimensionniin which the firstKielements store the distinct value ofηi0andm=mi1,mi2,…,miKi,0,…,0be anni-vector to represent the numbers ofηij∗occurring inηi0.(3)
Now assume that we have gotηib−1and the corresponding auxiliary vectorsc,mandη∗. The sampled valuesηibare generated by the following two-step procedure:Step 1. Forj=1,2,…,ni,(1) setmicj=micj−1, Ki−j=Ki−1ifmicj=0Kiotherwise;(2) draw a new value forcj from 1,2,…,ni according to the distribution lawqc=Prcj=c|Xi,c∖cj,ηi1∗,…,ηiKi−j∗=ℬmic−α^rfXij|ηic∗for1≤c≤ni,mic≠0b/ni−Ki−jθ^r+α^rKi−jhrjfor1≤c≤ni,mic=0and then, in casemicj=0, draw a value for ηicj from Gη⋅|Xij with g0 replaced by g^0r, where ℬ is the normalization constant subject to ∑c=0niqc=1; and then(3) setmicj=micj+1.Step 2. For allc∈c1,c2,…,cni which are obtained by above Step 1, draw ηb from the conditional densityfηic|c,Xi∝g^0rηic∏j:cj=cfXij|ηic.(2) Estimation. In phase estimation, with ηb,b=1,2,…,B, in hand, the estimates of α,θ, and g0 are obtained according to (10).By these newly obtained estimates, repeat the above two phases for the next iteration until the sequence of estimates is stable. The above computation process on the estimation ofα,θ,g0 under the framework in (3) can be briefly arranged as in Algorithm 3.Algorithm 3: Gibbs sampling under multigroup.
Step 1. Initializeα^0,θ^0andg^00;Step 2. At roundr+1,r=1,2,…, for alli=1,2,…,I,(i) arbitrarily selection initial values forηi0 and(ii) updateηi with Algorithm 2 to produce ηb,b=1,2,…,B;Step 3. Obtainα^r+1,θ^r+1andg^0r+1according to (10).We now conclude this section by the following remark on one of Neal’s [11] algorithms.Remark 1.
Algorithm2 is similar to, but with two minor differences, to Algorithm 8 of Neal [11] at Step 1. First, we only draw a single set of d i.i.d. values from g^0r in each iteration, while Neal [11] suggested drawing d i.i.d. values for each j=1,…,ni, hence needing ni−1 times more than the sampling computation here. The algorithm here hence may cause a significant reduction in the total computation load. Second, we draw a value for ηicj from Gη⋅|Xij if a new cluster for cj is produced, while Neal [11] simply drew a value from g0r. According to Theorem 1, we think Algorithm 2 draws a correct value.


## 4. Simulation

In this section, some simulations are reported to show the performance of estimates provided above. The data structure of the PDP mixture model is according to equation (3) and the simulations were conducted under the two settings.(i)
We considered two settings: (1)I=100 with 20 individuals in every group: n1=n2=⋯=n100=20 and (2) I=1 with 200 individuals in a single group.(ii)
The amount of Gibbs sampling isB=20 and total iteration of estimation is 200.(iii)
The combinations ofα,θ,G0: (a) standard normal distribution 0.3,10,N0,1, (b) t distribution with 3 degrees of freedom, and (c) the distribution of −logX12. α took values of 0.3,0.5,0.8, respectively. θ=10. The initial values α^0,θ^0 were (0.5, 5) and (0.1, 1).(iv)
The conditional density functionf was normal with mean ηij and standard deviation was 0.1⋅SDG0, where SDG0 is the standard deviation of G0.In the three base probabilities we chose, the first two are symmetric,t3 has heavy tails, and the last has a positive skew. We tried a few distinct initial values of α,θ in all cases; the simulations showed that the estimations were insensible to the initial values. For every combination of the values, α, θ, I, and g0, we produced 1000 replicates of the datasets under the setting described above. The simulation results are concisely discussed below. .(1)
The estimates of the base measuresg0. In the simulation, g0 was estimated more accurately under the setting I=100 than under the setting I=1. The results from some replicates out of the 1000 are presented in Figures 3–5 under some parameters α,θ, and g0, in terms of the iterations r=0,5, and 200 under I=1 and I=200, in which the dashed lines represent the curves of estimates and the solids are the true density curves. The figures show the two following clear features: the dashed lines are more close to the solids curves with the increase of iteration times in both data structures and apparently the performance of estimated curves under multigroup data is better than that of single-group data. The center of dashed curves with single data is a little deviation at round 200 in Figure 3, the estimated curve is more concentrated than real density with heavy tails in Figure 4, and the curves of g^0 are more fluctuating and skewed in Figure 5 under the single-group data. Overall, in both the single-group and multigroup environments, our estimates are still near the true density curves. The histograms of L1 distance between g^0 and true g0 with 1000 replicates are shown in Figure 6 under single group and multiple groups. Figure 7 delineates the total variation distance (TVD) between G0 and G^0. The TVD of two probability distributions P and Q is SupAPA−QA, with the supremum taken over measurable sets A. Each panel plots the TVD between G^0 and G0 as the rounds progress. As can be seen, the value of TVD gradually reaches a steady value after about 15 iterations with the bases G0 of N0,1 and t3, and the other needs to go through about 25 iterations.(2)
The estimates ofα,θ with the 1000 replicates are summarized in Table 1 in terms of the mean and standard deviation and their histograms are presented in Figure 8. The comparisons of estimates are listed with means, SD, and 95% credible intervals and the estimates of multigroup structure are better than those of the single-group structure no matter what kind of bases in three parametric combinations. It is more stable on the estimations of α,θ in the standard deviations criterion when I=100.Figure 3
α,θ=0.3,10 and g0=N0,1 with the upper panels for I = 1 and the lower for I = 100.Figure 4
α,θ=0.5,10 and g0=t3 with the upper panels for I = 1 and the lower for I = 100.Figure 5
α,θ=0.8,10 and g0=−logχ12 with the upper panels for I=1 and the lower for I=100.Figure 6
Histograms of L1 distance with different bases and groups. The left plots are forI = 1 and the bases from top to bottom are N0,1t3 and −logX12.Figure 7
Progression ofTVDG0,G^0 with the base G0s of N0,1,t3 and −logX12.Figure 8
Histograms ofα^ and θ^ with the plots in rows 1 and 3 for I = 1 and the others for I = 100. The bases from left to right are N0,1t3 and −logX12.Table 1
Performance of estimates forα,θ.
BasesN0,1t3−logX12α,θ(0.3, 10)(0.5, 10)(0.8, 10)I = 1α^,θ^(0.355, 10.057)(0.533, 9.585)(0.793, 8.605)S.D.(0.076, 2.026)(0.058, 2.264)(0.020, 2.444)C.I.α^(0.231, 0.512)(0.426, 0.645)(0.755, 0.831)C.I.θ^(5.716, 17.408)(5.377, 13.992)(3.609, 12.970)I = 100α^,θ^(0.283, 9.978)(0.521, 9.662)(0.811, 9.547)S.D.(0.017, 0.430)(0.021, 0.752)(0.010, 1.351)C.I.α^(0.225, 0.338)(0.475, 0.563)(0.784, 0.835)C.I.θ^(8.942, 10.975)(8.256, 11.226)(6.597, 12.014)


## 5. Applications to Real Data

In this section, the PDP mixture model is employed to analyze the three following well-known real-world datasets.The first is the galaxies data [8, 25] consisting of velocities of 82 galaxies from six well-separated conic sections in Corona Borealis region. The model with I=1 was applied and the algorithm provided in Section 3 was run 200 rounds, giving rise to the estimate α^,θ^=0.429,3.988 and g^0 that is depicted in Figure 9, simply showing that the base distribution has multimodality corresponding to identifying different clumped galaxies.Figure 9
Galaxyg^0: km/s (thousands).The second dataset consists of records of stamp thickness [26] in millimeters of 485 unwatermarked Hidalgo stamps dating from 1872 to 1874. Figure 10 presents estimates g^0 of stamp thickness at round 200 of the procedure mentioned above on the stamp data, which exhibits the close bimodality of the density and one is more steep than another. The curve of g^0 reflects the high probability around the bimodality of ηi∗s corresponding to the number of types of paper used. The stamp data depicts the concentration of shared variables on the modes. After 200 rounds’ iterations, we obtain the estimates α^,θ^=0.425,6.810.Figure 10
Thickness of stampg^0: mm (hundreds).The last is on Danish fire losses [27] which are heterogeneous as reflected in multimodality, skewness, and heavy tail distributions, collected by Copenhagen Reinsurance and consist of 2492 records over the period of 1980–1990. The losses are in millions of Danish Krone and not adjusted for inflation over time. The data is divided into multiple groups based on time and the observations between groups are independent as the conventional assumption in the industry of nonlife insurance. The estimated density g^0 is clearly described with two intervals of abscissa in Figure 11 at round 200 under the Poisson-Dirichlet process mixture models, having one concentrated modality and heavy tail. We get α^,θ^=0.481,0.793.Figure 11
Danish fire loss datag^0: DKK (millions).


## 6. Conclusions

In this article, we have shown the posterior distribution in the PDP mixture models and the estimates cannot be directly obtained because of unobserved latent variablesηijs. The unbiased “pseudoestimates” are constructed and the corresponding algorithm is depicted in detail based on nonparametric empirical Bayes. Some simulations are presented with both data structures, single-group and multigroup, on the pseudoestimates of base distribution G0, discount parameter α, and concentration parameter θ. By the numerical simulation, the performances of estimations for parameters under multigroup data are better than those of the single-group data with the same total size of individuals. Also, we applied the model to the well-known real-world datasets.



---

*Source: 6618548-2021-10-21.xml*