(C) Physica-Verlag/Springer 2004

Masahiro Kuroda Key words : Data augmentation algorithm, local computation, graphical mod- els, computational efficiency, missing data. COMPSTAT 2004 section: MCMC. In this paper, we discuss an efficient Bayesian computational Abstract: method when observed data are incomplete in discrete graphical models. The data augmentation (DA) algorithm of Tanner and Wong [8] is applied to finding the posterior distribution. Utilizing the idea of local computation, it is possible to improve the DA algorithm. We propose a local computation DA (LC-DA) algorithm and evaluate its computational efficiency.
1
Introduction
A graphical model is characterized by conditional independence relationships among variables of a statistical model.  Graphical models are broadly used in various fields to describe complex statistical models and to specify the multivariate distributions, see Whittaker [9] and Edwards [4]. For a large graphical model, it is rare to o obtain complete observed data. When observed data are incomplete, it is extremely difficult to obtain the exact posterior distribution for a graphical model and the calculation may take a long time when observed data are moderately large. To overcome this computational difficulty, various algorithms related to graph structures are proposed, see Cowell et al. [1]. In this paper, we apply the data augmenta- tion (DA) algorithm of Tanner and Wong [8] to approximating the posterior distribution of a graphical model.  Then, incorporating the idea of local computation into the DA algorithm, it is possible to reduce the computational effort. We propose a local computation DA (LC-DA) algorithm and evaluate its computational efficiency. In Section 2, we show the exact Bayesian computation to find the posterior distribution for a discrete graphical model with missing data. In Section 3, instead of doing the infeasible exact computation, we give the DA algorithm to approximate the posterior distribution. In Section 4, we present the LC- DA algorithm.  Section 5 discusses the efficiency of the LC-DA algorithm from the viewpoint of computational complexity.
Masahiro Kuroda
2
Graphical model with missing data and exact
Bayesian computation
Let V denote the set of vertices in a graph and X V = { X i | i ∈ V } be the set of discrete random variables. Associated with each vertex i ∈ V , a random variable X i takes values in a sample space Ω i . For a subset A ⊆ V , we write X A for {X i | i ∈ A} and Ω A = Π i ∈ A Ω i .  Let t he joint probability of X V denote
for every xv ∈ Ωv = Π i ∈ v Ω i and let θv = {pv(xv) | xv ∈ Ωv}. The marginal probability of X A for A ⊂ V can be written as
for every x A ∈ Ω A = Π i ∈ A Ω i and θ A = { P A ( X A ) | X A ∈ Ω A }. The symbol "\" denotes the operator of a difference set. The conditional probability of X A given X B = x b is defined as PA|B(XA|XB) = Pr(XA = XA | XB = XB) = PAUB(XAUB)/PB(xB), providing AU B ⊂ V and A ∩ B = ∅ where ∅ denotes the empty set, and also θ A|B = {PA|B(XA|XB) | XA E ΩA,XB E ΩB}. In this paper, we assume that the graph of Xv has the global independence that, for a triplet ( A , B , C ) of mutually disjoint subsets of V  and V = AUBU C , each vertex of A is s separated from each vertex of B given the subset C . Then, under the global independence structure, X A is independent of X B given X C . Thus we have PV (XV) = PC (XC)PA|C (XA|XC)PB|C (XB|XC), so that {θc,θ λ c,θ B | C} are mutually independent. Suppose that observed data can be classified into three groups such that one is complete data and the others are incomplete data with X B and X A missing. The observed data patterns are indicated by T = {t 0,t 1,t 2} = {V, A U C, B U C}. In addition, we assume missingness at random in the sense of Rubin [7]. Each of observed data is denoted by n 0 = { n t 0 ( x t 0 ) | x t 0 ∈ Ω t 0 }, n 1 = { n t 1 ( x t 1 ) | x t 1 ∈ Ω t } and n 2 = { n t 2 ( x t 2 ) | x t 2 ∈ Ω t 2 } . The sizes of the incomplete data n 1 and n 2 are considerably larger than the size of the complete data n 0 . Assuming that observed data n = ( n 0 , n 1 , n 2 ) have a multinomial distribution with θ ν , the likelihood L(n|θν) is given by

