【ML】EM（期望最大）算法

  EM算法（Expectation Maximization Algorithm，期望极大算法）是一种解决优化问题的迭代算法，用于求解含有隐变量的概率模型参数的极大似然估计（MLE）或极大后验概率估计（MAP）。EM算法是一种比较通用的参数估计算法，被广泛用于支持向量机（SMO算法）、朴素贝叶斯、GMM（高斯混合模型）、K-means（K均值聚类）和HMM（隐马尔可夫模型）的参数估计。

理解EM算法（例子）

  在统计学中，概率用于在已知一些参数的情况下，预测接下来的所得到的结果；而似然性则是用于在已知某些观测所得到的结果时，对有关事物的性质的参数进行估计。
  EM算法和极大似然估计的前提是一样的，都要假设数据总体的分布，如果不知道数据分布，是无法使用EM算法的。

三硬币模型

  假设有3枚硬币A，B，C，这些硬币正面出现的概率分别是$\pi$，$p$和$q$。进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，正面记作1，反面记作0；独立重复n此试验，观测结果：
$$1,1,0,1,0,0,1,0,1,1$$
假设只能观测到掷硬币的结果，不能观测掷硬币的过程，如何估计三硬币出现正面的概率？
  该模型的生成模型可写作
$$\begin{array}{rl}P(x|\theta )& =\sum _{z}P(x,z|\theta )=\sum _{z}P(z|\theta )P(x|z,\theta )\\ & =\pi p^x(1-p{)}^{1-x}+(1-\pi )q^x(1-q{)}^{1-x}\end{array}$$
其中，随机变量$X$是观测变量，表示依次试验观测结果是1或0；随机变量$Z$是隐变量，表示未观测到的掷硬币A的结果；$\theta =(\pi ,p,q)$是模型参数。

  设观测变量数据$X=\left(x_1,x_2,\dots ,x_n\right)$，隐变量数据$Z=\left(z_1,z_2,\dots ,z_n\right)$，$(X,Z)$称为完全数据。

1. 极大似然估计

  定义：
$$L(\theta )=\prod _{i=1}^{n}p\left(x_i|\theta \right)=\prod _{i=1}^n\sum _{z}P(z|\theta )P(x_i|z,\theta )$$
  $L(\theta )$反映的是在不同的参数$\theta$取值下，取得当前这个样本集的可能性，称为参数$\theta$相对于样本集X的似然函数。对似然函数取对数，将其变成连加的，称为对数似然函数：
$$\log L(\theta )=\log \prod _{i=1}^{n}p\left(x_i|\theta \right)=\sum _{i=1}^n\log p\left(x_i|\theta \right)$$

1. 取对数之后累积变为累和，求导更加方便；
2. 概率累积会出现数值非常小的情况，比如1e-30，由于计算机的精度是有限的，无法识别这一类数据，取对数之后，更易于计算机的识别。

  极大似然估计是，已经知道了结果，然后寻求使该结果出现的可能性极大的条件，以此作为估计值。即：
$$\hat{\theta }=\mathrm{argmax}\log L(\theta )$$
$\hat{\theta }$为$\theta$的极大似然估计量。

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

2.EM算法

(1)当我们知道了掷硬币$A$的结果，则可以确定下一步掷硬币B还是硬币C，利用极大似然很容易$(\pi ,p,q)$进行估计。
(2)当我们知道了$(\pi ,p,q)$，才有可能推断出每一个样本点对应的隐变量数据。

  EM 算法的思路是使用启发式的迭代方法，先猜想隐含参数（EM 算法的 E 步），接着基于观察数据和猜测的隐含参数一起来极大化对数似然，求解模型参数（EM算法的M步)。由于我们之前的隐含参数是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。我们基于当前得到的模型参数，继续猜测隐含参数（EM算法的 E 步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

EM算法流程

