从EM算法到GMM

从EM算法到GMM

1. 从一个简单的实例开始

假设你要对某大学所有学生的身高进行统计并且给出分析结果，你基本会自然而然地可能想到不妨设假设所有学生的身高服从高斯分布但参数未知即:
$$Y\sim N(\mu ,{\sigma }^2)$$
其中$Y$代表学生的身高，你这寻思：“瞧谁不起，直接从全体学生中抽个样然后，对样本求个样本均值和样本方差，这问题不就简简单单就解决了吗”。确实，这么做也属实好使，但是不完美。你想啊，那男生的身高和女生的身高它能服从一个分布吗？你一个男的长个1米8跟玩似的，要是换个女生长个这个头得费多大劲，所以你这个分析结果太粗糙不太行。你又想了：“多大点事儿，那我把男生和女生分开抽不就完了吗，一遍先抽他个1000个样本然后一算”。你说的很对，为了把这个实例推向一般化咱们加点难度，也就说啥意思呢：学校的全体学生对于你来说是个黑箱，你只能从里面获得一个一个的样本也就是身高，但是呢你不能加个过滤器，比如“我就要女生，男的不测”，那不行，只能来一个测一个。 这样一来就存在一个隐含的随机变量$Z$代表性别，比如男性$Z=1$，女性$Z=0$。那这样一来我们的未知参数就变成了$\theta =\left({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2\right)$。这次我们的模型就变成了：
$$Y=p{Y}_1+q{Y}_2$$
其中$P$表示个体是男生的概率，$q$表示个体是女生的概率。显然有$p+q=1$。两个高斯分布的叠加仍然是高斯分布只不过$Y_1和Y_2$如果不独立的话$Y$服从一个二维高斯分布。所以说之前我们直接假设整体服从一个高斯分布还是有那么亿点合理。
模型建立好了下面我们就可以开始我们的参数估计了。注意我要开始操作了。

2. 极大似然估计

对于上面的问题，老频率学派必然得极大似然估计一下，先写个似然函数：
$$l(\theta )=\log P(Y|\theta )$$
这里的$\log$是以$e$为底的。“那为你啥不写$\ln$呢？”，“爷就乐意写$\log$”。其实主要是想和李航老师的《统计学习方法》里面的推导保持一致。上面的似然函数不太完整，因为我们没有加入隐变量，加入隐变量之后：
$$l(\theta )=\log \sum _{Z}P(Y,Z|\theta )$$
上式也就简单的将边缘分布$P(Y|\theta )$变成了联合分布${\sum }_{Z}P(Y,Z|\theta )$。按正常来说接下来的工作也就简单了，极大化嘛。${\theta }^{\ast }=argmaxl(\theta )$。脑子想的很简单，手一算就完犊子了。首先隐变量$Z$是位置的，因为黑箱嘛。其次对数里面有求和符号这还是离散情况，连续情况里面就是一个积分即：$\log {\int }_{Z}P(Y,Z|\theta )$你以为求偏导然后令偏导等于0一解美滋滋？你去算吧，算不死你。“那咋办啊，极大似然都不行了，Fisher拉垮了！”。别急，解析解没有咱们还有迭代方法。

3. EM

我看网上大部分资料大概都是这么个意思：“EM算法分为两步E步求期望，M步求极大”。我个人觉得这话说的和放屁没啥区别，屁好歹还有个响，但是也不能说这句话的是错的因为EM算法的EM确实是这两个意思。
既然我们无法一步取得$l(\theta )$的最大值，我们可以使用迭代的方法来逐步逼近$l(\theta )$的最大值。把这句话翻译翻译就是，假设现在我已经取得了第$i$轮迭代的参数${\theta }^{(i)}$,只需在第$i+1$轮使得$l({\theta }^{(i+1)})>l({\theta }^{(i)})$即可，这样我们可以构造如下的函数：
$$B(\theta ,{\theta }^{(i)})=l(\theta )-l({\theta }^{(i)})=l(\theta )=\log \sum _{Z}P(Y,Z|\theta )-\log P(Y|{\theta }^{(i)})$$
注意这里的${\theta }^{(i)}$是一个已知数,$\theta$是一个未知数。显然这一顿操作好像并没有变得简单反而更复杂了，没事儿你看我接着变形：
$$\begin{array}{rlr}l(\theta )-l({\theta }^{(i)})& =\log \sum _{Z}P(Y|Z,\theta )P(Z|\theta )-\log P(Y|{\theta }^{(i)})& \\ & =\log \sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)})& \\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)})& =B(\theta ,{\theta }^{(i)})\end{array}$$
上面最后一步的放缩利用了jensen不等式具体过程在Jensen不等式初步理解及证明中有讲解。到了这一步概率论学得不错哥儿们应该已经能开出来苗头了，观察对数里面的内容${\sum }_{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$这是啥啊？这不就是期望吗？看出来的兄弟评论扣个1，没看出来的扣眼珠子。其实我当时也没看出来，上网一顿查啊，就是找不着。今天我细细地给你们掰扯一下。首先离散型随机变量的数学期望为：
E ( x ) = ∑ i P { X = x } x x ∈ A E(x)=\sum_{i}{P\{X=x\}x} \quad x\in A E(x)=i∑​P{X=x}xx∈A
A是随机变量X的样本空间，说白了就是x的所有可能取值。如果我想求$f(X)$的期望呢。随机变量函数的期望只要进行简单替换就完了：
E ( f ( x ) ) = ∑ i P { X = x } f ( x ) x ∈ A E(f(x))=\sum_{i}{P\{X=x\}f(x)} \quad x\in A E(f(x))=i∑​P{X=x}f(x)x∈A
这次看明白了吧，对数里面就是在某个条件下的期望，所以说EM算法的E步自然得就出现了。进一步变换我们可以得到：
$$\begin{array}{rl}l(\theta )& \ge B(\theta ,{\theta }^{(i)})+l({\theta }^{(i)}\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\end{array}$$
我们极大化$l(\theta )$的下界得到${\theta }^{(i+1)}$
$$\begin{array}{rl}{\theta }^{(i+1)}& =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y|Z,\theta )P(Z|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})\end{array}$$
至此，我们的EM算法已经差不多结束了。总的来说，E步就是求得当前的$Q$函数，M步就是对$Q$函数求极大获得下一轮的$\theta$。注意啊，EM虽然牛逼但是和大多数的迭代算法一样，不能保证全局最优！

