The EM Algorithm and Extensions

G.J.McLACHLAN and T.KRISHNAN

ii CONTENTS

To

Beryl, Jonathan, and Robbie

C o n t e n t s

Preface iv

1 GENERAL INTRODUCTION 1 1.1 Introduction 1 1.2 ML Estimation 3 1.3 Newton-Type Methods 5

1.3.1 Introduction 5 1.3.2 Newton-Raphson Method 6 1.3.3 Quasi-Newton Methods 7 1.3.4 Modified Newton Methods 8

1.4 Introductory Examples 9 1.4.1 Introduction 9 1.4.2 Example 1.1: A Multinomial Example 9 1.4.3 Example 1.2: Estimation of Mixing Proportions 16

1.5 Formulation of the EM Algorithm 21 1.5.1 EM Algorithm 21 1.5.2 Example 1.3: Censored Exponentially Distributed Sur-

vival Times 23 1.5.3 E- and M-Steps for the Regular Exponential Family 26 1.5.4 Example 1.4: Censored Exponentially Distributed Sur-

vival Times (Example 1.3 Continued) 27 1.5.5 Generalized EM Algorithm 28 1.5.6 GEM Algorithm Based on One Newton-Raphson Step 28 1.5.7 EM Gradient Algorithm 30 1.5.8 EM Mapping 30

1.6 EM Algorithm for MAP and MPL Estimation 31 1.6.1 Maximum a Posteriori Estimation 31 1.6.2 Example 1.5: Multinomial Example (Example 1.1 Con-

tinued) 31 1.6.3 Maximum Penalized Estimation 32

1.7 Brief Summary of the Properties of the EM Algorithm 33 1.8 History of the EM Algorithm 34

1.8.1 Work Before Dempster, Laird, and Rubin (1977) 34 1.8.2 EM Examples and Applications Since Dempster, Laird,

and Rubin (1977) 37 1.8.3 Two Interpretations of EM 38 1.8.4 Developments in EM Theory, Methodology, and Appli-

cations 39

1.9 Overview of the Book 41 1.10 Notation 42

2 EXAMPLES OF THE EM ALGORITHM 45 2.1 Introduction 45 2.2 Multivariate Data with Missing Values 45

2.2.1 Example 2.1: Bivariate Normal Data with Missing Values 45 2.2.2 Numerical Illustration 49 2.2.3 Multivariate Data: Buck’s Method 50

2.3 Least Squares with Missing Data 51 2.3.1 Healy–Westmacott Procedure 51 2.3.2 Example 2.2: Linear Regression with Missing Depen-

dent Values 52 2.3.3 Example 2.3: Missing Values in a Latin Square Design 53 2.3.4 Healy–Westmacott Procedure as an EM Algorithm 55

2.4 Example 2.4: Multinomial with Complex Cell Structure 56 2.5 Example 2.5: Analysis of PET and SPECT Data 59 2.6 Example 2.6: Multivariate t-Distribution (Known D.F.) 63

2.6.1 ML Estimation of Multivariate t-Distribution 63 2.6.2 Numerical Example: Stack Loss Data 66

2.7 Finite Normal Mixtures 67 2.7.1 Example 2.7: Univariate Component Densities 67 2.7.2 Example 2.8: Multivariate Component Densities 70 2.7.3 Numerical Example: Red Blood Cell Volume Data 71

2.8 Example 2.9: Grouped and Truncated Data 72 2.8.1 Introduction 72 2.8.2 Specification of Complete Data 73 2.8.3 E-Step 76 2.8.4 M-Step 77 2.8.5 Confirmation of Incomplete-Data Score Statistic 77 2.8.6 M-Step for Grouped Normal Data 78 2.8.7 Numerical Example: Grouped Log Normal Data 79

3 BASIC THEORY OF THE EM ALGORITHM 82 3.1 Introduction 82 3.2 Monotonicity of the EM Algorithm 83 3.3 Monotonicity of a Generalized EM Algorithm 84 3.4 Convergence of an EM Sequence to a Stationary Value 85

3.4.1 Introduction 85 3.4.2 Regularity Conditions of Wu (1983) 86 3.4.3 Main Convergence Theorem for a Generalized EM Se-

quence 87 3.4.4 A Convergence Theorem for an EM Sequence 88

3.5 Convergence of an EM Sequence of Iterates 88 3.5.1 Introduction 88 3.5.2 Two Convergence Theorems of Wu (1983) 89 3.5.3 Convergence of an EM Sequence to a Unique Maximum

