统计学习方法读书笔记第九章：EM算法及其推广

EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximization）。所以这一算法称为期望极大算法，简称EM算法。

EM算法的引入

概率模型有时既含有观测变量，又含有隐变量或潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

• EM算法
  将观测数据表示为$Y=(Y_1,Y_2,\cdots,Y_n)^T$Y=(Y1​,Y2​,⋯,Yn​)T，未观测数据表示为$Z=(Z_1,Z_2,\cdots,Z_n)^T$Z=(Z1​,Z2​,⋯,Zn​)T，则观测数据的似然函数为
  $P(Y|\theta)=\sum_ZP(Z|\theta)P(Y|Z,\theta) \tag{1}$P(Y∣θ)=Z∑​P(Z∣θ)P(Y∣Z,θ)(1)
  考虑求模型参数$\theta=(\pi,p,q)$θ=(π,p,q)的极大似然估计，即
  $\hat\theta=arg\max_\theta logP(Y|\theta) \tag{2}$θ^=argθmax​logP(Y∣θ)(2)
  这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。
  一般地，用$Y$Y表示观测随机变量的数据，$Z$Z表示隐随机变量的数据。$Y$Y和$Z$Z连在一起称为完全数据，观测数据$Y$Y又称为不完全数据。假设给定观测数据$Y$Y，其概率分布是$P(Y|\theta)$P(Y∣θ)，其中$\theta$θ是需要估计的模型参数，那么不完全数据$Y$Y的似然函数是$P(Y|\theta)$P(Y∣θ)，对数似然函数$L(\theta)=logP(Y|\theta)$L(θ)=logP(Y∣θ)；假设$Y$Y和$Z$Z的联合概率分布是$P(Y,Z|\theta)$P(Y,Z∣θ)，那么完全数据的对数似然函数是$logP(Y,Z|\theta)$logP(Y,Z∣θ)。
  EM算法通过迭代求$L(\theta)=logP(Y|\theta)$L(θ)=logP(Y∣θ)的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。下面来介绍EM算法。
  EM算法
  输入：观测变量数据$Y$Y，隐变量数据$Z$Z，联合分布$P(Y,Z|\theta)$P(Y,Z∣θ)，条件分布$P(Z|Y,\theta)$P(Z∣Y,θ)；
  输出：模型参数$\theta$θ。
  (1) 选择参数的初值$\theta^{(0)}$θ(0)，开始迭代；
  (2) E步：记$\theta^{(i)}$θ(i)为第$i$i次迭代参数$\theta$θ的估计值，在第$i+1$i+1次迭代的E步，计算
  $\begin{aligned} Q(\theta,\theta^{(i)})&=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}] \\ &=\sum_ZlogP(Y,Z|\theta)P(Z|Y,\theta^{(i)}) \tag{3} \end{aligned}$Q(θ,θ(i))​=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))​(3)
  这里，$P(Z|Y,\theta^{(i)})$P(Z∣Y,θ(i))实在给定观测数据$Y$Y和当前的参数估计$\theta^{(i)}$θ(i)下隐变量数据$Z$Z的条件概率分布；
  (3) M步：求使$Q(\theta,\theta^{(i)})$Q(θ,θ(i))极大化的$\theta$θ，确定第$i+1$i+1次迭代的参数的估计值$\theta^{(i+1)}$θ(i+1)
  $\theta^{(i+1)}=arg\max_\theta Q(\theta,\theta^{(i)}) \tag{4}$θ(i+1)=argθmax​Q(θ,θ(i))(4)
  (4) 重复第(2)步和第(3)步，直到收敛。
  第(2)步中的函数$Q(\theta,\theta^{(i)})$Q(θ,θ(i))是EM算法的核心，称为$Q$Q函数。
  定义1（Q函数） 完全数据的对数似然函数$logP(Y,Z|\theta)$logP(Y,Z∣θ)关于在给定观测数据$Y$Y和当前参数$\theta^{(i)}$θ(i)下对未观测数据$Z$Z的条件概率分布$P(Z|Y,\theta^{(i)})$P(Z∣Y,θ(i))的期望称为$Q$Q函数，即
  $Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}] \tag{5}$Q(θ,θ(i))=EZ​[logP(Y,Z∣θ)∣Y,θ(i)](5)
  下面关于EM算法作几点说明：
  步骤(1) 参数的初值可以任意选择，但需注意EM算法对初值是敏感的。
  步骤(2) E步求$Q(\theta,\theta^{(i)})$Q(θ,θ(i))。$Q$Q函数式中$Z$Z是未观测数据，$Y$Y是观测数据。注意，$Q(\theta,\theta^{(i)})$Q(θ,θ(i))的第１个变元表示要极大化的参数，第２个变元表示参数的当前估计值。每次迭代实际在求$Q$Q函数及其极大。
  步骤(3) M步求Q(\theta,\theta^{(i)})的极大化，得到$\theta^{(i+1)}$θ(i+1)，完成一次迭代$\theta^{(i)}\to\theta^{(i+1)}$θ(i)→θ(i+1)。后面将证明每次迭代使似然函数增大或达到局部极值。
  步骤(4) 给出停止迭代的条件，一般是对较小的正数$\varepsilon_1$ε1​，$\varepsilon_2$ε2​，若满足
  $||\theta^{(i+1)}-\theta^{(i)}||&lt;\varepsilon_1 或||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})||&lt;\varepsilon_2$∣∣θ(i+1)−θ(i)∣∣<ε1​或∣∣Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∣∣<ε2​
  则停止迭代。
