统计学习方法 EM 算法

学习方法 EM 算法"/>

统计学习方法 EM 算法

文章目录

• 统计学习方法 EM 算法
• • 引入
  • EM 算法
  • EM 算法的导出
  • EM 算法的收敛性
  • 三硬币模型

统计学习方法 EM 算法

学习李航《统计学习方法》时关于 EM 算法的笔记

引入

概率模型中有时候同时包含观测变量（observable variable）和隐变量（潜在变量，latent variable）。

• 如果只有观测变量的话，那我们利用观测得到的数据，使用参数估计的方法（如极大似然估计法、矩估计法、贝叶斯估计法），就可以估计参数；
• 但如果存在隐变量的话，就要使用 EM 算法，相当于含隐变量的极大似然估计法，或极大后验概率估计法；

首先给出一个抛硬币的例子，之后算法的解释会更形象一些。

例（三硬币模型）：假设有三枚硬币，分别记为 A、B 和 C，其抛一次出现正面的概率分别为 π \pi π 、 p p p 和 q q q 。进行如下抛硬币实验：

• 首先抛硬币 A，若是正面则选择硬币 B，否则选择硬币 C；
• 接着抛被选中的硬币，出现正面记作 1，出现反面记作 0；

独立地重复 n n n 次该实验，假设这里 n = 10 n=10 n=10 ，假设某次实验出现如下观测结果：
1 , 1 , 0 , 1 , 0 , 0 , 1 , 0 , 1 , 1 1,\,1,\,0,\,1,\,0,\,0,\,1,\,0,\,1,\,1 1,1,0,1,0,0,1,0,1,1
假设只能观测到实验结果（即观测变量），而不能观测抛硬币的过程（这里硬币 A 的结果就是隐变量），问如何估计三硬币正面出现的概率 π \pi π 、 p p p 和 q q q （即模型的参数）

解：三硬币模型可以写作：
P ( y ∣ θ ) = ∑ z P ( y , z ∣ θ ) = ∑ z P ( z ∣ θ ) P ( y ∣ z , θ ) = π p y ( 1 − p ) 1 − y + ( 1 − π ) q y ( 1 − q ) 1 − y \begin{aligned} P(y|\theta)=&\,\sum\limits_{z}P(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} \end{aligned} P(y∣θ)==​z∑​P(y,z∣θ)=z∑​P(z∣θ)P(y∣z,θ)πpy(1−p)1−y+(1−π)qy(1−q)1−y​

• 随机变量 y y y 指的是某一次实验的观测结果，为 0 0 0 或 1 1 1 ；
• 随机变量 z z z 指的是某一次实验的硬币 A 的结果，也可以记为 0 0 0 或 1 1 1 ；
• θ \theta θ 为模型参数，即 θ = ( π , p , q ) \theta=(\pi,\,p,\,q) θ=(π,p,q) ；

以向量形式表示 n n n 次实验结果，观测数据 Y = ( Y 1 , Y 2 , ⋯ , Y n ) T Y=(Y_1,\,Y_2,\,\cdots,\,Y_n)^T Y=(Y1​,Y2​,⋯,Yn​)T ，隐变量 Z = ( Z 1 , Z 2 , ⋯ , Z n ) T Z=(Z_1,\,Z_2,\,\cdots,\,Z_n)^T Z=(Z1​,Z2​,⋯,Zn​)T ；

• 我们将 Y Y Y 和 Z Z Z 连在一起称为完全数据 ， Y Y Y 称为不完全数据；

则观测数据的似然函数为：
P ( Y ∣ θ ) = ∑ Z P ( Z ∣ θ ) P ( Y ∣ Z , θ ) P(Y|\theta)=\sum_{Z}P(Z|\theta)P(Y|Z,\,\theta) \\ P(Y∣θ)=Z∑​P(Z∣θ)P(Y∣Z,θ)
这里 Y = ( 1 , 1 , 0 , 1 , 0 , 0 , 1 , 0 , 1 , 1 ) T Y=(1,\,1,\,0,\,1,\,0,\,0,\,1,\,0,\,1,\,1)^T Y=(1,1,0,1,0,0,1,0,1,1)T ，每个 Z i Z_i Zi​ 都可能为 0 或 1，对 Z Z Z 求和就相当于遍历 Z Z Z 的排列组合的所有可能的情况，一共有 2 n 2^n 2n 种。但这个向量形式的式子其实不好列，我们还是对向量的每一项列出式子来，或者说按照乘法原理对每一次实验列式子出来，即：
P ( Y ∣ θ ) = ∏ i = 1 n ( ∑ z i P ( z i ∣ θ ) P ( y i ∣ z i , θ ) ) P(Y|\theta)=\prod_{i=1}^{n}\left(\sum_{z_{i}}P(z_i|\theta)P(y_i|z_i,\,\theta)\right) P(Y∣θ)=i=1∏n​(zi​∑​P(zi​∣θ)P(yi​∣zi​,θ))
这相当于将得到每个 y i y_i yi​ 的概率乘在一起，得到：
P ( Y ∣ θ ) = ∏ i = 1 n ( π p y i ( 1 − p ) 1 − y i + ( 1 − π ) q y i ( 1 − q ) 1 − y i ) P(Y|\theta)=\prod_{i=1}^{n}\left(\pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i}\right) P(Y∣θ)=i=1∏n​(πpyi​(1−p)1−yi​+(1−π)qyi​(1−q)1−yi​)
我们考虑求模型参数 θ = ( π , p , q ) \theta=(\pi,\,p,\,q) θ=(π,p,q) 的极大似然估计，即：
θ ^ = arg ⁡ max ⁡ θ log ⁡ P ( Y ∣ θ ) \hat{\theta}=\arg \max_{\theta} \log P(Y|\theta) θ^=argθmax​logP(Y∣θ)
书上说这个问题没有解析解（咱也不知道为啥），只能通过迭代的方式求解。这里给出这个问题的 EM 算法，具体推导在后边有：

首先选取参数的初始值， θ ( 0 ) = ( π ( 0 ) , p ( 0 ) , q ( 0 ) ) \theta^{(0)}=(\pi^{(0)},\,p^{(0)},\,q^{(0)}) θ(0)=(π(0),p(0),q(0)) ，然后反复执行以下两步，直到参数收敛；设第 i i i 次迭代得到的参数为 θ ( i ) = ( π ( i ) , p ( i ) , q ( i ) ) \theta^{(i)}=(\pi^{(i)},\,p^{(i)},\,q^{(i)}) θ(i)=(π(i),p(i),q(i)) ，第 i + 1 i+1 i+1 次迭代为：

• E 步：计算在模型参数 θ ( i ) \theta^{(i)} θ(i) 下，观测数据 y j y_j yj​ 来自 B 的概率为：（有点像贝叶斯公式）

μ j ( i + 1 ) = π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j + ( 1 − π ( i ) ) ( q ( i ) ) y j ( 1 − q ( i ) ) 1 − y j \mu_j^{(i+1)}= \frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}} {\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j} + (1-\pi^{(i)}) (q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}} μj(i+1)​=π(i)(p(i))yj​(1−p(i))1−yj​+(1−π(i))(q(i))yj​(1−q(i))1−yj​π(i)(p(i))yj​(1−p(i))1−yj​​

