EM算法详解

从最大似然到EM算法

思考这样一个问题：假设我们要调查学校男生女生的体重分布，常见的作法是，在校园里抽样调查100个男生和100个女生的身高，这些男生身高（或女生身高）之间相互独立且都服从高斯分布：
$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$$

但我们不知道高斯分布的参数$(\mu ,{\sigma }^2)$，模型已知而参数未知让我们想到极大似然估计法，极大似然的原理说简单点就是在一次观测中发生了的事件其概率应该大，所以我们构造似然函数：
$$L(\theta )=\prod _{i=1}^{n}p(x_i;\mu ,{\sigma }^2),$$
为了方便分析我们定义对数似然函数：
$$H(\mu ,\sigma )=InL(\mu ,\sigma )=In\prod _{i=1}^{n}p(x_i;\mu ,{\sigma }^2)=\sum _{i=1}^{n}Inp(x_i;\mu ,{\sigma }^2)$$

要求 $\mu ,\sigma$，只需要使对数似然函数极大化，然后极大值对应的就 $\mu ,\sigma$是我们的估计。怎么求一个函数的最值？当然是求导，然后让导数为 0，解得 $\mu ,\sigma$。

补充：求极大似然估计值的一般步骤：

求最大似然函数估计值的一般步骤：

（1）写出似然函数；

（2）对似然函数取对数

（3）求导数，令导数为0，得到似然方程；

（4）解似然方程，得到的参数即为所求；

刚刚的背景是分男生、女生两个群体抽样，所以可以用似然估计法求出分布，现在假设我们我们随机抽样200人的身高，不去记录样本的性别。即我们只有一组身高数据，不知道这些样本来自哪个分布。这个时候样本之间独立不同分布，无法使用极大似然法。这个时候EM算法出场，先随便猜一下男生（身高）的正态分布的参数：如均值和方差是多少。例如男生的均值是1米7，方差是0.1米，然后计算出每个人更可能属于第一个还是第二个正态分布中的（如：某个人的身高是1米8，他最大可能属于男生的那个分布），这个是属于EM算法的第一步。有了每个人的归属，或者说我们已经大概地按上面的方法将这200个人分为男生和女生两部分，我们就可以根据之前说的最大似然那样，通过这些被大概分为男生的n个人来重新估计第一个分布的参数，女生的那个分布同样方法重新估计。这个是EM算法的第二步。然后，当我们更新了这两个分布的时候，每一个属于这两个分布的概率又变了，那么我们就再需要调整E步……如此往复，直到参数基本不再发生变化为止。这就是EM算法的基本思想。

总结一下：EM算法即期望(Expection)极大化(Maximization)算法，分两步走，E步：求期望，M步：求极大。用于含有隐变量的概率模型参数的极大似然估计。上述例子的隐变量是样本来自哪个群体未知。

Jensen不等式

主要讲解Jensen不等式，主要是如下的结论：

• $f(x)$是凸函数
  $$f(E(X))\le E(f(x))$$
  从函数的角度来看：若存在${\lambda }_i\le 0,{\sum }_i{\lambda }_i=1$，则有：
  $$f(\sum _{i=1}^M{\lambda }_i{x}_i)\le \sum _{i=1}^{M}f({\lambda }_i{x}_i)$$

简记：变量x期望的函数值小于等于函数的期望值

• $f(x)$是凹函数

$$f(E(X))\ge E(f(x))$$

简记：变量x期望的函数值大于等于函数的期望值

本文的重点在于EM算法，因此对Jensen不等式证明略，读者结合函数的凹凸性可轻易得出结论。

EM算法

这一部分的思路是先给出EM算法的基本流程，然后再证明这个算法的合理性。

算法流程

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

• (1)选择参数的初值${\theta }^0$，开始迭代

• (2) E步：记${\theta }^i$为第i次迭代参数$\theta$的估计值，在第i+1次迭代的E步，计算
  $$\begin{array}{rl}Q(\theta ,{\theta }^i)& =E_Z[logP(Y,Z|\theta )|Y,{\theta }^i]\\ & =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^i)\end{array}$$
  这里，$P(Z|Y,{\theta }^i)$是在给定观测数据Y和当前的参数估计${\theta }^i$下隐变量数据Z的条件概率分布；