• EM算法的导出
  上面叙述了EM算法。为什么EM算法能近似实现对观测数据的极大似然估计呢？下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法，由此可以清楚地看出EM算法的作用。
  我们面对一个含有隐含量的概率模型，目标是极大化观测数据（不完全数据）$Y$Y关于参数$\theta$θ的对数似然函数，即极大化
  $\begin{aligned} L(\theta)&amp;=logP(Y|\theta)=log\sum_ZP(Y,Z|\theta) \\ &amp;=log\bigg(\sum_ZP(Y|Z,\theta)P(Z|\theta)\bigg) \tag{6} \end{aligned}$L(θ)​=logP(Y∣θ)=logZ∑​P(Y,Z∣θ)=log(Z∑​P(Y∣Z,θ)P(Z∣θ))​(6)
  注意到这一极大化的主要困难是上式中有未观测数据并有包含和（或积分）的对数。
  事实上，EM算法是通过迭代逐步近似极大化$L(\theta)$L(θ)的。假设在第$i$i次迭代后$\theta$θ的估计值是$\theta^{(i)}$θ(i)。我们希望新估计值$\theta$θ能使$L(\theta)$L(θ)增加，即$L(\theta)&gt;L(\theta^{(i)})$L(θ)>L(θ(i))，并逐步达到极大值。为此，考虑两者的差：
  $L(\theta)-L(\theta^{(i)})=log\bigg(\sum_ZP(Y|Z,\theta)P(Z|\theta)\bigg)-logP(Y|\theta^{(i)})$L(θ)−L(θ(i))=log(Z∑​P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))
  利用Jensen不等式得到其下界：
  $\begin{aligned} L(\theta)-L(\theta^{(i)})&amp;=log\bigg(\sum_ZP(Y|Z,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta^{(i)})}\bigg)-logP(Y|\theta^{(i)}) \\ &amp;\geq\sum_ZP(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-logP(Y|\theta^{(i)}) \\ &amp;=\sum_ZP(Z|Y,\theta^{(i)})log\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(Y∣Z,θ(i))P(Y∣Z,θ(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∣θ(i))P(Y∣Z,θ)P(Z∣θ)​​
  令
  $B(\theta,\theta^{(i)})\hat{=}L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)} {P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \tag{7}$B(θ,θ(i))=^L(θ(i))+Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)​(7)
  则
  $L(\theta)\geq B(\theta,\theta^{(i)}) \tag{8}$L(θ)≥B(θ,θ(i))(8)
  即函数$B(\theta,\theta^{(i)})$B(θ,θ(i))是$L(\theta)$L(θ)的一个下界，且由上式可知，
  $L(\theta^{(i)})=B(\theta^{(i)},\theta^{(i)}) \tag{9}$L(θ(i))=B(θ(i),θ(i))(9)
  因此，任何可以使$B(\theta,\theta^{(i)})$B(θ,θ(i))增大的$\theta$θ，也可以使$L(\theta)$L(θ)增大。为了使$L(\theta)$L(θ)有尽可能大的增长，选择$\theta^{(i+1)}$θ(i+1)使$B(\theta,\theta^{(i)})$B(θ,θ(i))达到极大，即
  $\theta^{(i+1)}=arg\max_\theta B(\theta,\theta^{(i)}) \tag{10}$θ(i+1)=argθmax​B(θ,θ(i))(10)
  现在求$\theta^{(i+1)}$θ(i+1)的表达式。省去对$\theta$θ的极大化而言是常数的项，则有
  $\begin{aligned} \theta^{(i+1)}&amp;=arg\max_\theta\bigg(L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\bigg) \\ &amp;=arg\max_\theta\bigg(\sum_ZP(Z|Y,\theta^{(i)})log(P(Y|Z,\theta)P(Z|\theta))\bigg) \\ &amp;=arg\max_\theta\bigg(\sum_ZP(Z|Y,\theta^{(i)})logP(Y,Z|\theta)\bigg) \\ &amp;=arg\max_\theta Q(\theta,\theta^{(i)}) \tag{11} \end{aligned}$θ(i+1)​=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))log(P(Y∣Z,θ)P(Z∣θ)))=argθmax​(Z∑​P(Z∣Y,θ(i))logP(Y,Z∣θ))=argθmax​Q(θ,θ(i))​(11)
  上式等价于EM算法的一次迭代，即求$Q$Q函数及其极大化。$EM$EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
  下图给出EM算法的直观解释。途中上方曲线为$L(\theta)$L(θ)，下方曲线为$B(\theta,\theta^{(i)})$B(θ,θ(i))。$B(\theta,\theta^{(i)})$B(θ,θ(i))为对数似然函数$L(\theta)$L(θ)的下界。同时，两个函数在点$\theta=\theta^{(i)}$θ=θ(i)处相等。则EM算法找到下一个点$\theta^{(i+1)}$θ(i+1)使函数$B(\theta,\theta^{(i)})$B(θ,θ(i))极大化，也使函数$Q(\theta,\theta^{(i)})$Q(θ,θ(i))极大化。这时由于$L(\theta)\geq B(\theta,\theta^{(i)})$L(θ)≥B(θ,θ(i))，函数$B(\theta,\theta^{(i)})$B(θ,θ(i))的增加，保证对数似然函数$L(\theta)$L(θ)在每次迭代中也是增加的。EM算法在点$\theta^{(i+1)}$θ(i+1)重新计算$Q$Q函数值，进行下一次迭代。在这个过程中，对数似然函数$L(\theta)$L(θ)不断增大。从图可以推断出EM算法不能保证找到全局最优解。
  
