EM算法

问题引入

先思考这样一个问题：我们知道，人群中人的身高大致服从一个正态分布。那么现在，如果说我拿到了一个班的学生（就姑且假设是100人吧！）的身高，我想请你帮我估计一下，这个正态分布的参数$\theta:N(\mu,\sigma)$。如何估计？

好简单。应用极大似然估计的思想，把每一个样本拿出来相乘，求解得到概率最大的那个参数，即为我们想要的参数$\theta$

好，现在我们将问题增加一点点难度，倘若我想问，这个班的学生其实可以分为男生和女生。同样地，男生和女生的分布也大致服从正态分布。那么我想问，你能根据这些样本，算出这两个正态分布的参数$\theta_{1}:N(\mu_{1},\sigma_{1})$,$\theta_{2}:N(\mu_{2},\sigma_{2})$吗？

此时，你心里一定会有个小嘀咕。这问题看起来比较棘手，因为你告诉我的数据信息太“少”了。

倘若现在我告诉你，班里的100个人里，我已经自动地根据男女给你分好了集合。集合A全都是男生，集合B全都是女生，那么男生女生的正态分布参数你会做吗？
显然，很简单，就是两个分别做一下极大似然估计参数即可。

倘若现在反过来，我告诉你，我把男生服从的正态分布的参数$\theta_{1}:N(\mu_{1},\sigma_{1})$告诉你，把女生服从的正态分布参数$\theta_{2}:N(\mu_{2},\sigma_{2})$也告诉你，那么，你能告诉我对于每个身高，TA应该是男生还是女生吗？也很简单，我算一下这个身高属于男生和属于女生的概率，哪个大，我就认为TA属于哪一边儿。

但是，你有没有发现这个问题的好玩儿之处：我有两个缺失的条件A,B,你告诉了我A，我可以告诉你B是什么。或者你告诉我B，我也可以把A求出来告诉你。可是恼人的是，偏偏A和B，我们都不知道。这时我们该怎么办呢？

其实，刚刚提出的那个问题，就是一个典型的高斯混合模型的问题，接下来要讲的EM算法，就是专门来解决这个问题的一种算法。

EM算法初探

EM算法，英文全称是Expection Maximization，又叫做最大期望算法。
是在概率（probabilistic）模型中寻找参数最大似然估计的算法，最早在1950年被 Ceppellini等人用于基因频率估计，后来在隐马尔科夫模型（也常被称作Baum-Welch算法）方面被Hartley和Baum等人进一步推广分析。其中概率模型依赖于一些无法观测的隐藏变量（Latent Variable）。通过迭代，来计算含有隐变量的概率模型的参数的极大似然估计。

一句话总结EM算法的不同之处：EM算法可以利用不完整的信息实现概率模型的参数化估计。

EM算法的应用范畴：聚类，高斯混合模型，pLSA,隐马尔科夫模型（与B-W算法十分相似）

预备知识

预备知识1：极大似然估计

在我们做模型参数的预测时，常用的办法，就是利用极大似然估计的思想

$$L(\theta)=\prod_{i=1}^nf(x^{(i)};\theta)$$
之后，对似然函数取log，让其变为求和形式
$$log(L(\theta))=log(\prod_{i=1}^nf(x^{(i)};\theta))=\sum_{i}logf(x^{(i)};\theta)$$
然后取偏导数，另偏导数为零，即可求解参数。

预备知识2：凸函数

设f是定义域为实数的函数，如果对于所有的实数x。如果对于所有的实数x，f(x)的二次导数大于等于0，那么f是凸函数。

当x是向量时，如果其hessian矩阵H是半正定的，那么f是凸函数。如果只大于0，不等于0，那么称f是严格凸函数。

预备知识3：Jensen不等式

Jensen不等式表述如下：

如果f是凸函数，X是随机变量，那么：

$$E[f(X)]>=f(E[X])$$

特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。

EM算法的推导

1.似然函数的变换

第一步，我们同样地，来求似然函数的极大值。
只不过，这次的问题在于，我们的似然函数，包含了隐变量z。

$$\sum_{i}logp(x^{(i)};\theta) = \sum_{i}log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)$$
这一步很好理解，这是一个全概率公式，把隐变量zi露出来，并对z的所有可能加和

我们发现，如果用传统的求偏导，并令偏导数等于0的话，会无法求解，因为$log(f(x_1)+f(x_2)+...+f(x_n))$ 没有人愿意去算这个log里还带加和的公式的导数。