输入：观察数据$x=\left(x^{(1)},x^{(2)},\dots x^{(m)}\right)$，联合分布$p(x,z|\theta )$，条件分布$p(z|x,\theta )$，极大迭代次数$J$，其中z为隐变量。
输出：模型参数$\theta$。

1. 记第j次迭代后参数$\theta$的估计值为${\theta }^{(j)}$，随机初始化模型参数$\theta$的初值${\theta }^{(0)}$；
2. E步（求期望。固定模型参数的值，优化隐含数据的分布）
   第j+1次迭代：
   计算隐变量数据z的条件概率分布：
   $$Q_j\left(z^{(i)}\right):=P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right),z^{(i)}\in \mathcal{Z},i=1,2,...,m$$
   其中$\mathcal{Z}$为隐变量z的取值空间。
   计算期望：
   $$Q(\theta ,{\theta }^{(j)})=\sum _{i=1}^m\sum _{z^{(i)}\in \mathcal{Z}}{Q}_j\left(z^{(i)}\right)\log P\left(x^{(i)},z^{(i)}|\theta \right)$$
3. M步（取极大。固定隐含数据分布，优化模型参数）
   极大化$Q(\theta ,{\theta }^{(j)})$得到$\theta$：
   $${\theta }^{(j+1)}:=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(j)})$$
4. 重复E步、M步直至$\theta$收敛，停止迭代。

停止迭代条件：对较小的正数${ϵ}_1$、${ϵ}_2$满足$\Vert {\theta }^{(j+1)}-{\theta }^{(j)}\Vert <{ϵ}_1$或$\Vert L\left({\theta }^{(j+1)}\right)-L\left({\theta }^{(j)}\right)\Vert <{ϵ}_2$。

$Q$函数： 完全数据的对数似然函数$\log P(X,Z|\theta )$关于在给定观测数据$X$和当前参数${\theta }^{(j)}$下对隐变量数据$Z$的条件概率分布$P(Z|X,{\theta }^{(j)})$的期望称为$Q$函数，即
$$Q(\theta ,{\theta }^{(j)})=E_Z[\log P(X,Z|\theta )|X,{\theta }^{(j)}]$$

  总结来说，要计算$\theta$极大化对数似然函数$\log p(X|\theta )$，EM算法的思路是：考虑完全数据集$X,Z$的对数似然函数$\log p(X,Z|\theta )$在隐变量的后验分布$p(Z|X,{\theta }^{old})$下的期望，反复迭代极大化这个期望，即Q函数，最终得到模型的参数值。

EM算法的推导

  对于$m$个相互独立的样本$x=\left(x^{(1)},x^{(2)},\dots x^{(m)}\right)$，极大化对数函数找到模型参数$\theta$：
$$\theta =\arg \underset{\theta }{\max }\sum _{i=1}^m\log P\left(x^{(i)}|\theta \right)$$
  如果得到的观测数据有未观测到的隐变量数据$z=\left(z^{(1)},z^{(2)},\dots z^{(m)}\right)$，我们的目标变成了找到合适的$\theta$让对数似然函数极大：
$$\theta =\arg \underset{\theta }{\max }\sum _{i=1}^m\log P\left(x^{(i)}|\theta \right)=\arg \underset{\theta }{\max }\sum _{i=1}^m\log \sum _{z^{(i)}\in \mathcal{Z}}P\left(x^{(i)},z^{(i)}|\theta \right)$$

  如果对分别对未知的$\theta$和z分别求偏导，求导后形式会非常复杂，所以很难求解得到$\theta$和z。考虑引入一个与$z^{(i)}$、$\theta$相关的分布$Q_{\theta }\left(z^{(i)}\right)$满足：
$$\sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }(z^{(i)})=1,0\le Q_{\theta }(z^{(i)})\le 1$$
应用Jensen不等式对这个式子进行放缩：
$$\begin{array}{rl}\sum _{i=1}^m\log \sum _{z^{(i)}\in \mathcal{Z}}P\left(x^{(i)},z^{(i)}|\theta \right)& =\sum _{i=1}^m\log \sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}\\ & \ge \sum _{i=1}^m\sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)\log \frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}\end{array}$$

  对数函数是凹函数，满足Jensen不等式： $$log(E[y])\ge E[log(y)]$$
