EM算法

一、EM算法的引入

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

1. 为什么要使用EM算法

首先我们来看一个使用EM算法的例子。

例1（三硬币模型）    假设有3枚硬币，分别记作$A$，$B$，$C$。这些硬币正面出现的概率分别是$\pi$，$p$和$q$。进行如下掷硬币试验：先掷硬币$A$，根据其结果选出硬币$B$或硬币$C$，正面选硬币$B$，反面选硬币$C$；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复$n$次实验（这里$n=10$）观测结果如下：

$$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三枚硬币正面出现的概率，即三硬币模型的参数。
解    三硬币模型可写作

$$\begin{array}{rl}P(y\mid \mathbf{\theta })& =\sum _{z}P(y,z\mid \mathbf{\theta })=\sum _{z}P(z\mid \mathbf{\theta })P(y\mid z,\mathbf{\theta })\\ & =\pi p^y{(1-p)}^{1-y}+(1-\pi )q^y{(1-q)}^{1-y}\end{array}\text{\hspace{1em}(1.1)}$$

这里，随机变量$y$是观测变量，表示一次试验观测的结果是$1$和$0$；随机变量$z$是隐变量，表示未观测到的掷硬币$A$的结果；$\mathbf{\theta }=(\pi ,p,q)$是模型参数。这一模型是以上数据的生成模型。注意，随机变量$y$的数据可以观测，随机变量$z$的数据不可观测。

将观测数据表示为$\mathbf{Y}={\left(\begin{array}{cccc}{Y}_1& Y_2& \cdots & Y_n\end{array}\right)}^T$，未观测数据表示为$\mathbf{Z}={\left(\begin{array}{cccc}{Z}_1& Z_2& \cdots & Z_n\end{array}\right)}^T$，则观测数据的似然函数为

$$P(\mathbf{Y}\mid \mathbf{\theta })=\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{\theta })P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })\text{\hspace{1em}(1.2)}$$

即

$$P(\mathbf{Y}\mid \mathbf{\theta })=\prod _{j=1}^n\left[\pi p^{y_i}{(1-p)}^{1-y_i}+(1-\pi )q^{y_i}{(1-q)}^{1-y_i}\right]\text{\hspace{1em}(1.3)}$$

考虑求模型参数$\mathbf{\theta }=(\pi ,p,q)$的极大似然估计，即

$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }}{\max }\log P(\mathbf{Y}\mid \mathbf{\theta })\text{\hspace{1em}(1.4)}$$

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。下面我们就给出上面问题的EM算法。

EM算法首先呢选取参数的初值，记作${\mathbf{\theta }}^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，然后通过下面的步骤迭代计算参数的估计值，直至收敛为止。第$i$次迭代参数的估计值为${\mathbf{\theta }}^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$。EM算法的第$i+1$次迭代如下。

E步：计算在模型参数${\pi }^{(i)}$，$p^{(i)}$，$q^{(i)}$下观测数据$y_j$来自掷硬币$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-q^{(i)})}^{1-y_i}}\text{\hspace{1em}(1.5)}$$

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

$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}\text{\hspace{1em}(1.6)}$$

$$p^{(i+1)}=\frac{\sum _{j=1}^n{\mu }_j^{(i+1)}{y}_j}{\sum _{j=1}^n{\mu }_j^{(i+1)}}\text{\hspace{1em}(1.7)}$$

$$q^{(i+1)}=\frac{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{\sum _{j=1}^n(1-{\mu }_j^{(i+1)})}\text{\hspace{1em}(1.8)}$$

进行数值计算。假设模型参数的初值取为

$${\pi }^{(0)}=0.5,p^{(0)}=0.5,q^{(0)}=0.5$$

由式$(1.5)$，对$y_j=1$和$y_j=0$均有${\mu }_j^{(1)}=0.5$.

利用迭代公式$(1.6)$~公式$(1.8)$，得到

$${\pi }^{(1)}=0.5,p^{(1)}=0.6,q^{(1)}=0.6$$

由式$(1.5)$，

${\mu }_j^{(2)}=0.5,j=1,2,\cdots ,10$

继续迭代，得

$${\pi }^{(2)}=0.5,p^{(2)}=0.6,q^{(2)}=0.6$$

于是得到模型参数$\mathbf{\theta }$的极大似然估计：

$$\hat{\pi }=0.5,\hat{p}=0.6,\hat{q}=0.6$$

$\pi =0.5$表示硬币$A$是均匀的，这一结果容易理解。

如果取初值${\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算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般地，用$\mathbf{Y}$表示观测随机变量的数据，$\mathbf{Z}$表示隐随机变量的数据。$\mathbf{Y}$和$\mathbf{Z}$连在一起称为完全数据（complete-data），观测数据$\mathbf{Y}$又称为不完全数据（incomplete-data）。假设给定观测数据$\mathbf{Y}$，其概率分布是$P(\mathbf{Y}\mid \mathbf{\theta })$，其中$\mathbf{\theta }$是需要估计的模型参数，那么不完全数据$Y$的似然函数是$P(\mathbf{Y}\mid \mathbf{\theta })$，对数似然函数$L(\mathbf{\theta })=\log P(\mathbf{Y}\mid \mathbf{\theta })$；假设$\mathbf{Y}$和$\mathbf{Z}$的联合概率分布是$P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })$，那么完全数据的对数似然函数是$\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })$。

