1 简介

EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的迭代由两步组成：E步，求期望，M步，求极大。
概率模型有时既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯法估计模型参数。但是，当模型含有隐变量时，就不能用简单地使用这些估计方法。
是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

2 投硬币问题

假设有三枚硬币，分别记A、B、C，这些硬币正面出现的概率分别是 π \pi π,p和q。进行如下掷硬币试验：先掷硬币A，正面继续掷硬币B，反面掷硬币C；然后掷选出来的硬币B或C，出现正面记为1，反面记为0；独立重复n次试验，观测结果如下
1 1 0 1 0 0 1 0 1 1
假设只能观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币正面出现的概率，即三硬币模型的参数。
解：
三个硬币的模型可写作
P ( y ∣ θ ) = ∑ z ( y , z ∣ θ ) = ∑ z P ( z ∣ θ ) P ( y ∣ z , θ ) = π p y ( 1 − p ) 1 − y + ( 1 − π ) q y ( 1 − q ) 1 − y P(y|\theta) = \sum_z(y,z|\theta) = \sum_z P(z|\theta) P(y|z,\theta)\\ =\pi p_y(1-p)^{1-y}+(1-\pi)q_y(1-q)^{1-y} P(y∣θ)=z∑​(y,z∣θ)=z∑​P(z∣θ)P(y∣z,θ)=πpy​(1−p)1−y+(1−π)qy​(1−q)1−y
随机变量y是观测变量，表示一次试验观测的结果是1或0；随机变量z是隐变量，表示未观测到的掷硬币的A的结果； θ = ( π , p , q ) \theta = (\pi,p,q) θ=(π,p,q)是模型参数，这一模型是以上数据的生成模型。注意，随机变量y的数据可以观测，随机变量z的数据是不可观测。
将观测的数据表示为 Y = （ Y 1 , Y 2 , . . . , Y n ） T Y= （Y_1,Y_2,...,Y_n）^T Y=（Y1​,Y2​,...,Yn​）T，未观测数据表示为 Z = （ Z 1 , Z 2 , . . . , Z n ） T Z= （Z_1,Z_2,...,Z_n）^T Z=（Z1​,Z2​,...,Zn​）T，则观测数据的似然函数为

P ( Y ∣ θ ) = ∑ Z P ( Z ∣ θ ) P ( Y ∣ Z , θ ) = ∏ j = 1 n [ π p y i ( 1 − p ) 1 − y i + ( 1 − π ) q y i ( 1 − q ) 1 − y i ] P(Y|\theta) = \sum_ZP(Z|\theta)P(Y|Z,\theta)\\ =\prod _{j=1}^n [\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}] P(Y∣θ)=Z∑​P(Z∣θ)P(Y∣Z,θ)=j=1∏n​[πpyi​(1−p)1−yi​+(1−π)qyi​(1−q)1−yi​]
考虑求模型参数 θ = ( π , p , q ) \theta = (\pi,p,q) θ=(π,p,q)的极大似然估计，即
θ ^ = a r g m a x θ l o g P ( Y ∣ θ ) \hat{\theta} = argmax_{\theta} logP(Y|\theta) θ^=argmaxθ​logP(Y∣θ)
以下通过属于迭代算法的EM算法求解
EM算法首先选取参数的初值，记作 θ 0 = ( π 0 , p 0 , q 0 ) \theta ^0 = (\pi^0,p^0,q^0) θ0=(π0,p0,q0)，然后通过下面的步骤迭代计算参数的估计值，直至收敛为止。第i次迭代参数的估计值为 θ i = ( π i , p i , q i ) \theta^i = (\pi^i,p^i,q^i) θi=(πi,pi,qi)。EM算法的第i+1次迭代如下：
E步：计算在模型参数 π i ， p i ， q i \pi^i，p^i，q^i πi，pi，qi下观测数据 y j y_j yj​来自掷硬币B的概率
μ j i + 1 = π i ( p i ) y i ( 1 − p i ) 1 − y i π i ( p i ) y i ( 1 − p i ) 1 − y i + ( 1 − π i ) ( q i ) y i ( 1 − q i ) 1 − y i \mu_j^{i+1} = \frac{\pi^i (p^i)^{y_i} (1-p^i)^{1-y_i}}{\pi^i(p^i)^{y_i}(1-p^i)^{1-y_i}+(1-\pi^i)(q^i)^{y_i}(1-q^i)^{1-y_i}} μji+1​=πi(pi)yi​(1−pi)1−yi​+(1−πi)(qi)yi​(1−qi)1−yi​πi(pi)yi​(1−pi)1−yi​​
M步：计算模型参数的新估计值
π i + 1 = 1 n ∑ j = 1 n μ j i + 1 \pi^{i+1} = \frac{1}{n}\sum_{j=1}^n \mu_j^{i+1} πi+1=n1​j=1∑n​μji+1​
p i + 1 = ∑ j = 1 n μ j i + 1 y j ∑ j = 1 n μ j i + 1 p^{i+1} = \frac{\sum_{j=1}^n\mu _j^{i+1}y_j}{\sum_{j=1}^n\mu_j^{i+1}} pi+1=∑j=1n​μji+1​∑j=1n​μji+1​yj​​
q i + 1 = ∑ j = 1 n ( 1 − μ j i + 1 y j ) ∑ j = 1 n ( 1 − μ j i + 1 ) q^{i+1} = \frac{\sum_{j=1}^n(1-\mu_j^{i+1}y_j)}{\sum_{j=1}^n(1-\mu_j^{i+1})} qi+1=∑j=1n​(1−μji+1​)∑j=1n​(1−μji+1​yj​)​