• M 步：计算模型参数的新估计值：

π ( i + 1 ) = 1 n ∑ j = 1 n μ j ( i + 1 ) p ( i + 1 ) = ∑ j = 1 n μ j ( i + 1 ) y j ∑ j = 1 n μ j ( i + 1 ) q ( i + 1 ) = ∑ j = 1 n ( 1 − μ j ) ( i + 1 ) y j ∑ j = 1 n ( 1 − μ j ) ( i + 1 ) \begin{aligned} \pi^{(i+1)}=&\,\frac{1}{n}\sum\limits_{j=1}^{n}\mu_j^{(i+1)} \\ p^{(i+1)}=&\,\frac{\sum\limits_{j=1}^{n}\mu_j^{(i+1)}y_j}{\sum\limits_{j=1}^{n}\mu_j^{(i+1)}} \\ q^{(i+1)}=&\,\frac{\sum\limits_{j=1}^{n}(1-\mu_j)^{(i+1)}y_j}{\sum\limits_{j=1}^{n}(1-\mu_j)^{(i+1)}} \\ \end{aligned} π(i+1)=p(i+1)=q(i+1)=​n1​j=1∑n​μj(i+1)​j=1∑n​μj(i+1)​j=1∑n​μj(i+1)​yj​​j=1∑n​(1−μj​)(i+1)j=1∑n​(1−μj​)(i+1)yj​​​

例如：

• 初始值 π ( 0 ) = 0.5 \pi^{(0)}=0.5 π(0)=0.5 ， p ( 0 ) = 0.5 p^{(0)}=0.5 p(0)=0.5 ， q ( 0 ) = 0.5 q^{(0)}=0.5 q(0)=0.5 ；
• 第 1 轮， μ j ( 1 ) = 0.5 \mu_j^{(1)}=0.5 μj(1)​=0.5 （ j = 1 , 2 , ⋯ , 10 j=1,\,2,\,\cdots,\,10 j=1,2,⋯,10） ，得到 π ( 1 ) = 0.5 \pi^{(1)}=0.5 π(1)=0.5 ， p ( 1 ) = 0.6 p^{(1)}=0.6 p(1)=0.6 ， q ( 1 ) = 0.6 q^{(1)}=0.6 q(1)=0.6 ；
• 第 2 轮， μ j ( 2 ) = 0.5 \mu_j^{(2)}=0.5 μj(2)​=0.5 （ j = 1 , 2 , ⋯ , 10 j=1,\,2,\,\cdots,\,10 j=1,2,⋯,10） ，得到 π ( 2 ) = 0.5 \pi^{(2)}=0.5 π(2)=0.5 ， p ( 2 ) = 0.6 p^{(2)}=0.6 p(2)=0.6 ， q ( 2 ) = 0.6 q^{(2)}=0.6 q(2)=0.6 ；

已经收敛了。但是，如果初值是 π ( 0 ) = 0.4 \pi^{(0)}=0.4 π(0)=0.4 ， p ( 0 ) = 0.6 p^{(0)}=0.6 p(0)=0.6 ， q ( 0 ) = 0.7 q^{(0)}=0.7 q(0)=0.7 ，则最终估计值为 π ^ = 0.4064 \hat\pi=0.4064 π^=0.4064 ， p ^ = 0.5368 \hat p=0.5368 p^​=0.5368 ， q ^ = 0.6432 \hat q=0.6432 q^​=0.6432 ，说明 EM 算法的与初值的选择有关。

EM 算法

EM 算法：通过迭代求 L ( θ ) = log ⁡ P ( Y ∣ θ ) L(\theta)=\log P(Y|\theta) L(θ)=logP(Y∣θ) 的极大似然估计；由于每次迭代包含两步：E 步求期望，M 步求极大化，因此称为 EM 算法。

• 输入：观测变量数据 Y Y Y ，隐变量数据 Z Z Z ，联合分布 P ( Y , Z ∣ θ ) P(Y,\,Z|\theta) P(Y,Z∣θ) ，条件分布 P ( Z ∣ Y , θ ) P(Z|Y,\,\theta) P(Z∣Y,θ) ；

• 输出：模型参数 θ \theta θ ；

① 选择模型参数的初始值 θ ( 0 ) \theta^{(0)} θ(0) ；

② E 步：记 θ ( i ) \theta^{(i)} θ(i) 为第 i i i 次迭代时参数 θ \theta θ 的估计值，则第 i + 1 i+1 i+1 次迭代时，计算：
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) \begin{aligned} Q(\theta,\,\theta^{(i)})=&\,E_Z[\log P(Y,\,Z|\theta)|Y,\,\theta^{(i)}] \\ =&\, \sum_Z\log P(Y,\,Z|\theta)P(Z|Y,\,\theta^{(i)}) \end{aligned} Q(θ,θ(i))==​EZ​[logP(Y,Z∣θ)∣Y,θ(i)]Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))​
这里 Y Y Y 和 θ ( i ) \theta^{(i)} θ(i) 是已知的， Q Q Q 函数是关于 θ \theta θ 的函数，相当于在 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\,\theta^{(i)}) P(Z∣Y,θ(i)) 的概率测度下，计算以 θ \theta θ 为参数的模型中 log ⁡ P ( Y , Z ∣ θ ) \log P(Y,\,Z|\theta) logP(Y,Z∣θ) 的期望（有一点点绕）。

③ M 步：求使得 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i)) 极大化的 θ \theta θ ，作为本次迭代的估计值 θ ( i + 1 ) \theta^{(i+1)} θ(i+1) ，即：
θ ( i + 1 ) = arg ⁡ max ⁡ θ Q ( θ , θ ( i + 1 ) ) \theta^{(i+1)}=\arg \max_{\theta} Q(\theta,\,\theta^{(i+1)}) θ(i+1)=argθmax​Q(θ,θ(i+1))
④ 重复迭代 ② 和 ③ ，直到参数收敛，可以是取较小的正数 ε 1 \varepsilon_1 ε1​ 和 ε 2 \varepsilon_2 ε2​ ，使得：
∣ ∣ θ ( i + 1 ) − θ ( i ) ∣ ∣ < ε 1 或 ∣ ∣ Q ( θ ( i + 1 ) , Q ( i ) ) − Q ( θ ( i ) , Q ( i ) ) ∣ ∣ < ε 2 ||\theta^{(i+1)}-\theta^{(i)}||\lt \varepsilon_1 \quad\text{或}\quad ||Q(\theta^{(i+1)},\,Q^{(i)})-Q(\theta^{(i)},\,Q^{(i)})||<\varepsilon_2 ∣∣θ(i+1)−θ(i)∣∣<ε1​或∣∣Q(θ(i+1),Q(i))−Q(θ(i),Q(i))∣∣<ε2​
Q Q Q 函数：是 EM 算法的核心，是完全数据的对数似然函数 log ⁡ P ( Y , Z ∣ θ ) \log P(Y,\,Z|\theta) logP(Y,Z∣θ) 关于在给定观测数据 Y Y Y 和当前参数 θ ( i ) \theta^{(i)} θ(i) 下对为观测数据 Z Z Z 的条件概率分布 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\,\theta^{(i)}) P(Z∣Y,θ(i)) 的期望称为 Q Q Q 函数，即：
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] Q(\theta,\,\theta^{(i)})=E_Z[\log P(Y,\,Z|\theta)|Y,\,\theta^{(i)}] Q(θ,θ(i))=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]