2. EM算法的推导

我们面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据）$Y$关于参数$\theta$的对数似然函数，即极大化

$$\begin{array}{rl}L(\mathbf{\theta })& =\log P(\mathbf{Y}\mid \mathbf{\theta })=\log \sum _{\mathbf{Z}}P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\\ & =\log \left(\sum _{\mathbf{Z}}P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right)\end{array}\text{\hspace{1em}(1.9)}$$

注意到这一极大化的主要困难是式$(1.9)$中有未观测数据并有包含和（或积分）的对数。

事实上，EM算法是通过迭代逐步近似极大化$L(\mathbf{\theta })$的。假设在第$i$次迭代后$\mathbf{\theta }$的估计值是${\mathbf{\theta }}^{(i)}$。我们希望新估计的$\mathbf{\theta }$能使$L(\mathbf{\theta })$增加，即$L(\mathbf{\theta })>L({\mathbf{\theta }}^{(i)})$，并逐步达到极大值。因此，我们考虑两者的差：

$$L(\mathbf{\theta })-L({\mathbf{\theta }}^{(i)})=\log \left(\sum _{\mathbf{Z}}P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right)-\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})\text{\hspace{1em}(1.10)}$$

引理（琴生不等式（Jensen inequality））  若$f(x)$是区间$\left(a,b\right)$上的凸函数，则对于任意$x_1,x_2,\cdots ,x_n\in \left(a,b\right)$，${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\in {\mathbb{R}}^{+}$且$\sum _{i=1}^n{\lambda }_i=1$，有

$$f\left(\sum _{i=1}^n{\lambda }_i{x}_i\right)⩽\sum _{i=1}^n{\lambda }_{i}f\left(x_i\right)\text{\hspace{1em}(1.11)}$$

若$f(x)$是区间$\left(a,b\right)$上的凹函数，则对于任意$x_1,x_2,\cdots ,x_n\in \left(a,b\right)$，${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\in {\mathbb{R}}^{+}$且$\sum _{i=1}^n{\lambda }_i=1$，有

$$f\left(\sum _{i=1}^n{\lambda }_i{x}_i\right)⩾\sum _{i=1}^n{\lambda }_{i}f\left(x_i\right)\text{\hspace{1em}(1.12)}$$

证明    以$f(x)$是凸函数为例，采用数学归纳法来证明：

当$n=1$时，${\lambda }_1=1$，所以有

$$f(x_1)⩽f(x_1)\text{\hspace{1em}(1.13)}$$

恒成立。

当$n=2$时，${\lambda }_1+{\lambda }_2=1$，构造函数

$$F(x)={\lambda }_{1}f(x_1)+{\lambda }_{2}f(x)-f\left({\lambda }_1{x}_1+{\lambda }_{2}x\right)\text{\hspace{1em}(1.14)}$$

显然有

$$\begin{array}{rl}F(x_1)& ={\lambda }_{1}f(x_1)+{\lambda }_{2}f(x_1)-f\left({\lambda }_1{x}_1+{\lambda }_2{x}_1\right)\\ & =\left({\lambda }_1+{\lambda }_2\right)f(x_1)-f\left(\left({\lambda }_1+{\lambda }_2\right)x_1\right)\\ & =f(x_1)-f(x_1)\\ & =0\end{array}\text{\hspace{1em}(1.15)}$$

对$F(x)$求导，有

$$\begin{array}{rl}{F}^{{}^{\prime }}(x)& ={\lambda }_2{f}^{{}^{\prime }}(x)-{\lambda }_2{f}^{{}^{\prime }}({\lambda }_1{x}_1+{\lambda }_{2}x)\\ & ={\lambda }_2\left[f^{{}^{\prime }}(x)-f^{{}^{\prime }}({\lambda }_1{x}_1+{\lambda }_{2}x)\right]\end{array}\text{\hspace{1em}(1.16)}$$

因为$f(x)$是凸函数，所以恒有$f^{{}^{\prime \prime }}(x)⩾0$，这说明$f^{{}^{\prime }}(x)$单调递增。

所以

$$x-({\lambda }_1{x}_1+{\lambda }_{2}x)=(1-{\lambda }_2)x-{\lambda }_1{x}_1={\lambda }_1(x-x_1)\text{\hspace{1em}(1.17)}$$

当$x>x_1$时，$x>({\lambda }_1{x}_1+{\lambda }_{2}x)$，则$f^{{}^{\prime }}(x)>f^{{}^{\prime }}({\lambda }_1{x}_1+{\lambda }_{2}x)$，$F^{{}^{\prime }}(x)>0$；当$x<x_1$时，$x<({\lambda }_1{x}_1+{\lambda }_{2}x)$，则$f^{{}^{\prime }}(x)<f^{{}^{\prime }}({\lambda }_1{x}_1+{\lambda }_{2}x)$，$F^{{}^{\prime }}(x)<0$。

所以$F(x)$在$(a,x_1)$上单调递减，在$(x_1,b)$上单调递增。那么恒有

$$F(x)⩾F(x_1)=0\text{\hspace{1em}(1.18)}$$

即

$$f\left({\lambda }_1{x}_1+{\lambda }_2{x}_2\right)⩾{\lambda }_{1}f\left(x_1\right)+{\lambda }_{2}f\left(x_2\right)\text{\hspace{1em}(1.19)}$$

成立。

现在我们假设当$n=k$时，不等式成立，即

$$f\left(\sum _{i=1}^k{\lambda }_i{x}_i\right)⩾\sum _{i=1}^k{\lambda }_{i}f\left(x_i\right)\text{\hspace{1em}(1.20)}$$

则当$n=k+1$时，由$\sum _{i=1}^{k+1}{\lambda }_i=1$可得$\sum _{i=1}^k{\lambda }_i=1-{\lambda }_{k+1}$，从而有

$$\frac{1}{1-{\lambda }_{k+1}}\sum _{i=1}^k{\lambda }_i=1\text{\hspace{1em}(1.21)}$$

那么

$$\begin{array}{rl}\sum _{i=1}^{k+1}{\lambda }_{i}f(x_i)& =\sum _{i=1}^k{\lambda }_{i}f(x_i)+{\lambda }_{k+1}f(x_{k+1})\\ & =(1-{\lambda }_{k+1})\frac{1}{1-{\lambda }_{k+1}}\sum _{i=1}^k{\lambda }_{i}f(x_i)+{\lambda }_{k+1}f(x_{k+1})\\ & ⩾(1-{\lambda }_{k+1})f\left(\frac{\sum _{i=1}^k{\lambda }_i{x}_i}{1-{\lambda }_{k+1}}\right)+{\lambda }_{k+1}f(x_{k+1})\\ & ⩾f\left((1-{\lambda }_{k+1})\frac{\sum _{i=1}^k{\lambda }_i{x}_i}{1-{\lambda }_{k+1}}+{\lambda }_{k+1}{x}_{k+1}\right)\\ & =f\left(\sum _{i=1}^{k+1}{\lambda }_i{x}_i\right)\end{array}\text{\hspace{1em}(1.22)}$$

综上所述，$f\left(\sum _{i=1}^n{\lambda }_i{x}_i\right)⩽\sum _{i=1}^n{\lambda }_{i}f\left(x_i\right)$成立。

同理，若$f(x)$是区间$\left(a,b\right)$上的凹函数，有$f\left(\sum _{i=1}^n{\lambda }_i{x}_i\right)⩾\sum _{i=1}^n{\lambda }_{i}f\left(x_i\right)$成立。

那么我们注意到$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})$实际上就代表的是隐变量在一组数据下取遍所有可能的概率，所以有

