EM算法

• 引言
• EM算法例子及解法
• EM算法步骤和说明

引言

EM 算法是一种选代算法，1977 年 Dempster 等人总结提出，用于含有隐变量(hidden variable)的概率模型参数的极大似然估计，或极大后验概率估计EM算法的每次选代由两步组成:E步，求期望 (expectation);M步，求极大(maximization)所以这一算法称为期望极大算法(expectation maximizationalgorithm)，简称EM算法。

EM算法例子及解法

三硬币模型: 假设有3 硬币，分别作$A，B，C$这些硬币正面出现的概率分别是$\pi$，$p$和$q$进行如下硬币试验:先硬币$A$，根据其结果选出硬币$B$或硬币$C$，正面选硬币$B$，反面选硬币$C$;然后掷选出的硬币，掷硬币的结果，出现正面记作 1，出现反面记作 0;独立地重复$n $次试验(这里，n=10)，观测结果如下:
$$1,1,0,1,0,0,1,0,1,1$$
假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。
模型表达是为：
P(y∣θ) = ∑ zP(y,z∣θ)= ∑ zP(z∣θ)P(y∣z,θ) =π p y(1−p ) y+(1−π) q y(1−q ) 1−y \begin{align} P(y|\theta) &= \sum_{z}P(y,z|\theta)=\sum_{z}P(z|\theta)P(y|z,\theta) \nonumber\\ &=\pi p^y (1-p)^y+(1-\pi)q^y(1-q)^{1-y} \nonumber \end{align} P(y∣θ)​=z∑​P(y,z∣θ)=z∑​P(z∣θ)P(y∣z,θ)=πpy(1−p)y+(1−π)qy(1−q)1−y​
这里，随机变量$y$是观测变量，表示一次试验观测的结果是 1或0;随机变量是隐变量，表示未观测到的掷硬币 $A$ 的结果;$\theta =(\pi ,p,q)$是模型参这一模型是以上数据的生成模型，注意，随机变量$y$的数据可以观测，随机变量$z$的数据不可观测。
将观察数据表示为 Y = ( Y1, Y1, . . . , Y1)T Y=(Y_1,Y_1,…,Y_1)^TY=(Y1​,Y1​,…,Y1​)T,未观察数据表示为 Z = ( Z1, Z2, . . . , Zn)T Z=(Z_1,Z_2,…,Z_n)^TZ=(Z1​,Z2​,…,Zn​)T,则观察数据的似然函数为
P(Y∣θ)= ∑ zP(Z∣θ)P(Y∣Z,θ)P(Y|\theta)=\sum_zP(Z|\theta)P(Y|Z,\theta) P(Y∣θ)=z∑​P(Z∣θ)P(Y∣Z,θ)
即
P(Y∣θ)= ∏ j=1n[π p yj (1−p ) yj +(1−π) q yj (1−q ) 1− y j ]P(Y|\theta)=\prod_{j=1}^{n}[\pi p^{y_j} (1-p)^{y_j}+(1-\pi)q^{y_j}(1-q)^{1-{y_j}}] P(Y∣θ)=j=1∏n​[πpyj​(1−p)yj​+(1−π)qyj​(1−q)1−yj​]
考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即
θ ^=arg max⁡θlogP(Y∣θ)\hat{\theta}=arg \max\limits_{\theta}logP(Y|\theta) θ^=argθmax​logP(Y∣θ)
这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法下面给出针对以上问题的EM算法其推导过程省略。
EM算法首先选取参数的初值，记作 θ ( 0 )= ( π ( 0 ), p ( 0 ), q ( 0 ))\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})θ(0)=(π(0),p(0),q(0)),然后通过下面的迭代参数的估计值，直至收敛为止。第$i$次迭代参数的估计值为 θ ( 0 )= ( π ( i ), p ( i ), q ( i ))\theta^{(0)}=(\pi^{(i)},p^{(i)},q^{(i)})θ(0)=(π(i),p(i),q(i))。EM 算法的第$i+1$次迭代如下。
E步：计算在模型参数 π ( i ), p ( i ), q ( i ) \pi^{(i)},p^{(i)},q^{(i)}π(i),p(i),q(i)下观测数据 yi y_iyi​来自掷硬币B的概率
μ i+1 = π( p (i)) yj (1−( p (i) ) ) yj π( p (i)) yj (1−( p (i) ) ) yj +(1−π)( q (i)) yj (1−( q (i) ) ) 1− y j\mu^{i+1}=\frac{\pi (p^{(i)})^{y_j} (1- (p^{(i)}))^{y_j}}{\pi (p^{(i)})^{y_j} (1- (p^{(i)}))^{y_j}+(1-\pi)(q^{(i)})^{y_j}(1-(q^{(i)}))^{1-{y_j}}} μi+1=π(p(i))yj​(1−(p(i)))yj​+(1−π)(q(i))yj​(1−(q(i)))1−yj​π(p(i))yj​(1−(p(i)))yj​​
M步：计算模型参数的新估计值
π (i+1) = 1 n ∑ j=1n μ j (i+1) \pi^{(i+1)}=\frac{1}{n} \sum_{j=1}^{n} \mu_j^{(i+1)} π(i+1)=n1​j=1∑n​μj(i+1)​
p (i+1) =∑ j=1n μ j (i+1)y j ∑ j=1n μ 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)}} p(i+1)=∑j=1n​μj(i+1)​∑j=1n​μj(i+1)​yj​​
q (i+1) =∑ j=1n(1− μ j (i+1) ) y j ∑ j=1n(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)})} q(i+1)=∑j=1n​(1−μj(i+1)​)∑j=1n​(1−μj(i+1)​)yj​​
进行数字计算，假设模型参数的初值为
π (0) =0.5, p (0) =0.5, q (0) =0.5\pi^{(0)}=0.5,p^{(0)}=0.5,q^{(0)}=0.5 π(0)=0.5,p(0)=0.5,q(0)=0.5
对 yj= 1y_j=1yj​=1与 yj= 0y_j=0yj​=0均有 μj ( 1 )= 0.5\mu_j^{(1)}=0.5μj(1)​=0.5。
根据M步计算得到
π (0) =0.5, p (0) =0.6, q (0) =0.6\pi^{(0)}=0.5,p^{(0)}=0.6,q^{(0)}=0.6 π(0)=0.5,p(0)=0.6,q(0)=0.6
根据E步，可以得到
μ j (2) =0.5,j=1,2,…,10\mu_j^{(2)}=0.5,j=1,2,…,10 μj(2)​=0.5,j=1,2,…,10
继续迭代，得
π (0) =0.5, p (0) =0.6, q (0) =0.6\pi^{(0)}=0.5,p^{(0)}=0.6,q^{(0)}=0.6 π(0)=0.5,p(0)=0.6,q(0)=0.6
于是得到模型参数$\theta$的极大似然估计：
π ^=0.5, p ^=0.6, q ^=0.6\hat{\pi}=0.5,\hat p=0.6,\hat q=0.6 π^=0.5,p^​=0.6,q^​=0.6
$\pi =0.5$表示硬币A是匀称的，这一结果容易理解。
如果取初始值 π ( 0 )= 0.4 , p ( 0 )= 0.6 , q ( 0 )= 0.7\pi^{(0)}=0.4,p^{(0)}=0.6,q^{(0)}=0.7π(0)=0.4,p(0)=0.6,q(0)=0.7,那么得到的模型参数的极大似然估计 π^= 0.4064 , p^= 0.5368 , q^= 0.6432\hat{\pi}=0.4064,\hat p=0.5368,\hat q=0.6432π^=0.4064,p^​=0.5368,q^​=0.6432。这就是说，EM算法与初始值的选择有关，选择不同的初值可能得到不同的参数估计值。