Data Augmentation algorithm for graphical models
For the multinomial model, we assume that the prior distribution of θ v  is a Dirichlet distribution which has the density function
where αν = {αν(xν) | xv ∈ Ων} is a hyper-parameter. Then, accord- ing to the mutually independence relationships among {θ ς,θ λ ς,θ β | ς} , it is possible to factorize π(θν|αν) into
where ας = {ας(xc) | xc ∈ Ως}, αλς = {αλς(xλυς) | xλυς ∈ Ωλυς} and describes conditional independence of prior distributions and is called hyper Dirichlet prior distribution by Dawid and Lauritzen [3]. From the equations (1) and (2), we can obtain the mixture posterior distribution with the density
where, for i = 1,2, and ΣΩ(nt,(xt,)) denotes the sum over all possible ñ t,(xV) for all xV ∈ ΩV under the conditions ñ t i (xv) ≥ 0 and ∑ xv x , ñ t i (xv) = n t i (x t i ), and Because of combinational explosion, the posterior density (4) has a very com- plicated function.  Therefore, it is extremely difficult to calculate exactly π(θν | n ) and these computation may take a long time when the observed data are moderately large. Instead of performing the infeasible Bayesian computation, we use the DA algorithm which imputes incomplete data and finds π(θν | n ) using the Monte Carlo method.
1387
Masahiro Kuroda
3
DA algorithm to graphical models
The DA algorithm is a type of Markov chain Monte Carlo.  In the case that the incomplete-data posterior density is complicated as the posterior distribution (4) and the complete-data posterior is relative easy to handle and draw from, the DA algorithm is very suitable. Each iteration of the DA algorithm consists of an Imputation-step and a Posterior-step. For this case, the exact posterior distribution π(θ c | n ) can be obtained without any iterations, since the complete marginal data for X c are calcu- lated from n . Then the DA algorithm for the graphical model is given by the following iterative scheme: Initialization: Set an initial distribution π(o)(θν|n) = π(θν|αν). Imputation-step: Repeat the following steps for l = 1,..., L to obtain the augmented data ñ = (n o,ñ 1,ñ ), where
1. Generate θ ¥, from the current approximation π(t-1)(θν|n). 2. Generate ñ 1(l) and ñ 2(l) from the predictive multinomial distribu- tions f (ñ 1 θ̃Ḅḅḅn 1 ) and f (ñ 2̣lθ̃Ạt 2 , ñ 2 ), where
Posterior-step: Update π ( ′ )(θν | n ) given { ñ ( ′ ) | 1 ≤ l ≤ L } using the Monte Carlo method:
Then

Data Augmentation algorithm for graphical models
Until the approximations π ( t ) ( θ A | C | n ) and π ( t ) ( θ B | C | n ) converge to stationary distributions, the Imputation- and Posterior-steps are alternated repeatedly. Achieving convergence of the DA algorithm, the true posterior distribution π(θγ n ) can be found. In the practical implementation of the DA algorithm, the selection of the number of imputation ( L ) is perhaps more crucial. When the proportion of missing data in incomplete data is high and the size of the incomplete data is large, L must be considerably large. However, it is difficult to determine L on theoretical bases. In order to assess the convergence, diagnostic techniques are applied to output from the DA iteration. Cowles and Carlin [2] provide the comparative reviews of many convergence diagnostic techniques.
Application of LC-DA algorithm
44
In this section, we present the LC-DA algorithm. The important property of the LC-DA algorithm is that the DA algorithm is applied to each of factorized posterior distributions according to a graph structure and each posterior distribution is computed independently.  Then it is possible to reduce the computational efforts from the viewpoint of computational complexity. We denote the marginal data for X AUC and X BUC as n AC = (n AC,n 1,n 2,) and n bc = (n bc,n c,n 2), where n AC   =  {nto(xAUC) | xAUC E SAUC}, nBC = {nto(xBUC) | xBUC E SBUC}, n ξ  =  {n t,(xc)|xc ∈ Ωc}, n ξ = {nt,(xc)|xc ∈ Ωc}. With local computation to find π(θν | n ), we can obtain the following the- orem. Theorem 4.1. Suppose that C separates A and B in V . If C ⊆ t for all t ∈ T , then the calculation of the posterior distributions of θ A C  and θ B C can be done independently. Theorem 1 guarantees that the DA algorithm can execute separately to obtain the posterior distributions of θ A | C and θ B | C given n A C and n B C . The condition of C ⊂ t for all t ∈ T is caled “lossless decomposition” by Geng and Li [5]. Then the LC-DA algorithm realizes the computation according to the following iterative scheme:
The DA iteration of π(θ A C | n AC)
Initialization: Set an initial distribution π (0)(θ λ c | n Αc) = π (θ λ | c | α ν ). Imputation-step: Repeat the following steps for l = 1,..., L to impute ñ ≥ , where ñ²