$$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})=1\text{\hspace{1em}(1.23)}$$

那么我们就可以利用琴生不等式（Jensen inequality）来得到式$(1.10)$的下界：

$$\begin{array}{rl}L(\mathbf{\theta })-L({\mathbf{\theta }}^{(i)})& =\log \left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}\right)-\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})\\ & ⩾\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}-\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})\\ & =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}-\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})\\ & =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})}\end{array}\text{\hspace{1em}(1.24)}$$

令

$$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})\widehat{=}L({\mathbf{\theta }}^{(i)})+\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})}\text{\hspace{1em}(1.25)}$$

则

$$L(\mathbf{\theta })⩾B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})\text{\hspace{1em}(1.26)}$$

即函数$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$是$L(\mathbf{\theta })$的一个下界，而且由式$(1.25)$可知，

$$L({\mathbf{\theta }}^{(i)})=B({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})\text{\hspace{1em}(1.27)}$$

因此，任何可以使$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$增大的$\mathbf{\theta }$，也可以使$L(\mathbf{\theta })$增大。为了使$L(\mathbf{\theta })$有尽可能大的增长，选择${\mathbf{\theta }}^{(i+1)}$使$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$达到极大，即

$${\mathbf{\theta }}^{(i+1)}=\arg \underset{\mathbf{\theta }}{\max }B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})\text{\hspace{1em}(1.28)}$$

