EM算法原理介绍

1. 引言

    EM算法是Dempster等人在1977年提出来的一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望；M步，求极大，因此，该算法也被称为期望极大算法，简称EM算法。

2. EM算法原理介绍

2.1 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)=\log P(Y | \theta)$L(θ)=logP(Y∣θ)，假设$Y$Y和$Z$Z的联合概率分布是$P(Y, Z | \theta)$P(Y,Z∣θ)，那么完全数据对数似然函数是$\log P(Y, Z | \theta)$logP(Y,Z∣θ)。

    EM算法就是通过极大化不完全数据$Y$Y的对数似然函数来对参数$\theta$θ进行估计，即极大化：
$\begin{aligned} 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{aligned}$L(θ)​=logP(Y∣θ)=logZ∑​P(Y,Z∣θ)=log(Z∑​P(Y∣Z,θ)P(Z∣θ))​由于上式中含有未观测的数据和求和的对数，因此，没法直接对参数进行极大化估计。事实上，EM算法是通过迭代逐步近似极大化$L(\theta)$L(θ)，假设在第$i$i次迭代后$\theta$θ的估计值是$\theta^{(i)}$θ(i)，我们希望估计值$\theta$θ能使$L(\theta)$L(θ)增加，即$L(\theta)>L\left(\theta^{(i)}\right)$L(θ)>L(θ(i))，并逐步达到极大值，因此，可以直接考虑两者的差：
$L(\theta)-L\left(\theta^{(i)}\right)=\log \left(\sum_{Z} P(Y | Z, \theta) P(Z | \theta)\right)-\log P\left(Y | \theta^{(i)}\right)$L(θ)−L(θ(i))=log(Z∑​P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))利用Jensen不等式可以得到其下界：
$\begin{aligned} L(\theta)-L\left(\theta^{(i)}\right) &=\log \left(\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}\right)-\log P\left(Y | \theta^{(i)}\right) \\ & \geqslant \sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right)}-\log P\left(Y | \theta^{(i)}\right) \\ &=\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right) P\left(Y | \theta^{(i)}\right)} \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∣θ(i))P(Y∣Z,θ)P(Z∣θ)​​令
$B\left(\theta, \theta^{(i)}\right) \hat{=} L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right) P\left(Y | \theta^{(i)}\right)}$B(θ,θ(i))=^L(θ(i))+Z∑​P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)​则
$L(\theta) \geqslant B\left(\theta, \theta^{(i)}\right)$L(θ)⩾B(θ,θ(i))即函数$B\left(\theta, \theta^{(i)}\right)$B(θ,θ(i))是$L(\theta)$L(θ)的一个下界，并且有：
$L\left(\theta^{(i)}\right)=B\left(\theta^{(i)}, \theta^{(i)}\right)$L(θ(i))=B(θ(i),θ(i))因此，任何可以使$B\left(\theta, \theta^{(i)}\right)$B(θ,θ(i))增大的$\theta$θ也可以使$L(\theta)$L(θ)增大，因此，每次迭代时可以直接对$B\left(\theta, \theta^{(i)}\right)$B(θ,θ(i))进行极大化更新$\theta$θ：
$\theta^{(i+1)}=\arg \max _{\theta} B\left(\theta, \theta^{(i)}\right)$θ(i+1)=argθmax​B(θ,θ(i))对其求$\theta$θ偏导，有：
$\begin{aligned} \theta^{(i+1)} &=\arg \max _{\theta}\left(L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log \frac{P(Y | Z, \theta) P(Z | \theta)}{P\left(Z | Y, \theta^{(i)}\right) P\left(Y | \theta^{(i)}\right)}\right) \\ &=\arg \max _{\theta}\left(\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log (P(Y | Z, \theta) P(Z | \theta))\right) \\ &=\arg \max _{\theta}\left(\sum_{Z} P\left(Z | Y, \theta^{(i)}\right) \log P(Y, Z | \theta)\right) \\ &=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right) \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))​上式等价于EM算法的一次迭代，即求$Q$Q函数及其极大化，EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。最终，EM算法可以归纳如下：

    EM算法：