Likelihood Estimate 90

3.6 Examples of Nontypical Behavior of an EM (GEM) Sequence 91 3.6.1 Example 3.1: Convergence to a Saddle Point 91 3.6.2 Example 3.2: Convergence to a Local Minimum 94 3.6.3 Example 3.3: Nonconvergence of a Generalized EM Se-

quence 97

3.7 Score Statistic 99 3.8 Missing Information 100

3.8.1 Missing Information Principle 100 3.8.2 Example 3.4. (Example 1.3 Continued) 102

3.9 Rate of Convergence of the EM Algorithm 105 3.9.1 Rate Matrix for Linear Convergence 105 3.9.2 Measuring the Linear Rate of Converegence 106 3.9.3 Rate Matrix in Terms of Information Matrices 106 3.9.4 Rate Matrix for Maximum a Posteriori Estimation 107 3.9.5 Derivation of Rate Matrix in terms of Information Ma-

trices 108 3.9.6 Example 3.5. (Example 3.3 Continued) 109 4 STANDARD ERRORS AND SPEEDING UP CONVER-

GENCE 110 4.1 Introduction 110 4.2 Observed Information Matrix 111

4.2.1 Direct Evaluation 111 4.2.2 Extraction of Observed Information Matrix in Terms

of the Complete-Data Log Likelihood 111 4.2.3 Regular Case 113 4.2.4 Evaluation of the Conditional Expected Complete-Data

Information Matrix 114 4.2.5 Examples 115

4.3 Approximations to Observed Information Matrix: i.i.d. Case 120 4.4 Observed Information Matrix for Grouped Data 122

4.4.1 Approximation Based on Empirical Information 122 4.4.2 Example 4.3: Grouped Data from an Exponential Dis-

tribution 124

4.5 Supplemented EM Algorithm 128 4.5.1 Definition 128 4.5.2 Calculation of J(Ψ̂) via Numerical Differentiation 130 4.5.3 Stability 131 4.5.4 Monitoring Convergence 131 4.5.5 Example 4.4: Univariate Contaminated Normal Data 132 4.5.6 Example 4.5: Bivariate Normal Data with Missing Values136

4.6 Bootstrap Approach to Standard Error Approximation 139 4.7 Acceleration of the EM Algorithm via Aitken’s Method 140

4.7.1 Aitken’s Acceleration Method 140 4.7.2 Louis’ Method 141 4.7.3 Example 4.6: Multinomial Data 142 4.7.4 Example 4.7: Geometric Mixture 143 4.7.5 Grouped and Truncated Data. (Example 2.8 continued) 147

4.8 An Aitken Acceleration-Based Stopping Criterion 148 4.9 Conjugate Gradient Acceleration of EM Algorithm 149

4.9.1 Conjugate Gradient Method 149 4.9.2 A Generalized Conjugate Gradient Algorithm 150 4.9.3 Accelerating the EM Algorithm 150

4.10 Hybrid Methods for Finding the MLE 151 4.10.1Introduction 151 4.10.2Combined EM and Modified Newton-Raphson Algorithm152

4.11 A GEM Algorithm Based on One Newton-Raphson Step 154 4.11.1Derivation of a Condition to be a Generalized EM Se-

quence 154 4.11.2Simulation Experiment 155

4.12 EM Gradient Algorithm 155 4.13 A Quasi-Newton Acceleration of the EM Algorithm 158

4.13.1The Method 158 4.13.2Example 4.8: Dirichlet Distribution 161

5 EXTENSIONS OF THE EM ALGORITHM 165 5.1 Introduction 165 5.2 ECM Algorithm 166

5.2.1 Motivation 166 5.2.2 Formal Definition 166 5.2.3 Convergence Properties 168 5.2.4 Speed of Convergence 169 5.2.5 Discussion 170

5.3 Multicycle ECM Algorithm 170 5.4 Example 5.1: Normal Mixtures with Equal Correlations 172

5.4.1 Normal Components with Equal Correlations 172 5.4.2 Application of ECM Algorithm 172 5.4.3 Fisher’s Iris Data 174

5.5 Example 5.2: Mixture Models for Survival Data 174 5.5.1 Competing Risks in Survival Analysis 174 5.5.2 A Two-Component Mixture Regression Model 175 5.5.3 Observed Data 176 5.5.4 Application of EM Algorithm 177 5.5.5 M-Step for Gompertz Components 177 5.5.6 Application of a Multicycle ECM Algorithm 179

