写在前面

EM（Expectation Maximization 期望最大化）算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。其每次迭代由E、M两步构成。下面首先给出一般EM算法的求解过程（怎么做），然后结合一个例子来理解，然后讲为什么这么求解，即推导，最后讲述EM算法在高斯混合模型中的应用及小结。

EM算法

一般用Y表示观测随机变量，Z表示隐随机变量，Y和Z在一起称为完全数据，Y称为不完全数据。由于含有隐变量，我们不能直接由$P(Y|\theta)=\sum_{z}P(Z|\theta)P(Y|Z,\theta)$的最大似然估计来得到模型参数$\theta$。EM算法就是在给定Y和Z的联合概率分布为$P(Y,Z|\theta)$的情况下，通过迭代求解$L(\theta)=lnP(Y|\theta)$的极大似然估计来估算模型参数的算法。求解步骤如下：

1 选择一个初始的参数$\theta^{old}$。

2 E Step 估计$P(Z|Y,\theta^{old})$。

3 M Step 估计$\theta^{new}$

$$\theta^{new}=arg\max_{\theta}Q(\theta,\theta^{old})$$

其中

$$Q(\theta,\theta^{old})=\sum_{z}P(Z|Y,\theta^{old})logP(Y,Z|\theta)$$

4 检查是否到达停止迭代条件，一般是对较小的正数$\varepsilon_{1},\varepsilon_{2}$，若满足：

$$\lVert \theta^{new}-\theta^{old}\rVert<\varepsilon_{1} 或 \lVert Q(\theta^{new}-\theta^{old})-Q(\theta^{old}-\theta^{old})\rVert<\varepsilon_{2}$$

则停止迭代，否则$\theta^{old} \leftarrow \theta^{new}$转到步骤2继续迭代。

一个栗子

下面结合一个《统计学习方法》中的例子来加深下理解：

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

假设只能观察到投掷硬币的结果，而不知其过程，问如何估计三硬币正面出现的概率，即三硬币的模型参数$\pi$,$p$和$q$。

解答：我们现在只可以看到硬币最终的结果1111110000，即观测变量Y，至于结果来自于B还是C无从得知，我们设隐变量Z来表示来自于哪个变量，令$\theta=\{\pi,p,q\}$。观测数据的似然函数可以表示为：

$$P(Y|\theta)=\sum_{z}P(Z|\theta)P(Y|Z,\theta)=\prod_{j=1}^{n}[\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi)q^{y_{i}}(1-q)^{1-y_{j}}]$$

则模型参数$\theta$的最大log似然估计为：

$$\hat{\theta}=arg\max_{\theta}logP(Y|\theta)$$

由于隐变量的存在，使得观测数据的最大似然函数里log里有带有加和($\pi..+(1-\pi)..$)，导致上式是没有解析解的，只能通过迭代的方法求解，EM算法就是用于求解此类问题的一种迭代算法。

根据第二部分EM算法求解过程，假设第$i$次迭代参数的估计值为$\theta^{(i)}=(\pi^{(i)},p^{(i)},q^{(i)})$，EM算法的第i+1次迭代如下：

E步：估计$P(Z|Y,\theta^{(i)})$,即在模型参数$\theta^{(i)}$下观测数据$y_{j}$来自硬币B的概率

$$\mu^{(i+1)}=P(Z|Y,\theta^{(i)})=\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步：估计$\theta^{(i+1)}$,即：

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

其中：

\begin{equation} \begin{split} Q(\theta,\theta^{(i)})& =E_{z}[logP(Y,Z|\theta)|Y,\theta^{(i)}]\\ & =\sum_{z}P(Z|Y,\theta^{(i)})logP(Y,Z|\theta)\\ & =\sum_{j=1}^{n}[\mu^{(i+1)}log\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\mu^{(i+1)})log(1-\pi)q^{y_{i}}(1-q)^{1-y_{j}}] \end{split} \end{equation}

$\mu^{(i+1)}$是E步得到常数，可通过分别对参数$\pi$,$p$和$q$求偏导使其为零来最大化上式，获得$\pi$,$p$和$q$的新的估计值：

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

在选定参数初始值$\theta^{(0)}$后，根据E步，M步循环迭代，直至满足迭代停止条件，即可得到参数$\theta$的极大似然估计。

由上述EM计算过程和结合例子的应用，相信大家都会用EM算法解决问题了，即知道怎么做了，下面一节主要来讲述为什么这样做，即为什么这样就可以解决此类含有隐变量的最大似然估计。

EM算法的导出

面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据），即极大化

$$L(\theta)=logP(Y|\theta)=log\sum_{z}P(Y,Z|\theta)=log[\sum_{z}P(Z|\theta)P(Y|Z,\theta)]$$

极大化的主要困难是上式中含有未观测数据log里有包含和，EM算法则是通过迭代逐步近似极大化$L(\theta)$的。假设在第i次迭代后$\theta$的估计值是$\theta^{(i)}$，我们希望新估计的值$\theta$能使$L(\theta)$增加，即$L(\theta)>L(\theta^{(i)})$，并逐步达到极大值，考虑两者的差：

\begin{equation} \begin{split} L(\theta)-L(\theta^{(i)})& =log[\sum_{z}P(Z|\theta)P(Y|Z,\theta)]-logP(Y|\theta^{(i)})\\ & =log[\sum_{z}P(Z|Y,\theta^{(i)})\frac{P(Z|\theta)P(Y|Z,\theta)}{P(Z|Y,\theta^{(i)})}]-logP(Y|\theta^{(i)})\\ & \ge\sum_{z}P(Z|Y,\theta^{(i)})log\frac{P(Z|\theta)P(Y|Z,\theta)}{P(Z|Y,\theta^{(i)})}-logP(Y|\theta^{(i)})\\ &=\sum_{z}P(Z|Y,\theta^{(i)})log\frac{P(Z|\theta)P(Y|Z,\theta)}{P(Y|\theta^{(i)})P(Z|Y,\theta^{(i)})} \end{split} \end{equation}