如果直接理解 Q Q Q 函数的意义有点困难的话，我们可以先来看这个式子：
E Z [ Z ∣ Y , θ ( i ) ] E_Z[Z|Y,\,\theta^{(i)}] EZ​[Z∣Y,θ(i)]
这个式子相当于给定观测结果 Y Y Y 和模型当前的参数 θ ( i ) \theta^{(i)} θ(i) ，对于某一种隐变量的情况 Z = z Z=z Z=z ，我们可以算出它的概率 P ( Z = z ∣ Y , θ ( i ) ) P(Z=z|Y,\,\theta^{(i)}) P(Z=z∣Y,θ(i)) （比如说已知抛硬币的最终结果和三种硬币得到正面的概率，我们可以算出中间过程究竟是选择硬币 B 还是硬币 C 的概率；可能会用到贝叶斯公式）。既然存在这样的概率，那我们就可以算出 Z Z Z 的期望，即上面的式子：
E Z [ Z ∣ Y , θ ( i ) ] = ∑ z z P ( Z = z ∣ Y , θ ( i ) ) （假设是离散的情况） E_Z[Z|Y,\,\theta^{(i)}] = \sum_{z} zP(Z=z|Y,\,\theta^{(i)}) \quad\text{（假设是离散的情况）} EZ​[Z∣Y,θ(i)]=z∑​zP(Z=z∣Y,θ(i))（假设是离散的情况）
接着， log ⁡ P ( Y , Z ∣ θ ) \log P(Y,\,Z|\theta) logP(Y,Z∣θ) 相当于是已知观测数据 Y Y Y 和把模型参数 θ \theta θ 当作已知（ θ \theta θ 是个变量，实际上未知）的情况下，得到完全数据 Y Y Y 和 Z Z Z 的结果的概率。即对于某一种隐变量的情况 Z = z Z=z Z=z ， log ⁡ P ( Y , Z ∣ θ ) \log P(Y,\,Z|\theta) logP(Y,Z∣θ) 是一个关于 θ \theta θ 的函数。同时，我们这里把模型参数 θ \theta θ 当作已知，则 log ⁡ P ( Y , Z ∣ θ ) \log P(Y,\,Z|\theta) logP(Y,Z∣θ) 也可以看成是变量 Z Z Z 的函数。我们对这个 Z Z Z 的函数求期望，就得到了：
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ z log ⁡ P ( Y , Z = z ∣ θ ) P ( Z = z ∣ Y , θ ( i ) ) （对这个 Z 的函数求期望） = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) （等价于这种写法） \begin{aligned} Q(\theta,\,\theta^{(i)})=&\,E_Z[\log P(Y,\,Z|\theta)|Y,\,\theta^{(i)}] \\ =&\,\sum_z\log P(Y,\,Z=z|\theta)P(Z=z|Y,\,\theta^{(i)}) \quad\text{（对这个$Z$的函数求期望）}\\ =&\,\sum_Z\log P(Y,\,Z|\theta)P(Z|Y,\,\theta^{(i)}) \quad\text{（等价于这种写法）} \end{aligned} Q(θ,θ(i))===​EZ​[logP(Y,Z∣θ)∣Y,θ(i)]z∑​logP(Y,Z=z∣θ)P(Z=z∣Y,θ(i))（对这个Z的函数求期望）Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))（等价于这种写法）​
这与上面关于 EM 算法的描述中，对 Q Q Q 函数的定义是等价的。

EM 算法的导出

面对一个含有隐变量的概率模型，我们的目标是找到一个参数 θ \theta θ （当然， θ \theta θ 是一个向量，可能包含多个参数），使得从该模型得到观测数据 Y Y Y 的概率极大化，等同于极大化观测数据 Y Y Y 关于参数 θ \theta θ 的对数似然函数：
L ( θ ) = log ⁡ P ( Y ∣ θ ) = log ⁡ ∑ Z P ( Y , Z ∣ θ ) = log ⁡ ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) \begin{aligned} L(\theta)=&\,\log P(Y|\theta)=\log \sum_ZP(Y,\,Z|\theta) \\ =&\,\log \left( \sum_Z P(Y|Z,\,\theta)P(Z|\theta) \right) \end{aligned} L(θ)==​logP(Y∣θ)=logZ∑​P(Y,Z∣θ)log(Z∑​P(Y∣Z,θ)P(Z∣θ))​
这个函数相当于要遍历 Z Z Z 的各种情况，然后求和，所以比较困难。EM 算法则是通过迭代逐步近似最大化 L ( θ ) L(\theta) L(θ) 的。假设第 i i i 次迭代后得到的估计值是 θ ( i ) \theta^{(i)} θ(i) ，我们希望重新估计一个 θ \theta θ 使得 L ( θ ) L(\theta) L(θ) 比前一个估计值更大，即 L ( θ ) > L ( θ ( i ) ) L(\theta) \gt L(\theta^{(i)}) L(θ)>L(θ(i)) ，因此，考虑二者的差：
L ( θ ) − L ( θ ( i ) ) = log ⁡ ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) − log ⁡ P ( Y ∣ θ ( i ) ) L(\theta)-L(\theta^{(i)})=\log \left( \sum_Z P(Y|Z,\,\theta)P(Z|\theta) \right)-\log P(Y|\theta^{(i)}) L(θ)−L(θ(i))=log(Z∑​P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))
根据 Jensen 不等式，因为 log ⁡ \log log 是个 concave function，因此有：
L ( θ ) − L ( θ ( i ) ) = log ⁡ ( ∑ Z P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) ) − log ⁡ P ( Y ∣ θ ( i ) ) ≥ ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) − log ⁡ P ( Y ∣ θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) − ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) \begin{aligned} L(\theta)-L(\theta^{(i)}) =&\,\log \left( \sum_Z P(Z|Y,\,\theta^{(i)}) \frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})} \right)-\log P(Y|\theta^{(i)}) \\ \ge &\, \sum_Z P(Z|Y,\,\theta^{(i)})\log \frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})} -\log P(Y|\theta^{(i)}) \\ =&\, \sum_Z P(Z|Y,\,\theta^{(i)})\log \frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})} - \sum_Z P(Z|Y,\,\theta^{(i)})\log P(Y|\theta^{(i)}) \\ =&\,\sum_Z P(Z|Y,\,\theta^{(i)})\frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})P(Y|\theta^{(i)})} \end{aligned} L(θ)−L(θ(i))=≥==​log(Z∑​P(Z∣Y,θ(i))P(Z∣Y,θ(i))P(Y∣Z,θ)P(Z∣θ)​)−logP(Y∣θ(i))Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣Z,θ)P(Z∣θ)​−logP(Y∣θ(i))Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣Z,θ)P(Z∣θ)​−Z∑​P(Z∣Y,θ(i))logP(Y∣θ(i))Z∑​P(Z∣Y,θ(i))P(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)​​