• (3) M步：求使$Q(\theta ,{\theta }^i)$极大化的$\theta$，确定第i+1次迭代的参数的估计值${\theta }^{i+1}$，
  $${\theta }^{i+1}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)$$
  $Q(\theta ,{\theta }^i)$是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。

这里对Q函数做一些说明：完全数据的对数似然函数$logP(Y,Z|\theta )$ 关于在给定观测数据Y和当前参数${\theta }^{(i)}$ 下对未观测数据$Z$ 的条件概率分布$P(Z|Y,{\theta }^{(i)})$ 的期望称为Q函数。

• (4) 重复第(2)步和第(3)步，直到收敛，收敛条件：
  $$||{\theta }^{i+1}-{\theta }^i||<{\epsilon }_1$$
  或者：
  $$||Q({\theta }^{i+1},{\theta }^i)-Q({\theta }^i,{\theta }^i)||<{\epsilon }_2$$
  收敛迭代就结束了。我们来拆解一下这个M步骤，

EM算法的证明

为什么通过极大化这个$Q$ 函数对于EM算法来说是合理的？

下面我们从近似求解观测数据的对数似然函数的极大化问题来导出EM算法，由此可以清楚的看清EM算法的作用。 推导出Em算法可以近似实现对观测数据的极大似然估计的办法是找到E步骤的下界，让下界最大，通过逼近的方式实现对观测数据的最大似然估计。统计学习基础中采用的是相减方式，我们来看下具体的步骤。

（1）极大化观测数据$Y$关于参数 $\theta$ 的对数似然函数，并增加隐藏变量
$$\begin{array}{rl}L(\theta )& =logP(Y|\theta )=log\sum _{Z}P(Y,Z|\theta )\\ & =log\sum _{Z}P(Y|Z,\theta )P(Z,\theta )\end{array}$$
似然函数中有未观测数据，无法直接计算，因此EM算法通过迭代逐步实现最大化，假设第$i$ 次迭代后$\theta$ 的估计值是${\theta }^{(i)}$我们的目标是更新$\theta$，使得$L(\theta )$ 增加，即$L(\theta )>L({\theta }^i)$,并逐步达到最大值。我们考虑两者之差为：
$$\begin{array}{rl}L(\theta )-L({\theta }^{(i)})& =log(\sum _{Z}P(Y|Z,\theta )P(Z,\theta )-logP(Y|{\theta }^{(i)})\\ & =log(\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Z|Y,{\theta }^{(i)})})-logP(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)})})-logP(Y|{\theta }^{(i)})(1)\\ & =\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)})logP(Y|{\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)})}\end{array}$$
(1) 处 $\ge$这一个步骤就是采用了凹函数的Jensen不等式做转换。

令：
$$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 )$ 的一个下界，取$\theta ={\theta }^{(i)}$ ,根据链式法则：
$$P(Y,Z|\theta )=P(Y|Z,\theta )P(Z,\theta )=P(Z|Y,\theta )P(Y|\theta )$$
则有：
$$\begin{array}{rl}B({\theta }^{(i)},{\theta }^{(i)})& =L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,{\theta }^{(i)})P(Z,{\theta }^{(i)})}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\\ & =L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|{\theta }^{(i)})}{P(Y,Z|{\theta }^{(i)})}\\ & =L({\theta }^{(i)})\end{array}$$
所以，任何可以使$B(\theta ,{\theta }^{(i)})$ 增大的$\theta$ 都可以使$L(\theta )$ 增大。所以为了使$L(\theta )$ 尽可能增大，我们选择${\theta }^{(i+1)}$ 使得$B(\theta ,{\theta }^i)$ 尽可能大。

即最大化如下：

$$\begin{array}{rl}{\theta }^{(i+1)}& =\arg \underset{\theta }{\max }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)})}\\ & =\arg \underset{\theta }{\max }\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)})})\\ & =\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z,\theta ))-\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Z|Y,{\theta }^i)P(Y|{\theta }^{(i)}))\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z,\theta ))\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y,Z|\theta ))\\ & =arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}$$
所以EM算法的一-次迭代， 即求Q函数及其极大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
下图给出EM算法的直观解释。图中上方曲线为$L(\theta )$,下方曲线为$B(\theta ,{\theta }^{(i)})$。$B(\theta ,{\theta }^{(i)})$为对数似然函数$L(\theta )$的下界。