EM算法步骤和说明

一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据Y和Z连在一起称为完全数据 (complete-data)，观测数据Y又称为不完全数据(incomplete-data)。假设给定观测数据Y，其概率分布是$P(Y|\theta )$，其中是需要估计的模型参数，那么不完全数据Y的似然函数是 $P(Y|\theta )$，对数似然函数$L(0)=logP(Y|\theta )$;假设Y和Z的联合概率分布是$P(YZ|\theta )$，那完全的对数似然函数是$logP(Y,Z|\theta )$。
EM算法通过选代求$L(\theta )=logP(Y|\theta )$的极大似然估计每次代包含两步:E 步，求期望:M步，求极大化，下面来介绍 EM 算法。
输入:观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$;
输出:模型参数$\theta$。
(1)选择参数的初值 θ ( 0 ) \theta^{(0)}θ(0)，开始迭代;
(2)E步:记 θ ( i ) \theta^{(i)}θ(i)为第$i$次代参数的估计值，在第$i+1$次选代的E步，计算
Q(θ, θ (i) ) = E Z[logP(Y,Z∣θ)∣Y, θ (i) ] = ∑ ZlogP(Y,Z∣θ)P(Z∣Y, θ (i) )\begin{align} Q(\theta,\theta^{(i)}) &=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}] \nonumber\\ &=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)}) \nonumber\\ \end{align} Q(θ,θ(i))​=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))​​
这里，$P(Z|Y,\theta )$在给定观测数据$Y$和当前的参数估计 θ ( i ) \theta^{(i)}θ(i)下隐变量数据$Z$的条件概率分布；
(3)M步：求使 Q ( θ , θ ( i ))Q(\theta,\theta^{(i)})Q(θ,θ(i))极大化的$\theta$,确定第$i+1$次迭代的参数的估计值 θ ( i + 1 ) \theta^{(i+1)}θ(i+1)
θ (i+1) =arg⁡ max⁡θQ(θ, θ (i) )\theta^{(i+1)}=\arg \max \limits_{\theta}Q(\theta,\theta^{(i)}) θ(i+1)=argθmax​Q(θ,θ(i))
(4)重复第（2）步和第（3）步，直到收敛。
函数 Q ( θ , θ ( i ))Q(\theta,\theta^{(i)})Q(θ,θ(i))是EM算法的核心，称为Q函数（Q function）。
下面关于 EM 算法作几点说明:
步骤(1)参数的初值可以任意选择，但需注意 EM算法对初值是敏感的步骤(2)E步求 Q ( θ , θ ( i ))Q(\theta,\theta^{(i)})Q(θ,θ(i))。Q函数式中Z是未观察数据，Y是观测数据注意， Q ( θ , θ ( i ))Q(\theta,\theta^{(i)})Q(θ,θ(i))的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求$Q$函数及其极大。
步骤(3)M步求 Q ( θ , θ ( i ))Q(\theta,\theta^{(i)})Q(θ,θ(i))的极大化，得到 θ ( i + 1 ) \theta^{(i+1)}θ(i+1)，完成一次代 θ ( i )− > θ ( i + 1 ) \theta^{(i)}->\theta^{(i+1)}θ(i)−>θ(i+1)。后面将证明每次迭代使似然函数增大或达到局部极值。
步骤(4)给出停止迭代的条件，一般是对较小的正数 ϵ1, ϵ2 \epsilon_1,\epsilon_2ϵ1​,ϵ2​，若满足
∣∣ θ (i+1) − θ (i) ∣∣< ϵ 1或∣∣Q( θ (i+1) , θ (i) )−Q( θ (i) , θ (i) )∣∣< ϵ 2|| \theta^{(i+1)}- \theta^{(i)}||<\epsilon_1 或||Q(\theta^{(i+1)},\theta^{(i)})- Q(\theta^{(i)},\theta^{(i)})||<\epsilon_2 ∣∣θ(i+1)−θ(i)∣∣<ϵ1​或∣∣Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∣∣<ϵ2​
则停止迭代