• 第一行到第二行是 Jensen 不等式，注意到 ∑ Z P ( Z ∣ Y , θ ( i ) ) = 1 \sum_Z P(Z|Y,\,\theta^{(i)})=1 ∑Z​P(Z∣Y,θ(i))=1 ，所以我们可以用这个不等式：

log ⁡ ∑ j λ j y j ≥ ∑ j λ j log ⁡ y i 其中  λ j ≥ 0 , ∑ j λ j = 1 \log\sum_{j}\lambda_jy_j\geq \sum_j \lambda_j\log y_i \quad\text{其中 }\lambda_j\geq0,\,\sum_j\lambda_j=1 logj∑​λj​yj​≥j∑​λj​logyi​其中 λj​≥0,j∑​λj​=1

• 第二行到第三行也是因为 ∑ Z P ( Z ∣ Y , θ ( i ) ) = 1 \sum_Z P(Z|Y,\,\theta^{(i)})=1 ∑Z​P(Z∣Y,θ(i))=1 ，这样我们就可以提取公因数

记一个新的函数为：
B ( θ , θ ( i ) ) = L ( θ ( i ) ) + ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) B(\theta,\,\theta^{(i)})=L(\theta^{(i)})+\sum_Z P(Z|Y,\,\theta^{(i)})\log \frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})P(Y|\theta^{(i)})} B(θ,θ(i))=L(θ(i))+Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)​
则：
L ( θ ) ≥ B ( θ , θ ( i ) ) L(\theta) \ge B(\theta,\,\theta^{(i)}) L(θ)≥B(θ,θ(i))
即 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) 是 L ( θ ) L(\theta) L(θ) 的一个下界。并且对于 B ( θ ( i ) , θ ( i ) ) B(\theta^{(i)},\,\theta^{(i)}) B(θ(i),θ(i)) ，有：（这是条件概率的公式）
{ P ( Y ∣ Z , θ ( i ) ) P ( Z ∣ θ ( i ) ) = P ( Y , Z ∣ θ ( i ) ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) = P ( Y , Z ∣ θ ( i ) ) \left\{ \begin{array}{l} P(Y|Z,\,\theta^{(i)})P(Z|\theta^{(i)})=P(Y,\,Z|\theta^{(i)}) \\ P(Z|Y,\,\theta^{(i)})P(Y|\theta^{(i)})=P(Y,\,Z|\theta^{(i)}) \end{array} \right. {P(Y∣Z,θ(i))P(Z∣θ(i))=P(Y,Z∣θ(i))P(Z∣Y,θ(i))P(Y∣θ(i))=P(Y,Z∣θ(i))​
故 B ( θ ( i ) , θ ( i ) ) = L ( θ ( i ) ) B(\theta^{(i)},\,\theta^{(i)})=L(\theta^{(i)}) B(θ(i),θ(i))=L(θ(i)) 。因此，任何可以使得 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) 增大的 θ \theta θ ，都可以使得 L ( θ ) L(\theta) L(θ) 增大，我们选择 θ \theta θ 使得 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) 达到极大，即：
θ ( i + 1 ) = arg ⁡ max ⁡ θ B ( θ , θ ( i ) ) \theta^{(i+1)}=\arg\max_{\theta}B(\theta,\,\theta^{(i)}) θ(i+1)=argθmax​B(θ,θ(i))
有：
θ ( i + 1 ) = arg ⁡ max ⁡ θ B ( θ , θ ( i ) ) = arg ⁡ max ⁡ θ ( L ( θ ( i ) ) + ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) ) = arg ⁡ max ⁡ θ ( ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) = arg ⁡ max ⁡ θ ( ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y , Z ∣ θ ) ) = arg ⁡ max ⁡ θ Q ( θ , θ ( i ) ) \begin{aligned} \theta^{(i+1)} =&\, \arg\max_{\theta}B(\theta,\,\theta^{(i)}) \\ =&\, \arg\max_{\theta} \left( L(\theta^{(i)})+\sum_Z P(Z|Y,\,\theta^{(i)})\log \frac{P(Y|Z,\,\theta)P(Z|\theta)}{P(Z|Y,\,\theta^{(i)})P(Y|\theta^{(i)})} \right) \\ =&\, \arg\max_{\theta} \left( \sum_Z P(Z|Y,\,\theta^{(i)})\log P(Y|Z,\,\theta)P(Z|\theta) \right) \\ =&\, \arg\max_{\theta} \left( \sum_Z P(Z|Y,\,\theta^{(i)})\log P(Y,\,Z|\theta) \right) \\ =&\, \arg\max_{\theta} Q(\theta,\,\theta^{(i)}) \end{aligned} θ(i+1)=====​argθmax​B(θ,θ(i))argθmax​(L(θ(i))+Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)​)argθmax​(Z∑​P(Z∣Y,θ(i))logP(Y∣Z,θ)P(Z∣θ))argθmax​(Z∑​P(Z∣Y,θ(i))logP(Y,Z∣θ))argθmax​Q(θ,θ(i))​

• 第二行到第三行是因为，在第 i + 1 i+1 i+1 次迭代中， L ( θ ( i ) ) L(\theta^{(i)}) L(θ(i)) 已经是常数了， P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) P(Z|Y,\,\theta^{(i)})P(Y|\theta^{(i)}) P(Z∣Y,θ(i))P(Y∣θ(i)) 也是常数，可以提取公因子；

用图形进行直观解释：上方的曲线为 L ( θ ) L(\theta) L(θ) ，下方的曲线为 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) ，二者在点 θ = θ ( i ) \theta=\theta^{(i)} θ=θ(i) 处相等。当选择下一个点 θ ( i + 1 ) \theta^{(i+1)} θ(i+1) 使得 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) 极大化（也是使 Q ( θ , θ ( i ) ) Q(\theta,\,\theta^{(i)}) Q(θ,θ(i)) 极大化），这时由于 L ( θ ) ≥ B ( θ , θ ( i ) ) L(\theta) \ge B(\theta,\,\theta^{(i)}) L(θ)≥B(θ,θ(i)) ，因此 L ( θ ) L(\theta) L(θ) 在每次迭代中也是增加的（只要 B ( θ , θ ( i ) ) B(\theta,\,\theta^{(i)}) B(θ,θ(i)) 增加， L ( θ ) L(\theta) L(θ) 一定增加，因为二者在点 θ = θ ( i ) \theta=\theta^{(i)} θ=θ(i) 处相等）。

当然也可以看出，EM 算法不能保证找到全局最优值。

EM 算法的收敛性

这里有关于 EM 算法收敛性的两个定理。