通过图片我们总结ＥＭ算法的迭代过程如下：

(1) 初始状态：在$\theta ={\theta }^{(i)}$ 时，$L({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})$

(2) EM算法的更新过程：寻找${\theta }^{(i+1))}$ 使得$B(\theta ,{\theta }^{(i)})$极大化。因为$L(\theta )\ge B(\theta ,{\theta }^{(i)})$，所以得$B(\theta ,{\theta }^{(i)})$极大化也使得$L(\theta )$ 增加。

(3) 根据第二步找到的${\theta }^{(i+1))}$ 重新计算$B(\theta ,{\theta }^{(i)})$（或者说重新计算Q函数），进行下一次迭代。

上述过程暴露一个问题：对数似然函数不断增大，但EM算法不能保证找到全局最优解。

EM算法的收敛性

我们已经证明了EM算法的合理性，接下来我们来研究一下算法是否收敛。

要证EM算法的收敛性即证:
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$
因为：
$$P(Y|\theta )=\frac{P(Z,Y|\theta )}{P(Z|Y,\theta )}$$
两边同取对数：
$$logP(Y|\theta )=logP(Z,Y|\theta )-logP(Z|Y,\theta )$$
左右两边同时求期望，则有：
$$\begin{array}{rl}\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|\theta )& =logP(Y|\theta )\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z,Y|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )\end{array}$$
其中,${\sum }_{Z}P(Z|Y,{\theta }^{(i)})=1$,$logP(Y|\theta )$ 是与Z无关的变量，因此积分等于它本身。

根据Q函数的定义：
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
以及这里我们令
$$H(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )$$
所以：
$$\begin{array}{rl}logP(Y|\theta )& =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z,Y|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )\\ & =Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})\end{array}$$

$$\begin{array}{rl}logP(Y,{\theta }^{(i+1)})-logP(Y,{\theta }^{(i)})& =Q({\theta }^{(i+1)},{\theta }^{(i)}-H({\theta }^{(i+1)},{\theta }^{(i)})-(Q({\theta }^{(i)},{\theta }^{(i)})-H({\theta }^{(i)},{\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$$
所以，我们的目标是证明$H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})\le 0$ ，证明如下：

$$\begin{array}{rl}H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})& =\sum _Z(log(\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}))P(Z|Y,{\theta }^{(i)})\\ & \le log(\sum _Z\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)}))\\ & =logP(Z|Y,{\theta }^{(i+1)})=log1=0\end{array}$$
不等号由jensen不等式得到。

所以：
$$\begin{gathered} logP(Y,{\theta }^{(i+1)})-logP(Y,{\theta }^{(i)})\ge 0 \\ P(Y,{\theta }^{(i+1)})-P(Y,{\theta }^{(i)})\ge 0 \\ P(Y,{\theta }^{(i+1)})\ge P(Y,{\theta }^{(i)}) \end{gathered}$$
证明了EM算法的收敛性。但不能保证是全局最优，只能保证局部最优。

三硬币模型的EM解答

再来看一个例子，假设有3枚硬币，分别记作A，B和C。 这些硬币正面向上的概率分别是 $\pi$、$p$和$q$ 。进行如下抛硬币试验：先抛硬币A, 根据其结果选出硬币B或者硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，抛硬币的结果，出现正面记作1，出现反面记作0；独立重复10次试验。结果如下：

1，1，0，1，0，0，1，0，1，1

假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率。