因此，第二步，我们考虑对该公式做一些数学方法上的“调整”，以期望可以变成我们好解决的形式

$$\begin{eqnarray*} \sum_{i}logp(x^{(i)};\theta) =&\sum_{i}log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)\\ = & \sum_{i}log\sum_{z^{(i)}}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})} \end{eqnarray*}$$

我们添加了一个公式$Q_{i}(z^{(i)})$，乘了一个$Q_{i}(z^{(i)})$又除了一个$Q_{i}(z^{(i)})$

那么这个$Q_{i}(z^{(i)})$是什么呢？
对于每一个样例i，让$Q_{i}(z^{(i)})$表示该样例隐含变量z的某种分布，既然是分布函数，那么$Q_{i}(z^{(i)})$满足的条件是$\sum_{z}Q_{i}(z)=1$，且$Q_{i}(z)>0$ 至于如何选择$Q_{i}(z^{(i)})$，接下来会继续谈到

第三步，利用Jensen不等式

$$\sum_{i}log\sum_{z^{(i)}}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}>=\sum_{i}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$$

我们可以发现，经过我们Jensen不等式的处理，不等号右边的这一项，求和号已经都跑到log外面去了，虽然等号退化成为了不等号，但是好消息是，我们貌似已经有能力算右边了。

2.什么是E，什么是M

此时此刻，为了防止你不理解，我们停一下，捋一捋我们要做什么：
Q1:我们的目标是什么？
A1:是最大化不等号左边的项【即原似然函数的最大化】。那么
Q2:如何最大化左边的项？
A2：争取找到取等号的条件，并不断取右边项的极大。【因为，当我最优化我的下界时，我的似然函数也会得到更大的值。】

好了。回答了这两个问题，我们便有了问题解决的大方向：可以不断地建立似然函数【左边项】的下界（E步），然后最大化下界（M步）。不停迭代，直到我们的下界的极值逼近了函数的极值。如果你还是一头雾水，也许你会想问，不是说好的期望最大化（EM）算法么？期望在哪儿呢？最大化在哪儿呢？不要急，我们接下来就能回答这个问题

我们在第二步，“莫名其妙”地创造了一个z的分布函数$Q_{i}(z^{(i)})$，而且还费尽口舌地定义了这个分布函数的许多性质【其实也没有许多啦，我们仅仅要求每个$Q_{i}(z)>0$且加和为1】。现在我们就来解释下，为什么要这么定义。

答案很简单，这样定义的话，我们实际上相当于定义了一个期望。
如果看到这里，没有太多数学底子的你心里比较虚。没关系，我用最快的方式帮你回顾一下期望。
如果$X$是一个随机变量，那么$X$的期望$E(X)$为:

$$E(X)=\sum_{i}x_{i}*p(x_{i})$$
而且，根据期望的性质，离散型随机变量函数 $f(x)$的期望为：
$$E(f(X))=\sum_{i}f(x_{i})*p(x_{i})$$
其中， $\sum_{i}p(x_{i})=1$
若X是连续型随机变量，其实是一样的，只需要把求和号变成积分号即可。

好了，现在我们重看第二步，我们可以看到，$log\sum_{z^{(i)}}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$,实际上是这个$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$的期望,即

$$log\sum_{z^{(i)}}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}=logE(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})$$
在第三步我们可以看到 $$\sum_{i}log\sum_{z^{(i)}}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}>=\sum_{i}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$$
通过Jensen不等式，我们可以得到，对数似然函数的下界就是要计算 $log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$的，因此，在第二步里，求 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$的期望，这实际上是一个相对于Zi的分布Qi的一个期望，我们可以近似理解成“求下界的期望”，在第三步里，我们“求下界的最大化”。

这就是所谓的期望最大化。现在整理一下：

E-step：我们求对数似然下界函数的期望，即$logE(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})$

M-step: 我们寻找参数$\theta$,将下界$\sum_{i}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$最大化

好，解决了“什么是E，什么是M”的问题，我们下一步进入“如何求解EM”这一块

3.如何求解EM

还记得上面在A2处留了一个未解决的问题，即“Jensen不等式什么时候可以取等号？”

因为只有Jensen不等式取等号，才能让我的下界尽可能逼近原似然函数，我的下界才是一个尽可能“紧”的下界。也恰恰是遵循着这个原则，我要选择一个合适的分布函数Q，能让我这个等式取等号。这就解决了如何选择Q的问题。

