# Estimation of Poisson-Dirichlet Parameters with Monotone Missing Data

Authors: Xueqin Zhou; Jinlong Huang; Xianyi Wu
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This article considers the estimation of the unknown numerical parameters and the density of the base measure in a Poisson-Dirichlet process prior with grouped monotone missing data. The numerical parameters are estimated by the method of maximum likelihood estimates and the density function is estimated by kernel method. A set of simulations was conducted, which shows that the estimates perform well.


---

## Body

## 1. Introduction

As a young but fast-growing field of statistics, Bayesian nonparametrics (abbreviated as BNP below) focuses on Bayesian solutions of nonparametric and other infinite-dimensional statistical models. Compared with frequentist statistics and classical Bayesian statistics, BNP provides highly flexible and robust models for infinite-dimensional parameter spaces. The most extensively investigated priors in BNP include the famous Dirichlet process prior [1] and Polya tree process prior [2–5] which have played fundamental roles in the development of Bayesian nonparametrics.Dirichlet processes are also referred to as one-parameter Poisson-Dirichlet processes by Kingman [6]. As the first and foremost generalization of Dirichlet processes, two-parameter Poisson-Dirichlet processes (abbreviated as Poisson-Dirichlet process below) were first discussed by Pitman and Yor [7] and from then have made huge success in Bayesian nonparametric modeling in language, images, ecology, biology, genomics, and so on. Remarkable examples include Goldwater et al. [8] who used Poisson-Dirichlet process as an adaptor to justify the appearance of type frequencies in formal analyses of natural language and improved the performance of an earlier model for unsupervised learning of morphology, Sudderth and Jordan [9] who modeled the object frequencies and segment sizes by Poisson-Dirichlet processes and developed a statistical framework for the simultaneous, unsupervised segmentation and discovery of visual object categories from image databases, Favaro et al. [10] who used a Poisson-Dirichlet model to deal with the issue of prediction within species sampling problems, and Hoshino [11] who studied the microdata disclosure risk assessment with Pitman’s sampling formula, clarified some of its theoretical implications, and compared various models based on the Akaike Information Criterion by applying them to real data sets. For more references related to the application of Poisson-Dirichlet process in the area of language learning, one can be referred to Johnson et al. [12], Wood and Teh [13], and Wallach et al. [14].While the exact Bayesian methods take the assumption that prior distributions are completely specified, empirical Bayesian methods deal with the situations where prior distributions are at most partially specified and thus need to be estimated. Empirical Bayesian methods for parametric and semiparametric models have been investigated in a huge volume of literature. However, the study of empirical Bayesian methods in the framework of Bayesian nonparametrics is quite limited. A recently published paper is that by Yang and Wu [15] who studied the problem of estimation of the priors with monotonically missing data when the prior is a Dirichlet process.In this paper, we aim at estimating the unknown numerical parameters and the density of the base measure with independent and identically distributed (i.i.d.) groups of observations with a Poisson-Dirichlet process prior. Because the Dirichlet process prior is a special case of the Poisson-Dirichlet process prior, we in fact have extended the methodologies of Yang and Wu [15] to a bigger model. The estimation of the unknown parameters is carried out in two different methods, the maximum likelihood method and a naive method proposed by Carlton [16], of which performances are compared by a simulation study.Because there are two numerical parameters in Poisson-Dirichlet process priors, the maximum likelihood estimates (MLE) for the unknown parameters are discussed under three different settings (see the next section for the definition of the parametersα and θ and the density function): (i) the discount parameter α is unknown but the concentration parameter θ is known; (ii) the concentration parameter θ is unknown but the discount parameter α is known; and (iii) both α and θ are unknown. Favaro et al. [10] gave the empirical Bayes estimates when both α and θ were unknown on complete data without missing. The comparisons between the estimates of Favaro et al. [10] and ours are presented with the same sample size in terms of bias, standard deviations (SD), and mean squared errors (MSEs).The remainder of this paper is structured as follows. In Section2, we review the basic model and the definition of the Poisson-Dirichlet processes priors. Data structure and model assumptions are also described in this section. In Section 3, the MLEs of the prior parameters are discussed in detail under three different aforementioned situations. A naive estimate for the discount parameter α is also discussed. Section 4 discusses the estimates of base distribution density by kernel method. Section 5 presents a small simulation study to show the performance of the estimates discussed in Section 3.


## 2. The Data and Model