• EM算法在非监督学习中的应用
  监督学习是由训练数据$\{(x_1,y_1),(x_2,y_2),\cdots,(x_N,y_N)\}${(x1​,y1​),(x2​,y2​),⋯,(xN​,yN​)}学习条件概率分布$P(Y|X)$P(Y∣X)或决策函数$Y=f(X)$Y=f(X)作为模型，用于分类、回归、标注等任务。这时训练数据中的每个样本点由输入和输出对组成。
  有时训练数据只有输入没有对应的输出$\{(x_1,\cdot),(x_2,\cdot),\cdots,(x_N,\cdot)\}${(x1​,⋅),(x2​,⋅),⋯,(xN​,⋅)}，从这样的数据学习模型称为非监督学习问题。EM算法可以用于生成模型的非监督学习。生成模型由联合概率分布$P(X,Y)$P(X,Y)表示，可以认为非监督学习训练数据是联合概率分布产生的数据。$X$X为观测数据，$Y$Y为未观测数据。

EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。EM算法的最大优点是简单性和普适性。我们很自然地要问：EM算法得到的估计序列是否收敛？如果收敛，是否收敛到全局最大值或局部最大值？下面给出关于EM算法收敛性的两个定理。

• 定理１ 设$P(Y|\theta)$P(Y∣θ)为观测数据的似然函数，$\theta^{(i)}(i=1,2,\cdots)$θ(i)(i=1,2,⋯)为EM算法得到的参数估计序列，$P(Y|\theta^{(i)})(i=1,2,\cdots)$P(Y∣θ(i))(i=1,2,⋯)为对应的似然函数序列，则$P(Y|\theta^{(i)})$P(Y∣θ(i))是单调递增的，即
  $P(Y|\theta^{(i+1)})\geq P(Y|\theta^{(i)}) \tag{12}$P(Y∣θ(i+1))≥P(Y∣θ(i))(12)
  证明 由于
  $P(Y|\theta)=\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}$P(Y∣θ)=P(Z∣Y,θ)P(Y,Z∣θ)​
  取对数有
  $logP(Y|\theta)=logP(Y,Z|\theta)-logP(Z|Y,\theta)$logP(Y∣θ)=logP(Y,Z∣θ)−logP(Z∣Y,θ)
  由Q函数
  $Q(\theta,\theta^{(i)})=\sum_ZlogP(Y,Z|\theta)P(Z|Y,\theta^{(i)})$Q(θ,θ(i))=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))
  令
  $H(\theta,\theta^{(i)})=\sum_ZlogP(Z|Y,\theta)P(Z|Y,\theta^{(i)}) \tag{13}$H(θ,θ(i))=Z∑​logP(Z∣Y,θ)P(Z∣Y,θ(i))(13)
  于是对数似然函数可以写成
  $logP(Y|\theta)=Q(\theta,\theta^{(i)})-H(\theta,\theta^{(i)}) \tag{14}$logP(Y∣θ)=Q(θ,θ(i))−H(θ,θ(i))(14)
  在上式中分别取$\theta$θ为$\theta^{(i)}$θ(i)和$\theta^{i+1}$θi+1并相减，有
  $\begin{aligned} &amp;logP(Y|\theta^{(i+1)})-logP(Y|\theta^{(i)}) \\ &amp;=[Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})]-[H(\theta^{(i+1)})] \tag{15} \end{aligned}$​logP(Y∣θ(i+1))−logP(Y∣θ(i))=[Q(θ(i+1),θ(i))−Q(θ(i),θ(i))]−[H(θ(i+1))]​(15)
  为证式(12)，只需证式(15)右端是非负的。式(15)右端的第１项，由于$\theta^{(i+1)}$θ(i+1)使$Q(\theta,\theta^{(i)})$Q(θ,θ(i))达到极大，所以有
  $Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})\geq 0 \tag{16}$Q(θ(i+1),θ(i))−Q(θ(i),θ(i))≥0(16)
  其第２项，由式(13)可得：
  $\begin{aligned} H(\theta^{(i+1)},&amp;\theta^{(i)})-H(\theta^{(i)},\theta^{(i)}) \\ &amp;=\sum_Z\bigg(log\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}\bigg)P(Z|Y,\theta^{(i)}) \\ &amp;\leq log\bigg(\sum_Z\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}P(Z|Y,\theta^{(i)})\bigg) \\ &amp;=logP(Z|Y,\theta^{(i+1)})=0 \tag{17} \end{aligned}$H(θ(i+1),​θ(i))−H(θ(i),θ(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)))=logP(Z∣Y,θ(i+1))=0​(17)
  这里的不等号由Jensen不等式得到。
  由式(16)和式(17)即知式(15)右端是非负的。
  定理２ 设$L(\theta)=logP(Y|\theta)$L(θ)=logP(Y∣θ)为观测数据的对数似然函数，$\theta^{(i)}(i=1,2,\cdots)$θ(i)(i=1,2,⋯)为EM算法得到的参数估计序列，$L(\theta^{(i)})(i=1,2,\cdots)$L(θ(i))(i=1,2,⋯)为对应的对数似然函数序列。
  (1) 如果$P(Y|\theta)$P(Y∣θ)有上界，则$L(\theta^{(i)})=logP(Y|\theta^{(i)})$L(θ(i))=logP(Y∣θ(i))收敛到某一值$L^{*}$L∗；
  (2) 在函数$Q(\theta,\theta&#x27;)$Q(θ,θ′)与$L(\theta)$L(θ)满足一定条件下，由EM算法得到的参数估计序列$\theta^{(i)}$θ(i)的收敛值$\theta^{*}$θ∗是$L(\theta)$L(θ)的稳定点。
  证明 (1)由$L(\theta)=logP(Y|\theta^{(i)})$L(θ)=logP(Y∣θ(i))的单调性及$P(Y|\theta)$P(Y∣θ)的有界性立即得到。
  (2) 证明从略。
  定理２关于函数$Q(\theta,\theta&#x27;)$Q(θ,θ′)与$L(\theta)$L(θ)的条件在大多数情况下都是满足的。EM算法的收敛性包含关于对数似然函数序列$L(\theta^{(i)})$L(θ(i))的收敛性和关于参数估计序列$\theta^{(i)}$θ(i)的收敛性两层意思，前者并不蕴含后者。此外，定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点。所以在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

EM算法在高斯混合模型学习中的应用

EM算法的一个重要应用是高斯混合模型的参数估计。高斯混合模型应用广泛，在许多情况下，EM算法是学习高斯混合模型的有效方法。

• 高斯混合模型
  定义2（高斯混合模型） 高斯混合模型是指具有如下形式的概率分布模型：
  $P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k) \tag{18}$P(y∣θ)=k=1∑K​αk​ϕ(y∣θk​)(18)
  其中，$\alpha_k$αk​是系数，$\alpha_k\geq0$αk​≥0，$\sum_{k=1}^K\alpha_k=1$∑k=1K​αk​=1；$\phi(y|\theta_k)$ϕ(y∣θk​)是高斯分布密度，$\theta_k=(\mu_k,\sigma_k^2)$θk​=(μk​,σk2​)，
  $\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp\bigg(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\bigg)　\tag{19}$ϕ(y∣θk​)=2π​σk​1​exp(−2σk2​(y−μk​)2​)　(19)
  称为第$k$k个分模型。
  一般混合模型可以由任意概率分布密度代替式(19)中的高斯分布密度，我们只介绍最常用的高斯混合模型。

• 高斯混合模型参数估计的EM算法
  假设观测数据$y_1,y_2,\cdots,y_N$y1​,y2​,⋯,yN​由高斯混合模型生成，
  $P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta) \tag{20}$P(y∣θ)=k=1∑K​αk​ϕ(y∣θ)(20)
  其中，$\theta=(\alpha_1,\alpha_2,\cdots,\alpha_K;\theta_1,\theta_2,\cdots,\theta_K)$θ=(α1​,α2​,⋯,αK​;θ1​,θ2​,⋯,θK​)。我们用EM算法估计高斯混合模型的参数$\theta$θ。