我们就从这个点切入，首先思考一下不等式取等号的条件：不等式的等号成立，当且仅当这个值是一个常量。【因为常量的期望还是自己，所以那个E可以取得到等号】

因此，我们想要

$$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}=constant$$
这样，我们可以看到，其实 $Q_{i}(z^{(i)}$是可以用 ${p(x^{(i)},z^{(i)};\theta)}$来表示的！

我们继续推导，由于$\sum_{z}Q_{i}(z^{(i)})=1$，那么有$\sum_{z}p(x^{(i)},z^{(i)};\theta)=c$,那么

$$\begin{eqnarray*} Q_{i}(z^{(i)})=&\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_{z}p(x^{(i)},z^{(i)};\theta)}\\ =&\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}\\ =&p(z^{(i)}|x^{(i)};\theta) \end{eqnarray*}$$

我们可以看到，如果我们想让Jensen不等式取等号【换句话说，我们想得到一个似然函数非常紧的下界】，我们实际上所需要的那个分布函数$Q_{i}(z^{(i)})$，其实就是给定样本和参数的时候，隐变量$z^{(i)}$的条件概率

我们现在再重新看E过程：由于我们想让$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$为常数，即$E(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})=\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$因此,由于是常数了，我们确定了$Q_{i}(z^{(i)})$，也就确定了下界$logE(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})$

因此，我们的E过程可以修正为：确定概率分布$Q_{i}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)$

这比之前的E过程“确定一个下界$logE(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})$的期望”要简单容易多了。因此，我们得到最终的EM算法的步骤：

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

（2）E-step：对于每一个样本i，计算

$$Q_{i}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)$$

（3）M-step：计算

$$\theta^{(i+1)}:=argmax_{\theta}\sum_{i}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(i)})}{Q_{i}(z^{(i)})}$$

（4）重复步骤2,3，直至收敛【收敛条件可自己指定，如对于一个较小的正数$\varepsilon$，有$||\theta^{(i+1)}-\theta^{(i)}||<\varepsilon$

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

讲完了推导，我们来看几个实际应用的例子。现介绍一个EM算法在高斯混合模型中的应用，我们将使用EM算法求出高斯混合模型（Gaussian Mixture Model，以下简称GMM）的隐藏参数。

其实，回到我们最开始的那个问题，其实那个问题就是一个最简单的高斯混合模型：在一个数据集中，分布着多个服从不同$\mu$,$\sigma$的正态分布，他们混合在一起，不知道各自的比例是多少，也不知道各自的参数是多少，我们能看到的，只有共同生成的一群数据。我们既不知道每个数据属于哪个高斯分布的概率，也不知道每个高斯分布的参数的确切值。用一个链条化的方法来表示就是

所有数据里包括多个小数据集——有多少个小数据集是一个多项分布（$Multinomial$）——每一个小数据集里的数据，服从一个参数未知的高斯分布

没听懂？没关系，我反过来再说一遍，相信你就明白了。

我们假设数据是这样产生的：

从一个服从多项分布的类别里面抽取一个出来，比如选择了$Z_{i}$——再根据$Z_{i}$里的数据服从正态分布，根据$Z_{i}$自己的$\mu_{i}$,$\sigma_{i}$来生成一个数据。

所以，我们最开始的问题，类别只有两个：男生，女生，因此这个多项分布退化为了二项分布，然后，对于每一个类别男女，有各自自己的参数。那么，这里的隐变量Z，就是Z1=男，Z2=女。

好了。理解了一个最简单的高斯混合模型之后，我们来看一般情况并求解一般情况。多项分布的如果听懂了，那么二项分布简直so easy。

高斯混合模型问题重述：

随机变量X是由K个高斯分布混合而成，那么GMM的分布可以写成：高斯分布的线性叠加的形式

$$p(x)=\sum_{k=1}^K\pi_kN(x|\mu_k,\sigma_k)$$

其中，$z_k$代表属于哪类，$\pi_k$代表属于$z_k$这个类的概率

这时，我们得到了一个等价公式，讲潜在变量显示地表出。
接下来，引入所谓的$\gamma(z_k)$,即“给定x的条件下，隐变量z的条件概率”【其实这不就是求我们最上面的那个Q嘛！】

求法也很简单，最简单的贝叶斯公式：