The data are observed inI time periods and accordingly organized in groups i=1,2,…,I, with i representing the calendar time point the individuals in this group begin to be observed. Each group i contains ni individuals. Denote by (i,k) the kth individual in group i which is represented by a random vector Xik with Xikj being the jth coordinate of individual (i,k). The observations are presented by J-dimensional (J≤I) vectors for which the coordinates are observed sequentially in time, so that the observations are subject to monotone missing. A clearer picture of the data structure is exhibited in Table 1.Table 1
Data structure.
Groups
1
2
⋯
I
-
i
+
1
⋯
J
-
1
J
1
X
111
X
112
⋯
X
11
(
I
-
i
+
1
)
⋯
X
11
(
J
-
1
)
X
11
J
X
121
X
122
⋯
X
12
(
I
-
i
+
1
)
⋯
X
12
(
J
-
1
)
X
12
J
⋮
⋮
⋯
⋮
⋯
⋮
⋮
X
1
n
1
1
X
1
n
1
2
⋯
X
1
n
1
(
I
-
i
+
1
)
⋯
X
1
n
1
(
J
-
1
)
X
1
n
1
J
2
X
211
X
212
⋯
X
21
(
I
-
i
+
1
)
⋯
X
21
(
J
-
1
)
X
21
J
X
221
X
222
⋯
X
22
(
I
-
i
+
1
)
⋯
X
22
(
J
-
1
)
X
22
J
⋮
⋮
⋯
⋮
⋯
⋮
⋮
X
2
n
2
1
X
2
n
2
2
⋯
X
2
n
2
(
I
-
i
+
1
)
⋯
X
2
n
2
(
J
-
1
)
X
2
n
2
J
⋮
⋮
⋮
⋯
⋮
⋯
⋮
⋮
I
-
J
+
1
X
(
I
-
J
+
1
)
11
X
(
I
-
J
+
1
)
12
⋯
X
(
I
-
J
+
1
)
1
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
1
)
1
(
J
-
1
)
X
(
I
-
J
+
1
)
1
J
X
(
I
-
J
+
1
)
21
X
(
I
-
J
+
1
)
22
⋯
X
(
I
-
J
+
1
)
2
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
1
)
2
(
J
-
1
)
X
(
I
-
J
+
1
)
2
J
⋮
⋮
⋯
⋮
⋯
⋮
⋮
X
(
I
-
J
+
1
)
n
(
I
-
J
+
1
)
1
X
(
I
-
J
+
1
)
n
(
I
-
J
+
1
)
2
⋯
X
(
I
-
J
+
1
)
n
(
I
-
J
+
1
)
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
1
)
n
(
I
-
J
+
1
)
(
J
-
1
)
X
(
I
-
J
+
1
)
n
(
I
-
J
+
1
)
J
I
-
J
+
2
X
(
I
-
J
+
2
)
11
X
(
I
-
J
+
2
)
12
⋯
X
(
I
-
J
+
2
)
1
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
2
)
1
(
J
-
1
)
X
(
I
-
J
+
2
)
21
X
(
I
-
J
+
2
)
22
⋯
X
(
I
-
J
+
2
)
2
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
2
)
2
(
J
-
1
)
⋮
⋮
⋯
⋮
⋯
⋮
X
(
I
-
J
+
2
)
n
(
I
-
J
+
2
)
1
X
(
I
-
J
+
2
)
n
(
I
-
J
+
2
)
2
⋯
X
(
I
-
J
+
2
)
n
(
I
-
J
+
2
)
(
I
-
i
+
1
)
⋯
X
(
I
-
J
+
2
)
n
(
I
-
J
+
2
)
(
J
-
1
)
⋮
⋮
⋮
⋯
⋮
i
X
i
11
X
i
12
⋯
X
i
1
(
I
-
i
+
1
)
X
i
21
X
i
22
⋯
X
i
2
(
I
-
i
+
1
)
⋮
⋮
⋯
⋮
X
i
n
i
1
X
i
n
i
2
⋯
X
i
n
i
(
I
-
i
+
1
)
⋮
⋮
⋮
I
-
1
X
(
I
-
1
)
11
X
(
I
-
1
)
12
X
(
I
-
1
)
21
X
(
I
-
1
)
22
⋮
⋮
X
(
I
-
1
)
n
I
-
1
1
X
(
I
-
1
)
n
I
-
1
2
I
X
I
11
X
I
21
⋮
X
I
n
I
1Hence, the observed data up to timeI are monotone missing; that is, the components of an individual (i,k) are ordered in such a way that if an observation of a variable Xikj for individual (i,k) is missing, then so are the observations of all subsequent variables Xikl,l≥j (if any) for the same individual. Clearly, Table 1 indicates that variables for individuals in the first I-J+1 groups are completely observed and only the first I-i+1 components are observable for individuals in group i=I-J+2,…,I.Data structured as in Table1 is frequently encountered in real world. A typical example is loss data for claims used for the purpose of loss reserving in the industry of non-life insurance (see, e.g., [15, 17, 18]). Assume that the evaluation time of loss reserving is accident year I; each claim made at accident year i(i≤I) will be paid at the end of the subsequent J years after the claim (the first is the one paid at the end of the year the claim is made) so that the payments of an individual (i,k) correspond to a J-vector (Xik1,…,Xikj,…,XikJ). Then observing at the end of the evaluation year I, the payments for a claim in accident year i have been paid only with js such that j≤I-i+1; that is, the observable payments of claim (i,k) are just the subvector Xiko=(Xik1,Xik2,…,Xik,I-i+1), where the superscript “o” means “observable.” Other example data of this structure can also be found in Marini et al. [19] who discussed the maximum likelihood estimation in panel studies for the first time, Hao and Krishnamoorthy [20] who considered testing and estimation problems, and Raats et al. [21] who discussed the estimation and testing for multivariate regression model, among a large number of others, with monotone missing data.Before presenting the probabilistic characteristics of the data under a BNP framework in Assumption2, we recall the definition of Poisson-Dirichlet processes in terms of stick-breaking. Let (α,θ) be a pair of real numbers satisfying α∈[0,1) and θ>-α, which is always respected and hence will not be mentioned everywhere.Definition 1.
(
1
) Let (Vn)n≥1 be a sequence of mutually independent random variables with Vn~Beta(1-α,θ+nα) and p~1=V1 and p~k=Vk∏i=1k-1(1-Vi),k≥2, so that ∑k=1∞p~k=1. Then the ranked sequence p=(p1,p2,…) (i.e., p1≥p2≥⋯) of p~=(p~1,p~2,…) is said to be Poisson-Dirichlet distributed, writing p~PD(α,θ), where α and θ are called discount and concentration parameter, respectively.
(
2
) Let (Yn)n≥1 be a sequence of random variables independent and identically distributed as a probability measure μ (calledbase distribution below), taking values in a measurable space (X,B), and p=(pn)n≥1~PD(α,θ). Then the random probability measure(1)Q=∑k=1∞pkδYk·on (X,B) is referred to as a Poisson-Dirichlet process (indexed by B) with parameters α,θ and base distribution μ(·) or, in symbol, writing Q~PD(μ;α,θ).Clearly, for givenp and (Yn)n≥1, one has a realization of the Poisson-Dirichlet process Q that is a discrete distribution with mass pn at Yn,n=1,2,… in the domain X. Pitman and Yor [7] discussed in detail Poisson-Dirichlet processes in which a number of well-known properties of Dirichlet processes were generalized. Carlton [16] gave methods to estimate the parameters (α,θ) with completely observed data.Assumption 2.
(a)
For everyi=1,2,…,I, the data (Xi1,Xi2,…,Xini) have a Poisson-Dirichlet structure: Pi~PD(μ;α,θ) and, given Pi, the data Xi1,Xi2,…,Xini~i.i.d.Pi.
(b)
The random elements(Xi1,Xi2,…,Xini;Pi) are mutually independent over i=1,…,I.
(c)
The probability measureμ has a continuous density function, denoted also by μ(x), with respect to the J-dimensional Lebesgue measure on RJ.The intuition behind the assumption is as what follows. In many situations, the data are recorded chronologically so that the first subscript ofXikj reflects the calendar time when the data vector begins to be recorded and the third subscript represents the calendar time that component is recorded; hence the dependence between the data beginning at a same calendar time and independent over calendar times is reasonable. This structure is the typical case of loss reserving in general insurance with individual data (see, e.g., Huang et al. [17, 18]).Given the above data structure and model assumption, the objective of this paper is to estimate the unknown parametersα and θ and the density μ(x) of the base measure.


## 3. Estimation of Parameters(α,θ)

The unknown parameters(α,θ) of the distribution of a Poisson-Dirichlet process can generally be estimated by means of maximum likelihood as shown in Section 3.2. To do it under the data structure of Table 1, we first introduce some necessary preliminary results for Poisson-Dirichlet processes in Section 3.1


### 3.1. Preliminaries

In the following, writeXikj=(Xik1,Xik2,…,Xikj),j=1,2,…,J, for the first j components of Xik and Pij and μj, respectively, for the marginal distributions of Pi and μ on the first j-dimensional subspace Rj of RJ. The following lemma shows that if the observations Xi1,Xi2,…,Xini have a Poisson-Dirichlet structure, then so are their subvectors.Lemma 3.
If the observationsXi1,Xi2,…,Xini∣Pi~i.i.d.Pi and Pi~PD(μ; α,θ), then, for every j=1,2,…,J, given Pij, the partial observations Xi1j,Xi2j,…,Xinij~i.i.d.Pij and Pij~PD(μj;α,θ).Proof.
The assertion thatPij~PD(μj;α,θ) follows immediately from the definition of PD(μ;α,θ). For any Bj-measurable set Bkj, k=1,2,…,ni, we have (2)PrXi1j∈B1j,…,Xinij∈Bnij∣Pij=EPrXi1j∈B1j,…,Xinij∈Bnij∣Pi∣Pij=E∏k=1niPiBkj×RJ-j∣Pij=∏k=1niPijBkj.It indicates the desired result Xi1j,Xi2j,…,Xinij∣Pij~i.i.d.Pij.The next lemma shows that, with a Poisson-Dirichlet process structure, observations of any two individuals are identical with probability of1 if and only if so are their corresponding components. This provides an easy and convenient way to judge whether any two observations are the same only by their any (e.g., the first) component.Lemma 4.
For anys,t∈{1,2,…,ni}, the equality Xis1=Xit1 implies that Xisj=Xitj for all j=1,2,…,J almost surely.Proof.
We verify the casej=2 and the proof applies for any j>2. By the definition of Poisson-Dirichlet process, we take Pi=∑n=1∞pnδYn with an i.i.d. sequence Yii≥1 drawn from μ. Because μ is continuous, then we have Pr(Yi=Yj for some i,j∈{1,2,…})=0. Therefore, (3)PrX1=X2=EPrX1=X2∣Pi=E∑n=1∞PrX1=X2=Yn∣Pi=E∑n=1∞pn2.Corollary 3 of Pitman and Yor [7] implies that(4)PrX1=X2=EW2W=EW=1-α1+θ,where W~Beta(1-α,θ+α). On the other hand, because Lemma 3 implies that X11 and X21 are two observations from a PD structure with the parameters α and θ and μ1(t1)=∫μ(t1,t2,…,tJ)dt2⋯dtJ, similar discussion as above implies that Pr(X11=X21)=(1-α)/(1+θ). Therefore,(5)PrX1=X2∣X11=X21=PrX1=X2PrX11=X21=1.The proof is thus completed.