定理 9.1：设 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) 为观测数据的似然函数， θ ( i ) ( i = 1 , 2 , ⋯ ) \theta^{(i)}\,(i=1,\,2,\,\cdots) θ(i)(i=1,2,⋯) 为 EM 算法得到的参数估计序列，而 P ( Y ∣ θ ( i ) ) ( i = 1 , 2 , ⋯ ) P(Y|\theta^{(i)})\,(i=1,\,2,\,\cdots) P(Y∣θ(i))(i=1,2,⋯) 为对应的似然函数序列，则 P ( Y ∣ θ ( i ) ) P(Y|\theta^{(i)}) P(Y∣θ(i)) 是单调递增的，即：
P ( Y ∣ θ ( i + 1 ) ) ≥ P ( Y ∣ θ ( i ) ) P(Y|\theta^{(i+1)}) \ge P(Y|\theta^{(i)}) P(Y∣θ(i+1))≥P(Y∣θ(i))
证明：根据条件概率的公式有：
P ( Y ∣ θ ) = P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ) P(Y|\theta)=\frac{P(Y,\,Z|\theta)}{P(Z|Y,\,\theta)} P(Y∣θ)=P(Z∣Y,θ)P(Y,Z∣θ)​
取对数有：
log ⁡ P ( Y ∣ θ ) = log ⁡ P ( Y , Z ∣ θ ) − log ⁡ P ( Z ∣ Y , θ ) \log P(Y|\theta) = \log P(Y,\,Z|\theta) - \log P(Z|Y,\,\theta) logP(Y∣θ)=logP(Y,Z∣θ)−logP(Z∣Y,θ)
注意这里 Z Z Z 是一个具体的值，所以不需要对所有 Z Z Z 的取值情况求和；如果理解不了，可以写成以下形式：
P ( Y = y ∣ θ ) = P ( Y = y , Z = z ∣ θ ) P ( Z = z ∣ Y = y , θ ) P(Y=y|\theta)=\frac{P(Y=y,\,Z=z|\theta)}{P(Z=z|Y=y,\,\theta)} P(Y=y∣θ)=P(Z=z∣Y=y,θ)P(Y=y,Z=z∣θ)​
Q Q Q 函数的定义为：
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) Q(\theta,\,\theta^{(i)})=E_Z[\log P(Y,\,Z|\theta)|Y,\,\theta^{(i)}]=\sum_Z\log P(Y,\,Z|\theta)P(Z|Y,\,\theta^{(i)}) Q(θ,θ(i))=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))
再定义另一个函数：
H ( θ , θ ( i ) ) = ∑ Z log ⁡ P ( Z ∣ Y , θ ) P ( Z ∣ Y , θ ( i ) ) H(\theta,\,\theta^{(i)})=\sum_Z\log P(Z|Y,\,\theta)P(Z|Y,\,\theta^{(i)}) H(θ,θ(i))=Z∑​logP(Z∣Y,θ)P(Z∣Y,θ(i))
则对数似然函数可以写成：
log ⁡ P ( Y ∣ θ ) = Q ( θ , θ ( i ) ) − H ( θ , θ ( i ) ) \log P(Y|\theta) = Q(\theta,\,\theta^{(i)})-H(\theta,\,\theta^{(i)}) logP(Y∣θ)=Q(θ,θ(i))−H(θ,θ(i))
是不是很神奇，跟 i i i 一点关系没有。这个式子可以这样理解：
Q ( θ , θ ( i ) ) − H ( θ , θ ( i ) ) = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) − ∑ Z log ⁡ P ( Z ∣ Y , θ ) P ( Z ∣ Y , θ ( i ) ) = ∑ Z ( log ⁡ P ( Y , Z ∣ θ ) − log ⁡ P ( Z ∣ Y , θ ) ) P ( Z ∣ Y , θ ( i ) ) \begin{aligned} &\, Q(\theta,\,\theta^{(i)})-H(\theta,\,\theta^{(i)}) \\ =&\, \sum_Z\log P(Y,\,Z|\theta)P(Z|Y,\,\theta^{(i)})-\sum_Z\log P(Z|Y,\,\theta)P(Z|Y,\,\theta^{(i)}) \\ =&\, \sum_Z(\log P(Y,\,Z|\theta)-\log P(Z|Y,\,\theta))P(Z|Y,\,\theta^{(i)}) \end{aligned} ==​Q(θ,θ(i))−H(θ,θ(i))Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))−Z∑​logP(Z∣Y,θ)P(Z∣Y,θ(i))Z∑​(logP(Y,Z∣θ)−logP(Z∣Y,θ))P(Z∣Y,θ(i))​
前面说了，对于前面的取了对数的条件概率公式， Z Z Z 是可以任意取的，因此：
∑ Z ( log ⁡ P ( Y , Z ∣ θ ) − log ⁡ P ( Z ∣ Y , θ ) ) P ( Z ∣ Y , θ ( i ) ) = ∑ Z log ⁡ P ( Y ∣ θ ) P ( Z ∣ Y , θ ( i ) ) = log ⁡ P ( Y ∣ θ ) ∑ Z P ( Z ∣ Y , θ ( i ) ) = log ⁡ P ( Y ∣ θ ) \begin{aligned} &\, \sum_Z(\log P(Y,\,Z|\theta)-\log P(Z|Y,\,\theta))P(Z|Y,\,\theta^{(i)}) \\ =&\, \sum_Z\log P(Y|\theta)P(Z|Y,\,\theta^{(i)}) \\ =&\, \log P(Y|\theta)\sum_ZP(Z|Y,\,\theta^{(i)}) \\ =&\, \log P(Y|\theta) \end{aligned} ===​Z∑​(logP(Y,Z∣θ)−logP(Z∣Y,θ))P(Z∣Y,θ(i))Z∑​logP(Y∣θ)P(Z∣Y,θ(i))logP(Y∣θ)Z∑​P(Z∣Y,θ(i))logP(Y∣θ)​
接着我们要证明 log ⁡ P ( Y ∣ θ ( i ) ) \log P(Y|\theta^{(i)}) logP(Y∣θ(i)) 是单调递增的，有：
log ⁡ P ( Y ∣ θ ( i + 1 ) ) − log ⁡ P ( Y ∣ θ ( i ) ) = [ Q ( θ ( i + 1 ) , θ ( i ) ) − Q ( θ ( i ) , θ ( i ) ) ] − [ H ( θ ( i + 1 ) , θ ( i ) ) − H ( θ ( i ) , θ ( i ) ) ] \begin{aligned} &\, \log P(Y|\theta^{(i+1)})-\log P(Y|\theta^{(i)}) \\ =&\, [Q(\theta^{(i+1)},\,\theta^{(i)})-Q(\theta^{(i)},\,\theta^{(i)})]-[H(\theta^{(i+1)},\,\theta^{(i)})-H(\theta^{(i)},\,\theta^{(i)})] \end{aligned} =​logP(Y∣θ(i+1))−logP(Y∣θ(i))[Q(θ(i+1),θ(i))−Q(θ(i),θ(i))]−[H(θ(i+1),θ(i))−H(θ(i),θ(i))]​
前面证明了， θ ( i + 1 ) \theta^{(i+1)} θ(i+1) 使得 Q ( θ , θ ( i ) ) Q(\theta,\,\theta^{(i)}) Q(θ,θ(i)) 达到极大，因此右边第一项：
Q ( θ ( i + 1 ) , θ ( i ) ) − Q ( θ ( i ) , θ ( i ) ) ≥ 0 Q(\theta^{(i+1)},\,\theta^{(i)})-Q(\theta^{(i)},\,\theta^{(i)}) \geq 0 Q(θ(i+1),θ(i))−Q(θ(i),θ(i))≥0
右边第二项：
H ( θ ( i + 1 ) , θ ( i ) ) − H ( θ ( i ) , θ ( i ) ) = ∑ Z log ⁡ P ( Z ∣ Y , θ ( i + 1 ) ) P ( Z ∣ Y , θ ( i ) ) − ∑ Z log ⁡ P ( Z ∣ Y , θ ( i ) ) P ( Z ∣ Y , θ ( i ) ) = ∑ Z ( log ⁡ P ( Z ∣ Y , θ ( i + 1 ) ) P ( Z ∣ Y , θ ( i ) ) ) P ( Z ∣ Y , θ ( i ) ) ≤ log ⁡ ( ∑ Z P ( Z ∣ Y , θ ( i + 1 ) ) P ( Z ∣ Y , θ ( i ) ) P ( Z ∣ Y , θ ( i ) ) ) = 0 \begin{aligned} &\, H(\theta^{(i+1)},\,\theta^{(i)})-H(\theta^{(i)},\,\theta^{(i)}) \\ =&\, \sum_Z\log P(Z|Y,\,\theta^{(i+1)})P(Z|Y,\,\theta^{(i)}) - \sum_Z\log P(Z|Y,\,\theta^{(i)})P(Z|Y,\,\theta^{(i)}) \\ =&\, \sum_Z \left( \log\frac{P(Z|Y,\,\theta^{(i+1)})}{P(Z|Y,\,\theta^{(i)})} \right)P(Z|Y,\,\theta^{(i)}) \\ \leq&\, \log \left( \sum_Z \frac{P(Z|Y,\,\theta^{(i+1)})}{P(Z|Y,\,\theta^{(i)})}P(Z|Y,\,\theta^{(i)}) \right)=0 \end{aligned} ==≤​H(θ(i+1),θ(i))−H(θ(i),θ(i))Z∑​logP(Z∣Y,θ(i+1))P(Z∣Y,θ(i))−Z∑​logP(Z∣Y,θ(i))P(Z∣Y,θ(i))Z∑​(logP(Z∣Y,θ(i))P(Z∣Y,θ(i+1))​)P(Z∣Y,θ(i))log(Z∑​P(Z∣Y,θ(i))P(Z∣Y,θ(i+1))​P(Z∣Y,θ(i)))=0​

