Introduction
1.
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 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 augmentation (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 ©ICS AS CR 2005
343
Neural Network World 4/05, 343-350
of doing the infeasible exact computation, we give the DA algorithm to approxi- mate 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.
2.
Graphical model with missing data and exact
Bayesian computation
Let V denote the set of vertices in a graph and Xv = {Xi | i E V} be the set of discrete random variables. Associated with each vertex i E V, a random variable Xi takes values in a sample space Ωi. For a subset A C V, we write XA for {Xi | i E A} and Ω A = ΠiEA Ωi. Let the joint probability of Xv denote pv(xv) = Pr(Xv = xv), for every xv E Ωv = Π i ∈ v Ω i and let θν = {pv(xv) | xv ∈ Ωv}. The marginal probability of X A for A C V can be written as
for every xa E RA = fliEA Ri and BA = {PA(xA) | xA E RA}. The symbol "\" denotes the operator of a difference set. The conditional probability of X A given XB = xB is defined as PAIB(XAIXB) = Pr(XA = XA | XB = XB) = PAUB(XAUB)/PB(XB), providing AU B C V and A N B = 0 where 0 denotes the empty set, and also OA|B = {PA|B(XA|XB) | XA E SA,XB E SIB}. 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 = AUBUC, each vertex of A is 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 {0c,0alc,0blc} 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 XB and XA missing. The observed data patterns are indicated by T = {to, t1, t2} = {V, AUC, BUC}. In addition, we assume missingness at random in the sense of Rubin [7]. Each of observed data is denoted by n 0 = {nto(Xto) | Xto E Rto}, n 1 = {nt,(Xt,) | Xt,, E Rt,} and n 2 = {nt,,(Xt,) | Xt, E Rt,}.
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(nlθν) is given by
For the multinomial model, we assume that the prior distribution of θν is a Dirichlet distribution which has the density function
Then, according to where αν = {αν(xν) | xν ∈ Ων} is a hyper-parameter. the mutually independence relationships among {0c,0a|c,0b|c}, it is possible to factorize π(θν|αν) into
where ac = {ac(xc) | xc E IC}, aac = {aac(xauc) | xauc E Rauc} and abc = {abc(xbuc) | xbuc E {lbuc}. The Dirichlet prior distribution (3) de- scribes 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 Ln(n,,(x,,)) denotes the sum over all possible ñ,,(xv) for all xv E Rv under the conditions ñt,(xv) ≥ 0 and ∑zv\, ñt,(xv) = πt,(xt,), and âv(xv) = av(xv) + nto(xv) + ñt,(xv) + ñt,(xv). Because of combinational explosion, the posterior density (4) has a very compli- cated 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 al- gorithm which imputes incomplete data and finds π(θν | n ) using the Monte Carlo method.
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)
345
Neural Network World 4/05, 343-350
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 π(θς|n) can be obtained without any iterations, since the complete marginal data for Xc are calculated 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 aug- mented data ñ = (n°,ñ¹,ñ²), where
1. Generate θ i from the current approximation π(t-1)(θν|π). 2. Generate ñ 1 (1) and ñ 2 (1) from the predictive multinomial distributions f(ñ 1 | θ B|t,, n 1 ) and f(ñ 2 | θ ⋆ |t,, n 2 ), where
Posterior-step: Update 7r(t)(Ov ln) given {n(t) | 1 < l < L} using the Monte Carlo method:
Then
Until the approximations π ( ε ) (θ λ | ς | π ) and π ( ε ) (θ β | ς | π ) converge to stationary dis- tributions, the Imputation- and Posterior-steps are alternated repeatedly. Achiev- ing convergence of the DA algorithm, the true posterior distribution π(θν | n ) can be found. 346 In the practical implementation of the DA algorithm, the selection of the num- ber 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
4.
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 poste- rior 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 XAUC and XBUC as nAC = (nAC,n1,n2) and nbc = (n bc, n b, n 2 ), where n λς   =  {nto(xauc) | xauc ∈ Ωλυς}, πβς = {nto(xbuc) | xbuc ∈ Ωલυς}, n ξ   =   {n t, (xc) | xc ∈ Ωc},  n ξ = {n t, (xc) | xc ∈ Ωc}. With local computation to find π(θν | n ), we can obtain the following theorem. Theorem 4.1 Suppose that C separates A and B in V. If C C t for all t E T, then the calculation of the posterior distributions of θ A 1 C and θ B 1 C can be done independently. Theorem 4.1 guarantees that the DA algorithm can execute separately to obtain the posterior distributions of θ λ ι ς and θ β ι ς given n λ ς and n b ς . The condition of C C t for all t E T is called "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 π(θ λ | C | n λ C )
Initialization: Set an initial distribution π (0)(0A|C|nAC) = π(0A|C|0AC). Imputation-step: Repeat the following steps for l = 1,..., L to impute ñ2, where
1. Generate BAIC from the current approximation π(t-1)(BAICInAC). 2. Generate the imputed data ñ 2 (l) from the predictive multinomial distri- bution f(ñ²|θ A | C , n ² ).
347
Posterior-step: Update π (t-1)(θ A | C | n AC ) by the Monte Carlo method:
where ā')'C(XAUC) =  α AC (XAUC) + nto (XAUC) + nto (XAUC) + ñto (XAUC) The DA iteration of π(θ B | C | n BC ) The DA algorithm to obtain π(θ b | c | n bc ) are similar to the DA iteration of π(θ λ IC | n ΛC):  The Imputation-step generates {ñ l(l) | 1 ≤ l ≤ l ≤ L}, where
ñ λ = {ñ ι ( x buc) | x buc ∈ Ω buc , The Posterior-step finds π(t)(θ B I C I n BC ) using n BC , n 2 and {ñ I (t) | 1 ≤ l ≤ L}. When each of the approximations π ( έ ) ( θ λ | ς | π λ | ( θ b | ς | n b ς ) converges to a stationary distribution, the true posterior distribution π(θν | n ) can be calculated.
5.
Computational efficiency of LC-DA algorithm
We now evaluate the computational efficiency of the LC-DA algorithm from the viewpoint of computational complexity. Here we introduce two quantities: ● ||Ων|| = the number of all possible values in Ων • ||Ω 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 Next consider the time complexity under the exceed max(||SAUC||,||SBUC||). worst-case assumption. The time complexity of the DA algorithm can be ex- pressed by O(llflv ll). The implementation of the LC-DA algorithm can be done in 0(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 computational efforts. Finally, we briefly describe the convergence speed of the LC-DA algorithm. 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. [3] Dawid A. P., Lauritzen S. L.: Hyper Dirichlet laws in the statistical analysis of decomposable graphical models. The Annals of Statistics, 21, 1993, pp. 1272-1317. [4] Edwards D.: Introduction to graphical modeling. Second edition. Springer-Verlag, New York, 2000. [5] Geng Z., Li K.: Factorization of posteriors and partial imputation algorithm for graphical models with missing data. Statistics & Probability Letters, 64, 2003, pp. 369-379. [6] Liu J.: The collapsed Gibbs sampler in Bayesian computations with applications to a general regulation problem. Journal of the American Statistical Association, 89, 1994, pp. 958-966. [7] Rubin D. B.: Inference and missing data. Biometrika, 63, 1976, pp. 581-592. [8] Tanner M. A., Wong W. H.: The calculation of posterior distributions by data augmentation. Journal of the American Statistical Association, 82, 1987, pp. 528-540. [9] Whittaker J.: Graphical models in applied multivariate statistics. John Wiley & Sons, Chich- ester, 1990. Proof of Theorem 1 Since C ⊆ t for all t ∈ T, we have t ∩ C = C and pt(xt)nt(xt) x e n =
Then for any s C V, t N s = s, it holds nt(xs) = ñt(xs) and
Therefore it is possible to factorize the likelihood (1) as follows: L(θν|π) = L(θς|πς)L(θα|ς|πας)L(θв|ς|πας).
From the prior distribution (3) and the likelihood (5), we can obtain
Since it also holds the mutual independence among the posterior distributions, we can compute each posterior distribution independently. l