EM算法

问题引入

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

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

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

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

倘若现在我告诉你，班里的100个人里，我已经自动地根据男女给你分好了集合。集合A全都是男生，集合B全都是女生，那么男生女生的正态分布参数你会做吗？

显然，很简单，就是两个分别做一下极大似然估计参数即可。

倘若现在反过来，我告诉你，我把男生服从的正态分布的参数θ1:N(μ1,σ1)\theta_{1}:N(\mu_{1},\sigma_{1})告诉你，把女生服从的正态分布参数θ2:N(μ2,σ2)\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(θ)=∏i=1nf(x(i);θ)L(\theta)=\prod_{i=1}^nf(x^{(i)};\theta)

之后，对似然函数取log，让其变为求和形式

log(L(θ))=log(∏i=1nf(x(i);θ))=∑ilogf(x(i);θ)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])E[f(X)]>=f(E[X])

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

EM算法的推导

1.似然函数的变换

第一步，我们同样地，来求似然函数的极大值。

只不过，这次的问题在于，我们的似然函数，包含了隐变量z。

∑ilogp(x(i);θ)=∑ilog∑z(i)p(x(i),z(i);θ)\sum_{i}logp(x^{(i)};\theta) = \sum_{i}log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)

这一步很好理解，这是一个全概率公式，把隐变量zi露出来，并对z的所有可能加和

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

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

∑ilogp(x(i);θ)==∑ilog∑z(i)p(x(i),z(i);θ)∑ilog∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))\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*}

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

那么这个Qi(z(i))Q_{i}(z^{(i)})是什么呢？

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

第三步，利用Jensen不等式

∑ilog∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))>=∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(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外面去了，虽然等号退化成为了不等号，但是好消息是，我们貌似已经有能力算右边了。

2.什么是E，什么是M

此时此刻，为了防止你不理解，我们停一下，捋一捋我们要做什么：

Q1:我们的目标是什么？

A1:是最大化不等号左边的项【即原似然函数的最大化】。那么

Q2:如何最大化左边的项？

A2：争取找到取等号的条件，并不断取右边项的极大。【因为，当我最优化我的下界时，我的似然函数也会得到更大的值。】

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

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

答案很简单，这样定义的话，我们实际上相当于定义了一个期望。

如果看到这里，没有太多数学底子的你心里比较虚。没关系，我用最快的方式帮你回顾一下期望。

如果XX是一个随机变量，那么XX的期望E(X)E(X)为:

E(X)=∑ixi∗p(xi)E(X)=\sum_{i}x_{i}*p(x_{i})

而且，根据期望的性质，离散型随机变量函数f(x)f(x)的期望为：

E(f(X))=∑if(xi)∗p(xi)E(f(X))=\sum_{i}f(x_{i})*p(x_{i})

其中，∑ip(xi)=1\sum_{i}p(x_{i})=1

若X是连续型随机变量，其实是一样的，只需要把求和号变成积分号即可。

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

log∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))=logE(p(x(i),z(i);θ)Qi(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)})})

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

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

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

M-step: 我们寻找参数θ\theta,将下界∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(z(i))\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可以取得到等号】

因此，我们想要p(x(i),z(i);θ)Qi(z(i))=constant\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}=constant

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

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

Qi(z(i))===p(x(i),z(i);θ)∑zp(x(i),z(i);θ)p(x(i),z(i);θ)p(x(i);θ)p(z(i)|x(i);θ)\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不等式取等号【换句话说，我们想得到一个似然函数非常紧的下界】，我们实际上所需要的那个分布函数Qi(z(i))Q_{i}(z^{(i)})，其实就是给定样本和参数的时候，隐变量z(i)z^{(i)}的条件概率

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

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

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

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

（2）E-step：对于每一个样本i，计算Qi(z(i))=p(z(i)|x(i);θ)Q_{i}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)

（3）M-step：计算θ(i+1):=argmaxθ∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ(i))Qi(z(i))\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，有||θ(i+1)−θ(i)||<ε||\theta^{(i+1)}-\theta^{(i)}||<\varepsilon

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

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

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

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

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

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

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

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

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

高斯混合模型问题重述：

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

p(x)=∑k=1KπkN(x|μk,σk) p(x)=\sum_{k=1}^K\pi_kN(x|\mu_k,\sigma_k)

其中，zkz_k代表属于哪类，πk\pi_k代表属于zkz_k这个类的概率

这时，我们得到了一个等价公式，讲潜在变量显示地表出。

接下来，引入所谓的γ(zk)\gamma(z_k),即“给定x的条件下，隐变量z的条件概率”【其实这不就是求我们最上面的那个Q嘛！】

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

γ(zk)=p(zk=1|x)==p(zk=1)p(x|zk=1)∑Kj=1p(zj=1)p(x|zj=1)πkN(x|μk,σk)∑Kk=1πkN(x|μk,σk)\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*}

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

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

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

logp(X|π,μ,σ)=∑n=1Nlog∑k=1KπkN(xn|μk,∑)logp(X|\pi,\mu,\sigma)=\sum_{n=1}^{N}log\sum_{k=1}^{K}\pi_kN(x_n|\mu_k,\sum)

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

E-step：计算γ(zk)=p(zk=1|x)=p(zk=1)p(x|zk=1)∑Kj=1p(zj=1)p(x|zj=1)\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：计算θ(i+1):=argmaxθ∑i∑z(i)γ(zk)logp(x(i),z(i);θ(i))γ(zk)\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)}

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

μk=1Nk∑Ni=1γ(zk)\mu_k=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_k)

∑k=1Nk∑Ni=1γ(zk)(xi−μ)(xi−μ)T\sum_k=\frac{1}{N_k}\sum_{i=1}^N\gamma(z_k)(x_i-\mu)(x_i-\mu)^T

πk=1Nk∑Ni=1\pi_k=\frac{1}{N_k}\sum_{i=1}^N

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

总结

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

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

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

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

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

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

这样做的数学原理：

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

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

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

KL(q||p)=−∑Zq(Z)lnp(Z|X,θ)q(Z)KL(q||p) = - \sum_{Z}q(Z)ln{\frac{p(Z|X,\theta)}{q(Z)}}

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

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

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

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

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



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