5.6 Example 5.3: Contingency Tables with Incomplete Data 180 5.7 ECME Algorithm 182 5.8 MLE of t-Distribution with Unknown D.F. 182

5.8.1 Application of the EM Algorithm 182 5.8.2 M-Step 184 5.8.3 Application of ECM Algorithm 184 5.8.4 Application of ECME Algorithm 185 5.8.5 Some Standard Results 186 5.8.6 Missing Data 186 5.8.7 Numerical Examples 188 5.8.8 Theoretical Results on the Rate of Convergence 189 5.9 Example 5.5: Variance Components 190

5.9.1 A Variance Components Model 190 5.9.2 E-Step 191 5.9.3 M-Step 192 5.9.4 Application of Two Versions of ECME Algorithm 193 5.9.5 Numerical Example 194

5.10 Example 5.6: Factor Analysis 195 5.10.1ECM Algorithm for Factor Analysis 195 5.10.2Numerical Example 198

5.11 Efficient Data Augmentation 199 5.11.1Motivation 199 5.11.2Maximum Likelihood Estimation of t-Distribution 199 5.11.3Variance Components Model 204

5.12 Alternating ECM Algorithm 205 5.13 EMS Algorithm 207 5.14 One-Step-Late Algorithm 207 5.15 Variance Estimation for Penalized EM and OSL Algorithms 208

5.15.1Penalized EM Algorithm 208 5.15.2OSL Algorithm 209 5.15.3Example 5.7. (Example 4.1 continued) 210

5.16 Linear Inverse Problems 211

6 MISCELLANEOUS TOPICS 213 6.1 Iterative Simulation Algorithms 213 6.2 Monte Carlo E-Step 213 6.3 Stochastic EM Algorithm 215 6.4 Data Augmentation Algorithm 218 6.5 Multiple Imputation 219 6.6 Sampling–Importance Resampling 220 6.7 Gibbs Sampler 222

6.7.1 The Algorithm 222 6.7.2 Example 6.1: Gibbs Sampler for the Mixture Problem 223 6.7.3 Gibbs Sampler Analogs of ECM and ECME Algorithms 224 6.7.4 Some Other Analogs of the Gibbs Sampler 226

6.8 Missing-Data Mechanism and Ignorability 226 6.9 Competing Methods and Some Comparisons with EM 227

6.9.1 Introduction 227 6.9.2 Simulated Annealing 228 6.9.3 Comparison of SA and EM Algorithm for Normal Mix-

tures 229

6.10 The Delta Algorithm 230 6.11 Image Space Reconstruction Algorithm 231 6.12 Further Applications of the EM Algorithm 232

6.12.1Introduction 232 6.12.2Hidden Markov Models 233 6.12.3AIDS Epidemiology 236 6.12.4Neural Networks 238

REFERENCES 242

P r e f a c e

This book deals with the Expectation–Maximization algorithm, popularly known as the EM algorithm. This is a general-purpose algorithm for maximum likelihood estimation in a wide variety of situations best described as incomplete-data problems. The name EM algorithm was given by Dempster, Laird, and Rubin in a celebrated paper read before the Royal Statistical Society in 1976 and published in its journal in 1977. In this paper, a general formulation of the EM algorithm was presented, its basic properties established, and many examples and applications of it provided. The idea behind the EM algorithm is intuitive and natural and so algorithms like it were formulated and applied in a variety of problems even before this paper. However, it was in this seminal paper that the ideas in the earlier papers were synthesized, a general formulation and a theory developed, and a host of traditional and non-traditional applications indicated. Since then, the EM algorithm has become a standard piece in the statistician’s repertoire. The incomplete-data situations where the EM algorithm has been successfully applied include not only evidently incomplete-data situations, where there are missing data, truncated distributions, censored or grouped observations, but also a whole variety of situations where the incompleteness of the data is not natural or evident. Thus, in some situations, it requires a certain amount of ingenuity on the part of the statistician to formulate the incompleteness in a suitable manner to facilitate the application of the EM algorithm in a computationally profitable manner. Following the paper of Dempster, Laird, and Rubin (1977), a spate of applications of the algorithm have appeared in the literature.