### 3.2. Maximum Likelihood Estimation of the Unknown Parameters(α,θ)

Fori=1,2,…,I, let Ki be the number of distinct observations among the ni individuals of the ith group. Each distinct observation represents a species or cluster in biological perspective [7]. Let Aji be the number of those clusters, each of which appears j times in group i. For example, if Aji=3 for some j∈{1,2,…,ni}, then there are total 3 different species, each of which appears just j times in group i. Clearly,(
1
) there might exist many j∈{1,2,⋯,ni} with Aji=0,(
2
) generally,(6)Aji≥0forj=1,2,…,I,Ki=∑j=1niAji,∑j=1nijAji=ni.Write(7)Ai=A1i,A2i,…,Anii,i=1,2,…,I.Similar to Carlton [16] who used Ai to estimate the unknown parameters by maximum likelihood method when individuals are completely observed, we here discuss the MLE of (α,θ) in terms of the random variables Ai,i=1,2,…,I. Under Assumption 2, the log-likelihood function for (α,θ) given the observations Ai=ai in all the I groups is(8)lα,θ=∑i=1Ilog⁡Pr⁡Ai=ai=∑i=1Ilog⁡Nai-∑i=1I∑l=1ni-1log⁡θ+l+∑i=1I∑l=1Ki-1log⁡θ+lα+∑i=1I∑j=2niaji∑l=1j-1log⁡l-α,where N(ai)=ni!/∏j=1ni(j!)aji(aji)! is irrelevant to the parameter (α,θ).The estimation is analyzed in three different situations as what has been done in Carlton [16]:(i)α is unknown and θ is known.(ii)α is known and θ is unknown.(iii) Bothα and θ are unknown.Remark 5.
In the following, we assume thatKi≠1 and Ki≠ni for some i=1,2,…,I. Note that, in group i, Ki=1 means that all observations are the same and Ki=ni means that all the observations are different from one another. Thus, if Ki=1 or Ki=ni for all i=1,2,…,I, the MLE will be attained at the boundary of the parameter spaces. However, the joint events {Ki=ni,i=1,2,…,I} and {Ki=1,i=1,2,…,I} will become increasingly unlikely as n=∑i=1Ini tends to infinity; see Chapter 3 of Carlton [16].


#### 3.2.1. Estimateα with Known θ

In this case, note that(9)lαα,θ≔∂lα,θ∂α=∑i=1I∑l=1Ki-1lθ+lα-∑i=1I∑j=2niaji∑l=1j-11l-αis a decreasing function of α and lα(α,θ)→-∞ as α→1. The equality lα(α,θ)=0 has a unique solution if lα(0,θ)≥0, as stated in the theorem below.Theorem 6.
The MLEα^(θ) of α is unique and(10)α^θ>0ifθ<∑i=1IKiKi-12∑i=1I∑j=2niajiln⁡j-1+C,α^θ=0otherwise,where C is Euler’s constant.Proof.
Note that(11)lα0,θ=1θ∑i=1I∑l=1Ki-1l-∑i=1I∑j=2niaji∑l=1j-11l=12θ∑i=1IKiKi-1-∑i=1I∑j=2niajilnj-1+C.The proof is thus completed.


#### 3.2.2. Estimateθ with Known α

Ifα=0, then the Poisson-Dirichlet process reduces to a Dirichlet process, for which the estimation problem of θ can be referred to Yang and Wu [15]. We thus consider only the case of α∈(0,1). By the log-likelihood function l(α,θ) given in (8), we have(12)lθα,θ≔∂lα,θ∂θ=-∑i=1I∑l=1ni-11θ+l+∑i=1I∑l=1Ki-11θ+lα.A condition under which the MLE is unique is presented in the following theorem.Theorem 7.
Setkmax=max⁡{K1,K2,…,KI} and kmin=min⁡{K1,K2,…,KI}. Then the solution to the equation lθ(α,θ)=0 is unique if kmax-2kmin≤0.Proof.
For givenα∈(0,1), notice that the variable for lθ(α,θ) is only θ; we write(13)gθ=lθα,θ=∑i=1I∑l=1Ki-11θ+lα-1θ+l-∑l=Kini-11θ+l.Further recall that Ki>1 for i=1,2,…,I; it follows that(14)gθ=Iθ+α+∑i=1I∑l=2Ki-11θ+lα-∑l=1Ki-11θ+l-∑l=Kini-11θ+l,and consequently g(θ)→∞ as θ→-α for the second term is limited in the upper equation. With some algebraic transformation to (13), we obtain(15)gθ=∑i=1I1-α∑l=1Ki-1lθ+lαθ+l-∑l=Kini-11θ+l,which implies that, for sufficiently large θ,(16)gθ≤∑i=1IKi21-α-ni-Ki-1θ/2θ2<0.Then, by the continuity of g(θ), the Intermediate Value Theorem implies that c∈(-α,∞) such that g(c)=0.
For the uniqueness, it suffices to show thatg(θ)<0 whenever g′(θ)≥0. By (13), we have (17)g′θ=∑i=1I∑l=1ni-11θ+l2-∑l=1Ki-11θ+lα2=∑i=1Iα-1∑l=1Ki-1llα+l+2θθ+l2θ+lα2+∑l=Kini-11θ+l2,which implies that (18)g′θ<∑i=1I∑l=1Ki-1lα-1θ+l2θ+lα+∑l=Kini-11θ+l2.By the definition of kmax and kmin, (19)g′θ<∑i=1I1θ+Ki-1∑l=1Ki-12lα-1θ+lθ+lα+1θ+Ki∑l=Kini-11θ+l<2θ+kmax-1∑i=1I∑l=1Ki-1lα-1θ+lθ+lα+1θ+kmin∑i=1I∑l=Kini-11θ+l.The inequality is obtained when g′(θ)≥0:(20)∑i=1I∑l=1Ki-1lα-1θ+lθ+lα<θ+kmax-12θ+kmin∑i=1I∑l=Kini-11θ+l.With (15) and (20), we have (21)gθ<α-1+kmax-2kmin2θ+kmin∑i=1I∑l=Kini-11θ+l.When kmax-2kmin≤0, g(θ)<0. Particularly, there is a unique solution when I=1 for kmax=kmin=k1.The solution tolθ(α,θ)=0 is the MLE of θ and we denote it by θ^mle. When I=1, by Pitman and Yor [7], the PD(α,θ) distributions are mutually absolutely continuous when α∈(0,1) is fixed and θ is varying; hence θ^mle is not a consistent estimator. Some simulation results of θ^mle are reported in Section 5 when I>1; theoretically I can be ∞.


#### 3.2.3. Estimate Bothα and θ

