Unformatted text preview: 9 c A 9 c e 6 V 6 7 9 c e 6 V 7 6 6 9 e A 9 7 a GQ y s w C u t s E pQ E h V e e c E a Y R V T I R Q EI G E C A 9 7 @[email protected]@b8if8f@@[email protected] bbx a bFv8E rqHFigfbdPI [email protected] A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models
Jeff A. Bilmes ([email protected]) International Computer Science Institute Berkeley CA, 94704 and Computer Science Division Department of Electrical Engineering and Computer Science U.C. Berkeley TR97021 April 1998 Abstract We describe the maximumlikelihood parameter estimation problem and how the ExpectationMaximization (EM) algorithm can be used for its solution. We first describe the abstract form of the EM algorithm as it is often given in the literature. We then develop the EM parameter estimation procedure for two applications: 1) finding the parameters of a mixture of Gaussian densities, and 2) finding the parameters of a hidden Markov model (HMM) (i.e., the BaumWelch algorithm) for both discrete and Gaussian mixture observation models. We derive the update equations in fairly explicit detail but we do not prove any convergence properties. We try to emphasize intuition rather than mathematical rigor. % ) '543 32(10) % # ('&$" ! ii 1 Maximumlikelihood
Recall the definition of the maximumlikelihood estimation problem. We have a density function that is governed by the set of parameters (e.g., might be a set of Gaussians and could be the means and covariances). We also have a data set of size , supposedly drawn from this distribution, i.e., . That is, we assume that these data vectors are independent and identically distributed (i.i.d.) with distribution . Therefore, the resulting density for the samples is This function is called the likelihood of the parameters given the data, or just the likelihood function. The likelihood is thought of as a function of the parameters where the data is fixed. In the maximum likelihood problem, our goal is to find the that maximizes . That is, we wish to find where argmax Often we maximize instead because it is analytically easier. Depending on the form of this problem can be easy or hard. For example, if is simply a single Gaussian distribution where , then we can set the derivative of to zero, and solve directly for and (this, in fact, results in the standard formulas for the mean and variance of a data set). For many problems, however, it is not possible to find such analytical expressions, and we must resort to more elaborate techniques. 2 Basic EM
The EM algorithm is one such elaborate technique. The EM algorithm [ALR77, RW84, GJ95, JJ94, Bis95, Wu83] is a general method of finding the maximumlikelihood estimate of the parameters of an underlying distribution from a given data set when the data is incomplete or has missing values. There are two main applications of the EM algorithm. The first occurs when the data indeed has missing values, due to problems with or limitations of the observation process. The second occurs when optimizing the likelihood function is analytically intractable but when the likelihood function can be simplified by assuming the existence of and values for additional but missing (or hidden) parameters. The later application is more common in the computational pattern recognition community. As before, we assume that data is observed and is generated by some distribution. We call the incomplete data. We assume that a complete data set exists and also assume (or specify) a joint density function: . Where does this joint density come from? Often it "arises" from the marginal density function and the assumption of hidden variables and parameter value guesses (e.g., our two examples, Mixturedensities and BaumWelch). In other cases (e.g., missing data values in samples of a distribution), we must assume a joint relationship between the missing and observed values. 1 3 T &F SR X & 21baW `2 Y"W V2 X U ) 0( 7 65 4 ( 21 ' & 3 I H G P"2 7 @9 3 I H G # 8 DB F 65 3 ECA % # $"!!! 65 3 DB F @9 3 QCA 8 Where are the current parameters estimates that we used to evaluate the expectation and are the new parameters that we optimize to increase . This expression probably requires some explanation. The key thing to understand is that and are constants, is a normal variable that we wish to adjust, and is a random . The right side of Equation 1 can therefore be variable governed by the distribution rewritten as: Note that is the marginal distribution of the unobserved data and is dependent on both the observed data and on the current parameters, and is the space of values can take on. In the best of cases, this marginal distribution is a simple analytical expression of the assumed parameters and perhaps the data. In the worst of cases, this density might be very hard to obtain. Sometimes, in fact, the density actually used is but this doesn't effect subsequent steps since the extra factor, is not dependent on . As an analogy, suppose we have a function of two variables. Consider where is a constant and is a random variable governed by some distribution . Then is now a deterministic function that could be maximized if desired. The evaluation of this expectation is called the Estep of the algorithm. Notice the meaning of the two arguments in the function . The first argument corresponds to the parameters that ultimately will be optimized in an attempt to maximize the likelihood. The second argument corresponds to the parameters that we use to evaluate the expectation. The second step (the Mstep) of the EM algorithm is to maximize the expectation we computed in the first step. That is, we find: argmax These two steps are repeated as necessary. Each iteration is guaranteed to increase the loglikelihood and the algorithm is guaranteed to converge to a local maximum of the likelihood function. There are many rateofconvergence papers (e.g., [ALR77, RW84, Wu83, JX96, XJ96]) but we will not discuss them here.
different density functions since argument usage should should disambiguate different ones. v d f v` x wd v` t h r p h f d b` W 6"p ace y!'caY uI s'qig 9ecaY XV Recall that . In the following discussion, we drop the subscripts from 2 T DB F ECA With this new density function, we can define a new likelihood function, , called the completedata likelihood. Note that this function is in fact a random variable since the missing information is unknown, random, and presumably governed by an underlying distribution. That is, we can think of for some function where and are constant and is a random variable. The original likelihood is referred to as the incompletedata likelihood function. The EM algorithm first finds the expected value of the completedata loglikelihood with respect to the unknown data given the observed data and the current parameter estimates. That is, we define: (1) &F @9 4 R 65 3 T 3 (2) I X " Q# " P H 9 ! FD 2 B @ X 7 G 8 7 E @ & C 7 A# B @ 8 7 9 7! 8 6 # ( # ( # ( a `2 5 ( 2 # X X ( 4 X ( a Y2 # X D X 3 ( a Y # Y" 1 EB A 1)'% ! ( a & 6F 1 EB A 2 X X D 0 ( &$ T T @9 3 DB T QCA ( a 5 ! ( a P F a ( a 5 T &2 &F @9 3 T ( a `2 # X R 0Sa 5 ( T T T "( ( T 6F U R T A modified form of the Mstep is to, instead of maximizing , we find some such that . This form of the algorithm is called Generalized EM (GEM) and is also guaranteed to converge. As presented above, it's not clear how exactly to "code up" the algorithm. This is the way, however, that the algorithm is presented in its most general form. The details of the steps required to compute the given quantities are very dependent on the particular application so they are not discussed when the algorithm is presented in this abstract form. 3 Finding Maximum Likelihood Mixture Densities Parameters via EM
The mixturedensity parameter estimation problem is probably one of the most widely used applications of the EM algorithm in the computational pattern recognition community. In this case, we assume the following probabilistic model: where the parameters are such that and each is a density function parameterized by . In other words, we assume we have component densities mixed together with mixing coefficients . The incompletedata loglikelihood expression for this density from the data is given by: which is difficult to optimize because it contains the log of the sum. If we consider as incomplete, however, and posit the existence of unobserved data items whose values inform us which component density "generated" each data item, the likelihood expression is significantly simplified. That is, we assume that for each , and if the sample was mixture component. If we know the values of , the likelihood becomes: generated by the which, given a particular form of the component densities, can be optimized using a variety of techniques. The problem, of course, is that we do not know the values of . If we assume is a random vector, however, we can proceed. We first must derive an expression for the distribution of the unobserved data. Let's first guess at parameters for the mixture density, i.e., we guess that are the . Given , we can easily compute appropriate parameters for the likelihood for each and . In addition, the mixing parameters, can be though of as prior probabilities of each mixture component, that is component j . Therefore, using Bayes's rule, we can compute: 3 ( 4 7 ( ( F 7 ( 2 30 4 4 7 !!! 4 7 4 !!! 4 5 4 &$ '%! ) ) 2 20 30 31 EB A 0( 0( DB D B ) " " ( B ( EB A " F @ EB 4 F&F 65 3 QCA ) D T ) D A T T ( 0( ) ) ) D a DB ) QCA 0( ( 0( EB A F 65 3 QCA DB 7 ( ' # # ( a 5 " # ( B ) 30 4 ( 30 30 4 a ( ( " 84 7 3( 0 8 84 230 30 8 2 2 2 2 2 B 477 ( 4 477 ( 4 ) 0#( % ( B T # T T ! ( 7 2 ( ( ( (7 7 !!! 7 !!! T &F 4 9 3 ) 0( & 2 !!! ( B # ( a 5 "( a ( 9 & ($ " 6 ! ) 0( ) D B 4 a ( 1 bF 7 ( DEB A 4 a ( 0 QCA # D B 4 a ( F 7 ( QCA 0 4 a ! 1 0( ) ) 4 a ( 1 4 a ( 4 a B 4 a ( 1 4 a " B ( a ) ) ) # ' 0 !! )U 2 30 )U 2 30 ( a ) ) ' #0 )U
0 0 !! 0 4 a B ) 2 a 30 ) ) )U ' !! # !!! 4 a 9 0 0 0 0 ) )U 4 a " B ) 2 )a 30 ) ) )U D EB A 0( 4 a B" ' !! " 7 ( # 0 # 0 ) )a ' " 7 ( # 4 a " B ) )a ' # 4 0 0 )U ) ( C( )& # DB ( QCA )& ) ) 2 D B 30 C( ) !! ECA # 0 ) 230 2 2 30 30 D B C( ) " 7 ( QCA !! # ) 2 a 30 2 2 30 30 D B a " ' " 7 ( QCA B # 4 a Y21 "" @9 3 X X 4 a 5 where is an instance of the unobserved data independently drawn. When we now look at Equation 2, we see that in this case we have obtained the desired marginal density by assuming the existence of the hidden variables and making a guess at the initial parameters of their distribution. In this case, Equation 1 takes the form: B &# 9 !!! B X 4 a ( ( " B ) 0( X ' 4 a Y21 # for since and To maximize this expression, we can maximize the term containing independently since they are not related. In this form, , looks fairly daunting, yet it can be greatly simplified. We first note that . Using Equation 4, we can write Equation 3 as: 4 7 and the term containing (5) (4) (3) ) C( ) # ) C( ) # 4 a 5 To find the expression for , we introduce the Lagrange multiplier , and solve the following equation: For some distributions, it is possible to get an analytical expressions for as functions of everything else. For example, if we assume dimensional Gaussian component distributions with mean and then covariance matrix , i.e., To derive the update equations for this distribution, we need to recall some results from matrix algebra. The trace of a square matrix tr is equal to the sum of 's diagonal elements. The trace of a scalar equals that scalar. Also, tr tr tr , and tr tr which implies that tr where . Also note that indicates the determinant of a . matrix, and that We'll need to take derivatives of a function of a matrix with respect to elements of that matrix. Therefore, we define to be the matrix with entry where is the entry of . The definition also applies taking derivatives with respect to a vector. First, . Second, it can be shown that when is a symmetric matrix: by the definition of the inverse of a matrix. Finally, it can be shown that: Diag 5 7 6 8 B U 8 8 tr 5 d c a `6 b! 6 ! if if diag 7 5 65 ( 65 C ( C Y 5 8 e65 `W Y D B 5 QCA 5 Y where is the cofactor of . Given the above, we see that: 6 ! 6 ! a b if if G ( S 5 8 &A7 G & G QR2 IPE F 0HEE )' 0( &% U ) 421 0( &% U 3 )' 7 8 @75 65 & ($ 6 ! 6 # 5 5 4 a ( 7 8 5 # $I ! I " ! S ( Y`WX ( 5 5 &7 ( ) 0( Y # B( 4 ( ( 8 9 7 8 5 6 5 Summing both sizes over , we get that resulting in: 4 a ( ) 0( # or (6) with the constraint that 4 a ( 1b ) D B 0( ) ECA G & 2 7 G 5 G E 0E F 2 # 65 CD 5 G W 8 @75 & '$ 6 ! BV5U 7 5 5 ( 5 4 ( B( ( E % TG % 1 % E & '$ 6 ! ( . Setting the derivative to zero, i.e., V 4 a ( . This gives ) 0( # where diag ( 4 a ( 0( I ) # W # " ( ( ) 0( V 4 a ( " ( ( $ 4 a ( 1 # " ( B G ( G ( 4 a ( 0( ) # G ) 4 a ( ( 0#( G 4 a ( 0 G ( ) 0( # G 4 a ( 1 G T( B G ( 4 a ( 1 F 0)( ) B a DECA # ) ) D C( G ( EB A # To find with which we can easily solve for Taking the log of Equation 6, ignoring any constant terms (since they disappear after taking derivatives), and substituting into the right side of Equation 5, we get: Taking the derivative of Equation 7 with respect to , note that we can write Equation 7 as: to obtain: and setting it equal to zero, we get: tr tr P 4 a ( 1 ) 0( # # 4 a ( 1 4 a ( ) 0( 4 a ( where . Taking the derivative with respect to , we get: diag ) 0( $ 4 a ( 1 ) 0( DB QCA a # ) 0( DB QCA a # and where , implies that ( ( ) 0( # ! G B G T( ( ( ) ) or B G 4 a ( 1 C( ) ) 4 a ( 0#( ( ! G T( ! 4 a ( 1 C( ( 4 a ( 1 C( # )# )# diag 6 diag diag (7) Summarizing, the estimates of the new parameters in terms of the old parameters are as follows: Note that the above equations perform both the expectation step and the maximization step simultaneously. The algorithm proceeds by using the newly derived parameters as the guess for the next iteration. 4 Learning the parameters of an HMM, EM, and the BaumWelch algorithm
A Hidden Markov Model is a probabilistic model of the joint probability of a collection of random variables . The variables are either continuous or discrete observations and the variables are "hidden" and discrete. Under an HMM, there are two conditional independence assumptions made about these random variables that make associated algorithms tractable. These independence assumptions are 1), the hidden variable, given the hidden variable, is independent of previous variables, or: In this section, we derive the EM algorithm for finding the maximumlikelihood estimate of the parameters of a hidden Markov model given a set of observed feature vectors. This algorithm is also known as the BaumWelch algorithm. is a discrete random variable with possible values . We further assume that the underlying "hidden" Markov chain defined by is timehomogeneous (i.e., is independent of the time ). Therefore, we can represent as a timeindependent stochastic transition matrix . The special case of time is described by the initial state distribution, . We say that we are in state at time if .A particular sequence of states is described by where is the state at time . A particular observation sequence is described as . The probability of a particular observation vector at a particular time for state is described by: . The complete collection of parameters for all observation distributions is represented by . There are two forms of output distributions we will consider. The first is a discrete observation assumption where we assume that an observation is one of possible observation symbols 7 $ 6 $ 7 $ B B $ $ $ $ $ $ @) !!! P!!! $ % !! C Y !!! 6 6 % !! B C !!!P C C ! $ 1$ ( S ! 6 % ( 0 $ $ $ $ ) @) & '$ % 8 $ $ 6 B B B B 5 $ & '$ and 2), the observation, given the hidden variable, is independent of other variables, or: $ B a $ 4 a ( 0( ) # G T( ! G ( 4 a ( 1 0#( ) ) 4 a ( 1 C( # G ) 4 a ( 1 ( C#( ) 0( $ $ $ $ @) !!! @) 4 a ( & ($ # $ % B !!! B P!!! $ $ 1 $ @)
$ 1. Find for some . We use the forward (or the backward) procedure for this since it is much more efficient than direct evaluation. 2. Given some and some , find the best state sequence that explains The Viterbi algorithm solves this problem but we won't discuss it in this paper. HMMs) algorithm solves this problem, and we will develop it presently. In subsequent sections, we will consider only the first and third problems. The second is addressed in [RJ93]. 4.1 Efficient Calculation of Desired Quantities
One of the advantages of HMMs is that relatively efficient algorithms can be derived for the three problems mentioned above. Before we derive the EM algorithm directly using the function, we review these efficient procedures. Recall the forward procedure. We define 1. 2. 3. The backward procedure is similar: 1. 2. 3. 8 ! $ which is the probability of the ending partial sequence at time . We can efficiently define as: given that we started at state ! B !!! $ ! B $ which is the probability of seeing the partial sequence We can efficiently define recursively as: and ending up in state at time . !!!P $ $ $ ! !!!P B !!! $ 7 1 ( $ $ ( ( 1 ( ( ( )0#( $ ) a ( # ) ( ( 0( S ( ( ( # 3. Find argmax . The BaumWelch (also called forwardbackward or EM for 1 8 6 5 B C !!! C C . We describe the complete set of HMM parameters for a given model by: are three basic problems associated with HMMs: . In this case, if , then of probably distribution we consider is a mixture of . The second form multivariate Gaussians for each state where . There $ 6 C 8 $ $ ) G $ $ $ 1 8 B P!!! 1 7 ( ( S ) 0( 1 # ( % ( 8 P!!! $ ) . 6 ! $ 6! $ $ ! ! $ 6 a ) $ @) ! $ # $ ! 1 ! 1 which is the probability of being in state at time and being in state at time be expanded as: where is an indicator random variable that is when we are in state at time , and a random variable that is when we move from state to state after time . ! $ $ G H 6 ! E H 6 ! $ E G $ $ $ $ $ G 3 ! E H ! $ E G ( $ ( ( S ( $ $ ( ! $ B !! $ $ $ ! 6 B !! 0 ! 1 $ ( S ( )a C( ) ( $ ! $ # $ 6 ( 6 ) a# ( ( ( B P!!! ( ( $ $ $ 1b ! !!! ( ( which is the probability of being in state at time for the state sequence ! We now define $ ! ( so we can define things in terms of Also note that because of Markovian conditional independence We also define $ $ 6 ! 1 ( is the expected number of transitions from state to state for 6 ! ( ) $ B is the expected number of times in state and therefore is the expected number of transitions away from state for . Similarly, ! ! ) ( B S # ( ! or as: and If we sum these quantities across time, we can get some useful values. I.e., the expression and 9 $ $ & ( $ ! as . These follow from the fact that . Note that: . This can also is Jumping the gun a bit, our goal in forming an EM algorithm to estimate new parameters for the HMM by using the old parameters and the data. Intuitively, we can do this simply using relative frequencies. I.e., we can define update rules as follows: The quantity (8) is the expected relative frequency spent in state at time 1. The quantity is the expected number of transitions from state to state relative to the expected total number of transitions away from state . And, for discrete distributions, the quantity is the expected number of times the output observations have been equal to while in state relative to the expected total number of times in state . component of the mixture For Gaussian mixtures, we define the probability that the generated observation as where is a random variable indicating the mixture component at time for state . From the previous section on Gaussian Mixtures, we might guess that the update equations for this case are: 10 ( $ ( ( ( ( ) B $ ) ) ( G B ) ) $ B $ ) ) ( B ) ) $ & '$ # ( When there are come: observation sequences the being of length , the update equations be ! &$ '%! ! 8 B ( G & ($ $ ( ! $ 1 $ ( ( ( ( ( ) B $ "( ( ) B ( ) ( $ $ $B ( G ! ( ) B $ ( )B ( G $ $ ( ) B $ ( ) B$ ( ( )B ( ( 6 ( ! ) B$ ( ) B $ $ $ ! ! ( $ (9) S ( ! (10) $ $ ( These relatively intuitive equations are in fact the EM algorithm (or BalmWelch) for HMM parameter estimation. We derive these using the more typical EM notation in the next section. 4.2 Estimation formula using the The function then becomes: (11) Since the parameters we wish to optimize are now independently split into the three terms in the sum, we can optimize each term individually. The first term in Equation 11 becomes since by selecting all , we are simply repeatedly selecting the values of , so the right hand side is just the marginal expression for time . Adding the Lagrange multiplier , using the constraint that , and setting the derivative equal to zero, we get: For the remainder of the discussion any primed parameters are assumed to be the initial, guessed, or previous parameters whereas the unprimed parameters are being optimized. Note here that we assume the initial distribution starts at instead of for notational convenience. The basic results are the same however. 11 R C @) $ QCA DB where are our initial (or guessed, previous, etc.) estimates of the parameters and where space of all state sequences of length . Given a particular state sequence , representing is quite easy. I.e., C B $ ( ) ( D B 0( R ! C ( ECA R C @) EB A D # h I D B ( R C @) C ) ECA R R C U ( $ 0R C @) ) 0( # W U S h S R ! B' ) $ D EB A ) ) C 1 ( EB A C( D B C @) $ ( C D B ( R C ) QCA R # is: is the C C @) ( C ( ( We consider to be the observed data and the underlying state sequence to be hidden or unobserved. The incompletedata likelihood function is given by whereas the completedata likelihood function is . The function therefore ( B $ ) ) ( ( ) B ) function. and B ( G $ ( ) ) ( $ $ $B ( G ! ( B ) ) $ S B P!!! @) B C P!!! C R Taking the derivative, summing over to get , and solving for The second term in Equation 11 becomes: because for this term, we are, for each time looking over all transitions from to and weighting that by the corresponding probability the right hand side is just sum of the jointmarginal for time and . In a similar way, we can use a Lagrange multiplier with the constraint to get: The third term in Equation 11 becomes: because for this term, we are, for each time , looking at the emissions for all states and weighting each possible emission by the corresponding probability the right hand side is just the sum of the marginal for time . For discrete distributions, we can, again, use use a Lagrange multiplier but this time with the constraint . Only the observations that are equal to contribute to the probability value, so we get: For Gaussian Mixtures, the form of the function is slightly different, i.e., the hidden variables must include not only the hidden state sequence, but also a variable indicating the mixture component for each state at each time. Therefore, we can write as: where is the vector that indicates the mixture component for each state at each time. If we expand this as in Equation 11, the first and second terms are unchanged because the parameters are independent of which is thus marginalized away by the sum. The third term in Equation 11 becomes: This equation is almost identical to Equation 5, except for an addition sum component over the hidden state variables. We can optimize this in an exactly analogous way as we did in Section 3, and we get: 12 R ( a ) $ ! & '$ " S $ $ R 6 C ! C ) ( ) $ $C D ) ) 0( 1b" ( ( EB A R B # # ) $ $ D ) 0( R ! C 0 ( EB A B R # $ R ! $ C $ ) R 6 C ! C 6 ! R $ $ R $ ! C @$ ) ) CR ! C @) 8 C @) ( $ $ R ! C $ ) ) B $ " ( R ! C @) ) B , we get: S R
$ DEB A ) )a )0( R C ! R C@) ( @) $ C @) QCA DB B $ % # B 1 # !!! $ )B $ ( C @) ) ) ( $ $ C D ) @) EB A B I U $ )B ( $ )B B ( ! $ R S U D EB A ) S 6 ( a ) B DB ECA )
$ ( B $ ( $ ( ( As can be seen, these are the same set of update equations as given in the previous section. The update equations for HMMs with multiple observation sequences can similarly be derived and are addressed in [RJ93]. References
[ALR77] A.P.Dempster, N.M. Laird, and D.B. Rubin. Maximumlikelihood from incomplete data via the em algorithm. J. Royal Statist. Soc. Ser. B., 39, 1977. [Bis95] [GJ95] [JJ94] [JX96] [RJ93] C. Bishop. Neural Networks for Pattern Recognition. Clarendon Press, Oxford, 1995. Z. Ghahramami and M. Jordan. Learning from incomplete data. Technical Report AI Lab Memo No. 1509, CBCL Paper No. 108, MIT AI Lab, August 1995. M. Jordan and R. Jacobs. Hierarchical mixtures of experts and the em algorithm. Neural Computation, 6:181214, 1994. M. Jordon and L. Xu. Convergence results for the em approach to mixtures of experts architectures. Neural Networks, 8:14091431, 1996. L. Rabiner and B.H. Juang. Fundamentals of Speech Recognition. Prentice Hall Signal Processing Series, 1993. [RW84] R. Redner and H. Walker. Mixture densities, maximum likelihood and the em algorithm. SIAM Review, 26(2), 1984. [Wu83] [XJ96] C.F.J. Wu. On the convergence properties of the em algorithm. The Annals of Statistics, 11(1):95103, 1983. L. Xu and M.I. Jordan. On convergence properties of the em algorithm for gaussian mixtures. Neural Computation, 8:129151, 1996. 13 R $ $ $ R ! C @) ) B $ ! $ C @) B ( G $ ! ( G $ ) B$ ( and $ $ $ R ! C ) $ $ @$ ) B $ ( G R ! C @) ) B ...
View
Full
Document
This note was uploaded on 02/07/2012 for the course CSCI 5512 taught by Professor Staff during the Spring '08 term at Minnesota.
 Spring '08
 Staff

Click to edit the document details