EM algorithm and applications Lecture #9

1. EM algorithm and applicationsLecture #9 Background Readings: Chapters 11.2, 11.6 in the text book, Biological Sequence Analysis, Durbin et al., 2001. .

2. The EM algorithm This lecture plan: • Presentation and Correctness Proof of the EM algorithm. • Examples of Implementations

3. Model, Parameters, ML A “model with parameters θ ” is a probabilistic space M, in which each simple event y is determined by values of random variables (dice). The parameters θ are the probabilities associated with the random variables. (In HMM of length L, the simple events are HMM-sequences of length L, and the parameters are the transition probabilities mkland the emission probabilities ek(b)). An “observed data” is a non empty subset xM. (In HMM, it is usually all the simple events which fit with a given output sequence). Given observed data x, the ML method seeks parameters θ* which maximize the likelihood of the data p(x|θ)=∑yp(x,y|θ). (In HMM, x can be the transmitted letters ,and y the hidden states) Finding such θ* is easy when the observed data is a simple event, but hard in general.

4. The EM algorithm Assume a model with parameters as in the previous slide. Given observed data x, the likelihood of x under model parameters θis given by p(x|θ)=∑yp(x,y|θ). (The pairs (x,y)are the simple events which comprise x. Informally, y denotes the possible values of the“hidden data”). The EM algorithm receives x and parameters θ, and returns new parameters * s.t. p(x|*) ≥ p(x|θ), with equality only if λ*=θ. i.e., the new parameters increase the likelihood of the observed data.

5. The EM algorithm • The graphs below are the logarithms of the likelihood functions log P(x| λ) Log(Lθ)= E [log P(x,y|λ)] θ λ* λ EM uses the current parameters θ to construct a simpler ML problem Lθ: Guarantee: if Lθ(λ)>Lθ(θ), than P(x| λ)>P(x| θ).

6. Derivation of the EM Algorithm Let x be the observed data. Let {(x,y1),…,(x,yk)} be the set of (simple) events which comprise x. Our goal is to find parameters θ* which maximize the sum As this is hard, we start with some parameters θ, and only find λ*s.t. if λ*≠θ then: Finding λ* is obtained via “virtual sampling”, defined next.

7. For given parameters θ, Let pi=p(yi|x,θ). (note that p1+…+pk=1). We use the pi’sto define “virtual” sampling, in which: y1 occurs p1times, y2 occurs p2times, … yk occurs pktimes

8. The EM algorithm In each iteration the EM algorithm does the following. • (E step): Given θ, compute the function • (M step): Find* which maximizes Lθ() (Next iteration sets  * and repeat). Comment: At the M-step we only need that Lθ(*)>Lθ(θ). This change yields the so called Generalized EM algorithm. It is used when it is hard to find the optimal *. Usually, the computations use the function:

9. Correctness Theorem for the EM Algorithm

10. Correctness proof of EM

11. Correctness proof of EM (end)

12. Example: Baum Welsh = EM for HMM The Baum-Welsh algorithm is the EM algorithm for HMM: • E step for HMM: where λ are the new parameters {mkl,ek(b)}. • M step for HMM: look for λ which maximizes Lθ().

13. Baum Welsh = EM for HMM (cont) Mkl Ek(b)

14. A simple example: EM for 2 coin tosses Consider the following experiment: Given a coin with two possible outcomes: H (head) and T (tail), with probabilities θH, θT= 1- θH. The coin is tossed twice, but only the 1st outcome, T, is seen. So the data is x = (T,*). We wish to apply the EM algorithm to get parameters that increase the likelihood of the data. Let the initial parameters be θ = (θH, θT) = ( ¼, ¾ ).

15. EM for 2 coin tosses (cont) The “hidden data” which produce x are the sequences y1= (T,H); y2=(T,T); Hence the likelihood of x with parameters (θH, θT), is p(x| θ) = P(x,y1|) + P(x,y2|) = qHqT+qT2 For the initial parameters θ = ( ¼, ¾ ), we have: p(x| θ) = ¼ ∙ ¾ + ¾ ∙ ¾ = ¾ Note that in this case P(x,yi|) = P(yi|), for i = 1,2. we can always define y so that (x,y) = y (otherwise we set y’(x,y) and replace the “ y ”s by “ y’ ”s).