• 小于等于由 Jensen 不等式得到，和前面使用 Jensen 不等式导出 EM 算法的技巧是一样的

所以 log ⁡ P ( Y ∣ θ ( i ) ) \log P(Y|\theta^{(i)}) logP(Y∣θ(i)) 是单调递增的。

定理 9.2：设 L ( θ ) = log ⁡ P ( Y ∣ θ ) L(\theta)=\log P(Y|\theta) L(θ)=logP(Y∣θ) 为观测数据的对数似然函数，， θ ( i ) ( i = 1 , 2 , ⋯ ) \theta^{(i)}\,(i=1,\,2,\,\cdots) θ(i)(i=1,2,⋯) 为 EM 算法得到的参数估计序列，而 L ( θ ( i ) ) ( i = 1 , 2 , ⋯ ) L(\theta^{(i)})\,(i=1,\,2,\,\cdots) L(θ(i))(i=1,2,⋯) 为对应的似然函数序列：

1. 若 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) 有上界，则 L ( θ ( i ) ) = log ⁡ P ( Y ∣ θ ( i ) ) L(\theta^{(i)})=\log P(Y|\theta^{(i)}) L(θ(i))=logP(Y∣θ(i)) 收敛到某一值 L ∗ L^\ast L∗ ；
2. 在函数 Q ( θ , θ ( i ) ) Q(\theta,\,\theta^{(i)}) Q(θ,θ(i)) 与 L ( θ ) L(\theta) L(θ) 满足一定条件下，由 EM 算法得到的参数估计序列 θ ( i ) \theta^{(i)} θ(i) 的收敛值 θ ∗ \theta^\ast θ∗ 是 L ( θ ) L(\theta) L(θ) 的稳定点。

第一点由定理 9.1 得到的 L ( θ ) L(\theta) L(θ) 的单调性和有界性可以得到；第二点参考茆的高级数理统计（

三硬币模型

现在重新来看三硬币模型，这个模型的隐变量是硬币 A 的结果 Z = ( z 1 , z 2 , ⋯ , z n ) Z=(z_1,\,z_2,\,\cdots,\,z_n) Z=(z1​,z2​,⋯,zn​) 。