我们先来构建这样一个三硬币模型：
$$\begin{array}{rl}P(y|\theta )& =\sum _{z}P(y,z|\theta )=\sum _{z}P(y|z,\theta )P(z|\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$

• 若$y=1$，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则$P(1|\theta )=\pi p+(1-\pi )q$
• 若$y=0$，则$P(0|\theta )=\pi (1-p)+(1-\pi )(1-q)$

y是观测变量，表示一次观测结果是1或0，z是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，$\theta =(\pi ,p,q)$表示模型参数，将观测数据表示为$Y=(Y_1,Y_2,...,Y_n{)}^T$，未观测的数据表示为$Z=(Z_1,Z_2,...,Z_n{)}^T$，则观测函数的似然函数是：
$$\begin{array}{rl}P(Y|\theta )& =\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\\ & =\prod _{i=0}(\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i})\end{array}$$
考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即：
$$\hat{\theta }=\arg \underset{\theta }{\max }logP(Y|\theta )=\arg \underset{\theta }{\max }\sum _{i=1}^{10}log(\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}$$
我们发现直接对上式求偏导等于0很难求得参数的值。

这个问题含有隐变量(不知道观测样本值是来自哪个硬币，B还是C)，只有通过迭代方法来求解，我们借助EM算法进行求解。

计算在模型参数${\pi }^i，p^i，q^i$下观测变量$y_i$来源于硬币B的概率：
$${\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-p^i{)}^{1-y_i}}$$
注：这个公式的分母是$P(Y|\theta )$，表示观测变量$y_i$来源于B硬币和C硬币的概率和。分子表示观测变量$y_i$来源于B硬币的概率。

（1）E步：一共n个样本，构造Q函数：
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =\sum _{j=1}^n\sum _{Z}P(Z|y_j,{\theta }^{(i)})\log P(y_j,Z|\theta )\\ & =\sum _{j=1}^n\left\{P(Z=1|y_j,{\theta }^{(i)})\mathrm{l}\mathrm{o}\mathrm{g}P(y_j,Z=1|\theta )+P(Z=0|y_j,{\theta }^{(i)})\mathrm{l}\mathrm{o}\mathrm{g}P(y_j,Z=0|\theta )\right\}\\ & =\sum _{j=1}^n\left\{{\mu }_j^{(i+1)}logP(y_j,Z=1|\theta )+(1-{\mu }_j^{(i+1)})\mathrm{l}\mathrm{o}\mathrm{g}P(y_j,Z=0|\theta )\right\}\\ & =\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}$$
（2）M步：极大化Q函数
$$\begin{array}{r}{\theta }^{(i+1)}={\mathrm{argmax}}_{\theta }Q(\theta ,{\theta }^{(i)})={\mathrm{argmax}}_{\theta }\left\{\sum _{j=1}^n\sum _{Z}P(Z|y_j,{\theta }^{(i)})\log P(y_j,Z|\theta )\right\}\end{array}$$
极大值问题：即令导函数为0，即：
∂ Q ( θ , θ ( i ) ) ∂ π = ∑ j = 1 N { μ j ( i + 1 ) ln ⁡ [ π p y j ( 1 − p ) 1 − y j ] + ( 1 − μ j ( i + 1 ) ) ln ⁡ [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] ∂ π } = ∑ j = 1 N { μ j ( i + 1 ) p y j ( 1 − p ) 1 − y j π p y j ( 1 − p ) 1 − y j + ( 1 − μ j ( i + 1 ) ) − q y j ( 1 − q ) 1 − y j ( 1 − π ) q y j ( 1 − q ) 1 − y j } = ∑ j = 1 N { μ j ( i + 1 ) − π π ( 1 − π ) } = ( ∑ j = 1 N μ j ( i + 1 ) ) − n π π ( 1 − π ) \begin{aligned} \frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi}&=\sum_{j=1}^N\{\frac{\mu_{j}^{(i+1)}\ln [\pi p^{y_j}(1-p)^{1-y_j}]+(1-\mu_{j}^{(i+1)})\ln [(1-\pi )q^{y_j}(1-q)^{1-y_j}] }{\partial \pi}\}\\&=\sum_{j=1}^N\{ \mu_{j}^{(i+1)}\frac{p^{y_j}(1-p)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j}}+(1-\mu_{j}^{(i+1)})\frac{-q^{y_j}(1-q)^{1-y_j}}{(1-\pi )q^{y_j}(1-q)^{1-y_j}} \}\\&=\sum_{j=1}^N\{ \frac{\mu_{j}^{(i+1)}-\pi }{\pi (1-\pi)}\}\\&=\frac{(\sum_{j=1}^N\mu_{j}^{(i+1)})-n\pi }{\pi (1-\pi)} \end{aligned} ∂π∂Q(θ,θ(i))​​=j=1∑N​{∂πμj(i+1)​ln[πpyj​(1−p)1−yj​]+(1−μj(i+1)​)ln[(1−π)qyj​(1−q)1−yj​]​}=j=1∑N​{μj(i+1)​πpyj​(1−p)1−yj​pyj​(1−p)1−yj​​+(1−μj(i+1)​)(1−π)qyj​(1−q)1−yj​−qyj​(1−q)1−yj​​}=j=1∑N​{π(1−π)μj(i+1)​−π​}=π(1−π)(∑j=1N​μj(i+1)​)−nπ​​