As discussed above, the MLEs forα and θ when both are unknown can be obtained by solving the equations lα(α,θ)=lθ(α,θ)=0. To guarantee the existence of the MLEs for α and θ, we need (α^,θ^) of lα(α,θ)=lθ(α,θ)=0 to be such that(22)detHα,θ=lααlθθ-lθα2>0,where H(α,θ)≔lααlθαlαθlθθ is the Hessian matrix; that is, (23)lαα≔∂lα∂α=-∑i=1I∑l=1Ki-1l2θ+lα2-∑i=1I∑j=2niaji∑l=1j-11l-α2,lαθ≔∂lα∂θ=-∑i=1I∑l=1Ki-1lθ+lα2=∂lθ∂α≕lθα,lθθ≔∂lθ∂α=∑i=1I∑l=1ni-11θ+l2-∑i=1I∑l=1Ki-11θ+lα2.


#### 3.2.4. Remarks

While we have so far proved the existence of MLEs under the three different cases, analytical expressions of the estimates under none of the cases are given and, thus, the MLEs can only be obtained with numerical methods, for example, theNewton-Raphson iteration method ordichotomy method.In addition, by Carlton [16], a naive estimate of α in terms of merely the observations in group i is α^i=logKi/logni, which is consistent as ni tends to ∞. Now with the I groups of data, we can estimate α by(24)α~=1I∑i=1Iα^i=1I∑i=1Ilog⁡Kilog⁡niorα^w=∑i=1Ilog⁡Ki∑log⁡ni=∑i=1Ilog⁡Kilog⁡ni·log⁡ni∑i=1Ilog⁡ni=∑i=1Iωi·log⁡Kilog⁡ni,where ωi=log⁡ni/∑i=1Ilog⁡ni. It is then straightforward from Lemma 5.2 in Carlton [16] that both α~ and α^w are strongly consistent with the true value of α as max1≤i≤Ini→∞.


## 3.1. Preliminaries

In the following, writeXikj=(Xik1,Xik2,…,Xikj),j=1,2,…,J, for the first j components of Xik and Pij and μj, respectively, for the marginal distributions of Pi and μ on the first j-dimensional subspace Rj of RJ. The following lemma shows that if the observations Xi1,Xi2,…,Xini have a Poisson-Dirichlet structure, then so are their subvectors.Lemma 3.
If the observationsXi1,Xi2,…,Xini∣Pi~i.i.d.Pi and Pi~PD(μ; α,θ), then, for every j=1,2,…,J, given Pij, the partial observations Xi1j,Xi2j,…,Xinij~i.i.d.Pij and Pij~PD(μj;α,θ).Proof.
The assertion thatPij~PD(μj;α,θ) follows immediately from the definition of PD(μ;α,θ). For any Bj-measurable set Bkj, k=1,2,…,ni, we have (2)PrXi1j∈B1j,…,Xinij∈Bnij∣Pij=EPrXi1j∈B1j,…,Xinij∈Bnij∣Pi∣Pij=E∏k=1niPiBkj×RJ-j∣Pij=∏k=1niPijBkj.It indicates the desired result Xi1j,Xi2j,…,Xinij∣Pij~i.i.d.Pij.The next lemma shows that, with a Poisson-Dirichlet process structure, observations of any two individuals are identical with probability of1 if and only if so are their corresponding components. This provides an easy and convenient way to judge whether any two observations are the same only by their any (e.g., the first) component.Lemma 4.
For anys,t∈{1,2,…,ni}, the equality Xis1=Xit1 implies that Xisj=Xitj for all j=1,2,…,J almost surely.Proof.
We verify the casej=2 and the proof applies for any j>2. By the definition of Poisson-Dirichlet process, we take Pi=∑n=1∞pnδYn with an i.i.d. sequence Yii≥1 drawn from μ. Because μ is continuous, then we have Pr(Yi=Yj for some i,j∈{1,2,…})=0. Therefore, (3)PrX1=X2=EPrX1=X2∣Pi=E∑n=1∞PrX1=X2=Yn∣Pi=E∑n=1∞pn2.Corollary 3 of Pitman and Yor [7] implies that(4)PrX1=X2=EW2W=EW=1-α1+θ,where W~Beta(1-α,θ+α). On the other hand, because Lemma 3 implies that X11 and X21 are two observations from a PD structure with the parameters α and θ and μ1(t1)=∫μ(t1,t2,…,tJ)dt2⋯dtJ, similar discussion as above implies that Pr(X11=X21)=(1-α)/(1+θ). Therefore,(5)PrX1=X2∣X11=X21=PrX1=X2PrX11=X21=1.The proof is thus completed.


## 3.2. Maximum Likelihood Estimation of the Unknown Parameters(α,θ)

Fori=1,2,…,I, let Ki be the number of distinct observations among the ni individuals of the ith group. Each distinct observation represents a species or cluster in biological perspective [7]. Let Aji be the number of those clusters, each of which appears j times in group i. For example, if Aji=3 for some j∈{1,2,…,ni}, then there are total 3 different species, each of which appears just j times in group i. Clearly,(
1
) there might exist many j∈{1,2,⋯,ni} with Aji=0,(
2
) generally,(6)Aji≥0forj=1,2,…,I,Ki=∑j=1niAji,∑j=1nijAji=ni.Write(7)Ai=A1i,A2i,…,Anii,i=1,2,…,I.Similar to Carlton [16] who used Ai to estimate the unknown parameters by maximum likelihood method when individuals are completely observed, we here discuss the MLE of (α,θ) in terms of the random variables Ai,i=1,2,…,I. Under Assumption 2, the log-likelihood function for (α,θ) given the observations Ai=ai in all the I groups is(8)lα,θ=∑i=1Ilog⁡Pr⁡Ai=ai=∑i=1Ilog⁡Nai-∑i=1I∑l=1ni-1log⁡θ+l+∑i=1I∑l=1Ki-1log⁡θ+lα+∑i=1I∑j=2niaji∑l=1j-1log⁡l-α,where N(ai)=ni!/∏j=1ni(j!)aji(aji)! is irrelevant to the parameter (α,θ).The estimation is analyzed in three different situations as what has been done in Carlton [16]:(i)α is unknown and θ is known.(ii)α is known and θ is unknown.(iii) Bothα and θ are unknown.Remark 5.
In the following, we assume thatKi≠1 and Ki≠ni for some i=1,2,…,I. Note that, in group i, Ki=1 means that all observations are the same and Ki=ni means that all the observations are different from one another. Thus, if Ki=1 or Ki=ni for all i=1,2,…,I, the MLE will be attained at the boundary of the parameter spaces. However, the joint events {Ki=ni,i=1,2,…,I} and {Ki=1,i=1,2,…,I} will become increasingly unlikely as n=∑i=1Ini tends to infinity; see Chapter 3 of Carlton [16].


### 3.2.1. Estimateα with Known θ

In this case, note that(9)lαα,θ≔∂lα,θ∂α=∑i=1I∑l=1Ki-1lθ+lα-∑i=1I∑j=2niaji∑l=1j-11l-αis a decreasing function of α and lα(α,θ)→-∞ as α→1. The equality lα(α,θ)=0 has a unique solution if lα(0,θ)≥0, as stated in the theorem below.Theorem 6.
The MLEα^(θ) of α is unique and(10)α^θ>0ifθ<∑i=1IKiKi-12∑i=1I∑j=2niajiln⁡j-1+C,α^θ=0otherwise,where C is Euler’s constant.Proof.
Note that(11)lα0,θ=1θ∑i=1I∑l=1Ki-1l-∑i=1I∑j=2niaji∑l=1j-11l=12θ∑i=1IKiKi-1-∑i=1I∑j=2niajilnj-1+C.The proof is thus completed.


### 3.2.2. Estimateθ with Known α