将式$(1.25)$代入，有

$$\begin{array}{rl}{\mathbf{\theta }}^{(i+1)}& =\arg \underset{\mathbf{\theta }}{\max }\left(L({\mathbf{\theta }}^{(i)})+\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})}\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})+\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})}\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})+\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})}\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \frac{P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \left[P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right]-\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\right)\end{array}\text{\hspace{1em}(1.29)}$$

注意到，$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})$对$\mathbf{\theta }$而言是常数，换言之，想让$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \left[P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right]-\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})$取最大值，其实就是让$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \left[P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right]$取最大值，即$(1.29)$可变换为

$$\begin{array}{rl}{\mathbf{\theta }}^{(i+1)}& =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \left[P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right]-\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log \left[P(\mathbf{Y}\mid \mathbf{Z},\mathbf{\theta })P(\mathbf{Z}\mid \mathbf{\theta })\right]\right)\\ & =\arg \underset{\mathbf{\theta }}{\max }\left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\right)\end{array}\text{\hspace{1em}(1.30)}$$

我们记

$$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})=\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\text{\hspace{1em}(1.31)}$$

如此，我们如下定义EM算法的核心：$Q$函数。

定义9.1（$\mathbf{Q}$函数）  完全数据的对数似然函数$\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })$关于在给定观测数据$\mathbf{Y}$和当前参数${\mathbf{\theta }}^{(i)}$下对未观测数据$\mathbf{Z}$的条件概率分布$\log P(\mathbf{Y},\mathbf{Z}\mid {\mathbf{\theta }}^{(i)})$的期望称为$Q$函数，即

$$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})=E_{\mathbf{Z}}\left[\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\mid \mathbf{Y},{\mathbf{\theta }}^{(i)}\right]\text{\hspace{1em}(1.32)}$$

于是式$(1.30)$等价于EM算法的一次迭代，即求$Q$函数极其极大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

图1.1给出了EM算法的直观解释。图中上方曲线为$L(\mathbf{\theta })$，下方曲线为$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$。由式$(1.26)$，两个函数在点$\mathbf{\theta }={\mathbf{\theta }}^{(i)}$处相等。由式$(1.28)$，式$(1.29)$和式$(1.30)$，EM算法找到下一个点${\mathbf{\theta }}^{(i+1)}$使函数$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$极大化，也使函数$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$极大化。这时由于$L(\mathbf{\theta })⩾B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$，函数$B(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$的增加，保证对数似然函数$L(\mathbf{\theta })$在每次迭代中也是增加的。EM算法在点${\mathbf{\theta }}^{(i+1)}$重新计算$Q$函数值，进行下一次迭代。在这个过程中，对数似然函数$L(\mathbf{\theta })$不断增大。从图中可以推断出EM算法不能保证找到全局最优值。

图1.1  EM算法的解释

二、EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。EM算法的最大优点是简单性和普适性。那么EM算法得到的估计序列是否收敛？如果收敛，是否收敛到全局最大值或局部极大值？

定理2.1    设$P(\mathbf{Y}\mid \mathbf{\theta })$为观测数据的似然函数，${\mathbf{\theta }}^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})(i=1,2,\cdots )$为对应的似然函数序列，则$P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})$是单调递增的，即

$$P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i+1)})⩾P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})\text{\hspace{1em}(2.1)}$$

证明    由于

$$P(\mathbf{Y}\mid \mathbf{\theta })=\frac{P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })}{P(\mathbf{Z}\mid \mathbf{Y},\mathbf{\theta })}$$

取对数有

$$\log P(\mathbf{Y}\mid \mathbf{\theta })=\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })-\log P(\mathbf{Z}\mid \mathbf{Y},\mathbf{\theta })$$

由式$(1.32)$

$$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})=\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })$$

令

$$H(\mathbf{\theta },{\mathbf{\theta }}^{(i)})=\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},\mathbf{\theta })\text{\hspace{1em}(2.2)}$$

注意到$\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})=1$，且对数似然函数不含有隐变量$\mathbf{Z}$，于是

$$\begin{array}{rl}\log P(\mathbf{Y}\mid \mathbf{\theta })& =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y}\mid \mathbf{\theta })\\ & =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\left(\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })-\log P(\mathbf{Z}\mid \mathbf{Y},\mathbf{\theta })\right)\\ & =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })-\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Z}\mid \mathbf{Y},\mathbf{\theta })\\ & =Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})-H(\mathbf{\theta },{\mathbf{\theta }}^{(i)})\end{array}\text{\hspace{1em}(2.3)}$$