其中$$\{\begin{array}{l}E[y]={\sum }_i{p}_i{y}_i,p_i\ge 0,{\sum }_i{p}_i=1\\ p_i=Q_{\theta }\left(z^{(i)}\right)\\ \\ y_i=\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}\end{array}$$

  这样，对于任意一种分布$Q_{\theta }(z^{(i)})$，都给出了目标函数的一个下界，这个下界是$log\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}$的加权平均，也是我们所说的期望，这就是Expectation的来历。接下来的思路是：选择一个合适的$Q_{\theta }(z^{(i)})$，使其给出的下界与目标函数值最接近，然后就有理由通过极大化这个下界来极大化目标函数值。

$Q_{\theta }(z^{(i)})$的选择：假设$\theta$已经给定

  此时，目标函数的下界可由$Q_{\theta }\left(z^{(i)}\right)$和$P\left(x^{(i)},z^{(i)}|\theta \right)$完全确定。我们可以通过调整$Q_{\theta }\left(z^{(i)}\right)$使下界逼近目标函数的真实值，当不等式变成等式时，说明我们调整后的下界就是目标函数真实值了。由 Jensen不等式可知，等式成立的条件是随机变量是常数，则有：
$$\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}=c$$
其中c为不依赖$z^{(i)}$的常数。于是对于任意$i$我们得到：
$$P\left(x^{(i)},z^{(i)}|\theta \right)=c{Q}_{\theta }\left(z^{(i)}\right)$$
方程两边同时累加和：
$$\sum _{z^{(i)}\in \mathcal{Z}}P\left(x^{(i)},z^{(i)}|\theta \right)=c\sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)$$
由于${\sum }_{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)=1$，有：
$$\sum _{z^{(i)}\in \mathcal{Z}}P\left(x^{(i)},z^{(i)}|\theta \right)=c$$
从而，
$$Q_{\theta }\left(z^{(i)}\right)=\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{c}=\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{{\sum }_{z^{(i)}\in \mathcal{Z}}P\left(x^{(i)},z^{(i)}|\theta \right)}=\frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{P\left(x^{(i)}|\theta \right)}=P\left(z^{(i)}|x^{(i)},\theta \right)$$
上式中，$Q_{\theta }(z^{(i)})$是已知观测变量$x^{(i)}$和模型参数$\theta$下的隐变量$z^{(i)}$分布。

  经过上述推导，我们可以得出结论：当参数$\theta$给定时，选择$Q_{\theta }\left(z^{(i)}\right)=P\left(z^{(i)}|x^{(i)},\theta \right)$，这样给出的下界与目标函数值最为接近（就等于目标函数值）。

极大化$Q_{\theta }\left(z^{(i)}\right)$给出的下界来极大化目标函数

    至此，我们推出了在固定参数$\theta$后分布$Q_{\theta }\left(z^{(i)}\right)$的选择问题，从而建立了目标函数的下界，并且这个下界与目标函数值及其接近，所以，极大化对数似然函数等价于极大化这个下界，接下来的就是固定$Q_{\theta }\left(z^{(i)}\right)$后，调整$\theta$，去极大化目标函数的下界，即：

$$\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)\log \frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{Q_{\theta }\left(z^{(i)}\right)}$$

  去掉上式中常数的部分：
$$\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}\in \mathcal{Z}}{Q}_{\theta }\left(z^{(i)}\right)\log P\left(x^{(i)},z^{(i)}|\theta \right)$$