4. GMM

前面磨叽那么多现在我要开始写快点了，高斯混合模型(Gaussian mixture model)首先它是一个聚类的算法。它的前提假设是总体为多个高斯分布的叠加即总体的概率密度为：
$$\begin{gathered} F(x)={\alpha }_1{\varphi }_1(x)+{\alpha }_2{\varphi }_2(x)+\cdots +{\alpha }_k{\varphi }_k(x) \\ \sum _{i=1}^k{\alpha }_i=1 \end{gathered}$$
其中$\varphi (x)$为高斯分布的概率密度函数。如果我们从总体中进行抽样那它就可能来自某个高斯分布。

上面这幅图是从三个不同的高斯分布所抽取的。好了，现在我们的任务就是把这些图的颜色去掉，然后将这些点进行聚类并得到模型的参数。对照EM算法，设我们的隐变量为${\gamma }_{jk}$。
$${\gamma }_{jk}=\{\begin{array}{ll}1& y_j\text{来自第}k个分布\\ 0& 其他\end{array}$$
表示$y_j$来自第$k$个高斯分布。完全数据的似然函数可以写为(这里有个技巧请看完)：
$$\begin{array}{rl}P(Y,\gamma |\theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jk}|\theta )\\ & =\prod _{j=1}^N\sum _{k=1}^K{\gamma }_{jk}{\alpha }_k\varphi (y_j|\theta )\\ & =\prod _{j=1}^N\sum _{k=1}^K{\gamma }_{jk}{\alpha }_k\left(\frac{1}{2\pi {\left|{\Sigma }_k\right|}^{\frac{1}{2}}}\text{exp}\left\{-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right\}\right)\end{array}$$
对完全数据去对数后得到似然函数：
$$l(\theta )=\sum _{j=1}^N\log \sum _{k=1}^K{\gamma }_{jk}{\alpha }_k\left(\frac{1}{2\pi {\left|{\Sigma }_k\right|}^{\frac{1}{2}}}\text{exp}\left\{-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right\}\right)$$
这时候我的心理咯噔了一下，这不完犊子了吗，怎么对数里面还有求和符号啊，一会求偏导的时候不完犊子了吗？我喝了一口水，仔细观察了一下推导的过程没错…突然发现了转机，${\sum }_k{\alpha }_k=1$可以用Jensen不等式进行缩小把求和符号拿到外面去！
$$l(\theta )\ge \sum _{j=1}^N\sum _{k=1}^K{\alpha }_k\log {\gamma }_{jk}\left(\frac{1}{2\pi {\left|{\Sigma }_k\right|}^{\frac{1}{2}}}\text{exp}\left\{-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right\}\right)=C(\theta ,{\theta }^k)$$
虽然$C(\theta ,{\theta }^k)\le Q(\theta ,{\theta }^k)$比之前的下界更小了，但这无伤大雅。But，球逗麻袋，就在我在纸上写完这个$C$函数之后我他妈马上就意识到了如果${\gamma }_{jk}=0$对数里面为$0$超出定义域了，这个方法只能放弃、拉倒了。其实问题在于我写似然函数的时候是完全自己写的，也确实没问题，李航老师的书上是乘法形式没有求和，刚开始我还有疑惑，为什么李航老师要这么写？他这么写肯定有理由，我这个形式行不行呢？带着这个问题我开始往下写，就遇到上述的种种问题然后才豁然开朗，李航老师都是使用乘法形式是为了取对数之后能够全部分离开，而隐变量在乘法形式中是放在$\varphi$函数的指数位置也就不存在对数的真数等于$0$的情况了。
下面给出李航老师版本的似然函数：
$$\begin{array}{rl}P(Y,\gamma |\theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jk}|\theta )\\ & =\prod _{k=1}^K\prod _{j=1}^N{\left[{\alpha }_k\varphi (y_j|\theta )\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left(\frac{1}{2\pi {\left|{\Sigma }_k\right|}^{\frac{1}{2}}}\text{exp}\left\{-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right\}\right)}^{{\gamma }_{jk}}\\ \text{其中}{n}_k=\sum _j{\gamma }_{jk},\sum _k{n}_k=N& \end{array}$$
对数似然为：
$$l(\theta ,{\theta }^k)=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{2\pi }\right)-\frac{1}{2}\log \left|{\Sigma }_k\right|-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right]\right\}$$
至此，我们得到了完全数据(完全是指包含隐变量$\gamma$在内的似然)的似然函数。根据EM算法，接下来我们要用E步来求得$Q$函数，把鼠标滚轮扒拉到EM算法那，然后可知：
$$\begin{array}{rl}Q(\theta ,{\theta }^k)& =\sum _{Z}P(Z|{\theta }^i,Y)\log P(Y|Z,\theta )P(Z|\theta )\\ & ={\mathbb{E}}_Z\left[\log P(Y|Z,\theta )P(Z|\theta )\right]\\ & ={\mathbb{E}}_Z\left[\log P(Y,Z|\theta )\right]\end{array}$$
有的兄弟可能会问了，既然最后求得还是完全数据后验的似然即：${\mathbb{E}}_Z\left[\log P(Y,Z|\theta )\right]$那之前为啥还用贝斯公式变换成：${\mathbb{E}}_Z\left[\log P(Y|Z,\theta )P(Z|\theta )\right]$？这个问题其实我在看书的时候也是挺费解的，可能推导的时候没有考虑那么多，也可能是为了引出$P(Z|{\theta }^i,Y)$，这其实很符合常理，因为一般联合概率都是很难求的，但是边缘分布好求，所以把$\log P(Y,Z|\theta )$拆开还是很自然的，有点牵强哈，懂得都懂，不懂的也尽量理解。
现在回到我们原来的问题上，我们只求求得$Q$函数之后再极大化就可以获得第$i+1$轮迭代的参数。$Q$函数为：
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i+1)})& ={\mathbb{E}}_{\gamma }\left[\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{2\pi }\right)-\frac{1}{2}\log \left|{\Sigma }_k\right|-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right]\right\}\right]\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^N{\mathbb{E}}_{\gamma }[{\gamma }_{jk}]\log {\alpha }_k+\sum _{j=1}^N{\mathbb{E}}_{\gamma }[{\gamma }_{jk}]\left[\log \left(\frac{1}{2\pi }\right)-\frac{1}{2}\log \left|{\Sigma }_k\right|-\frac{1}{2}(y_j-{\mu }_k{)}^T{\Sigma }_k^{-1}(y_j-{\mu }_k)\right]\right\}\end{array}$$
${\gamma }_{jk}$的期望为:
$$\begin{array}{rl}{\hat{\gamma }}_{jk}& =E({\gamma }_{jk}|y_j,\theta )=0\times P({\gamma }_{jk}=0|y_j,\theta )+1\times P({\gamma }_{jk}=1|y_j,\theta )\\ & =P({\gamma }_{jk}=1|y_j,\theta )\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{P(y_j|\theta )}\\ & =\frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{k}P({\gamma }_{jk}=1,y_j|{\theta }_k)}\\ & =\frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{k}P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}\\ & =\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_k{\alpha }_k\varphi (y_j|{\theta }_k)}k=1,2,\cdots ,K\end{array}$$
简简单单地把E步求了出来，接下来就是M步，求Q函数的极大。那M步怎么求？直接对参数求导就完了。用${\hat{\mu }}_k,{\hat{\Sigma }}_k,{\hat{\alpha }}_k$表示${\theta }^{(i+1)}$的各个参数。
$$\begin{array}{rl}{\hat{\mu }}_k& =\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}\\ {\hat{\Sigma }}_k& =\frac{{\sum }_{j=1}^N(y_j-{\mu }_k)(y_j-{\mu }_k{)}^T}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}\\ {\hat{\alpha }}_k& =\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N}\end{array}$$
其中$k=1,2,\dots ,K$，具体的推导过程太长了，需要一定的矩阵求导知识。我可能单独整一篇来讲求导内容。现在迭代公式ok了剩下的内容是什么傻子都能猜出来了，编程测试呗。