在式$(2.3)$中分别取$\mathbf{\theta }$为${\mathbf{\theta }}^{(i)}$和${\mathbf{\theta }}^{(i+1)}$并相减，有

$$\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i+1)})-\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})=\left[Q({\mathbf{\theta }}^{(i+1)},{\mathbf{\theta }}^{(i)})-Q({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})\right]-\left[H({\mathbf{\theta }}^{(i+1)},{\mathbf{\theta }}^{(i)})-H({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})\right]\text{\hspace{1em}(2.4)}$$

为证式$(2.1)$，只需证式$(2.4)$右端是非负的。式$(2.4)$右端的第$1$项，由于${\mathbf{\theta }}^{(i+1)}$使$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$达到极大，所以有

$$Q({\mathbf{\theta }}^{(i+1)},{\mathbf{\theta }}^{(i)})-Q({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})⩾0\text{\hspace{1em}(2.5)}$$

其第$2$项，由式$(2.2)$可得：

$$\begin{array}{rl}H({\mathbf{\theta }}^{(i+1)},{\mathbf{\theta }}^{(i)})-H({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})& =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\left(\log \frac{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i+1)})}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}\right)\\ & ⩽\log \left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\frac{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i+1)})}{P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})}\right)\\ & =\log \left(\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i+1)})\right)\\ & =0\end{array}\text{\hspace{1em}(2.6)}$$

这里的不等号由上面提到的琴生不等式得到。
由式$(2.5)$和式$(2.6)$即知式$(2.3)$右端是非负的。

定理2.2    设$L(\mathbf{\theta })=P(\mathbf{Y}\mid \mathbf{\theta })$为观测数据的对数似然函数，${\mathbf{\theta }}^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$L({\mathbf{\theta }}^{(i)})(i=1,2,\cdots )$为对应的对数似然函数序列。
（1）如果$P(\mathbf{Y}\mid \mathbf{\theta })$有上界，则$L({\mathbf{\theta }}^{(i)})=\log P(\mathbf{Y}\mid {\mathbf{\theta }}^{(i)})$收敛到某一值$L^{\ast }$；
（2）在函数$Q(\mathbf{\theta },{\mathbf{\theta }}^{{}^{\prime }})$与$L(\mathbf{\theta })$满足一定条件下，由EM算法得到的参数估计序列${\mathbf{\theta }}^{(i)}$的收敛值${\mathbf{\theta }}^{\ast }$是$L(\mathbf{\theta })$的稳定点。
证明    （1）由$L(\mathbf{\theta })=\log P(\mathbf{Y}\mid \mathbf{\theta })$的单调性及$P(\mathbf{Y}\mid \mathbf{\theta })$的有界性立即得到。
（2）参见论文《On the Convergence Properties of the EM Algorithm》。

定理2.2关于函数$Q(\mathbf{\theta },{\mathbf{\theta }}^{{}^{\prime }})$与$L(\mathbf{\theta })$的条件在大多数情况下都是满足的。EM算法的收敛性包含关于对数似然函数序列$L({\mathbf{\theta }}^{(i)})$的收敛性和关于参数估计序列${\mathbf{\theta }}^{(i)}$的收敛性两层意思，前者并不蕴涵后者。此外，定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点。所以在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

三、EM算法

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

输出：观测模型参数$\mathbf{\theta }$。

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

（2）E步：记${\mathbf{\theta }}^{(i)}$为第$i$次迭代参数$\mathbf{\theta }$的估计值，在第$i+1$次迭代的E步，计算

$$\begin{array}{rl}Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})& =E_{\mathbf{Z}}\left[\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\mid \mathbf{Y},{\mathbf{\theta }}^{(i)}\right]\\ & =\sum _{\mathbf{Z}}P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})\log P(\mathbf{Y},\mathbf{Z}\mid \mathbf{\theta })\end{array}\text{\hspace{1em}(3.1)}$$

这里，$P(\mathbf{Z}\mid \mathbf{Y},{\mathbf{\theta }}^{(i)})$是在给定观测数据$\mathbf{Y}$和当前的参数估计${\mathbf{\theta }}^{(i)}$下隐变量数据$\mathbf{Z}$的条件概率分布；

（3）M步：求使$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$极大化的$\mathbf{\theta }$，确定第$i+1$次迭代的参数的估计值${\mathbf{\theta }}^{(i+1)}$

$${\mathbf{\theta }}^{(i+1)}=\arg \underset{\mathbf{\theta }}{\max }Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})\text{\hspace{1em}(3.2)}$$

（4）重复第（2）步和第（3）步，直到收敛。

下面关于EM算法作几点说明：