EM算法的收敛性

  设$\mathrm{L}(\theta )=\log p(\mathrm{x}|\theta )$是观测数据$x$的对数似然函数，${\theta }^{(i)}(i=1,2,\dots ,n)$是EM算法得到的参数估计序列，$\mathrm{L}\left({\theta }^{(\mathrm{i})}\right)$为对应的对数似然函数序列，则$\mathrm{L}\left({\theta }^{(\mathrm{i})}\right)=\log p\left(x|{\theta }^{(\mathrm{i})}\right)$必定会收敛到某一值$L^{\ast }$。

  证明：要证明 EM 算法收敛，则我们需要证明我们的对数似然函数的值在迭代的过程中一直在增大，即：
$$\sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j+1)}\right)\ge \sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j)}\right)$$
由于
$$Q\left(\theta ,{\theta }^{(j)}\right)=\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log P\left(x^{(i)},z^{(i)}|\theta \right)$$
令
$$H\left(\theta ,{\theta }^{(j)}\right)=\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log P\left(z^{(i)}|x^{(i)},\theta \right)$$
上两式相减得到：
$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(j)}\right)-H\left(\theta ,{\theta }^{(j)}\right)& =\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log P\left(x^{(i)},z^{(i)}|\theta \right)\\ & -\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log P\left(z^{(i)}|x^{(i)},\theta \right)\\ & =\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log \frac{P\left(x^{(i)},z^{(i)}|\theta \right)}{P\left(z^{(i)}|x^{(i)},\theta \right)}\\ & =\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log P\left(x^{(i)}|\theta \right)\\ & =\sum _{i=1}^m\log P\left(x^{(i)}|\theta \right)\end{array}$$
  在上式中分别取$\theta$为${\theta }^{(j)}$和${\theta }^{(j+1)}$，并相减得到：
$$\sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j+1)}\right)-\sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j)}\right)=\left[L\left({\theta }^{(j+1)},{\theta }^{(j)}\right)-L\left({\theta }^{(j)},{\theta }^{(j)}\right)\right]-\left[H\left({\theta }^{(j+1)},{\theta }^{(j)}\right)-H\left({\theta }^{(j)},{\theta }^{(j)}\right)\right]$$
  要证明EM算法的收敛性，我们只需要证明上式的右边是非负的即可。
  由于${\theta }^{(j+1)}$使得$L(\theta ,{\theta }^{(j)})$极大，因此有：
$$Q\left({\theta }^{(j+1)},{\theta }^{(j)}\right)-Q\left({\theta }^{(j)},{\theta }^{(j)}\right)\ge 0$$
  而对于第二部分，我们有：
$$\begin{array}{rl}H\left({\theta }^{(j+1)},{\theta }^{(j)}\right)-H\left({\theta }^{(j)},{\theta }^{(j)}\right)& =\sum _{i=1}^m\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\log \frac{P\left(z^{(i)}|x^{(i)},{\theta }^{(j+1)}\right)}{P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)}\\ & \le \sum _{i=1}^m\log \left(\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)\frac{P\left(z^{(i)}|x^{(i)},{\theta }^{(j+1)}\right)}{P\left(z^{(i)}|x^{(i)},{\theta }^{(j)}\right)}\right)\\ & =\sum _{i=1}^m\log \left(\sum _{z^{(i)}}P\left(z^{(i)}|x^{(i)},{\theta }^{(j+1)}\right)\right)=0\end{array}$$
  上式中不等式用到了Jensen不等式。至此，得到了
$$\sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j+1)}\right)-\sum _{i=1}^m\log P\left(x^{(i)}|{\theta }^{(j)}\right)\ge 0$$
证明了EM算法的收敛性。

  从上面的推导可以看出，EM 算法可以保证收敛到一个稳定点，但是却不能保证收敛到全局的极大值点，因此它是局部最优的算法，解决方法是初始化几次不同的参数进行迭代，取结果最好的那次。当然，如果我们的优化目标$Q(\theta ,{\theta }^{(j)})$是凸的，则EM算法可以保证收敛到全局极大值，这点和梯度下降法这样的迭代算法相同。