1. 明确隐变量，写出完全数据的对数似然函数
   可以摄像观测数据$y_j,j=1,2,\cdots,N$yj​,j=1,2,⋯,N，是这样产生的：首先依概率$\alpha_k$αk​选择第$k$k个高斯分布模型$\phi(y|\theta_k)$ϕ(y∣θk​)；然后依第$k$k个分模型的概率分布$\phi(y|\theta_k)$ϕ(y∣θk​)生成观测数据$y_j$yj​。这时观测数据$y_j,j=1,2,\cdots,N$yj​,j=1,2,⋯,N，是已知的；反映观测数据$y_j$yj​来自第$k$k个分模型的数据是未知的，$k=1,2,\cdots,K$k=1,2,⋯,K，以隐变量$\gamma_{jk}$γjk​表示，定义如下：
   $\gamma_{jk}=\left\{ \begin{array}{ll} 1, &amp; 第j个观测来自第k个模型 \\ 0, &amp; 否则 \end{array}\right. \\ j=1,2,\cdots,N; k=1,2,\cdots,K \tag{21}$γjk​={1,0,​第j个观测来自第k个模型否则​j=1,2,⋯,N;k=1,2,⋯,K(21)
   $\gamma_{jk}$γjk​是0-1随机变量。
   有了观测数据$y_j$yj​及未观测数据$\gamma_{jk}$γjk​，那么完全数据是
   $(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}), j=1,2,\cdots,N$(yj​,γj1​,γj2​,⋯,γjK​),j=1,2,⋯,N
   于是，可以写出完全数据的似然函数：
   $\begin{aligned} P(y,\gamma|\theta)&amp;=\prod_{j=1}^NP(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}|\theta) \\ &amp;=\prod_{k=1}^K\prod_{j=1}^N[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &amp;=\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &amp;=\prod_{k=1}^K\alpha_k^{n_k}\prod_{k=1}^N\bigg[\frac{1}{\sqrt{2\pi}\sigma_k}\exp\bigg(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\bigg)\bigg]^{\gamma_{jk}} \end{aligned}$P(y,γ∣θ)​=j=1∏N​P(yj​,γj1​,γj2​,⋯,γjK​∣θ)=k=1∏K​j=1∏N​[αk​ϕ(yj​∣θk​)]γjk​=k=1∏K​αknk​​j=1∏N​[ϕ(yj​∣θk​)]γjk​=k=1∏K​αknk​​k=1∏N​[2π​σk​1​exp(−2σk2​(yj​−μk​)2​)]γjk​​
   式中，$n_k=\sum_{j=1}^N\gamma_{jk}$nk​=∑j=1N​γjk​，$\sum_{k=1}^Kn_k=N$∑k=1K​nk​=N。
   那么，完全数据的对数似然函数为
   $\log P(y,\gamma|\theta)=\sum_{k=1}^Kn_k\log\alpha_k+\sum_{j=1}^N\gamma_{jk}\bigg[\log\bigg(\frac{1}{\sqrt{2\pi}}\bigg)-\log\sigma_k-\frac{1}{2\sigma^2}(y_j-\mu_k)^2\bigg] \tag{22}$logP(y,γ∣θ)=k=1∑K​nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σ21​(yj​−μk​)2](22)