16. EM for 2 coin tosses - E step Calculate Lθ() = Lθ(λH,λT). Recall: λH,λTare the new parameters, which we need to optimize This is the “virtual sampling” p(y1|x,θ) = p(y1,x|θ)/p(x|θ) = (¾∙ ¼)/ (¾) = ¼ p(y2|x,θ) = p(y2,x|θ)/p(x|θ) = (¾∙ ¾)/ (¾) = ¾ Thus we have

17. EM for 2 coin tosses - E step For a sequence y of coin tosses, let NH(y) be the number of H’s in y, and NT(y) be the number of T’s in y. Then In our example: y1= (T,H); y2=(T,T), hence: NH(y1) = NT(y1)=1, NH(y2) =0, NT(y2)=2

18. Example: 2 coin tosses - E step Thus NT= 7/4 NH= ¼ And in general:

19. EM for 2 coin tosses - M step Find * which maximizes Lθ() And as we already saw, is maximized when: [The optimal parameters (0,1), will never be reached by the EM algorithm!]

20. EM for single random variable (dice) Now, the probability of each y(≡(x,y)) is given by a sequence of dice tosses. The dice has m outcomes, with probabilities λ1,..,λm. Let Nk(y) = #(outcome k occurs in y). Then Let Nk be the expected value of Nk(y), given x and θ: Nk=E(Nk|x,θ) = ∑y p(y|x,θ) Nk(y), Then we have:

21. L(λ) for one dice Nk

22. EM algorithm for n independent observations x1,…, xn : Expectation step It can be shown that, if the xjare independent, then:

23. The ABO locus has six possible genotypes {a/a, a/o, b/o, b/b, a/b, o/o}. The first two genotypes determine blood type A, the next two determine blood type B, then blood type AB, and finally blood type O. We wish to estimate the proportion in a population of the 6 genotypes. Example: The ABO locus A locus is a particular place on the chromosome. Each locus’ state (called genotype) consists of two alleles – one parental and one maternal. Some loci (plural of locus) determine distinguished features. The ABO locus, for example, determines blood type. Suppose we randomly sampled N individuals and found that Na/a have genotype a/a, Na/b have genotype a/b, etc. Then, the MLE is given by:

24. The ABO locus (Cont.) However, testing individuals for their genotype is a very expensive. Can we estimate the proportions of genotype using the common cheap blood test with outcome being one of the four blood types (A, B, AB, O) ? The problem is that among individuals measured to have blood type A, we don’t know how many have genotype a/a and how many have genotype a/o. So what can we do ?

25. The ABO locus (Cont.) The Hardy-Weinberg equilibrium rule states that in equilibrium the frequencies of the three alleles qa,qb,qoin the population determine the frequencies of the genotypes as follows: qa/b= 2qa qb, qa/o= 2qa qo, qb/o= 2qb qo, qa/a= [qa]2, qb/b= [qb]2, qo/o= [qo]2. In fact, Hardy-Weinberg equilibrium rule follows from modeling this problem as data x with hidden parameters y:

26. The ABO locus (Cont.) The dice’ outcome are the three possible alleles a, b and o. The observed data are the blood types A, B, AB or O. Each blood type isdetermined by two successive random sampling of alleles, which is an “ordered genotypes pair” – this is the hidden data. A ={(a,a), (a,o),(o,a)}; B={(b,b),(b,o),(o,b); AB={(a,b),(b,a)}; O={(o,o)}. So we have three parameters of one dice – qa,qb,qo - that we need to estimate. We start with parameters θ= (qa,qb,qo), and then use EM to improve them.

27. EM setting for the ABO locus The observed data x =(x1,..,xn) is a sequence of elements (blood types) from the set {A,B,AB,O}. eg: (B,A,B,B,O,A,B,A,O,B, AB) are observations (x1,…x11). The hidden data (ie the y’s)for each xjis the set of ordered pairs of alleles that generates it. For instance, for A it is the set {aa, ao, oa}. The parameters= {qa ,qb, qo} are the (current) probabilities of the alleles. The complete implementation of the EM algorithm for this problem will be given in the tutorial.