$$\begin{array}{rl}\because \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial \pi }=0& ⟹\pi =\frac{1}{n}\sum _{j=1}^N{\mu }_j^{(i+1)}\\ \therefore {\pi }^{(i+1)}& =\frac{1}{n}\sum _{j=1}^N{\mu }_j^{(i+1)}\end{array}$$

∂ Q ( θ , θ ( i ) ) ∂ p = ∑ j = 1 N { μ j ( i + 1 ) ln ⁡ [ π p y j ( 1 − p ) 1 − y j ] + ( 1 − μ j ( i + 1 ) ) ln ⁡ [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] ∂ p } = ∑ j = 1 N { μ j ( i + 1 ) π ( y j p y j − 1 ( 1 − p ) 1 − y j + p y j ( − 1 ) ( 1 − y j ) ( 1 − p ) 1 − y j − 1 ) π p y j ( 1 − p ) 1 − y j + 0 } = ∑ j = 1 N { μ j ( i + 1 ) ( y j − p ) p ( 1 − p ) } = ( ∑ j = 1 N μ j ( i + 1 ) y j ) − ( p ∑ j = 1 N μ j ( i + 1 ) ) p ( 1 − p ) \begin{aligned} \frac{\partial Q(\theta,\theta^{(i)})}{\partial p}&=\sum_{j=1}^N\{\frac{\mu_{j}^{(i+1)}\ln [\pi p^{y_j}(1-p)^{1-y_j}]+(1-\mu_{j}^{(i+1)})\ln [(1-\pi )q^{y_j}(1-q)^{1-y_j}] }{\partial p}\}\\&=\sum_{j=1}^N\{\mu_{j}^{(i+1)}\frac{\pi (y_jp^{y_j-1}(1-p)^{1-y_j}+p^{y_j}(-1)(1-y_j)(1-p)^{1-y_j-1})}{\pi p^{y_j}(1-p)^{1-y_j}}+0 \}\\&=\sum_{j=1}^N\{ \frac{\mu_{j}^{(i+1)}(y_j-p) }{p(1-p)}\}\\&=\frac{(\sum_{j=1}^N\mu_{j}^{(i+1)}y_j)-(p\sum_{j=1}^N\mu_{j}^{(i+1)}) }{p(1-p)} \end{aligned} ∂p∂Q(θ,θ(i))​​=j=1∑N​{∂pμj(i+1)​ln[πpyj​(1−p)1−yj​]+(1−μj(i+1)​)ln[(1−π)qyj​(1−q)1−yj​]​}=j=1∑N​{μj(i+1)​πpyj​(1−p)1−yj​π(yj​pyj​−1(1−p)1−yj​+pyj​(−1)(1−yj​)(1−p)1−yj​−1)​+0}=j=1∑N​{p(1−p)μj(i+1)​(yj​−p)​}=p(1−p)(∑j=1N​μj(i+1)​yj​)−(p∑j=1N​μj(i+1)​)​​

$$\begin{array}{rl}\because \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial p}=0& ⟹p=\frac{{\sum }_{j=1}^N{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^N{\mu }_j^{(i+1)}}\\ \therefore 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$，因为对$y_i=1$和$y_i=0$均有${\mu }_j^1=0.5$，利用迭代公式计算得到${\pi }^1=0.5，p^1=0.6，q^1=0.6$，继续迭代得到最终的参数：
$$\widehat{{\pi }^0}=0.5，\widehat{p^0}=0.6，\widehat{q^0}=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算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。