Ifα=0, then the Poisson-Dirichlet process reduces to a Dirichlet process, for which the estimation problem of θ can be referred to Yang and Wu [15]. We thus consider only the case of α∈(0,1). By the log-likelihood function l(α,θ) given in (8), we have(12)lθα,θ≔∂lα,θ∂θ=-∑i=1I∑l=1ni-11θ+l+∑i=1I∑l=1Ki-11θ+lα.A condition under which the MLE is unique is presented in the following theorem.Theorem 7.
Setkmax=max⁡{K1,K2,…,KI} and kmin=min⁡{K1,K2,…,KI}. Then the solution to the equation lθ(α,θ)=0 is unique if kmax-2kmin≤0.Proof.
For givenα∈(0,1), notice that the variable for lθ(α,θ) is only θ; we write(13)gθ=lθα,θ=∑i=1I∑l=1Ki-11θ+lα-1θ+l-∑l=Kini-11θ+l.Further recall that Ki>1 for i=1,2,…,I; it follows that(14)gθ=Iθ+α+∑i=1I∑l=2Ki-11θ+lα-∑l=1Ki-11θ+l-∑l=Kini-11θ+l,and consequently g(θ)→∞ as θ→-α for the second term is limited in the upper equation. With some algebraic transformation to (13), we obtain(15)gθ=∑i=1I1-α∑l=1Ki-1lθ+lαθ+l-∑l=Kini-11θ+l,which implies that, for sufficiently large θ,(16)gθ≤∑i=1IKi21-α-ni-Ki-1θ/2θ2<0.Then, by the continuity of g(θ), the Intermediate Value Theorem implies that c∈(-α,∞) such that g(c)=0.
For the uniqueness, it suffices to show thatg(θ)<0 whenever g′(θ)≥0. By (13), we have (17)g′θ=∑i=1I∑l=1ni-11θ+l2-∑l=1Ki-11θ+lα2=∑i=1Iα-1∑l=1Ki-1llα+l+2θθ+l2θ+lα2+∑l=Kini-11θ+l2,which implies that (18)g′θ<∑i=1I∑l=1Ki-1lα-1θ+l2θ+lα+∑l=Kini-11θ+l2.By the definition of kmax and kmin, (19)g′θ<∑i=1I1θ+Ki-1∑l=1Ki-12lα-1θ+lθ+lα+1θ+Ki∑l=Kini-11θ+l<2θ+kmax-1∑i=1I∑l=1Ki-1lα-1θ+lθ+lα+1θ+kmin∑i=1I∑l=Kini-11θ+l.The inequality is obtained when g′(θ)≥0:(20)∑i=1I∑l=1Ki-1lα-1θ+lθ+lα<θ+kmax-12θ+kmin∑i=1I∑l=Kini-11θ+l.With (15) and (20), we have (21)gθ<α-1+kmax-2kmin2θ+kmin∑i=1I∑l=Kini-11θ+l.When kmax-2kmin≤0, g(θ)<0. Particularly, there is a unique solution when I=1 for kmax=kmin=k1.The solution tolθ(α,θ)=0 is the MLE of θ and we denote it by θ^mle. When I=1, by Pitman and Yor [7], the PD(α,θ) distributions are mutually absolutely continuous when α∈(0,1) is fixed and θ is varying; hence θ^mle is not a consistent estimator. Some simulation results of θ^mle are reported in Section 5 when I>1; theoretically I can be ∞.


### 3.2.3. Estimate Bothα and θ

As discussed above, the MLEs forα and θ when both are unknown can be obtained by solving the equations lα(α,θ)=lθ(α,θ)=0. To guarantee the existence of the MLEs for α and θ, we need (α^,θ^) of lα(α,θ)=lθ(α,θ)=0 to be such that(22)detHα,θ=lααlθθ-lθα2>0,where H(α,θ)≔lααlθαlαθlθθ is the Hessian matrix; that is, (23)lαα≔∂lα∂α=-∑i=1I∑l=1Ki-1l2θ+lα2-∑i=1I∑j=2niaji∑l=1j-11l-α2,lαθ≔∂lα∂θ=-∑i=1I∑l=1Ki-1lθ+lα2=∂lθ∂α≕lθα,lθθ≔∂lθ∂α=∑i=1I∑l=1ni-11θ+l2-∑i=1I∑l=1Ki-11θ+lα2.


### 3.2.4. Remarks

While we have so far proved the existence of MLEs under the three different cases, analytical expressions of the estimates under none of the cases are given and, thus, the MLEs can only be obtained with numerical methods, for example, theNewton-Raphson iteration method ordichotomy method.In addition, by Carlton [16], a naive estimate of α in terms of merely the observations in group i is α^i=logKi/logni, which is consistent as ni tends to ∞. Now with the I groups of data, we can estimate α by(24)α~=1I∑i=1Iα^i=1I∑i=1Ilog⁡Kilog⁡niorα^w=∑i=1Ilog⁡Ki∑log⁡ni=∑i=1Ilog⁡Kilog⁡ni·log⁡ni∑i=1Ilog⁡ni=∑i=1Iωi·log⁡Kilog⁡ni,where ωi=log⁡ni/∑i=1Ilog⁡ni. It is then straightforward from Lemma 5.2 in Carlton [16] that both α~ and α^w are strongly consistent with the true value of α as max1≤i≤Ini→∞.


## 3.2.1. Estimateα with Known θ

In this case, note that(9)lαα,θ≔∂lα,θ∂α=∑i=1I∑l=1Ki-1lθ+lα-∑i=1I∑j=2niaji∑l=1j-11l-αis a decreasing function of α and lα(α,θ)→-∞ as α→1. The equality lα(α,θ)=0 has a unique solution if lα(0,θ)≥0, as stated in the theorem below.Theorem 6.
The MLEα^(θ) of α is unique and(10)α^θ>0ifθ<∑i=1IKiKi-12∑i=1I∑j=2niajiln⁡j-1+C,α^θ=0otherwise,where C is Euler’s constant.Proof.
Note that(11)lα0,θ=1θ∑i=1I∑l=1Ki-1l-∑i=1I∑j=2niaji∑l=1j-11l=12θ∑i=1IKiKi-1-∑i=1I∑j=2niajilnj-1+C.The proof is thus completed.


## 3.2.2. Estimateθ with Known α