• 输入：观测变量$Y$Y，隐变量数据$Z$Z，联合分布$P(Y, Z | \theta)$P(Y,Z∣θ)，条件分布$P(Z | Y, \theta)$P(Z∣Y,θ)；
• 输出：模型参数$\theta$θ
• 选择参数初始值$\theta^{(0)}$θ(0)，开始迭代；
• E步：记$\theta^{(i)}$θ(i)为第$i$i次迭代参数$\theta$θ的估计值，在第$i+1$i+1次迭代的E步，计算期望，即Q函数：$\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &=E_{z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right] \\ &=\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right) \end{aligned}$Q(θ,θ(i))​=Ez​[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))​这里，$P\left(Z | Y, \theta^{(i)}\right)$P(Z∣Y,θ(i))是在给定观测数据$Y$Y和当前的参数估计$\theta^{(i)}$θ(i)下的隐变量数据$Z$Z的条件概率分布；
• M步：求使$Q\left(\theta, \theta^{(i)}\right)$Q(θ,θ(i))极大化的$\theta$θ，确定第$i+1$i+1次迭代的参数的估计值：
  $\theta^{(i+1)}=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right)$θ(i+1)=argθmax​Q(θ,θ(i))
• 重复4、5步骤，直到收敛，即$\left\|\theta^{(i+1)}-\theta^{(i)}\right\|&lt;\varepsilon_{1}$∥∥​θ(i+1)−θ(i)∥∥​<ε1​或者$\left\|Q\left(\theta^{(i+1)}, \theta^{(i)}\right)-Q\left(\theta^{(i)}, \theta^{(i)}\right)\right\|&lt;\varepsilon_{2}$∥∥​Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∥∥​<ε2​，其中，$\varepsilon_{1}, \varepsilon_{2}$ε1​,ε2​为设定的阈值。

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

    EM算法的一个重要的应用是高斯混合模型的参数估计。首先，先看下什么是高斯混合模型。

    高斯混合模型： 高斯混合模型是指具有如下形式的概率分布模型：
$P(y | \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y | \theta_{k}\right)$P(y∣θ)=k=1∑K​αk​ϕ(y∣θk​)其中，$\alpha_{k}$αk​是系数，$\alpha_{k} \geqslant 0, \quad \sum_{k=1}^{K} \alpha_{k}=1$αk​⩾0,∑k=1K​αk​=1，$\phi\left(y | \theta_{k}\right)$ϕ(y∣θk​)是高斯分布函数，$\theta_{k}=\left(\mu_{k}, \sigma_{k}^{2}\right)$θk​=(μk​,σk2​)：
$\phi\left(y | \theta_{k}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)$ϕ(y∣θk​)=2π​σk​1​exp(−2σk2​(y−μk​)2​)称为第$k$k个分模型。

    假设观测数据$y_{1}, y_{2}, \cdots, y_{N}$y1​,y2​,⋯,yN​由高斯混合模型生成，
$P(y | \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y | \theta_{k}\right)$P(y∣θ)=k=1∑K​αk​ϕ(y∣θk​)其中，$\theta=\left(\alpha_{1}, \alpha_{2}, \cdots, \alpha_{K} ; \theta_{1}, \theta_{2}, \cdots, \theta_{K}\right)$θ=(α1​,α2​,⋯,αK​;θ1​,θ2​,⋯,θK​)，接下来用EM算法估计高斯混合模型的参数$\theta$θ。

    可以设想观测数据$y_{j}, \quad j=1,2, \cdots, N$yj​,j=1,2,⋯,N是这样产生的：首先依概率$\alpha_{k}$αk​选择第$k$k个高斯分布模型$\phi\left(y | \theta_{k}\right)$ϕ(y∣θk​)，然后依第$k$k个分模型的概率分布$\phi\left(y | \theta_{k}\right)$ϕ(y∣θk​)生成观测数据$y_{j}$yj​，这时，观测数据$y_{j}, \quad j=1,2, \cdots, N$yj​,j=1,2,⋯,N是已知的，反映观测数据$y_{j}$yj​来自第$k$k个分模型是未知的，用隐变量$\gamma_{j k}$γjk​表示，其定义如下：
