统计学习方法 EM 算法

统计学习方法 EM 算法

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

引入

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

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

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

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

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

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

解：三硬币模型可以写作：
$$\begin{array}{rl}P(y|\theta )=& \sum _{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{array}$$

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

以向量形式表示 $n$ 次实验结果，观测数据 $Y=(Y_1,Y_2,\cdots ,Y_n{)}^T$ ，隐变量 $Z=(Z_1,Z_2,\cdots ,Z_n{)}^T$ ；

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

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

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

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

$${\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}}$$

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

$$\begin{array}{rl}{\pi }^{(i+1)}=& \frac{1}{n}\sum _{j=1}^n{\mu }_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)}}\\ q^{(i+1)}=& \frac{\sum _{j=1}^n(1-{\mu }_j{)}^{(i+1)}{y}_j}{\sum _{j=1}^n(1-{\mu }_j{)}^{(i+1)}}\end{array}$$

例如：

• 初始值 ${\pi }^{(0)}=0.5$ ，$p^{(0)}=0.5$ ，$q^{(0)}=0.5$ ；
• 第 1 轮，${\mu }_j^{(1)}=0.5$ （$j=1,2,\cdots ,10$） ，得到 ${\pi }^{(1)}=0.5$ ，$p^{(1)}=0.6$ ，$q^{(1)}=0.6$ ；
• 第 2 轮，${\mu }_j^{(2)}=0.5$ （$j=1,2,\cdots ,10$） ，得到 ${\pi }^{(2)}=0.5$ ，$p^{(2)}=0.6$ ，$q^{(2)}=0.6$ ；

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

EM 算法

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

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

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

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

② E 步：记 ${\theta }^{(i)}$ 为第 $i$ 次迭代时参数 $\theta$ 的估计值，则第 $i+1$ 次迭代时，计算：
$$\begin{array}{rl}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{array}$$
这里 $Y$ 和 ${\theta }^{(i)}$ 是已知的， $Q$ 函数是关于 $\theta$ 的函数，相当于在 $P(Z|Y,{\theta }^{(i)})$ 的概率测度下，计算以 $\theta$ 为参数的模型中 $\log P(Y,Z|\theta )$ 的期望（有一点点绕）。

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

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

EM 算法的导出

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

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

$$\log \sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_j\log y_i\text{其中 }{\lambda }_j\ge 0,\sum _j{\lambda }_j=1$$

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

记一个新的函数为：
$$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)})}$$
则：
$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$
即 $B(\theta ,{\theta }^{(i)})$ 是 $L(\theta )$ 的一个下界。并且对于 $B({\theta }^{(i)},{\theta }^{(i)})$ ，有：（这是条件概率的公式）
$$\{\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}$$
故 $B({\theta }^{(i)},{\theta }^{(i)})=L({\theta }^{(i)})$ 。因此，任何可以使得 $B(\theta ,{\theta }^{(i)})$ 增大的 $\theta$ ，都可以使得 $L(\theta )$ 增大，我们选择 $\theta$ 使得 $B(\theta ,{\theta }^{(i)})$ 达到极大，即：
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})$$
有：
$$\begin{array}{rl}{\theta }^{(i+1)}=& \arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})\\ =& \arg \underset{\theta }{\max }\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 \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y|Z,\theta )P(Z|\theta )\right)\\ =& \arg \underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )\right)\\ =& \arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}$$

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

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

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

EM 算法的收敛性

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

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

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

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

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

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

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

三硬币模型

现在重新来看三硬币模型，这个模型的隐变量是硬币 A 的结果 $Z=(z_1,z_2,\cdots ,z_n)$ 。

其先求其对数似然函数：
$$\log P(Y,Z|\theta )=\log [P(Y|Z,\theta )P(Z|\theta )]$$
按照条件概率的公式，其实也可以写成：
$$\log P(Y,Z|\theta )=\log [P(Z|Y,\theta )P(Y|\theta )]$$
但是你想想抛硬币的过程，显然第一种是顺着抛硬币的过程写的，这样比较好算；有：
$$P(Z|\theta )=\prod _{j=1}^n{\pi }^{z_j}(1-\pi {)}^{(1-z_j)}$$
下面这个式子可能比较难理解，也是让我想了很久，可以按照 $z_i$ 为 0 或者为 1 来分类讨论：
$$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|\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=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(Z|Y,\theta )$ （因为需要 $P(Z|Y,{\theta }^{(i)})$ 来算期望），按照条件概率的公式有：
$$P(Z|Y,\theta )=\frac{P(Y,Z|\theta )}{P(Y|\theta )}$$
我们现在有 $P(Y,Z|\theta )$ 了，可以求一下 $P(Y|\theta )$ ，有：
$$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,Z|\theta )$ 求和得到。比较难理解，有一点像二项式展开，对每个 $z_j$ 取 0 和 1 的情况遍历，${\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)}$ 就像是从每一个括号中选出其中一个来乘：
$$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(Z|Y,\theta )$ 了：
$$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}}$$
我们引入记号，在模型参数 ${\theta }^{(i)}$ 下，观测数据 $y_j$ 来自 B 的概率为：（有点像贝叶斯公式）
$${\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}}$$

则：
$$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=z_j|Y=y_j,{\theta }^{(i)})=({\mu }_j^{(i+1)}{)}^{z_j}(1-{\mu }_j^{(i+1)}{)}^{(1-z_j)}$$
现在有了似然函数，也有了已知 $Y$ 和 ${\theta }^{(i)}$ 下出现 $Z$ 的概率，我们现在可以求 $Q$ 函数了：
$$\begin{array}{rl}& 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]\text{（可以交换求和顺序）}\\ =& \sum _{j=1}^n\sum _{z_j}\sum _{z_k,k\ne j}\left[\log P(y_j,z_j|\theta )P(Z|Y,{\theta }^{(i)})\right]\end{array}$$
发现当 $z_j$ 的值固定时，有 ${\sum }_{z_k,k\ne j}P(Z|Y,{\theta }^{(i)})=P(z_j|y_j,{\theta }^{(i)})$ ：
$$\begin{array}{rl}& \sum _{z_k,k\ne j}\log P(y_j,z_j|\theta )P(Z|Y,{\theta }^{(i)})\\ =& \log P(y_j,z_j|\theta )\sum _{z_k,k\ne j}P(Z|Y,{\theta }^{(i)})\\ =& P(z_j|y_j,{\theta }^{(i)})\log P(y_j,z_j|\theta )\end{array}$$
因此：
$$\begin{array}{rl}& 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{array}$$

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

$Q$ 对 $\pi$ 求导：
$$\begin{array}{rl}\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{array}$$
$Q$ 对 $p$ 求导：
$$\begin{array}{rl}\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 _{j=1}^n{\mu }_j^{(i+1)}{y}_j}{\sum _{j=1}^n{\mu }_j^{(i+1)}}\end{array}$$
$Q$ 对 $q$ 求导：
$$\begin{array}{rl}\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 _{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})}\end{array}$$