进行数字计算。假设模型参数的初值取
( π 0 , p 0 , q 0 ) = （ 0.5 , 0.5 , 0.5 ） (\pi^0,p^0,q^0) = （0.5,0.5,0.5） (π0,p0,q0)=（0.5,0.5,0.5）
由公式（3），对 y i = 1 ， y i = 0 y_i=1，y_i = 0 yi​=1，yi​=0均有 μ j 1 = 0.5 \mu_j^1 =0.5 μj1​=0.5，利用迭代公式（4）～（6）得到
( π 1 , p 1 , q 1 ) = （ 0.5 , 0.6 , 0.6 ） (\pi^1,p^1,q^1) = （0.5,0.6,0.6） (π1,p1,q1)=（0.5,0.6,0.6）
由于公式（2）得到 μ j 2 = 0.5 , j = 1 , 2 , . . , 10 \mu_j^{2} =0.5,j=1,2,..,10 μj2​=0.5,j=1,2,..,10
继续迭代，得到 ( π 2 , p 2 , q 2 ) = （ 0.5 , 0.6 , 0.6 ） (\pi^2,p^2,q^2) = （0.5,0.6,0.6） (π2,p2,q2)=（0.5,0.6,0.6）
于是得到模型参数 θ \theta θ的极大似然估计
( π ^ , p ^ , q ^ ) = （ 0.5 , 0.6 , 0.6 ） (\hat{\pi},\hat{p},\hat{q}) = （0.5,0.6,0.6） (π^,p^​,q^​)=（0.5,0.6,0.6）
如果 ( π 0 , p 0 , q 0 ) = （ 0.4 , 0.6 , 0.7 ） (\pi^0,p^0,q^0) = （0.4,0.6,0.7） (π0,p0,q0)=（0.4,0.6,0.7），那么得到的模型参数的极大似然估计 ( π ^ , p ^ , q ^ ) = （ 0.4064 , 0.5368 , 0.6432 ） (\hat{\pi},\hat{p},\hat{q}) = （0.4064,0.5368,0.6432） (π^,p^​,q^​)=（0.4064,0.5368,0.6432），也就是说，EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。
一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据，Y和Z连在一起称为完全数据，观测数据Y又称为不完全数据。假设给定观测数据Y，其概率分布是 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ)，其中 θ \theta θ是需要估计的模型参数，那么不完全数据Y的似然函数是 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ)，对数似然函数 L ( θ ) = l o g P ( Y ∣ θ ) L(\theta) = logP(Y|\theta) L(θ)=logP(Y∣θ)，假设Y和Z的联合概率分布是 P ( Y , Z ∣ θ ) P(Y,Z|\theta) P(Y,Z∣θ),那么完全数据的对数似然函数是 l o g P ( Y , Z ∣ θ ) logP(Y,Z|\theta) logP(Y,Z∣θ).