2. EM算法的E步：确定Q函数
   $\begin{aligned} Q(\theta,\theta^{(i)})&amp;=E[\log P(y,\gamma|\theta)|y,\theta^{(i)}] \\ &amp;=E\bigg\{\sum_{k=1}^Kn_k\log\alpha_k+\sum_{j=1}^N\gamma_{jk}\bigg[\log\bigg(\frac{1}{\sqrt{2\pi}}\bigg)-\log\sigma_k-\frac{1}{2\sigma_k^2(y_j-\mu_k)^2}\bigg]\bigg\} \\ &amp;=\sum_{k=1}^K\bigg\{\sum_{j=1}^N(E\gamma_{jk})\log\alpha_k+\sum_{j=1}^N(E\gamma_{jk})\bigg[\log\bigg(\frac{1}{\sqrt{2\pi}}\bigg)-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\bigg]\bigg\} \end{aligned}$Q(θ,θ(i))​=E[logP(y,γ∣θ)∣y,θ(i)]=E{k=1∑K​nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σk2​(yj​−μk​)21​]}=k=1∑K​{j=1∑N​(Eγjk​)logαk​+j=1∑N​(Eγjk​)[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}​
   这里需要计算$E(\gamma_{jk}|y,\theta)$E(γjk​∣y,θ)，记为$\hat\gamma_{jk}$γ^​jk​。
   $\begin{aligned} \hat\gamma_{jk}&amp;=E(\gamma_{jk}|y,\theta)=P(\gamma_{jk}=1|y,\theta) \\ &amp;=\frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^KP(\gamma_{jk}=1,y_j|\theta)} \\ &amp;=\frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)} \\ &amp;=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta)}, j=1,2,\cdots; k=1,2,\cdots,K \end{aligned}$γ^​jk​​=E(γjk​∣y,θ)=P(γjk​=1∣y,θ)=∑k=1K​P(γjk​=1,yj​∣θ)P(γjk​=1,yj​∣θ)​=∑k=1K​P(yj​∣γjk​=1,θ)P(γjk​=1∣θ)P(yj​∣γjk​=1,θ)P(γjk​=1∣theta)​=∑k=1K​αk​ϕ(yj​∣θ)αk​ϕ(yj​∣θk​)​,j=1,2,⋯;k=1,2,⋯,K​
   $\hat\gamma_{jk}$γ^​jk​是在当前模型参数下第$j$j个观测数据来自第$k$k个分模型的概率，称为分模型$k$k对观测数据$y_j$yj​的响应度。
   将$\hat\gamma_{jk}=E\gamma_{jk}$γ^​jk​=Eγjk​及$n_k=\sum_{j=1}^NE\gamma_{jk}$nk​=∑j=1N​Eγjk​代入式(22)即得
   $Q(\theta,\theta^{(i)})=\sum_{k=1}^Kn_k\log\alpha_k+\sum_{k=1}^N\hat\gamma_{jk}\bigg[\log\bigg(\frac{1}{\sqrt{2\pi}}\bigg)-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\bigg] \tag{23}$Q(θ,θ(i))=k=1∑K​nk​logαk​+k=1∑N​γ^​jk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2](23)