Ifα=0, then the Poisson-Dirichlet process reduces to a Dirichlet process, for which the estimation problem of θ can be referred to Yang and Wu [15]. We thus consider only the case of α∈(0,1). By the log-likelihood function l(α,θ) given in (8), we have(12)lθα,θ≔∂lα,θ∂θ=-∑i=1I∑l=1ni-11θ+l+∑i=1I∑l=1Ki-11θ+lα.A condition under which the MLE is unique is presented in the following theorem.Theorem 7.
Setkmax=max⁡{K1,K2,…,KI} and kmin=min⁡{K1,K2,…,KI}. Then the solution to the equation lθ(α,θ)=0 is unique if kmax-2kmin≤0.Proof.
For givenα∈(0,1), notice that the variable for lθ(α,θ) is only θ; we write(13)gθ=lθα,θ=∑i=1I∑l=1Ki-11θ+lα-1θ+l-∑l=Kini-11θ+l.Further recall that Ki>1 for i=1,2,…,I; it follows that(14)gθ=Iθ+α+∑i=1I∑l=2Ki-11θ+lα-∑l=1Ki-11θ+l-∑l=Kini-11θ+l,and consequently g(θ)→∞ as θ→-α for the second term is limited in the upper equation. With some algebraic transformation to (13), we obtain(15)gθ=∑i=1I1-α∑l=1Ki-1lθ+lαθ+l-∑l=Kini-11θ+l,which implies that, for sufficiently large θ,(16)gθ≤∑i=1IKi21-α-ni-Ki-1θ/2θ2<0.Then, by the continuity of g(θ), the Intermediate Value Theorem implies that c∈(-α,∞) such that g(c)=0.
For the uniqueness, it suffices to show thatg(θ)<0 whenever g′(θ)≥0. By (13), we have (17)g′θ=∑i=1I∑l=1ni-11θ+l2-∑l=1Ki-11θ+lα2=∑i=1Iα-1∑l=1Ki-1llα+l+2θθ+l2θ+lα2+∑l=Kini-11θ+l2,which implies that (18)g′θ<∑i=1I∑l=1Ki-1lα-1θ+l2θ+lα+∑l=Kini-11θ+l2.By the definition of kmax and kmin, (19)g′θ<∑i=1I1θ+Ki-1∑l=1Ki-12lα-1θ+lθ+lα+1θ+Ki∑l=Kini-11θ+l<2θ+kmax-1∑i=1I∑l=1Ki-1lα-1θ+lθ+lα+1θ+kmin∑i=1I∑l=Kini-11θ+l.The inequality is obtained when g′(θ)≥0:(20)∑i=1I∑l=1Ki-1lα-1θ+lθ+lα<θ+kmax-12θ+kmin∑i=1I∑l=Kini-11θ+l.With (15) and (20), we have (21)gθ<α-1+kmax-2kmin2θ+kmin∑i=1I∑l=Kini-11θ+l.When kmax-2kmin≤0, g(θ)<0. Particularly, there is a unique solution when I=1 for kmax=kmin=k1.The solution tolθ(α,θ)=0 is the MLE of θ and we denote it by θ^mle. When I=1, by Pitman and Yor [7], the PD(α,θ) distributions are mutually absolutely continuous when α∈(0,1) is fixed and θ is varying; hence θ^mle is not a consistent estimator. Some simulation results of θ^mle are reported in Section 5 when I>1; theoretically I can be ∞.


## 3.2.3. Estimate Bothα and θ

As discussed above, the MLEs forα and θ when both are unknown can be obtained by solving the equations lα(α,θ)=lθ(α,θ)=0. To guarantee the existence of the MLEs for α and θ, we need (α^,θ^) of lα(α,θ)=lθ(α,θ)=0 to be such that(22)detHα,θ=lααlθθ-lθα2>0,where H(α,θ)≔lααlθαlαθlθθ is the Hessian matrix; that is, (23)lαα≔∂lα∂α=-∑i=1I∑l=1Ki-1l2θ+lα2-∑i=1I∑j=2niaji∑l=1j-11l-α2,lαθ≔∂lα∂θ=-∑i=1I∑l=1Ki-1lθ+lα2=∂lθ∂α≕lθα,lθθ≔∂lθ∂α=∑i=1I∑l=1ni-11θ+l2-∑i=1I∑l=1Ki-11θ+lα2.


## 3.2.4. Remarks

While we have so far proved the existence of MLEs under the three different cases, analytical expressions of the estimates under none of the cases are given and, thus, the MLEs can only be obtained with numerical methods, for example, theNewton-Raphson iteration method ordichotomy method.In addition, by Carlton [16], a naive estimate of α in terms of merely the observations in group i is α^i=logKi/logni, which is consistent as ni tends to ∞. Now with the I groups of data, we can estimate α by(24)α~=1I∑i=1Iα^i=1I∑i=1Ilog⁡Kilog⁡niorα^w=∑i=1Ilog⁡Ki∑log⁡ni=∑i=1Ilog⁡Kilog⁡ni·log⁡ni∑i=1Ilog⁡ni=∑i=1Iωi·log⁡Kilog⁡ni,where ωi=log⁡ni/∑i=1Ilog⁡ni. It is then straightforward from Lemma 5.2 in Carlton [16] that both α~ and α^w are strongly consistent with the true value of α as max1≤i≤Ini→∞.


## 4. The Estimation of Base Distribution

In this section, we estimate the base distributionμ of PD(μ;α,θ) under the monotone missing data structure. In Section 4.1, we pick up a set of i.i.d. observations from μ and then estimate the density of μ by kernel estimation method. With the i.i.d. observations, we can estimate the base density μ by multivariate kernel density method that is similar to the one used by Yang and Wu [15] who studied the multivariate density estimation in a Dirichlet process prior under the same missing mechanism and showed that their estimate is superior to Titterington and Mill [22] under asymptotic mean integrated squared errors (AMISE) criterion (note that the classical versions of [22] for complete data can be found in, e.g., [23, 24]).


### 4.1. Deduce i.i.d. Observations fromμ

Denote byZi1,Zi2,…,Zi,Ki the distinct observations of Xi1,Xi2,…,Xi,ni. The following theorem shows that, given Ki, Zi1,Zi2,…,Zi,Ki are independent and identically μ(t1,t2,…,tJ)-distributed—note that, by Lemmas 3 and 4, under the structure of this model, for every group i, though the data are partly missing, the value of Ki can be easily worked out.Theorem 8.
Suppose thatμ(t1,t2,…,tJ) is continuous; then, given Ki, the distinct observations Zi1, Zi2,…,Zi,Ki are i.i.d. with distribution μ(t1,t2,…,tJ).Proof.
The proof can be found in Theorem2.5 in Korwar and Hollander [25].By Theorem8, there are Ki i.i.d. observations in group i and hence the number of i.i.d. observations in the whole I groups is K=∑i=1IKi. Note that the observed part of the K i.i.d. observations, namely,(25)Zik1,Zik2,…,Zik,I-i∧J:k=1,2,…,Ki,is still with a monotone missing structure. For every j=1,2,…,J, write Zikj=(Zik1,Zik2,…,Zikj) for the subvector of the first j components of Zik from individual (i,k),k=1,2,…,Ki, and(26)Sj=Zimj:1≤m≤Ki,1≤i≤I-j+1for the maximal set of the subvectors Zikj,k=1,2…,Ki, i=1,2,…,I, which are completely observed. Obviously, SJ is made of all the observations with no missing components and S1 is the set of the first components of all Zik in the I groups. Similarly, define(27)S~j=Zimj:1≤m≤Ki,1≤i≤I-j.Obviously, S~j⊂Sj, for all j.


### 4.2. Kernel Density Estimation

