Abstract

1. Introduction
The expectation and maximization (EM) algorithm proposed by Dempster et al. (1977) is often used to find maximization likelihood estimates (MLEs) for log-linear models with partially classified data because of its stability, flexibility and simplicity, see Fuchs (1982) and Chen et al. (1984). However the EM algorithm converges slowly when there is a relatively large proportion of partially classified data. Many algorithms have been proposed to speed up the convergence of the EM algorithm, see McLachlan and Krishnan (1997). Newton-type accelerators require the computation of a matrix inversion at each iteration of the EM algorithm. The computation becomes increasingly more complicated as the number of parameters increases, and it also becomes numerically unstable. In this paper, we apply the Aitken δ² method of Aitken (1926) to the EM and ECM algorithms to speed up their convergence. The Aitken δ 2 acceleration does not affect the simplicity, stability and flexibility of the EM algorithm. The Aitken δ2 accelerated method is then applied to log-linear models with partially classified data. In Section 2, we introduce log-linear models and the EM and ECM algorithms for partially classified data. In Section 3, we provide the Aitken δ² acceleration for the EM and ECM algorithms and present theoretical results concerning their convergence. Numerical experiments in Section 4 illustrate the performance of Aitken δ² acceleration for the EM and ECM algorithms. M. Kuroda et al. / Statistics and Probability Letters 78 (2008) 2332–2338
2.  Log-linear models and the EM algorithm
Let X V  = (X 1 , . . , X k ) be a k-dimensional discrete random vector indexed by V = {1, . . ,k}. We also denote a finite set of values of X i by Ω , for i ∈ V , and the space of possible values of X V as the Cartesian product Ωv = Π j ∈ v Ω j . For a subset A ⊆ V , we write X A for {X j |i ∈ A} and Ω A = Π j ∈ A Ω j Let p(xv) denote the cell probability that Xv = xv and let θ = {p(xv)|xv ∈ Ωv} be the set of cell probabilities. The marginal probability that X A = x A for A ⊂ V can be calculated by where the symbol ''\ '' denotes the set difference. In this paper, we consider hierarchical log-linear models for contingency tables with a multinomial distribution. A hierarchical log-linear model is represented by a generating class E = {e 1 , . . , e 1 , in maximal interaction terms, and e i ∈ E is called a generator. For example, consider a log-linear model with the generating class E = {{1, 2}, {1, 3}, {2, 3}}. Then the log-linear model is log p(xv) = u + u 1 (x(1)) + u 2 (x(2)) + u 3 (x(3)) + u 12 (x(1,2)) + u 12 (x(1,3)) + u 23 (x(2,3)).
For the situation of missing data, the completely classified data are represented in the full contingency tables and the partially classified data are represented in supplemental tables. These tables are indexed with an observing pattern T  = {V, t 1 , . . , t 5 }, where each configuration t is a set of observed variables for a supplemental table. The cell counts of the contingency table classified by the observed variables X t are denoted by n t = {n t (x i )| x t ∈ Ω t }. Assume that the missing data are missing at random in the sense of Rubin (1976). Below we show the EM algorithm for finding the MLE of the parameter vector θ for a log-linear model with a generating class E and an observing pattern T = {V, t 1 , . . , ts }. For log-linear models with closed-form MLEs, the EM algorithm repeats the following two steps until convergence to a desired accuracy is obtained: E-step:  Calculate the expected marginal cell counts ñ{r' = {ñc (xc)(r' xc  ∈ Ωc } for each generator e ∈ E :
(r+1)|x e  E  S2e} for each generator e ∈ E: M-step:  Calculate marginal probabilities {p(x, n ( x t ). Then find θ(r+1) from
where C  = (c 1 , . . ,c M –1) is a set of separators of E which may have the same elements, and v(c) counts the number of separators equal to c . C can be obtained by c i = e σ (j) ∩ (e σ (1) ∪ · ·· ∪ e σ (j − 1)), where the generating class E = {e 1, . . , e M} can be ordered as E = (eσ(1), eσ(2), . . , eσ(M)) in such a way that it has the running intersection property: eσ(j) ∩ (eσ(1) ∪ · ·· ∪ eσ(j−1)) ⊂ eσ(i) for some i < j, see Lauritzen (1995). M. Kuroda et al. / Statistics and Probability Letters 78 (2008) 2332–2338 When log-linear models do not have closed-form MLEs, such as in the case of the model (1), the ECM algorithm of Meng and Rubin (1993) may be applied. In the ECM algorithm, the M-step of the EM algorithm is replaced by several computationally simpler CM-steps: CM-step : Find θ(r+1) with the iterative proportional fitting procedure of Bishop et al. (1975) as follows:
for i = 1, 2, . . , M . At each CM-step, the iteration is implemented for every generator once only, and is not repeated until convergence.
3. The Aitken δ² acceleration
The Aitken δ2 method is a nonlinear method for accelerating the convergence of scalar sequences and it is particularly powerful for linear sequence convergence. The method is a simple and computationally inexpensive procedure. First, we describe the Aitken δ 2 method. Let {φ(r)}r≥0 be a scalar sequence which converges to φ*. For the scalar sequence {φ(r)}r>0, the Aitken δ2 method generates a sequence {φ(r)}r>0 by
(5) For convergence of the scalar sequence {φ (r)} r≥0, Traub (1964) provided the following lemma. Lemma 1. If {φ(r)}r≥0 converges to a stationary point φ* as r → ∞, then {φ(r)}r≥0 generated by Eq. (5) converges to the same stationary point φ* To compare the speed of convergence of the sequence {φ ( γ )} r>0 from the Aitken δ  method with that of the sequence {φ(r)}r≥0, we use the following notion of Brezinski and Zaglia (1991). Definition 1. Let {φ(γ)}r>0 be a scalar sequence obtained by applying an extrapolation method to {φ(γ)}r>0. Assume that lim f →∞ φ (r) = lim f →∞ φ (r) = φ *. If then we say that the sequence {φ(r)}r>0 converges to φ* faster than {φ(r)}r>0 or that the extrapolation method accelerates the convergence of {φ(r)}r>0 Traub (1964) proved that the Aitken δ 2 method accelerates the convergence of {φ(r)}r>0 in the sense that
Next we apply the Aitken δ2 acceleration to the EM and ECM algorithms. For simplicity, we write θ  = {pv(xv)){v  ∈ Ωv} = (θ i,θ i, . . ,θ i), where d is the number of cells in the contingency table and the probability vector sequence, we transform the parameter vector θ into distinct scalar parameters. A multinomial distribution with d cells can be factorized into d − 1 independently conditional binomial distributions whose parameters φ = φ(θ) are defined respectively as
(7) M. Kuroda et al. / Statistics and Probability Letters 78 (2008) 2332–2338 (see the Appendix). Since φ  = φ(θ) is a one–one function of θ, the stationary point of φ  = φ(θ) is φ(θ*) for the stationary point θ * of θ . In order to accelerate the convergence of the sequence {θ( γ } r≥0 obtained by the EM or ECM algorithms, we apply the Aitken δ² acceleration to generate the sequence { ∂( γ } r>0. The Aitken δ² acceleration for the EM and ECM algorithms is presented as follows. Let θ(0) = θ(0) denote the initial value. E-step:   Using θ(r) = {p(xv)(r)|xv ∈ Ωv} and the observed frequencies, calculate the expected marginal counts ñ(r) for each generator e ∈ E by Eq. (2). M (or CM)-step:  Find θ(r+1) = {p(xv)(r+1)|xv ∈ Ωv} by using Eq. (3) or (4). Aitken δ²acceleration:  Calculate
by Eq. (7) where (θ(r-1),θ(r),θ(r+1)) is obtained at t the previous M (or CM)-steps. Generate a vector
and check the convergence by where δ is the desired accuracy. Note that to find the expected frequencies at each E-step, we use θ(r) = {p(xy)(7)|xy ∈ Ωy} obtained at the previous M-step but not the ∂( ) obtained at t the previous Aitken δ² acceleration. Below we give the properties of the convergence of the Aitken δ² acceleration for the EM and ECM algorithms. Let θ* denote a stationary point of the sequence {θ(r)}r≥0 obtained by the EM or ECM algorithm. Theorem 1. For a given initial value θ(0), the sequence {θ(r)}r≥0 obtained by using the Aitken δ² acceleration converges to the same stationary point θ*. Proof. Since φ = φ(θ) of Eq. (7) is a one–one monotone function of θ, we need only prove convergence of the sequence {φ(r)}r>0 to the stationary point φ* of the sequence {φ(r)}r>0. If {φ(r)}r≥0 is a convergent vector sequence to the vector φ*, then ∥φ(r) − φ*‖∞ converges to zero as r → ∞, moreover |φ,{ } — φ, | converges to zero as r → ∞∞, where || · ||∞ = max ; {| · |}. Thus, from the above facts and Lemma 1, the sequence {φ(r)},r>0 converges to the stationary point φ* of the sequence {φ(r)},r>0 as r → ∞. In the following theorem, we show the speed of convergence of the Aitken δ 2 accelerated EM and ECM algorithms. Theorem 2. The Aitken 8² acceleration speeds up the convergence of the EM and ECM algorithms, that is, for all i, Proof. From Eq. (6), we have M. Kuroda et al. / Statistics and Probability Letters 78 (2008) 2332–2338
for all i . Since ὁ ( r ) of Eq. (8) is a one–one monotone function of φ ( r ), it follows that for all i . Thus the sequence {θ(r)}r>0 converges to θ* faster than {θ(r)}r>0 does. ■
4.  Numerical experiments
In this section, we illustrate how much faster the Aitken δ² accelerated EM algorithm converges compared to the EM algorithm using numerical experiments. Note that the convergence speed of the Aitken δ² acceleration does not depend on the structure of the log-linear models but on the EM and ECM sequences. Example 1. Consider a 2 × 2 × 2 contingency table concerning the survival of infants in Table 1. The completely classified data in Table 1(a) was previously analyzed in Bishop et al. (1975), and Table 1(b) and (c) give the artificial data. For the data of Table 1(b), Clinic (X 1 ) is missing; for the data of Table 1(c), Prenatal Care (X 2 ) is missing. The observed pattern is T  = {{1, 2, 3}, {1, 2}, {2, 3}}. Suppose that the data have a multinomial distribution with unknown parameter θ . Table 2 summarizes the MLEs for several log-linear models and the numbers of iterations for the EM and ECM algorithms and the Aitken δ² acceleration with the desired accuracy δ = 10−9. We applied the ECM algorithm to the no third-order-interaction term model with E = {{1, 2}, {1, 3}, {2, 3}}. As shown in Table 2, the Aitken δ² acceleration converges to the MLEs faster than the EM algorithm, and the convergence speed is about twice as fast as that of the EM algorithm for all of these log-linear models. Example 2. Consider a 2 × 2 contingency table with completely and partially classified observations. Let X 1 and X 2 be dichotomous variables. Suppose that the observed frequencies have a multinomial distribution with a parameter vector θ . The observed frequencies classified by X 1 and X 2 are n 12 = (n 12(1, 1), n 12(2, 1), n 12(1, 2), n 12(2, 2)) = (5, 2, 4, 1), and the frequencies partially classified by X 1 are n 1 = (n 1 (1), n 1 (2)) = (75, 25). Several data sets of the observed frequencies n 2 partially classified by X 2 are given in the first two columns of Table 3. Set the desired accuracy δ to be 10−9. From the third column of Table 3, we can see that for these data sets, the convergence of the EM algorithm is quite slow, and its convergence depends largely on the proportion of incompletely observed frequencies. From the fourth column of Table 3, we can see that Aitken δ² acceleration speeds up the convergence of the EM algorithm by a factor between 2 and 4.
Table 3

Acknowledgements
The authors are grateful to an anonymous referee and the executive editor Dr. Richard A. Johnson for valuable comments. This research was supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Young Scientists, No. 187002796. M. Kuroda et al. / Statistics and Probability Letters 78 (2008) 2332–2338 Appendix , Y d ) be a random vector from a multinomial distribution with a parameter vector θ Let – ( ∑ d –1 θ i ). That is, the distribution of Y is given by (θ 1 , , Y i . Using a chain of binomial random variables, the multinomial distribution can be factorized into where n a product of conditional binomial distributions of Y i 's as follows:
Thus we can obtain the parameterization ϕ = ϕ(θ) given by