其先求其对数似然函数：
log ⁡ P ( Y , Z ∣ θ ) = log ⁡ [ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ] \log P(Y,\,Z|\theta)=\log [P(Y|Z,\,\theta)P(Z|\theta) ] logP(Y,Z∣θ)=log[P(Y∣Z,θ)P(Z∣θ)]
按照条件概率的公式，其实也可以写成：
log ⁡ P ( Y , Z ∣ θ ) = log ⁡ [ P ( Z ∣ Y , θ ) P ( Y ∣ θ ) ] \log P(Y,\,Z|\theta)=\log [P(Z|Y,\,\theta)P(Y|\theta) ] logP(Y,Z∣θ)=log[P(Z∣Y,θ)P(Y∣θ)]
但是你想想抛硬币的过程，显然第一种是顺着抛硬币的过程写的，这样比较好算；有：
P ( Z ∣ θ ) = ∏ j = 1 n π z j ( 1 − π ) ( 1 − z j ) P(Z|\theta)=\prod_{j=1}^{n}\pi^{z_j}(1-\pi)^{(1-z_j)} P(Z∣θ)=j=1∏n​πzj​(1−π)(1−zj​)
下面这个式子可能比较难理解，也是让我想了很久，可以按照 z i z_i zi​ 为 0 或者为 1 来分类讨论：
P ( Y ∣ Z , θ ) = ( p y j ( 1 − p ) ( 1 − y j ) ) z j ( q y j ( 1 − q ) ( 1 − y j ) ) ( 1 − z j ) P(Y|Z,\,\theta)=\left( p^{y_j}(1-p)^{(1-y_j)} \right)^{z_j}\left( q^{y_j}(1-q)^{(1-y_j)} \right)^{(1-z_j)} P(Y∣Z,θ)=(pyj​(1−p)(1−yj​))zj​(qyj​(1−q)(1−yj​))(1−zj​)
则整理得到：
P ( Y , Z ∣ θ ) = ∏ j = 1 n ( π p y j ( 1 − p ) ( 1 − y j ) ) z j ( ( 1 − π ) q y j ( 1 − q ) ( 1 − y j ) ) ( 1 − z j ) P(Y,\,Z|\theta)=\prod_{j=1}^{n}\left( \pi p^{y_j}(1-p)^{(1-y_j)} \right)^{z_j}\left( (1-\pi)q^{y_j}(1-q)^{(1-y_j)} \right)^{(1-z_j)} P(Y,Z∣θ)=j=1∏n​(πpyj​(1−p)(1−yj​))zj​((1−π)qyj​(1−q)(1−yj​))(1−zj​)
其实可以发现：
P ( Y = y j , Z = z j ∣ θ ) = ( π p y j ( 1 − p ) ( 1 − y j ) ) z j ( ( 1 − π ) q y j ( 1 − q ) ( 1 − y j ) ) ( 1 − z j ) P(Y=y_j,\,Z=z_j|\theta)=\left( \pi p^{y_j}(1-p)^{(1-y_j)} \right)^{z_j}\left( (1-\pi)q^{y_j}(1-q)^{(1-y_j)} \right)^{(1-z_j)} P(Y=yj​,Z=zj​∣θ)=(πpyj​(1−p)(1−yj​))zj​((1−π)qyj​(1−q)(1−yj​))(1−zj​)
有了似然函数，我们还需要求出 P ( Z ∣ Y , θ ) P(Z|Y,\,\theta) P(Z∣Y,θ) （因为需要 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\,\theta^{(i)}) P(Z∣Y,θ(i)) 来算期望），按照条件概率的公式有：
P ( Z ∣ Y , θ ) = P ( Y , Z ∣ θ ) P ( Y ∣ θ ) P(Z|Y,\,\theta)=\frac{P(Y,\,Z|\theta)}{P(Y|\theta)} P(Z∣Y,θ)=P(Y∣θ)P(Y,Z∣θ)​
我们现在有 P ( Y , Z ∣ θ ) P(Y,\,Z|\theta) P(Y,Z∣θ) 了，可以求一下 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) ，有：
P ( Y ∣ θ ) = ∏ j = 1 n ( π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j ) P(Y|\theta)=\prod_{j=1}^{n}\left(\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi) q^{y_j}(1-q)^{1-y_j}\right) P(Y∣θ)=j=1∏n​(πpyj​(1−p)1−yj​+(1−π)qyj​(1−q)1−yj​)
这个式子可以直接按照含义列出来，但其实也可以通过对 P ( Y , Z ∣ θ ) P(Y,\,Z|\theta) P(Y,Z∣θ) 求和得到。比较难理解，有一点像二项式展开，对每个 z j z_j zj​ 取 0 和 1 的情况遍历， ( π p y i ( 1 − p ) ( 1 − y i ) ) z i ( ( 1 − π ) q y i ( 1 − q ) ( 1 − y i ) ) ( 1 − z i ) \left( \pi p^{y_i}(1-p)^{(1-y_i)} \right)^{z_i}\left( (1-\pi)q^{y_i}(1-q)^{(1-y_i)} \right)^{(1-z_i)} (πpyi​(1−p)(1−yi​))zi​((1−π)qyi​(1−q)(1−yi​))(1−zi​) 就像是从每一个括号中选出其中一个来乘：
P ( Y ∣ θ ) = ∑ Z P ( Y , Z ∣ θ ) = ∏ j = 1 n ( π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j ) P(Y|\theta)=\sum_Z P(Y,\,Z|\theta)=\prod_{j=1}^{n}\left(\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi) q^{y_j}(1-q)^{1-y_j}\right) P(Y∣θ)=Z∑​P(Y,Z∣θ)=j=1∏n​(πpyj​(1−p)1−yj​+(1−π)qyj​(1−q)1−yj​)
那我们就可以得到 P ( Z ∣ Y , θ ) P(Z|Y,\,\theta) P(Z∣Y,θ) 了：
P ( Z ∣ Y , θ ) = P ( Y , Z ∣ θ ) P ( Y ∣ θ ) = ∏ j = 1 n ( π p y j ( 1 − p ) ( 1 − y j ) ) z j ( ( 1 − π ) q y j ( 1 − q ) ( 1 − y j ) ) ( 1 − z j ) π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j P(Z|Y,\,\theta)=\frac{P(Y,\,Z|\theta)}{P(Y|\theta)}=\prod_{j=1}^{n} \frac {\left( \pi p^{y_j}(1-p)^{(1-y_j)} \right)^{z_j}\left( (1-\pi)q^{y_j}(1-q)^{(1-y_j)} \right)^{(1-z_j)}} {\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi) q^{y_j}(1-q)^{1-y_j}} P(Z∣Y,θ)=P(Y∣θ)P(Y,Z∣θ)​=j=1∏n​πpyj​(1−p)1−yj​+(1−π)qyj​(1−q)1−yj​(πpyj​(1−p)(1−yj​))zj​((1−π)qyj​(1−q)(1−yj​))(1−zj​)​
我们引入记号，在模型参数 θ ( i ) \theta^{(i)} θ(i) 下，观测数据 y j y_j yj​ 来自 B 的概率为：（有点像贝叶斯公式）
μ j ( i + 1 ) = π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j + ( 1 − π ( i ) ) ( q ( i ) ) y j ( 1 − q ( i ) ) 1 − y j \mu_j^{(i+1)}= \frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}} {\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j} + (1-\pi^{(i)}) (q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}} μj(i+1)​=π(i)(p(i))yj​(1−p(i))1−yj​+(1−π(i))(q(i))yj​(1−q(i))1−yj​π(i)(p(i))yj​(1−p(i))1−yj​​