There are two methods to estimate the joint densityμ of the random vector Zik, proposed by Yang and Wu [15] and Titterington and Mill [22], respectively. They are described below.(
1) Yang and Wu’s Method. Let {ωj,j=1,2,…,J} be a sequence ofconsistent kernel functions—a terminology used in Yang and Wu [15] to indicate that ωj are density functions on space Rj satisfying ωj-1(t1,t2,…,tj-1)=∫ωj(t1,t2,…,tj-1,tj)dtj for j=2,…,J. Then, under the monotone missing data structure, μj can be estimated in two different ways: (i) using the consistent property of kernel ωj+1 and the data Sj+1 and (ii) using the ordinary kernel density estimation with data Sj∖S~j. Formally, the two estimates (denoted by μ^j1 and μ^j2, resp.) are given by(28)μ^j1t1,t2,…,tj=∫1hj+1j+1∑i=1I-jKi∑i=1I-j∑k=1Kiωj+1t1-Zik1hj+1,…,tj+1-Zikj+1hj+1dtj+1=1hj+1j∑i=1I-jKi∑i=1I-j∑k=1Kiωjt1-Zik1hj+1,t2-Zik2hj+1,…,tj-Zikjhj+1,μ^j2t1,…,tj=1hjjKI-j+1∑k=1KI-j+1ωjt1-ZI-j+1k1hj,…,tj-Zi-j+1kjhj.Therefore, the density μ(t1,t2,…,tj) for j=1,2,…,J-1 can be estimated by(29)μ^t1,t2,…,tj=∑i=1I-jKiμ^j1t1,…,tj+KI-j+1μ^j2t1,…,tj∑i=1I-j+1Ki.And μ(t1,t2,…,tJ) can be estimated by μ^J2(t1,t2,…,tJ). Furthermore, the consistence of ωj also indicates that(30)μ^tj∣t1,t2,…,tj-1=μ^t1,t2,…,tj∫μ^t1,t2,…,tjdtjforj=2,3,…,J.By a decomposition of the joint densityμ as products of a series of conditional densities, that is, μ(t1,t2,…,tJ)=μ(t1)∏j=1Jμ(tj∣t1,t2,…,tj-1), one can estimate the density μ by(31)μ^mt1,t2,…,tJ=μ^t1∏j=2Jμ^tj∣t1,t2,…,tj-1,where the superscript “m" indicates that the estimate is obtained by the monotone missing data structure and μ^(t1) and μ^(tj∣t1,t2,…,tj-1) are, respectively, proper estimates for μ(t1) and μ(tj∣t1,t2,…,tj-1).(
2) Titterington and Mill’s Method. For a vector X=(Xo,Xm), Titterington and Mill [22] used ∫K(·∣(xo,xm),h)f^(xm∣xo)dxm instead of the usual kernel K(·∣xo) when estimating the density of Xo, where Xo and Xm are, respectively, the observed part and missing part of X and f^(xm∣xo) is a probability density function to estimate true f(xm∣xo).With the above idea, dataSj, and group (I-j+2), another estimate for μj is given by(32)μ^t1,t2,…,tj=∑i=1I-j+1Kiμ^j3t1,t2,…,tj+KI-j+2μ^j4t1,t2,…,tj∑i=1I-j+2Ki,where (33)μ^j3t1,t2,…,tj=1hjj∑i=1I-j+1Ki∑i=1I-j+1∑k=1Kiωjt1-Zik1hj,t2-Zik2hj,…,tj-Zikjhj,μ^j4t1,t2,…,tj=1hj+1jKI-j+2∑k=1KI-j+2∫ωj+1t1-ZI-j+2k1hj+1,…,tj+1-zI-j+2kj+1hj+1f^jzI-j+2kj+1∣zI-j+2k1,zI-j+2k2,…,zI-j+2kjdzI-j+2kj+1.We should note that f^j(tj+1∣t1,t2,…,tj) is an appropriate estimate of the density of Z(I-j+2)k(j+1) given (Z(I-j+2)k1, Z(I-j+2)k2,…,Z(I-j+2)kj).


## 4.1. Deduce i.i.d. Observations fromμ

Denote byZi1,Zi2,…,Zi,Ki the distinct observations of Xi1,Xi2,…,Xi,ni. The following theorem shows that, given Ki, Zi1,Zi2,…,Zi,Ki are independent and identically μ(t1,t2,…,tJ)-distributed—note that, by Lemmas 3 and 4, under the structure of this model, for every group i, though the data are partly missing, the value of Ki can be easily worked out.Theorem 8.
Suppose thatμ(t1,t2,…,tJ) is continuous; then, given Ki, the distinct observations Zi1, Zi2,…,Zi,Ki are i.i.d. with distribution μ(t1,t2,…,tJ).Proof.
The proof can be found in Theorem2.5 in Korwar and Hollander [25].By Theorem8, there are Ki i.i.d. observations in group i and hence the number of i.i.d. observations in the whole I groups is K=∑i=1IKi. Note that the observed part of the K i.i.d. observations, namely,(25)Zik1,Zik2,…,Zik,I-i∧J:k=1,2,…,Ki,is still with a monotone missing structure. For every j=1,2,…,J, write Zikj=(Zik1,Zik2,…,Zikj) for the subvector of the first j components of Zik from individual (i,k),k=1,2,…,Ki, and(26)Sj=Zimj:1≤m≤Ki,1≤i≤I-j+1for the maximal set of the subvectors Zikj,k=1,2…,Ki, i=1,2,…,I, which are completely observed. Obviously, SJ is made of all the observations with no missing components and S1 is the set of the first components of all Zik in the I groups. Similarly, define(27)S~j=Zimj:1≤m≤Ki,1≤i≤I-j.Obviously, S~j⊂Sj, for all j.


## 4.2. Kernel Density Estimation

There are two methods to estimate the joint densityμ of the random vector Zik, proposed by Yang and Wu [15] and Titterington and Mill [22], respectively. They are described below.(
1) Yang and Wu’s Method. Let {ωj,j=1,2,…,J} be a sequence ofconsistent kernel functions—a terminology used in Yang and Wu [15] to indicate that ωj are density functions on space Rj satisfying ωj-1(t1,t2,…,tj-1)=∫ωj(t1,t2,…,tj-1,tj)dtj for j=2,…,J. Then, under the monotone missing data structure, μj can be estimated in two different ways: (i) using the consistent property of kernel ωj+1 and the data Sj+1 and (ii) using the ordinary kernel density estimation with data Sj∖S~j. Formally, the two estimates (denoted by μ^j1 and μ^j2, resp.) are given by(28)μ^j1t1,t2,…,tj=∫1hj+1j+1∑i=1I-jKi∑i=1I-j∑k=1Kiωj+1t1-Zik1hj+1,…,tj+1-Zikj+1hj+1dtj+1=1hj+1j∑i=1I-jKi∑i=1I-j∑k=1Kiωjt1-Zik1hj+1,t2-Zik2hj+1,…,tj-Zikjhj+1,μ^j2t1,…,tj=1hjjKI-j+1∑k=1KI-j+1ωjt1-ZI-j+1k1hj,…,tj-Zi-j+1kjhj.Therefore, the density μ(t1,t2,…,tj) for j=1,2,…,J-1 can be estimated by(29)μ^t1,t2,…,tj=∑i=1I-jKiμ^j1t1,…,tj+KI-j+1μ^j2t1,…,tj∑i=1I-j+1Ki.And μ(t1,t2,…,tJ) can be estimated by μ^J2(t1,t2,…,tJ). Furthermore, the consistence of ωj also indicates that(30)μ^tj∣t1,t2,…,tj-1=μ^t1,t2,…,tj∫μ^t1,t2,…,tjdtjforj=2,3,…,J.By a decomposition of the joint densityμ as products of a series of conditional densities, that is, μ(t1,t2,…,tJ)=μ(t1)∏j=1Jμ(tj∣t1,t2,…,tj-1), one can estimate the density μ by(31)μ^mt1,t2,…,tJ=μ^t1∏j=2Jμ^tj∣t1,t2,…,tj-1,where the superscript “m" indicates that the estimate is obtained by the monotone missing data structure and μ^(t1) and μ^(tj∣t1,t2,…,tj-1) are, respectively, proper estimates for μ(t1) and μ(tj∣t1,t2,…,tj-1).(
2) Titterington and Mill’s Method. For a vector X=(Xo,Xm), Titterington and Mill [22] used ∫K(·∣(xo,xm),h)f^(xm∣xo)dxm instead of the usual kernel K(·∣xo) when estimating the density of Xo, where Xo and Xm are, respectively, the observed part and missing part of X and f^(xm∣xo) is a probability density function to estimate true f(xm∣xo).With the above idea, dataSj, and group (I-j+2), another estimate for μj is given by(32)μ^t1,t2,…,tj=∑i=1I-j+1Kiμ^j3t1,t2,…,tj+KI-j+2μ^j4t1,t2,…,tj∑i=1I-j+2Ki,where (33)μ^j3t1,t2,…,tj=1hjj∑i=1I-j+1Ki∑i=1I-j+1∑k=1Kiωjt1-Zik1hj,t2-Zik2hj,…,tj-Zikjhj,μ^j4t1,t2,…,tj=1hj+1jKI-j+2∑k=1KI-j+2∫ωj+1t1-ZI-j+2k1hj+1,…,tj+1-zI-j+2kj+1hj+1f^jzI-j+2kj+1∣zI-j+2k1,zI-j+2k2,…,zI-j+2kjdzI-j+2kj+1.We should note that f^j(tj+1∣t1,t2,…,tj) is an appropriate estimate of the density of Z(I-j+2)k(j+1) given (Z(I-j+2)k1, Z(I-j+2)k2,…,Z(I-j+2)kj).