上式中的不等号是由Jensen不等式得到，下面对Jesen不等式做一个简单回顾。

Jensen不等式：如果f是凸函数，X是随机变量，则$E[f(X)] \ge f(EX)$，特别地，如果f是严格凸函数，那么当$E[f(X)]=f(EX)$当且仅当$P(x=E[X])=1$，即X为常量。凹函数的性质和凸函数相反。

上式中$f(X)$为$log(X)$为凹函数，则$E[f(X)] \le f(EX)$，是上式中不等号的由来。

继续转回正题。令

$$B(\theta,\theta^{(i)}) \triangleq L(\theta^{(i)})+\sum_{z}P(Z|Y,\theta^{(i)})log\frac{P(Z|\theta)P(Y|Z,\theta)}{P(Z|\theta^{(i)})P(Z|Y,\theta^{(i)})}$$

则

$$L(\theta) \ge B(\theta,\theta^{(i)})$$

即$B(\theta,\theta^{(i)})$是$L(\theta)$的下界，且$L(\theta^{(i)})=B(\theta^{(i)},\theta^{(i)})$

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

$$\theta^{(i+1)}=arg\max_{\theta}B(\theta,\theta^{(i)})$$

省去对$\theta$极大化而言是常熟的项，有

\begin{equation} \begin{split} \theta^{(i+1)}& =arg\max_{\theta}[L(\theta^{(i)})+\sum_{z}P(Z|Y,\theta^{(i)})log\frac{P(Z|\theta)P(Y|Z,\theta)}{P(Z|\theta^{(i)})P(Z|Y,\theta^{(i)})}]\\ & =arg\max_{\theta}[\sum_{z}P(Z|Y,\theta^{(i)})logP(Z|\theta)P(Y|Z,\theta)]\\ & =arg\max_{\theta}[\sum_{z}P(Z|Y,\theta^{(i)})logP(Y,Z|\theta)]\\ & =arg\max_{\theta}Q(\theta,\theta^{(i)}) \end{split} \end{equation}

推导完毕。

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

EM算法可以用来估计高斯混合模型中的参数，例如：

假设观测数据$y_{1},y_{2},...,y_{N}$由高斯混合模型生成

$$P(y|\theta)=\sum_{k=1}^{K}\alpha_{k}\phi(y|\theta_{k})$$

其中$\theta=(\alpha_{1},\alpha_{2},...,\alpha_{k};\theta_{1},\theta_{2},...,\theta_{k})$。

可以设想观测数据$y_{j},j=1,2,...,N$是这样产生，先根据概率$\alpha_{k}$选择第$k$个高斯分布$\phi(y|\theta_{k})$，然后根据第$k$个模型的概率分布$\phi(y|\theta_{k})$生成观测数据$y_{j}$。此时观测数据$y_{j}$是已知的，反映观测数据$y_{j}$来自哪个分模型是未知的，看作隐变量。可以看出，可以用EM算法估计高斯混合模型（含有隐变量的概率模型参数）的参数$\theta$。

1 首先确定E步，估计$P(Z|Y,\theta^{(i)})$，即在已知第$i$次迭代参数的情况下，观测数据$y_{j}$来自计算分模型$k$的概率。

$$\gamma_{jk}^{(i+1)}=P(Z|Y,\theta^{(i)})=\frac{\alpha_{k}^{(i)}\phi(y|\theta_{k}^{(i)})}{\sum_{k=1}^{K}\alpha_{k}^{(i)}\phi(y|\theta_{k}^{(i)})}, \,\, j=1,2,...,N; \, k=1,2,...,K$$

2 M步，将新估计的Z的概率代进最大似然的公式，对待估计参数分别求偏导，以计算新一轮迭代模型参数（详细的推导不再赘述，感兴趣的可以自行推导）。

$$\mu_{k}^{(i+1)}=\frac{\sum_{j=1}^{N}\gamma_{jk}^{(i+1)}y_{j}}{\sum_{j=1}^{N}\gamma_{jk}^{(i+1)}}, \,\, k=1,2,...,K$$

$$(\sigma_{k}^{2})^{(i+1)}=\frac{\sum_{j=1}^{N}\gamma_{jk}^{(i+1)}(y_{j}-\mu_{k}^{(i+1)})^{2}}{\sum_{j=1}^{N}\gamma_{jk}^{(i+1)}}, \,\, k=1,2,...,K$$

$$\alpha_{k}^{(i+1)}=\frac{\sum_{j=1}^{N}\gamma_{jk}^{(i+1)}}{N}, \,\, k=1,2,...,K$$

重复E步，M步直至收敛。

小结

1 EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。含有隐变量的概率数据表示为$P(Y,Z|\theta)$，Y表示观测变量，Z是隐变量，$\theta$是模型参数。EM算法通过迭代求解观测数据的对数似然函数$L(\theta)=log(P|\theta)$的极大化来实现极大似然估计。每次迭代包含两步：

E步，求解$P(Z|Y,\theta^{old})$；

M步，估计$\theta^{new}$

$$\theta^{new}=arg\max_{\theta}Q(\theta,\theta^{old})$$

2 EM算法应用极其广泛，主要用于含有隐变量的概率模型的学习，但其对参数初始值比较敏感，而且不能保证收敛到全局最优。

参考文献

[1] 统计学习方法.李航

[2] Pattern Recognition And Machine Learning.Christopher M. Bishop

[3] The EM Algorithm,JerryLead's Blog

[4] 三硬币问题