3 EM算法过程

输入：观测变量数据Y，隐变量数据Z，联合分布 P ( Y ， Z ∣ θ ) P(Y，Z|\theta) P(Y，Z∣θ)，条件分布 P ( Z ∣ Y , θ ) P(Z|Y,\theta) P(Z∣Y,θ)
输出：模型参数 θ \theta θ
（1）选择参数的初值 θ 0 \theta ^0 θ0，开始迭代
（2）E步：记 θ i \theta ^i θi为迭代参数 θ \theta θ的估计值，在第i+1次迭代的E步，计算
Q ( θ , θ i ) = E Z [ l o g P ( Y , Z ∣ θ ) ∣ Y , θ i ] = ∑ Z l o g P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ i ) Q(\theta,\theta^i) = E_Z[logP(Y,Z|\theta)|Y,\theta^i]\\ =\sum_Z logP(Y,Z|\theta)P(Z|Y,\theta^i) Q(θ,θi)=EZ​[logP(Y,Z∣θ)∣Y,θi]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θi)
其中， P ( Z ∣ Y , θ i ) P(Z|Y,\theta^i) P(Z∣Y,θi)是在给定观测数据Y和当前的参数估计 θ i \theta ^i θi下隐变量数据Z的条件概率分布。 Q ( θ , θ i ) Q(\theta,\theta^i) Q(θ,θi)称为Q函数。定义是完成数据的对数似然函数 l o g P ( Y , Z ∣ θ ) logP(Y,Z|\theta) logP(Y,Z∣θ)关于在给定观测数据Y和当前参数 θ i \theta^i θi下对未观测数据Z的条件概率分布 P ( Z ∣ Y , θ i ) P(Z|Y,\theta^i) P(Z∣Y,θi)的期望称为Q函数。
（3）M步：求使 Q ( θ , θ i ) Q(\theta,\theta^i) Q(θ,θi)极大化的 θ \theta θ，确定第i+1次迭代的参数的估计值 θ i + 1 \theta ^{i+1} θi+1
θ i + 1 = a r g m a x θ Q ( θ , θ i ) \theta ^{i+1} = argmax_{\theta}Q(\theta,\theta^i) θi+1=argmaxθ​Q(θ,θi)
（4）重复第（2）和第（3）步，直到收敛
注意：

• 步骤（1）参数的初值可以任意选择，但需要注意EM算法对初值是敏感的
• 步骤（2）E步求 Q ( θ , θ i ) Q(\theta,\theta^i) Q(θ,θi)，Q函数中Z是未观测数据，Y 是观测数据。注意 Q ( θ , θ i ) Q(\theta,\theta^i) Q(θ,θi)的第一个变元表示要极大化的参数，第二个变元表示餐忽视的当前估计值。每次迭代实际在求Q函数及其极大。
• 步骤（3）M步求 Q ( θ , θ i ) Q(\theta,\theta^i) Q(θ,θi)的极大化，得到 θ i + 1 \theta^{i+1} θi+1，完成一次迭代 θ i → θ i + 1 \theta^i \rightarrow \theta^{i+1} θi→θi+1.
• 步骤（4）给出停止迭代的条件，一般是对较小的正数 ϵ 1 ， ϵ 2 \epsilon_1，\epsilon_2 ϵ1​，ϵ2​，若满足 ∣ ∣ θ i + 1 − θ i ∣ ∣ < ϵ ||\theta^{i+1}-\theta^i||< \epsilon ∣∣θi+1−θi∣∣