## 5. Simulation

In this section, we report a small numerical simulation regarding the performance of the estimates ofα and θ given in Section 3. The simulation was proceeded under the following settings:(i)
I=10 and because the estimates depend only on the observations of Kis which can be determined by the first component, we set J=1 for convenience.(ii)
The standard normal distributionN(0,1) was taken as the base measure μ.(iii)
The sizes of individuals(n1,n2,…,n10) took three levels of values:level 1:n1=50,n2=51,…,n9=58,n10=59,level 2:n1=500,n2=510,…,n9=580,n10=590,level 3:n1=100,n2=200,…,n9=900,n10=1000.The procedure of the simulation under those settings is simply described as follows.Algorithm 9.
​Step 1.
Generate the dataXi1,Xi2,…,Xini for every group i=1,2,…,I by means of a generalized Polya urn model as follows: (i)
Xi1∣α,θ,μ~μ(·)=N(0,1).(ii)
Xij,j=2,3,…,ni, is generated by the prediction rule (34)Xij+1∣Xi1,…,Xij,α,θ,μ~θ+αKijθ+jμ·+1θ+j∑l=1Kijml-αδZil·,whereZil is the lth distinct value among Xi1,…,Xij, ml=#{r:Xir=Zil,1≤r≤j} is the number of occurrences of Zil, and Kij is the number of distinct values of Xi1,…,Xij.Step 2.
Compute the values ofKi and aji from the data generated in Step 1.Step 3.
Compute the estimates ofα,θ with Newton-Raphson method.The simulation compared the estimatesα^mle and α^w of α for known θ in which we set θ=10 and α took three values: 0.1, 0.3, and 0.5; the case where α>0.5 can simply be obtained by replacing α with 1-α. For each combination of the parameters, the simulation was proceeded 1000 replications. The performances of the two estimates are shown in Table 2, from which it is clearly shown that the simulated biases of α^mle are all smaller than those of α^w. Although the simulated standard deviations of α^w are stable and smaller in three levels under three different values, the MSEs of α^mle are smaller than those of α^w. That is, α^mle performs better than α^w under the MSE criterion.Table 2
Performance of estimates forα with θ=10.
Size
α
^
Bias
SD
MSE
0.1
0.3
0.5
0.1
0.3
0.5
0.1
0.3
0.5
Level 1
α
^
m
l
e
0.091
0.001
0.078
0.043
0.014
0.036
0.010
0.000
0.007
α
^
w
0.530
0.507
0.150
0.021
0.004
0.036
0.282
0.257
0.024
Level 2
α
^
m
l
e
0.040
0.008
−0.001
0.005
0.005
0.018
0.002
0.000
0.000
α
^
w
0.518
0.391
0.104
0.002
0.003
0.024
0.269
0.153
0.011
Level 3
α
^
m
l
e
0.042
0.009
−0.001
0.006
0.005
0.004
0.002
0.000
0.090
α
^
w
0.525
0.397
0.276
0.003
0.002
0.003
0.276
0.157
0.331Next, the simulation studied the performance of the estimates ofθ with α=0.5 at four true values: θ=0.1,1,5, and 10. The simulation was also proceeded 1000 replications for each combination of the parameters. The simulation results are summarized in Table 3, from which we can see that the estimates θ^mle generally performed better with the increasing of the sample size.Table 3
Performance of the estimate ofθ with α=0.5.
Sample size
θ
0.1
1
5
10
Bias
Level 1
−0.018
0.067
0.107
0.047
Level 2
0.019
0.052
0.022
0.081
Level 3
0.012
0.012
0.114
0.032
SD
Level 1
0.230
0.429
0.971
1.653
Level 2
0.210
0.369
0.746
1.105
Level 3
0.202
0.369
0.779
1.187
MSE
Level 1
0.053
0.189
0.954
2.821
Level 2
0.044
0.138
0.556
1.228
Level 3
0.041
0.137
0.620
1.410Finally, we considered the MLEs forα and θ when neither of the parameters is known. The true values of the parameters and the performance of the estimates of the 1000 replications are shown in Table 4. The simulation used the Newton- Raphson iteration method.Table 4
Performance of the estimates for(α,θ) under level 1.
(
α
,
θ
)
(
α
^
,
θ
^
)
Bias
SD
MSEs
(0.1, 10)
(
α
^
,
θ
^
)
m
l
e
(0.064, −1.209)
(0.095, 2.098)
(0.013, 5.864)
(
α
^
,
θ
^
)
F
(0.153, −3.762)
(0.380, 3.644)
(0.038, 27.429)
(0.2, 4)
(
α
^
,
θ
^
)
m
l
e
(0.171, −1.621)
(0.100, 1.193)
(0.039, 4.052)
(
α
^
,
θ
^
)
F
(0.207, −2.200)
(0.106, 1.757)
(0.054, 7.925)
(0.25, 6)
(
α
^
,
θ
^
)
m
l
e
(0.049, −0.663)
(0.109, 2.080)
(0.014, 4.765)
(
α
^
,
θ
^
)
F
(0.105, −1.637)
(0.118, 3.457)
(0.025, 14.626)
(0.3, 7)
(
α
^
,
θ
^
)
(0.011, −0.094)
(0.100, 2.234)
(0.010, 5.001)
(
α
^
,
θ
^
)
F
(0.039, −0.121)
(0.0.117, 4.852)
(0.015, 23.557)
(0.4, 8)
(
α
^
,
θ
^
)
m
l
e
(−0.009, 0.429)
(0.095, 2.867)
(0.009, 8.405)
(
α
^
,
θ
^
)
F
(−0.020, 1.083)
(0.103, 5.242)
(0.011, 28.653)
(0.5, 5)
(
α
^
,
θ
^
)
m
l
e
(−0.001, 0.263)
(0.077, 2.160)
(0.006, 4.734)
(
α
^
,
θ
^
)
F
(−0.008, 1.434)
(0.084, 4.528)
(0.007, 22.562)Favaro et al. [10] adopted an empirical Bayes procedure for estimating (α,θ), writing (α^,θ^)F, when the data is not divided into groups. They investigated the issue of prediction within species sampling problems and applications to genomics based on the estimates (α^,θ^)F when both parameters are unknown. Favaro et al. [10] also studied the asymptotic behavior of the number of new species conditionally on the observed sample and derived asymptotic highest posterior density intervals for the estimates of their interest. The comparisons of simulation results between (α^,θ^)F and the MLEs which are under the multigroups data structure are presented in Table 4 with the total size of individuals in the perspective of bias, SD, and MSEs. From Table 4, we can see that the MSE of α^ is very small and the MSE of θ^ is a bit bigger. Apparently, the simulation shows that the MLEs generally are better than (α^,θ^)F of single group data under the critera of bias, SD, and MSEs.


## 6. Conclusion

In this paper, we studied the estimates of the unknown parameters as well as the density of the base measure in Poisson-Dirichlet process priors under monotonically missing data. The parameters are estimated by the method of MLEs and a set of simulations shows that the estimates perform well.



---

*Source: 7892507-2017-10-12.xml*