The EM algorithm is not without its limitations, many of which came to light in attempting to apply it in certain complex incomplete-data problems and some even in innocuously simple incomplete-data problems. However, a number of modifications and extensions of the algorithm has been developed to overcome some of these limitations. Thus there is a whole battery of EM-related algorithms and more are still being developed. The current developments are, however, in the direction of iterative simulation techniques or Markov Chain Monte Carlo methods, many of which can be looked upon as simulation-based versions of various EM-type algorithms.

Incomplete-data problems arise in all statistical contexts. Hence in these problems where maximum likelihood estimates usually have to be computed iteratively, there is the scope and need for an EM algorithm to tackle them. Further, even if there are no missing data or other forms of data incompleteness, it is often profitable to express the given problem as an incomplete-data one within an EM framework. For example, in some multiparameter problems like in random effects models, where an averaging over some parameters is to be carried out, an incomplete-data approach via the the EM algorithm and its variants has been found useful. No wonder then that the EM algorithm has become an ubiquitous statistical tool, is a part of the entire spectrum of statistical methods, and has found applications in almost all fields where statistical techniques have been applied. The EM algorithm and its variants have been applied in such fields as medical imaging, dairy science, correcting census undercount, and AIDS epidemiology, to mention a few. Articles containing applications of the EM algorithm and even some with some methodological content have appeared in a variety of journals on statistical theory and methodology, statistical computing, and statistical applications in engineering, biology, medicine, social sciences, etc. Meng and Pedlow (1992) list a bibliography of over 1000 items and now there are at least 1700 publications related to the EM algorithm.

It is surprising that despite the obvious importance of the technique and its ramifications, no book on the subject has so far appeared. Indeed, many modern books dealing with some aspect of statistical estimation have at least some EM algorithm content. The books by Little and Rubin (1987), Tanner (1991, 1993), and Schafer (1996) have substantial EM algorithm content. But still, there seems to be a need for a full-fledged book on the subject. In our experience of lecturing to audiences of professional statisticians and to users of statistics, it appears that there is a definite need for a unified and complete treatment of the theory and methodology of the EM algorithm and its extensions, and their applications. The purpose of our writing this book is to fulfill this need. The various extensions of the EM algorithm due to Rubin, Meng, C. Liu, and others that have appeared in the last few years, have made this need even greater. Many extensions of the EM algorithm in the direction of iterative simulation have also appeared in recent years. Inclusion of these techniques in this book may have resulted in a more evenhanded and comprehensive treatment of the EM algorithm and its extensions. However, we decided against it, since this methodology is still evolving and rapid developments in this area may make this material soon obsolete. So we have restricted this book to the EM algorithm and its variants and have only just touched upon the iterative simulation versions of it.

The book is aimed at theoreticians and practitioners of Statistics and its objective is to introduce to them the principles and methodology of the EM algorithm and its tremendous potential for applications. The main parts of the book describing the formulation of the EM algorithm, detailing its methodology, discussing aspects of its implementation, and illustrating its application in many simple statistical contexts, should be comprehensible to graduates with Statistics as their major subject. Throughout the book, the theory and methodology are illustrated with a number of examples. Where relevant, analytical examples are followed up with numerical examples. There are about thirty examples in the book. Some parts of the book, especially examples like factor analysis and variance components analysis, will need basic knowledge of these techniques to comprehend the full impact of the use of the EM algorithm. But our treatment of these examples is self-contained, although brief. However, these examples can be skipped without losing continuity.

Chapter 1 begins with a brief discussion of maximum likelihood (ML) estimation and standard approaches to the calculation of the maximum likelihood estimate (MLE) when it does not exist as a closed form solution of the likelihood equation. This is followed by a few examples of incomplete-data problems for which an heuristic derivation of the EM algorithm is given. The EM algorithm is then formulated and its basic terminology and notation established. The case of the regular exponential family (for the complete-data problem) for which the EM algorithm results in a particularly elegant solution, is specially treated. Throughout the treatment, the Bayesian perspective is also included by showing how the EM algorithm and its variants can be adapted to compute maximum a posteriori (MAP) estimates. The use of the EM algorithm and its variants in maximum penalized likelihood estimation (MPLE), a technique by which the MLE is smoothed, is also included.

Chapter 1 also includes a summary of the properties of the EM algorithm. Towards the end of this chapter, a comprehensive discussion of the history of the algorithm is presented, with a listing of the earlier ideas and examples upon which the general formulation is based. The chapter closes with a summary of the developments in the methodology since the Dempster et al. (1977) paper and with an indication of the range of applications of the algorithm.