1389
Masahiro Kuroda
1. Generate θ * | ς from the current approximation π (t−1)(θ λ c | n ac). 2. Generate the imputed data ñ2 (l) from the predictive multinomial distribution f(ñ²|θ*,ιc,n²). Posterior-step:  Update π (t-1) (θ λ c | n Α c) by the Monte Carlo method:
where The DA iteration of π(θ β [ C | π BC ) The DA algorithm to obtain π(θ β | C | π B c) are similar to the DA iteration of π(θ | C|n AC):  The Imputation-step generates {ñ 1 (l) | 1 ≤ l ≤ l ≤ l), nt,(xC),ñt,(xBuC) ≥ 0}.The Posterior-step finds π(t)(θB|C|nBc) using n Bc, n² and {ñ¹(l) | 1 ≤ l ≤ L}. When each of the approximations π ( t )(θ λ c | n λ c ) and π ( t )(θ β | c | n bc ) con- verges to a stationary distribution, the true posterior distribution π(θν n ) can be calculated.
Computational efficiency of LC-DA algorithm
5
We now evaluate the computational efficiency of the LC-DA algorithm from the viewpoint of computational complexity. Here we introduce two quantities: ● ||Ωv|| = the number of all possible values in Ωv ● ||Ω A || = the number of all possible values in Ω A where A ⊂ V As for the space complexity, the amount of the storage space required by the DA algorithm is ||Ωv||. Alternatively, in the LC-DA algorithm, it can not exceed max(||Ω AUC||, ||Ω BUC||). Next consider the time complexity under the worst-case assumption. The time complexity of the DA algorithm can be expressed by O (|| Ω v ||). The implementation of the LC-DA algorithm can be done in O (max(||Ω AUC||, ||Ω BUC||)). The LC-DA algorithm is more efficient than the DA algorithm from both aspects of the space and time complexities and then can reduce the compu- tational efforts. Finally, we briefly describe the convergence speed of the LC-DA algo- rithm.  The LC-DA algorithm is regarded as the collapsed Gibbs sampler of Liu [6]. Then, according to Liu's [6] result, the convergence speed of the LC-DA algorithm is faster than the speed of the DA algorithm.  We shall investigate its convergence speed in detail. Data Augmentation algorithm for graphical models
Appendix
Proof of Theorem 1 Since C ⊆ t for all t ∈ T , we have t ∩ C = C and P t ( x t ) π τ ( x t ) c e Cl
(5) L(θν|n) = L(θc|nc)L(θα|c|nac)L(θb|c|nbc). From the prior distribution (3) and the likelihood (5), we can obtain the posterior distribution π(θν n ) X {L(θc|πc)π(θc|αc)} × {L(θa|c|πac)π(θa|c|αac)} × { L(θB|C|nBc)π(θB|C|αBc)} = π(θc|nc)π(θa|c|nac)π(θb|c|nbc).
1391
Masahiro Kuroda
Since it also holds the mutual independence among the posterior distribu- tions, we can compute each posterior distribution independently. ■ [1] Cowell R.G., Dawid A.P., Lauritzen S.L., Spiegelhalter D.J. (1999). Probabilistic networks and expert systems. Springer-Verlag, New York. [2] Cowles M.K., Carlin D.B. (1996). Markov chain Monte Carlo conver- gence diagnostics: A comparative review. Journal of the American Sta- tistical Association 91, 883–904. [3] Dawid A.P., Lauritzen S.L. (1993). Hyper Dirichlet laws in the statistical analysis of decomposable graphical models. The Annals of Statistics 21, 1272–1317. [4] Edwards D. (2000). Introduction to graphical modeling. Second edition. Springer-Verlag, New York. [5] Geng Z., Li K. (2003). Factorization of posteriors and partial imputation algorithm for graphical models with missing data. Statistics & Probabil- ity Letters 64, 369–379. [6] Liu J. (1994). The collapsed Gibbs sampler in Bayesian computations with applications to a gene regulation problem.  Journal of the American Statistical Association 89, 958–966. [7] Rubin D.B. (1976). Inference and missing data. Biometrika 63, 581 – 592. [8] Tanner M.A., Wong W.H. (1987). The calculation of posterior distribu- tions by data augmentation. Journal of the American Statistical Associ- ation 82, 528 – 540. [9] Whittaker J. (1990). Graphical models in applied multivariate statistics. John Wiley & Sons, Chichester. Acknowledgement : The author would like to thank the editor and an anony- mous referee for providing valuable comments on the paper. Address : M. Kuroda, Department of Socio-Information, Okayama University of Science, 1-1 Ridai-cho, Okayama 700-0005, Japan E-mail: kuroda@soci.ous.ac.jp