3. 确定EM算法的M步
   迭代的M步是求函数$Q(\theta,\theta^{(i)})$Q(θ,θ(i))对$\theta$θ的极大值，即求新一轮迭代的模型参数：
   $\theta^{(i+1)}=\arg\max_\theta Q(\theta,\theta^{(i)})$θ(i+1)=argθmax​Q(θ,θ(i))
   用$\hat\mu_k$μ^​k​，$\hat\sigma_k^2$σ^k2​及$\hat\alpha_k$α^k​，$k=1,2,\cdots,K$k=1,2,⋯,K，表示$\theta^{(i+1)}$θ(i+1)的各参数。求$\hat\mu_k$μ^​k​，$\hat\sigma_k^2$σ^k2​只需将式(23)分别对$\mu_k$μk​，$\sigma_k^2$σk2​求偏导数并令其为０，即可得到；求$\hat\alpha_k$α^k​是在$\sum_{k=1}^K\alpha_k=1$∑k=1K​αk​=1条件下求偏导数并令其为０得到的。结果如下：
   $\hat\mu_k=\frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}}, k=1,2,\cdots,K \tag{24}$μ^​k​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​yj​​,k=1,2,⋯,K(24)
   $\hat\sigma_k^2=\frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}}, k=1,2,\cdots,K \tag{25}$σ^k2​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​(yj​−μk​)2​,k=1,2,⋯,K(25)
   $\hat\alpha_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N\hat\gamma_{jk}}{N}, k=1,2,\cdots,N \tag{26}$α^k​=Nnk​​=N∑j=1N​γ^​jk​​,k=1,2,⋯,N(26)
   重复以上计算，直到对数似然函数值不再有明显的变化为止。
   现将估计高斯混合模型参数的EM算法总结如下：
   算法2（高斯混合模型参数估计的EM算法）
   输入：观测数据$y_1,y_2,\cdots,y_N$y1​,y2​,⋯,yN​，高斯混合模型；
   输出：高斯混合模型参数。
   (1) 取参数的初始值开始迭代
   (2) E步：依据当前模型参数，计算分模型$k$k对观测数据$y_j$yj​的响应度
   $\hat\gamma_{jk}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta)}, j=1,2,\cdots; k=1,2,\cdots,K$γ^​jk​=∑k=1K​αk​ϕ(yj​∣θ)αk​ϕ(yj​∣θk​)​,j=1,2,⋯;k=1,2,⋯,K
   (3) M步：计算新一轮迭代的模型参数
   $\hat\mu_k=\frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}}, k=1,2,\cdots,K$μ^​k​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​yj​​,k=1,2,⋯,K
   $\hat\sigma_k^2=\frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}}, k=1,2,\cdots,K$σ^k2​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​(yj​−μk​)2​,k=1,2,⋯,K
   $\hat\alpha_k=\frac{\sum_{j=1}^N\hat\gamma_{jk}}{N}, k=1,2,\cdots,N$α^k​=N∑j=1N​γ^​jk​​,k=1,2,⋯,N
   (4) 重复第(2)步和第(3)步，直到收敛。

EM算法的推广

EM算法还可以解释为F函数的极大-极大算法，基于这个解释有若干变形与推广，如广义期望极大算法。下面予以介绍。