$$\begin{eqnarray*} \gamma(z_k)=p(z_k=1|x)=&\frac{p(z_k=1)p(x|z_k=1)}{\sum_{j=1}^Kp(z_j=1)p(x|z_j=1)}\\ =&\frac{\pi_kN(x|\mu_k,\sigma_k)}{\sum_{k=1}^K\pi_kN(x|\mu_k,\sigma_k)} \end{eqnarray*}$$

其实我们可以这样看，我们把$\pi_k$看作$z_k=1$的先验概率，然后得到的这个$\gamma(z_k)$看作观测到x之后，对应的“后验概率”，用通俗的语言文字解释就是“这个类，对这个样本付了多大责任”，或者说，“这个样本，有多少是由这个类解释的”。

这地方有点绕，不过请多读两遍，肯定会有豁然开朗的感觉。

接下来，理论层面搞懂了，下一步就是求解极大似然了。依照上面的方法，写出极大似然的函数：

$$logp(X|\pi,\mu,\sigma)=\sum_{n=1}^{N}log\sum_{k=1}^{K}\pi_kN(x_n|\mu_k,\sum)$$

接下来，EM的两个步骤变为：

E-step：计算$\gamma(z_k)=p(z_k=1|x)=\frac{p(z_k=1)p(x|z_k=1)}{\sum_{j=1}^Kp(z_j=1)p(x|z_j=1)}$

M-step：计算$\theta^{(i+1)}:=argmax_{\theta}\sum_{i}\sum_{z^{(i)}}\gamma(z_k)log\frac{p(x^{(i)},z^{(i)};\theta^{(i)})}{\gamma(z_k)}$

此时已经可以解析求解了，解析求得的解为：

$\mu_k=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_k)$

$\sum_k=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_k)(x_i-\mu)(x_i-\mu)^T$

$\pi_k=\frac{1}{N_k}\sum_{i=1}^N$

好了，此时我们便得到了高斯混合模型第K个高斯分布的参数$\theta_k=(\mu_k,\sum_k)$ 此时我们回到最初男女身高的那个问题，那么我们实际上就是求$\theta_1$和$\theta_2$啦

总结

EM算法是一种求解隐变量的数值计算方法，本文通过讲述如何在高斯混合模型中求解参数，借此推导了一遍EM算法。其实，EM算法不仅仅能应用在高斯混合模型，还可以应用在诸如隐马模型，LDA模型中，与变分，采样的结合，形成了更为强大的数值计算武器。

【附录：可选读】EM算法的另一种数学理解

最后，用一种数学的方法去理解EM算法的每一步，在数学上是怎么一回事儿

假如我们现在有一个对数似然函数,我们E步和M步分别做了什么呢？

E步，就是找蓝色的这条线，即“找一个十分紧的下界”

M步，就是把这个紧的下界求极大值，让新的theta new取在原来theta old函数【原下界】的最大值上，依次迭代。

这样做的数学原理：
我们定义了一个隐变量的分布q（Z），对于任意的q（Z），以下公式成立：

$$lnp(X|\theta) = \mathcal{L}(q,\theta)+KL(q||p)$$

其中$\mathcal{L}(q,\theta) = \sum_{Z}q(Z)ln{\frac{p(X,Z|\theta)}{q(Z)}}$

$KL(q||p) = - \sum_{Z}q(Z)ln{\frac{p(Z|X,\theta)}{q(Z)}}$

$\mathcal{L}(q,\theta)$是概率分布q(Z)的一个泛函，KL(q||p)是KL散度，它包含了给定X的条件下，Z的条件概率分布

由于KL散度的性质【此处不再继续延伸】，KL(q||p)>=0恒成立。当且仅当$q(Z) = p(Z|X,\theta)$时，等号成立。

因此，我们想要等号成立，实际上是我们想把KL散度趋于0。换句话说，$\mathcal{L}(q,\theta)$是原函数$lnp(X|\theta)$的一个下界。

在E步时，实际上是下界$\mathcal{L}(q,\theta^{旧})$关于q（Z）被最大化。因此，下界的最大化出现在KL散度为0的时候，即q(Z)等于它的后验概率分布的时候。

而M步的时候，q(Z)保持固定，我们把下界$\mathcal{L}(q,\theta^{旧})$关于$\theta$进行最大化，得到一个新的$\theta^{新}$，此时，下界L会继续增大。

这样也就证明了，EM算法每一步都是让下界L函数最大化的方法，也就保证了每一步都会趋向于最优值，也就间接地证明了，EM算法会最终得到一个收敛的解。