In Chapter 2, a variety of examples of the EM algorithm is presented, following the general formulation in Chapter 1. These examples include missing values (in the conventional sense) in various experimental designs, the multinomial distribution with complex cell structure as used in genetics, the multivariate t-distribution for the provision of a robust estimate of a location parameter, Poisson regression models in a computerized image reconstruction process such as SPECT/PET, and the fitting of normal mixture models to grouped and truncated data as in the modeling of the volume of red blood cells.

In Chapter 3, the basic theory of the EM algorithm is systematically presented, and the monotonicity of the algorithm, convergence, and rates of convergence properties are established. The Generalized EM (GEM) algorithm and its properties are also presented. The principles of Missing Information and Self-Consistency are discussed. In this chapter, attention is inevitably given to mathematical details. However, mathematical details and theoretical points are explained and illustrated with the help of earlier and new examples. Readers not interested in the more esoteric aspects of the EM algorithm may only study the examples in this chapter or omit the chapter altogether without losing continuity.

In Chapter 4, two issues which have led to some criticism of the EM algorithm are addressed. The first concerns the provision of standard errors, or the full covariance matrix in multivariate situations, of the MLE obtained via the EM algorithm. One initial criticism of the EM algorithm was that it does not automatically provide an estimate of the covariance matrix of the MLE, as do some other approaches such as Newton-type methods. Hence we consider a number of methods for assessing the covariance matrix of the MLE Ψ̂ of the parameter vector Ψ, obtained via the EM algorithm. Most of these methods are based on the observed information matrix. A coverage is given of methods such as the Supplemented EM algorithm that allow the observed information matrix to be calculated within the EM framework. The other common criticism that has been leveled at the EM algorithm is that its convergence can be quite slow. We therefore consider some methods that have been proposed for accelerating the convergence of the EM algorithm. They include methods based on Aitken’s acceleration procedure and the generalized conjugate gradient approach, and hybrid methods that switch from the EM algorithm after a few iterations to some Newton-type method. We consider also the use of the EM gradient algorithm as a basis of a quasi-Newton approach to accelerate convergence of the EM algorithm. This algorithm approximates the M-step of the EM algorithm by one Newton-Raphson step when the solution does not exist in closed form.

In Chapter 5, further modifications and extensions to the EM algorithm are discussed. The focus is on the Expectation–Conditional Maximization (ECM) algorithm and its extensions, including the the Expectation–Conditional Maximization Either (ECME) and Alternating ECM (AECM) algorithms. The ECM algorithm is a natural extension of the EM algorithm in situations where the maximization process on the M-step is relatively simple when conditional on some function of the parameters under estimation. The ECM algorithm therefore replaces the M-step of the EM algorithm by a number of computationally simpler conditional maximization (CM) steps. These extensions of the EM algorithm typically produce an appreciable reduction in total computer time. More importantly, they preserve the appealing convergence properties of the EM algorithm, such as its monotone convergence.

In Chapter 6, very brief overviews are presented of iterative simulation techniques such as the Monte Carlo E-step, Stochastic EM algorithm, Data Augmentation, and the Gibbs sampler and their connections with the various versions of the EM algorithm. Then, a few methods such as Simulated Annealing, which are considered competitors to the EM algorithm, are described and a few examples comparing the performance of the EM algorithm with these competing methods are presented. The book is concluded with a brief account of the applications of the EM algorithm in such topical and interesting areas as Hidden Markov Models, AIDS epidemiology, and Neural Networks.

One of the authors (GJM) would like to acknowledge gratefully financial support from the Australian Research Council. Work on the book by one of us (TK) was facilitated by two visits to the University of Queensland, Brisbane, and a visit to Curtin University of Technology, Perth. One of the visits to Brisbane was under the Ethel Raybould Fellowship scheme. Thanks are due to these two Universities for their hospitality. Thanks are also due to the authors and owners of copyrighted material for permission to reproduce tables and figures. The authors also wish to thank Ramen Kar and Amiya Das of the Indian Statistical Institute, Calcutta, and Pauline Wilson of the University of Queensland, Brisbane, for their help with LaTEX and word processing, and Rudy Blazek of Michigan State University and Matthew York and David Peel of the University of Queensland for assistance with the preparation of some of the figures.

G.J. McLachlan T. Krishnan Brisbane Calcutta

November 1995