• F函数的极大-极大算法
  首先引进F函数并讨论其性质。
  定义3（F函数） 假设隐变量数据$Z$Z的概率分布为$\tilde P(Z)$P~(Z)，定义分布$\tilde P$P~与参数$\theta$θ的函数$F(\tilde P, \theta)$F(P~,θ)如下：
  $F(\tilde P,\theta)=E_{\tilde P}[\log P(Y,Z|\theta)]+H(\tilde P) \tag{27}$F(P~,θ)=EP~​[logP(Y,Z∣θ)]+H(P~)(27)
  称为F函数。式中$H(\tilde P)=-E_{\tilde P}\log\tilde P(Z)$H(P~)=−EP~​logP~(Z)是分布$\tilde P(Z)$P~(Z)的熵。
  在定义3中，通常假设$P(Y,Z|\theta)$P(Y,Z∣θ)是$\theta$θ的连续函数，因而$F(\tilde P,\theta)$F(P~,θ)是$\tilde P$P~和$\theta$θ的连续函数。函数$F(\tilde P,\theta)$F(P~,θ)还有以下重要性质：
  引理1 对于固定的$\theta$θ，存在唯一的分布$\tilde P_\theta$P~θ​极大化$F(\tilde P,\theta)$F(P~,θ)，这时$\tilde P_\theta$P~θ​由下式给出：
  $\tilde P_\theta(Z)=P(Z|Y,\theta) \tag{28}$P~θ​(Z)=P(Z∣Y,θ)(28)
  并且$\tilde P_\theta$P~θ​随$\theta$θ连续变化。
  证明 对于固定的$\theta$θ，可以求得使$F(\tilde P,\theta)$F(P~,θ)达到极大的分布$\tilde P_\theta(Z)$P~θ​(Z)。为此，引进拉格朗日乘子$\lambda$λ，拉格朗日函数为
  $L=E_{\tilde P}\log P(Y,Z|\theta)-E_{\tilde P}\log \tilde P(Z)+\lambda\bigg(1-\sum_Z\tilde P(Z)\bigg) \tag{29}$L=EP~​logP(Y,Z∣θ)−EP~​logP~(Z)+λ(1−Z∑​P~(Z))(29)
  将其对$\tilde P$P~求偏导数：
  $\frac{\partial L}{\partial\tilde P(Z)}=\log P(Y,Z|\theta)-\log\tilde P(Z)-1-\lambda$∂P~(Z)∂L​=logP(Y,Z∣θ)−logP~(Z)−1−λ
  令偏导数等于0，得出
  $\lambda=\log P(Y,Z|\theta)-\log \tilde P_\theta(Z)-1$λ=logP(Y,Z∣θ)−logP~θ​(Z)−1
  由此推出$\tilde P_\theta(Z)$P~θ​(Z)与$P(Y,Z|\theta)$P(Y,Z∣θ)成比例
  $\frac{P(Y,Z|\theta)}{\tilde P_\theta(Z)}=e^{1+\lambda}$P~θ​(Z)P(Y,Z∣θ)​=e1+λ
  再从约束条件$\sum_Z\tilde P_\theta(Z)=1$∑Z​P~θ​(Z)=1得式(28)。
  由假设$P(Y,Z|\theta)$P(Y,Z∣θ)是$\theta$θ的连续函数，得到$\tilde P_\theta$P~θ​是$\theta$θ的连续函数。
  引理2 若$\tilde P_\theta(Z)=P(Z|Y,\theta)$P~θ​(Z)=P(Z∣Y,θ)，则
  $F(\tilde P,\theta)=\log P(Y|\theta) \tag{30}$F(P~,θ)=logP(Y∣θ)(30)
  由以上引理，可以得到关于EM算法用F函数的极大-极大算法的解释。
  定理3 设$L(\theta)=\log P(Y|\theta)$L(θ)=logP(Y∣θ)为观测数据的对数似然函数，$\theta^{(i)},i=1,2,\cdots$θ(i),i=1,2,⋯，为EM算法得到的参数估计序列，函数$F(\tilde P,\theta)$F(P~,θ)由式(27)定义。若果$F(\tilde P,\theta)$F(P~,θ)在$\tilde P^*$P~∗和$\theta^*$θ∗有局部极大值，那么$L(\theta)$L(θ)也在$\theta^*$θ∗有局部极大值。类似地，如果$F(\tilde P,\theta)$F(P~,θ)在$\tilde P^*$P~∗和$\theta^*$θ∗达到全局最大值，那么$L(\theta)$L(θ)也在$\theta^*$θ∗达到全局最大值。
  证明 由引理1和引理2可知，$L(\theta)=\log P(Y|\theta)=F(\tilde P_\theta,\theta)$L(θ)=logP(Y∣θ)=F(P~θ​,θ)对任意$\theta$θ成立。特别地，对于使$F(\tilde P,\theta)$F(P~,θ)达到极大的参数$\theta^*$θ∗，有
  $L(\theta^*)=F(\tilde P_{\theta^*},\theta^*)=F(\tilde P^*,\theta^*) \tag{31}$L(θ∗)=F(P~θ∗​,θ∗)=F(P~∗,θ∗)(31)
  为了证明$\theta^*$θ∗是$L(\theta)$L(θ)的极大点，需要证明不存在接近$\theta^*$θ∗的点$\theta^{**}$θ∗∗，使$L(\theta^{**})&gt;L(\theta^*)$L(θ∗∗)>L(θ∗)。加入存在这样的点$\theta^{**}$θ∗∗，那么应有$F(\tilde P^{**},\theta^{**})&gt;F(\tilde P^*,\theta^*)$F(P~∗∗,θ∗∗)>F(P~∗,θ∗)，这里$\tilde P^{**}=\tilde P_{\theta^{**}}$P~∗∗=P~θ∗∗​。但因$\tilde P_\theta$P~θ​是随$\theta$θ连续变化的，$\tilde P^{**}$P~∗∗应接近$\tilde P^*$P~∗，这与$\tilde P^*$P~∗和$\theta^*$θ∗是$F(\tilde P,\theta)$F(P~,θ)的局部极大点的假设矛盾。
  类似可以证明关于全局最大值的结论。
  定理4 EM算法的一次迭代可由F函数的极大-极大算法实现。
  设$\theta^{(i)}$θ(i)为第$i$i次迭代参数$\theta$θ的估计，$\tilde P^{(i)}$P~(i)为第$i$i次迭代函数$\tilde P$P~的估计。在第$i+1$i+1次迭代的两步为
  (1) 对固定的$\theta^{(i)}$θ(i)，求$\tilde P^{(i+1)}$P~(i+1)使$F(\tilde P,\theta^{(i)})$F(P~,θ(i))极大化；
  (2) 对固定的$\tilde P^{(i+1)}$P~(i+1)，求$\theta^{i+1}$θi+1使$F(\tilde P^{(i+1)},\theta)$F(P~(i+1),θ)极大化。
  证明 (1) 由引理1，对于固定的$\theta^{(i)}$θ(i)，
  $\tilde P^{(i+1)}(Z)=\tilde P_{\theta^{(i)}}(Z)=P(Z|Y,\theta^{(i)})$P~(i+1)(Z)=P~θ(i)​(Z)=P(Z∣Y,θ(i))
  使$F(\tilde P,\theta^{(i)})$F(P~,θ(i))极大化。此时
  $\begin{aligned} F(\tilde P^{(i+1)},\theta)&amp;=E_{\tilde P^{(i+1)}}[\log P(Y,Z|\theta)]+H(\tilde P^{(i+1)}) \\ &amp;=\sum_Z\log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})+H(\tilde P^{(i+1)}) \end{aligned}$F(P~(i+1),θ)​=EP~(i+1)​[logP(Y,Z∣θ)]+H(P~(i+1))=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))+H(P~(i+1))​
  由$Q(\theta,\theta^{(i)})$Q(θ,θ(i))的定义式(5)有
  $F(\tilde P^{(i+1)},\theta)=Q(\theta,\theta^{(i)})+H(\tilde P^{(i+1)})$F(P~(i+1),θ)=Q(θ,θ(i))+H(P~(i+1))
  (2) 固定$\tilde P^{(i+1)}$P~(i+1)，求$\theta^{(i+1)}$θ(i+1)使$F(\tilde P^{(i+1)},\theta)$F(P~(i+1),θ)极大化。得到
  $\theta^{(i+1)}=\arg\max_\theta F(\tilde P^{(i+1)},\theta)=\arg\max_\theta Q(\theta,\theta^{(i)})$θ(i+1)=argθmax​F(P~(i+1),θ)=argθmax​Q(θ,θ(i))
  通过以上两步完成了EM算法的一次迭代。由此可知，由EM算法与F函数的极大-极大算法得到的参数估计序列$\theta^{(i)},i=1,2,\cdots$θ(i),i=1,2,⋯，是一致的。
  这样，就有EM算法的推广。