步骤（1）    参数的初值 可以任意选择，但需注意EM算法对初值是敏感的。

步骤（2）    E步求$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$。$Q$函数式中$\mathbf{Z}$是未观测数据，$\mathbf{Y}$是观测数据。注意，$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求$Q$函数及其极大。

步骤（3）    M步求$Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$的极大化，得到${\mathbf{\theta }}^{(i+1)}$，完成一次迭代${\mathbf{\theta }}^{(i)}\to {\mathbf{\theta }}^{(i+1)}$。后面将证明每次迭代使似然函数增大或达到局部极值。

步骤（4）    给出停止迭代的条件，一般是对较小的正数${\epsilon }_1$，${\epsilon }_2$，若满足

$$\Vert {\mathbf{\theta }}^{(i+1)}-{\mathbf{\theta }}^{(i)}\Vert <{\epsilon }_1或\Vert Q({\mathbf{\theta }}^{(i+1)},{\mathbf{\theta }}^{(i)})-Q({\mathbf{\theta }}^{(i)},{\mathbf{\theta }}^{(i)})\Vert <{\epsilon }_2$$

则停止迭代。

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

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

1. 高斯混合模型

定义4.1（高斯混合模型）    高斯混合模型是指具有如下形式的概率分布模型：

$$P(y\mid \theta )=\sum _{k=1}^K{\alpha }_k\varphi \left(y\mid {\theta }_k\right)\text{\hspace{1em}(4.1)}$$

其中，${\alpha }_k$是系数，${\alpha }_k⩾0$，$\sum _{k=1}^K{\alpha }_k=1$；$\varphi \left(y\mid {\theta }_k\right)$是高斯分布密度，${\theta }_k=\left({\mu }_k,{\sigma }_k^2\right)$，