则：
P ( Z ∣ Y , θ ( i ) ) = ∏ j = 1 n ( μ j ( i + 1 ) ) z j ( 1 − μ j ( i + 1 ) ) ( 1 − z j ) P(Z|Y,\,\theta^{(i)})=\prod_{j=1}^{n}(\mu_j^{(i+1)})^{z_j}(1-\mu_j^{(i+1)})^{(1-z_j)} P(Z∣Y,θ(i))=j=1∏n​(μj(i+1)​)zj​(1−μj(i+1)​)(1−zj​)
其实也可以发现：
P ( Z = z j ∣ Y = y j , θ ( i ) ) = ( μ j ( i + 1 ) ) z j ( 1 − μ j ( i + 1 ) ) ( 1 − z j ) P(Z=z_j|Y=y_j,\,\theta^{(i)})=(\mu_j^{(i+1)})^{z_j}(1-\mu_j^{(i+1)})^{(1-z_j)} P(Z=zj​∣Y=yj​,θ(i))=(μj(i+1)​)zj​(1−μj(i+1)​)(1−zj​)
现在有了似然函数，也有了已知 Y Y Y 和 θ ( i ) \theta^{(i)} θ(i) 下出现 Z Z Z 的概率，我们现在可以求 Q Q Q 函数了：
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) = ∑ Z ∑ j = 1 n log ⁡ P ( y j , z j ∣ θ ) P ( Z ∣ Y , θ ( i ) ) = ∑ Z ∑ j = 1 n [ log ⁡ P ( y j , z j ∣ θ ) P ( Z ∣ Y , θ ( i ) ) ] = ∑ j = 1 n ∑ Z [ log ⁡ P ( y j , z j ∣ θ ) P ( Z ∣ Y , θ ( i ) ) ] （可以交换求和顺序） = ∑ j = 1 n ∑ z j ∑ z k , k ≠ j [ log ⁡ P ( y j , z j ∣ θ ) P ( Z ∣ Y , θ ( i ) ) ] \begin{aligned} &\, Q(\theta,\,\theta^{(i)}) \\ =&\, E_Z[\log P(Y,\,Z|\theta)|Y,\,\theta^{(i)}] \\ =&\, \sum_Z\log P(Y,\,Z|\theta)P(Z|Y,\,\theta^{(i)}) \\ =&\, \sum_Z \sum_{j=1}^{n}\log P(y_j,z_j|\theta) P(Z|Y,\,\theta^{(i)}) \\ =&\, \sum_Z \sum_{j=1}^{n} \left[ \log P(y_j,z_j|\theta) P(Z|Y,\,\theta^{(i)}) \right]\\ =&\, \sum_{j=1}^{n} \sum_Z \left[ \log P(y_j,z_j|\theta) P(Z|Y,\,\theta^{(i)}) \right] \quad\text{（可以交换求和顺序）}\\ =&\, \sum_{j=1}^{n} \sum_{z_j} \sum_{z_k,\,k\not=j} \left[ \log P(y_j,z_j|\theta) P(Z|Y,\,\theta^{(i)}) \right]\\ \end{aligned} ======​Q(θ,θ(i))EZ​[logP(Y,Z∣θ)∣Y,θ(i)]Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))Z∑​j=1∑n​logP(yj​,zj​∣θ)P(Z∣Y,θ(i))Z∑​j=1∑n​[logP(yj​,zj​∣θ)P(Z∣Y,θ(i))]j=1∑n​Z∑​[logP(yj​,zj​∣θ)P(Z∣Y,θ(i))]（可以交换求和顺序）j=1∑n​zj​∑​zk​,k=j∑​[logP(yj​,zj​∣θ)P(Z∣Y,θ(i))]​
发现当 z j z_j zj​ 的值固定时，有 ∑ z k , k ≠ j P ( Z ∣ Y , θ ( i ) ) = P ( z j ∣ y j , θ ( i ) ) \sum_{z_k,\,k\not=j} P(Z|Y,\,\theta^{(i)})=P(z_j|y_j,\,\theta^{(i)}) ∑zk​,k=j​P(Z∣Y,θ(i))=P(zj​∣yj​,θ(i)) ：
∑ z k , k ≠ j log ⁡ P ( y j , z j ∣ θ ) P ( Z ∣ Y , θ ( i ) ) = log ⁡ P ( y j , z j ∣ θ ) ∑ z k , k ≠ j P ( Z ∣ Y , θ ( i ) ) = P ( z j ∣ y j , θ ( i ) ) log ⁡ P ( y j , z j ∣ θ ) \begin{aligned} &\, \sum_{z_k,\,k\not=j} \log P(y_j,z_j|\theta) P(Z|Y,\,\theta^{(i)}) \\ =&\, \log P(y_j,z_j|\theta) \sum_{z_k,\,k\not=j} P(Z|Y,\,\theta^{(i)}) \\ =&\, P(z_j|y_j,\,\theta^{(i)})\log P(y_j,z_j|\theta) \end{aligned} ==​zk​,k=j∑​logP(yj​,zj​∣θ)P(Z∣Y,θ(i))logP(yj​,zj​∣θ)zk​,k=j∑​P(Z∣Y,θ(i))P(zj​∣yj​,θ(i))logP(yj​,zj​∣θ)​
因此：
Q ( θ , θ ( i ) ) = ∑ j = 1 n ∑ z j P ( z j ∣ y j , θ ( i ) ) log ⁡ P ( y j , z j ∣ θ ) = ∑ j = 1 n ( μ j ( i + 1 ) log ⁡ π p y j ( 1 − p ) ( 1 − y j ) + ( 1 − μ j ( i + 1 ) ) log ⁡ ( 1 − π ) q y j ( 1 − q ) ( 1 − y j ) ) \begin{aligned} &\, Q(\theta,\,\theta^{(i)}) \\ =&\, \sum_{j=1}^{n} \sum_{z_j}P(z_j|y_j,\,\theta^{(i)})\log P(y_j,z_j|\theta) \\ =&\, \sum_{j=1}^{n} \left( \mu_j^{(i+1)}\log \pi p^{y_j}(1-p)^{(1-y_j)} + (1-\mu_j^{(i+1)})\log (1-\pi)q^{y_j}(1-q)^{(1-y_j)} \right) \end{aligned} ==​Q(θ,θ(i))j=1∑n​zj​∑​P(zj​∣yj​,θ(i))logP(yj​,zj​∣θ)j=1∑n​(μj(i+1)​logπpyj​(1−p)(1−yj​)+(1−μj(i+1)​)log(1−π)qyj​(1−q)(1−yj​))​

接下来是 M 步，即求得 θ = ( π , p , q ) \theta=(\pi,\,p,\,q) θ=(π,p,q) 使得 Q ( θ , θ ( i ) ) Q(\theta,\,\theta^{(i)}) Q(θ,θ(i)) 极大化：

Q Q Q 对 π \pi π 求导：
∂ Q ∂ π = 1 π ∑ j = 1 n μ j ( i + 1 ) − 1 1 − π ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) = 0 ⇒ π ( i + 1 ) = 1 n ∑ j = 1 n μ j ( i + 1 ) \begin{aligned} \frac{\partial Q}{\partial \pi} =&\,\frac{1}{\pi}\sum_{j=1}^n\mu_{j}^{(i+1)}-\frac{1}{1-\pi}\sum_{j=1}^n(1-\mu_{j}^{(i+1)})=0 \\ \Rightarrow&\, \pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^{n}\mu_j^{(i+1)} \end{aligned} ∂π∂Q​=⇒​π1​j=1∑n​μj(i+1)​−1−π1​j=1∑n​(1−μj(i+1)​)=0π(i+1)=n1​j=1∑n​μj(i+1)​​
Q Q Q 对 p p p 求导：
∂ Q ∂ p = ∑ j = 1 n μ j ( i + 1 ) ( y j p − 1 − y j 1 − p ) = 0 ⇒ p ( i + 1 ) = ∑ j = 1 n μ j ( i + 1 ) y j ∑ j = 1 n μ j ( i + 1 ) \begin{aligned} \frac{\partial Q}{\partial p} =&\, \sum_{j=1}^{n}\mu_{j}^{(i+1)}\left( \frac{y_j}{p}-\frac{1-y_j}{1-p} \right) =0\\ \Rightarrow &\, p^{(i+1)}=\frac{\sum\limits_{j=1}^{n}\mu_j^{(i+1)}y_j}{\sum\limits_{j=1}^{n}\mu_j^{(i+1)}} \end{aligned} ∂p∂Q​=⇒​j=1∑n​μj(i+1)​(pyj​​−1−p1−yj​​)=0p(i+1)=j=1∑n​μj(i+1)​j=1∑n​μj(i+1)​yj​​​
Q Q Q 对 q q q 求导：
∂ Q ∂ q = ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) ( y j q − 1 − y j 1 − q ) = 0 ⇒ q ( i + 1 ) = ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) y j ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) \begin{aligned} \frac{\partial Q}{\partial q} =&\, \sum_{j=1}^{n}(1-\mu_{j}^{(i+1)})\left( \frac{y_j}{q}-\frac{1-y_j}{1-q} \right) =0\\ \Rightarrow &\, q^{(i+1)}=\frac{\sum\limits_{j=1}^{n}(1-\mu_j^{(i+1)})y_j}{\sum\limits_{j=1}^{n}(1-\mu_j^{(i+1)})} \end{aligned} ∂q∂Q​=⇒​j=1∑n​(1−μj(i+1)​)(qyj​​−1−q1−yj​​)=0q(i+1)=j=1∑n​(1−μj(i+1)​)j=1∑n​(1−μj(i+1)​)yj​​​