• GEM算法
  算法3（GEM算法1）
  输入：观测数据，F函数；
  输出：模型参数。
  (1) 初始化参数$\theta^{(0)}$θ(0)，开始迭代
  (2) 第$i+1$i+1次迭代，第1步：记$\theta^{(i)}$θ(i)为参数$\theta$θ的估计值，$\tilde P^{(i)}$P~(i)为函数$\tilde P$P~的估计。求$\tilde P^{(i+1)}$P~(i+1)使$\tilde P$P~极大化$F(\tilde P,\theta^{(i)})$F(P~,θ(i))
  (3) 第2步：求$\theta^{(i+1)}$θ(i+1)使$F(\tilde P^{(i+1)},\theta)$F(P~(i+1),θ)极大化
  (4) 重复(2)和(3)，直到收敛。
  在GEM算法1中，有时求$Q(\theta,\theta^{(i)})$Q(θ,θ(i))的极大化是很困难的。下面介绍的GEM算法2和GEM算法3并不是直接求$\theta^{(i+1)}$θ(i+1)使$Q(\theta,\theta^{(i)})$Q(θ,θ(i))达到极大的$\theta$θ，而是找一个$\theta^{(i+1)}$θ(i+1)使得$Q(\theta^{(i+1)},\theta^{(i)})&gt;Q(\theta^{(i)},\theta^{(i)})$Q(θ(i+1),θ(i))>Q(θ(i),θ(i))。
  算法4（GEM算法2）
  输入：观测数据，Q函数；
  输出：模型参数。
  (1) 初始化参数$\theta^{(0)}$θ(0)，开始迭代
  (2) 第$i+1$i+1次迭代，第1步：记$\theta^{(i)}$θ(i)为参数$\theta$θ的估计值，计算
  $\begin{aligned} Q(\theta,\theta^{(i)})&amp;=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}] \\ &amp;=\sum_ZP(Z|Y,\theta^{(i)})\log P(Y,Z|\theta) \end{aligned}$Q(θ,θ(i))​=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​P(Z∣Y,θ(i))logP(Y,Z∣θ)​
  (3) 第2步：求$\theta^{(i+1)}$θ(i+1)使
  KaTeX parse error: Expected 'EOF', got '\thata' at position 35: …theta^{(i)})>Q(\̲t̲h̲a̲t̲a̲{(i)},\theta^{(…
  (4) 重复(2)和(3)，直到收敛。
  当参数$\theta$θ的维数为$d(d\geq 2)$d(d≥2)时，可采用一种特殊的GEM算法，它将EM算法的M步分解为d次条件极大化，每次只改变参数向量的一个分量，其余分量不改变。
  算法5（GEM算法3）
  输入：观测数据，Q函数；
  输出：模型参数。
  (1) 初始化参数$\theta^{(0)}=(\theta_1^{(0)},\theta_2^{(0)},\cdots,\theta_d^{(0)})$θ(0)=(θ1(0)​,θ2(0)​,⋯,θd(0)​)，开始迭代
  (2) 第$i+1$i+1次迭代，第1步：记$\theta^{(i)}=(\theta_1^{(i)},\theta_2^{(i)},\cdots,\theta_d^{(i)})$θ(i)=(θ1(i)​,θ2(i)​,⋯,θd(i)​)为参数$\theta=(\theta_1,\theta_2,\cdots,\theta_d)$θ=(θ1​,θ2​,⋯,θd​)的估计值，计算
  $\begin{aligned} Q(\theta,\theta^{(i)})&amp;=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}] \\ &amp;=\sum_ZP(Z|Y,\theta^{(i)})\log P(Y,Z|\theta) \end{aligned}$Q(θ,θ(i))​=EZ​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​P(Z∣Y,θ(i))logP(Y,Z∣θ)​
  (3) 第2步：进行d次条件极大化：
  首先，在$\theta_2^{(i)},\cdots,\theta_k^{(i)}$θ2(i)​,⋯,θk(i)​保持不变的条件下求使$Q(\theta,\theta^{(i)})$Q(θ,θ(i))达到极大的$\theta_2^{(i+1)}$θ2(i+1)​；
  然后，在$\theta_1=\theta_1^{(i+1)},\theta_j=\theta_j^{(i)},j=2,3,\cdots,k$θ1​=θ1(i+1)​,θj​=θj(i)​,j=2,3,⋯,k的条件下求使$Q(\theta,\theta^{(i)})$Q(θ,θ(i))达到极大的$\theta_2^{(i+1)}$θ2(i+1)​；
  如此继续，经过d次条件极大化，得到$\theta^{(i+1)}=(\theta_1^{(i+1)},\theta_2^{(i+1)},\cdots,\theta_d^{(i+1)})$θ(i+1)=(θ1(i+1)​,θ2(i+1)​,⋯,θd(i+1)​)使得
  KaTeX parse error: Expected 'EOF', got '\thata' at position 35: …theta^{(i)})>Q(\̲t̲h̲a̲t̲a̲{(i)},\theta^{(…
  (4) 重复(2)和(3)，直到收敛。