$$\varphi \left(y\mid {\theta }_k\right)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{{\left(y-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)\text{\hspace{1em}(4.2)}$$

称为第$k$个分模型。

一般混合模型可以由任意概率分布密度代替式$(4.2)$中的高斯分布密度，例如：01分布，二项分布，泊松分布等等。

高斯混合模型本质上就是将多个高斯分布产生的数据放在一起，第$k$个模型产生的数据的数量占总数据数量的比例就是${\alpha }_k$，现在这些数据我们都不知道是由哪一个高斯分布产生的，EM算法的作用就是进行划分，分清楚每一个数据到底是由哪一个高斯分布产生的。

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

假设观测数据$y_1,y_2,\cdots ,y_N$由高斯混合模型生成，那么

$$P(y_j\mid \theta )=\sum _{k=1}^K{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right),j=1,2,\cdots ,N\text{\hspace{1em}(4.3)}$$

其中，$\theta =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K;{\theta }_1,{\theta }_2,\cdots ,{\theta }_K\right)$。我们用EM算法估计高斯混合模型的参数$\theta$。

1. 明确隐变量，写出完全数据的对数似然函数

可以设想观测数据$y_j$，$j=1,2,\cdots ,N$，是这样产生的：首先依概率${\alpha }_k$选择第$k$个高斯分布分模型$\varphi \left(y\mid {\theta }_k\right)$，然后依第$k$个分模型的概率分布$\varphi \left(y\mid {\theta }_k\right)$生成观测数据$y_j$。这时观测数据$y_j$，$j=1,2,\cdots ,N$是已知的；反映观测数据$y_j$来自第$k$个分模型的数据是未知的，$k=1,2,\cdots ,K$，以隐变量${\gamma }_{jk}$表示，其定义如下：

$$\begin{gathered} {\gamma }_{jk}=\{\begin{array}{ll}1& ,\text{第j个观测来自第k个分模量}\\ 0& ,\text{否则}\end{array} \\ j=1,2,\cdots ,N;k=1,2,\cdots ,K\text{\hspace{1em}(4.4)} \end{gathered}$$

${\gamma }_{jk}$是0-1随机变量。

有了观测数据$y_j$及未观测数据${\gamma }_{jk}$，那么完全数据是

$$\left(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}\right),j=1,2,\cdots ,N$$

于是，可以写出完全数据的似然函数：

$$\begin{array}{rl}P(y,\gamma \mid \theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}\mid \theta )\\ & =\prod _{k=1}^K\prod _{j=1}^N{\left[{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\varphi \left(y_j\mid {\theta }_k\right)\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{{\left(y_j-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\end{array}$$

式中，$n_k=\sum _{j=1}^N{\gamma }_{jk}$，由于对于每个数据，其实只能是由一个模型生成的，所以对于${\gamma }_{jk}$，相同的$j$下，只能有一个$k$使${\gamma }_{jk}=1$，其余的都为$0$，所以有$\sum _{k=1}^K{\gamma }_{jk}=1$，那么$\sum _{k=1}^K{n}_k=\sum _{k=1}^K\sum _{j=1}^N{\gamma }_{jk}=\sum _{j=1}^N\sum _{k=1}^K{\gamma }_{jk}=N$。

那么，完全数据的对数似然函数为

$$\log P(y,\gamma \mid \theta )=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\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\}$$

2. EM算法的E步：确定$Q$函数

$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E\left[\log P(y,\gamma \mid \theta )\mid y,{\theta }^{(i)}\right]\\ & =E\left\{\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\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\}\right\}\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^N\left(E_{{\gamma }_{jk}}\right)\log {\alpha }_k+\sum _{j=1}^N\left(E_{{\gamma }_{jk}}\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{array}\text{\hspace{1em}(4.5)}$$

这里需要计算$E\left({\gamma }_{jk}\mid y,\theta \right)$，记为${\hat{\gamma }}_{jk}$。

$$\begin{array}{rl}{\hat{\gamma }}_{jk}& =E\left({\gamma }_{jk}\mid y,\theta \right)=P\left({\gamma }_{jk}=1\mid y,\theta \right)\\ & =\frac{P\left({\gamma }_{jk}=1,y_j\mid \theta \right)}{\sum _{k=1}^{K}P\left({\gamma }_{jk}=1,y_j\mid \theta \right)}\\ & =\frac{P\left(y_j\mid {\gamma }_{jk}=1,\theta \right)P\left({\gamma }_{jk}=1\mid \theta \right)}{\sum _{k=1}^{K}P\left(y_j\mid {\gamma }_{jk}=1,\theta \right)P\left({\gamma }_{jk}=1\mid \theta \right)}\\ & =\frac{{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)}{\sum _{k=1}^K{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)},j=1,2,\cdots ,N;k=1,2,\cdots ,K\end{array}$$

${\hat{\gamma }}_{jk}$是在当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_j$的响应度。

将${\hat{\gamma }}_{jk}=E_{{\gamma }_{jk}}$及$n_k=\sum _{j=1}^N{E}_{{\gamma }_{jk}}$代入式$(4.5)$，即得

$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}\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\}\text{\hspace{1em}(4.6)}$$

3. 确定EM算法的M步

迭代的M步是求函数$Q(\theta ,{\theta }^{(i)})$对$\theta$的极大值，即求新一轮迭代的模型参数：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

用${\hat{\mu }}_k$，${\hat{\sigma }}_k^2$及${\hat{\alpha }}_k$，$k=1,2,\cdots ,K$，表示${\theta }^{(i+1)}$的各参数。求${\hat{\mu }}_k$，${\hat{\sigma }}_k^2$只需将式$(4.6)$分别对${\hat{\mu }}_k$，${\hat{\sigma }}_k^2$求偏导数并令其为$0$，即可得到；求${\hat{\alpha }}_k$是在$\sum _{k=1}^K{\alpha }_k=1$条件下求偏导数并令其为$0$得到的。j结果如下：

$${\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\text{\hspace{1em}(4.7)}$$

$${\hat{\sigma }}_k^2=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{\sum _{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\text{\hspace{1em}(4.8)}$$

$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K\text{\hspace{1em}(4.9)}$$

重复以上计算，直到对数似然函数值不再有明显的变化为止。

如此我们便得到了高斯混合模型参数估计的EM算法：

算法4.1（高斯混合模型参数估计的EM算法）

输入：观测数据$y_1,y_2,\cdots ,y_N$，高斯混合模型；

输出：高斯混合模型参数。

（1）取参数的初始值开始迭代；

（2）E步：依据当前模型参数，计算分模型$k$对观测数据$y_j$的响应度

$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)}{\sum _{k=1}^K{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)},j=1,2,\cdots ,N;k=1,2,\cdots ,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$$

$${\hat{\sigma }}_k^2=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{\sum _{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K$$

$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K$$

（4）重复第（2）步和第（3）步，直到收敛。

五、代码实现

该代码实现了高斯混合模型参数估计的EM算法，选取了三个二维高斯分布模型${\mathbf{\mu }}_1=\left(\begin{array}{c}2\\ 1\end{array}\right)$，${\mathbf{\sigma }}_1^2=\left(\begin{array}{cc}0.5& 0\\ 0& 0.2\end{array}\right)$；${\mathbf{\mu }}_2=\left(\begin{array}{c}5\\ 9\end{array}\right)$，${\mathbf{\sigma }}_2^2=\left(\begin{array}{cc}0.3& 0\\ 0& 0.3\end{array}\right)$；${\mathbf{\mu }}_3=\left(\begin{array}{c}10\\ 6\end{array}\right)$，${\mathbf{\sigma }}_3^2=\left(\begin{array}{cc}0.3& 0\\ 0& 0.4\end{array}\right)$，每个模型都随机生成100个点，然后用EM算法对这300个点进行划分和聚类，其中初值设定为

$$\alpha =\left(\begin{array}{c}0.3\\ 0.3\\ 0.4\end{array}\right),\mu =\left(\begin{array}{cc}2.0+random& 1.0+random\\ 5.0+random& 9.0+random\\ 10.0+random& 6.0+random\end{array}\right),{\sigma }^2=\left(\begin{array}{cc}0.5+random& 0.2+random\\ 0.3+random& 0.3+random\\ 0.3+random& 0.4+random\end{array}\right)$$

然后进行迭代，最终得到结果如下：

图5.1  迭代结果 图5.2  EM算法对高斯混合模型的实现

完整代码如下：

import numpy as np
import random
import math
import matplotlib.pyplot as plt


# 创造数据
def create_data(mu1, mu2, sigma1, sigma2, size):
    mean = np.array([mu1, mu2])
    cov = np.array([[sigma1, 0], [0, sigma2]])
    Y = np.random.multivariate_normal(mean, cov, size=size)
    return Y


# 高斯分布模型
def Gau(y, mu, sigma):
    tmp1 = -0.5 * ((y[0] - mu[0]) ** 2 / sigma[0] + (y[1] - mu[1]) ** 2 / sigma[1])
    tmp2 = 1.0 / (2 * np.pi * math.sqrt(sigma[0]) * math.sqrt(sigma[1]))
    res = tmp2 * math.exp(tmp1)
    return res


# EM算法
def EM(y, K):
    N = len(y)
    alpha = np.array([0.3, 0.3, 0.4])
    mu = np.array([[2.0 + random.random(), 1.0 + random.random()],
                   [5.0 + random.random(), 9.0 + random.random()],
                   [10.0 + random.random(), 6.0 + random.random()]])
    sigma = np.array([[0.5 + random.random(), 0.2 + random.random()],
                      [0.3 + random.random(), 0.3 + random.random()],
                      [0.3 + random.random(), 0.4 + random.random()]])
    gama = []
    for i in range(N):
        gama.append([0.0, 0.0, 0.0])
    gama = np.array(gama)
    gama1 = gama
    t = 1

    while True:
        # E步
        for j in range(N):
            tmp = 0.0
            for k in range(K):
                tmp = tmp + alpha[k] * Gau(y[j], mu[k], sigma[k])
            for k in range(K):
                gama[j, k] = alpha[k] * Gau(y[j], mu[k], sigma[k]) / tmp

        # 判断
        if t == 2:
            sum = 0.0
            for j in range(N):
                sum += np.linalg.norm(gama[j] - gama1[j])
            if sum < 0.0000001:
                break

        # M步
        tmp_gama = np.array([0.0, 0.0, 0.0])
        for k in range(K):
            for j in range(N):
                tmp_gama[k] = tmp_gama[k] + gama[j, k]

        for k in range(K):
            tmp = np.array([0.0, 0.0])
            for j in range(N):
                tmp = tmp + gama[j, k] * y[j]
            mu[k] = tmp / tmp_gama[k]

        for k in range(K):
            tmp = np.array([0.0, 0.0])
            for j in range(N):
                tmp = tmp + gama[j, k] * (y[j] - mu[k]) ** 2
            sigma[k] = tmp / tmp_gama[k]

        for k in range(K):
            alpha[k] = tmp_gama[k] / N

        gama1 = gama
        t = 2

    return N, alpha, mu, sigma, gama


if __name__ == '__main__':
    # K = 3
    # n = 100, N = 300
    s = np.array([3, 100, 2, 1, 0.5, 0.2, 5, 9, 0.3, 0.3, 10, 6, 0.3, 0.4])
    lens = len(s)
    for i in range(2, lens):
        s[i] = np.float_(s[i])
    K = np.int_(s[0])
    size = np.int_(s[1])
    Y = []
    for i in range(2, 14, 4):
        tmp_Y = create_data(s[i], s[i + 1], s[i + 2], s[i + 3], size)
        Y.extend(tmp_Y)
    Y = np.array(Y)

    N, alpha, mu, sigma, gama = EM(Y, K)

    print("α =", alpha)
    print("μ =", mu)
    print("σ^2 =", sigma)

    for j in range(N):
        index = 0
        big = gama[j, 0]
        if big < gama[j, 1]:
            index = 1
            big = gama[j, 1]
        if big < gama[j, 2]:
            index = 2

        if index == 0:
            c = 'r'
        elif index == 1:
            c = 'g'
        else:
            c = 'b'
        plt.scatter(Y[j, 0], Y[j, 1], c=c, marker='+')

    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()