其中，$\gamma_{j k}$γjk​是0-1随机变量。因此，可以得到完全数据的似然函数：
其中，$n_{k}=\sum_{j=1}^{N} \gamma_{j k}, \quad \sum_{k=1}^{K} n_{k}=N$nk​=∑j=1N​γjk​,∑k=1K​nk​=N。因此，完全数据的对数似然函数可以表达为：
$\log P(y, \gamma | \theta)=\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]$logP(y,γ∣θ)=k=1∑K​nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]

    接下来，可以计算Q函数：
$\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &amp;=E\left[\log P(y, \gamma | \theta) | y, \theta^{(i)}\right] \\ &amp;=E\left\{\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} \\ &amp;=\sum_{k=1}^{K}\left\{\sum_{j=1}^{N}\left(E \gamma_{j k}\right) \log \alpha_{k}+\sum_{j=1}^{N}\left(E \gamma_{j k}\right)\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} \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​1​(yj​−μk​)2]}=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\left(\gamma_{j k} | y, \theta\right)$E(γjk​∣y,θ)，记为$\hat{\gamma}_{j k}$γ^​jk​，其计算如下：

将$\hat{\gamma}_{j k}=E \gamma_{j k}$γ^​jk​=Eγjk​，$n_{k}=\sum_{j=1}^{N} E \gamma_{j k}$nk​=∑j=1N​Eγjk​代入Q函数得：
$Q\left(\theta, \theta^{(i)}\right)=\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \hat{\gamma}_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]$Q(θ,θ(i))=k=1∑K​nk​logαk​+j=1∑N​γ^​jk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]

    接下来是M步，计算Q函数的极大值，
$\theta^{(i+1)}=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right)$θ(i+1)=argθmax​Q(θ,θ(i))用$\hat{\mu}_{k}, \hat{\sigma}_{k}^{2}\hat{\alpha}_{k}, k=1,2, \cdots, K$μ^​k​,σ^k2​α^k​,k=1,2,⋯,K分别表示$\theta^{(i+1)}$θ(i+1)的各参数，分别对$\mu_{k}, \sigma_{k}^{2}, \alpha_{k}$μk​,σk2​,αk​求偏导并令其为0得：
$\begin{array}{c}{\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K} \\\\ {\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K} \\\\ {\hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K}\end{array}$μ^​k​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​yj​​,k=1,2,⋯,Kσ^k2​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​(yj​−μk​)2​,k=1,2,⋯,Kα^k​=Nnk​​=N∑j=1N​γ^​jk​​,k=1,2,⋯,K​ 重复以上过程，直到收敛为止。高斯混合模型参数估计的EM算法可以总结如下：

    高斯混合模型参数估计的EM算法：

1. 输入：观测数据$y_{1}, y_{2}, \cdots, y_{N}$y1​,y2​,⋯,yN​，高斯混合模型；
2. 输出：高斯混合模型的参数；
3. 取参数的初始值开始迭代
4. E步：依据当前的参数，计算分模型$k$k对观测数据$y_j$yj​的响应度：
   $\hat{\gamma}_{j k}=\frac{\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K$γ^​jk​=∑k=1K​αk​ϕ(yj​∣θk​)αk​ϕ(yj​∣θk​)​,j=1,2,⋯,N;k=1,2,⋯,K
5. M步：计算新一轮迭代的模型参数：
   $\begin{array}{c}{\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K} \\\\ {\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K} \\\\ {\hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K}\end{array}$μ^​k​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​yj​​,k=1,2,⋯,Kσ^k2​=∑j=1N​γ^​jk​∑j=1N​γ^​jk​(yj​−μk​)2​,k=1,2,⋯,Kα^k​=Nnk​​=N∑j=1N​γ^​jk​​,k=1,2,⋯,K​
6. 重复4、5步，直到收敛。

3. 总结

• EM算法对初始值比较敏感。
• EM算法由于是通过迭代的思想，采用下界不断逼近对数似然函数，因此，得到的参数估计可能是局部最优解。