机器学习笔记EM算法及实践（以混合高斯模型（GMM）为例来次完整的EM）

Posted 2020-09-06 王大宝的CD

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了机器学习笔记EM算法及实践（以混合高斯模型（GMM）为例来次完整的EM）相关的知识，希望对你有一定的参考价值。

今天要来讨论的是EM算法。第一眼看到EM我就想到了我大枫哥，EM Master，千里马，RUA!!!不知道看这个博客的人有没有懂这个梗的。好的，言归正传，今天要讲的EM算法，全称是Expectation maximization，期望最大化。怎么个意思呢，就是给你一堆观测样本，让你给出这个模型的参数估计。我靠，这套路我们前面讨论各种回归的时候不是已经用烂了吗？求期望，求对数期望，求导为0，得到参数估计值，这套路我懂啊，MLE！但问题在于，如果这个问题存在中间的隐变量呢？会不会把我们的套路给带崩呢，我们通过两个例子来认识一下这两种情况。

====================================================================

不存在中间变量的EM。

假设有一天人类消除了性别的差别，所有的人都是同一个性别。这个时候，我给了你一群人的身高让你给我做一个估计人身高的模型。

怎么办呢？感觉上身高应该是服从高斯分布吧，所以假设人的身高分布服从高斯分布N(Mu,Sigma^2)，现在我又有了N个人的身高的数据，我就可以照着上面的套路进行了。先求对数似然函数

接着对两个参数求偏导为0

这样就得到了我们的参数估计

喜闻乐见的结果，又好求又符合我们的直觉，那我们再来看看另一种情况。

====================================================================

存在中间变量的EM。

正如你所知，身高和人种的关系挺大的，而人类又有那么多种族，所以，再给你一堆人的数据，要做一个估计人身高的模型，那我们应该怎么做呢？

首先，现在分为不同种族若干类了，这些类别的概率肯定有个分布，其次，各种族当中身高是服从不同的分布的，那么这样身高的估计就变成了

Alpha代表了该样本属于某一人种的比例，其实就是隐藏的中间变量。Muk和sigmak^2为各类高斯分布的参数。按照我们上面的套路就是求对数似然概率再求导得到参数的估计，那么先来看看似然函数

这下尴尬的情况出现了，对数里面带加号，这下求导就变得复杂异常了，而且没法求解，事实上，这种式子确实没有解析解。不过憋灰心啊，假设我们随便猜一个alpha的分布为Q，那么对数似然函数可以写成

由于Q是alpha的一个分布，所以似然函数可以看成是一个log(E(x))，log是个凹函数啊，割线始终在函数图像下方，Jensen不等式反向应用一下，有log(E(x))>=E(log(x))，所以上面的对数似然有

冷静分析一下现在的情况，我们现在得到了一个对数似然函数的下界函数，我们采用曲线救国的战略，我们求解它的局部最大值，那么更新后的参数带入这个下界函数一定比之前的参数值大，而它本身又是对数似然函数的下界函数，所以参数更新后，我们的对数似然函数一定是变大了！所以，就利用这种方法进行迭代，最后就能得到比较好的参数估计。还有点晕吗，没事，我从百度扒个图给你来个形象的解释

红色那条线就是我们的对数似然函数，蓝色那条是我们在当前参数下找到的对数似然的下界函数，可以看到，我们找到它的局部极值那么参数更新成thetanew，此时对数似然函数的值也得到了上升，这样重复进行下去，是不是就可以收敛到对数似然函数的一个局部极值了嘛。对的，局部极值，并不能保证是全局最优，但它就是个估计嘛，你还要她怎样？！

到了这里，我们好像跳着先把第二步参数更新的工作做完了，那么还有一个事情是我们要注意的，Q呢，Q是啥，没Q你算啥极值，更新啥参数。我们已经知道Q是alpha的一个分布，然后我们肯定是希望这个下界函数尽量贴近原来的对数似然函数，这样我们才能更快地更新参数，那下界函数啥时最大呢，等号成立呗，等号成立说明你求期望的对象是个常数呀，所以log和Q谁在前后都无所谓，那么就有了

直观地可以理解成第i个样本来自第k个类别的可能性。好了，现在Q也确定了，我们根据上面所说的方法更新参数，再更新Q，再更新参数，迭代进行下去就可以了。

如果你能坚持看到这里，少侠我只能说你大功已成。因为其实我们已经把EM算法整个推导完了，也许你还是有点云里雾里，那我们再来仔细梳理一下这个流程

1 拿到所有的观测样本，根据先验或者喜好先给一个参数估计。

2 根据这个参数估计和样本计算类别分布Q，得到最贴近对数似然函数的下界函数。

3 对下界函数求极值，更新参数分布。

4 迭代计算，直至收敛。

说起来啊，EM算法据说是机器学习进阶的一个算法，但至少目前来看，它的思路还是很容易理解的嘛，整个过程中唯一一个可能初学者觉得有点绕的地方就是应用Jensen不等式的那一步，那我再啰嗦两句。所谓Jensen不等式，你可以这么理解，对于一个凸函数f而言，它的割线始终在函数图像上方你承认吧，我在上面任取两点x1，x2，参数theta介于0到1之间，那么theta*x1+(1-theta)*x2就是介于x1和x2之间的一点吧，在这点上过x1x2割线的值大于函数值吧，是不是就有了theta*f(x1)+(1-theta)*f(x2)>f(theta*x1+(1-theta)*x2)，根据这个结论再推广开来，就有E(f(x))>f(E(x))，在对数似然函数中，由于log是个凹函数，所以把它反过来用，老铁没毛病吧？！这一点想通了我觉得整个EM算法的流程还是蛮好懂的。

下面呢，我们还回到这个身高模型的预测，假设给了m个样本，有k个种族，每个种族的身高都是服从高斯分布的，那么这就变成了EM算法中最具代表性的一个例子，高斯混合模型GMM。

====================================================================高斯混合模型（GMM）

刚才已经讲了EM算法的套路了，假设我们现在处于某一步迭代中，那么我们该干嘛呢？

E-step 求最佳的类别分布

可以将其理解为第i个样本属于第J类的概率。

M-step 更新参数

求得了Q之后，我们就得到了最贴近原对数似然函数的下界函数，那我们对它求极值就可以得到更新后的参数，先看一下这个下界函数

Log函数里面全是乘积项这是我们最喜欢的形式，这样求导的时候但凡不相关的我们直接扔掉就行，待求参数mu,sigma^2,psi，依次求导为0就成。

对于psi的求解可能复杂一些，首先我们把下界函数中与psi不相关的项全部去掉，然后psi作为各类别的比例有一个天然的约束条件就是所有的psi之和为1，所以目标函数变成

这种带约束的优化前面在SVM的时候不知道用了多少回，拉格朗日乘子法

接下来对psi求导

两边同时再对j从1到k连加，psi那一项就没了，右式就变成样本数目m，这样就求得了beta，回代我们就可以求得psi参数的更新

至此，所有的参数更新工作就已完成，下面重复进行迭代就行了。

我们先把GMM的算法梳理一下

1 给参数取初始值，开始迭代。

2 求每个样本对每个类别的概率，科学的叫法叫求响应度

3 更新模型参数

4 重复23两步直至收敛。

我们再来看看这些参数的意义，其实未尝不符合我们的直觉认识。W(i,j)可以看做第i个样本属于第j类的概率，那么所有样本中属于第j类的个数就是w(i,j)之和，每个样本xi对应第j类的值就是W(i,j) xi，这样算的平均数就是第j类对应的mu，继续按照这个思路算的方差就是第j类的sigma^2，第j类的概率就是第j类的个数除以总样本数。所以，GMM模型虽然推导起来有点吓人，但仔细想想它最后的结果也是符合我们的直觉认识的，每个样本都是一部分属于某一类，所有样本中的某一类的部分构成了这一类的分布，perfect！！！

====================================================================

这样的话，理论部分我们就讲完了，接下来又是调包侠的时刻了，上次写完后我想到鸢尾花数据无监督算法也能做啊，不给标签我们强行给它分类看看效果如何。所以这里我们K-Means和GMM算法分别对鸢尾花进行处理，看看它们的聚类效果如何。

代码如下

import numpy as np
from sklearn import datasets
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
#读取数据
iris=datasets.load_iris()
x=iris.data[:,:2]
y=iris.target
mu = np.array([np.mean(x[y == i], axis=0) for i in range(3)])  
print '实际均值 = \\n', mu
#K-Means
kmeans=KMeans(n_clusters=3,init='k-means++',random_state=0)
y_hat1=kmeans.fit_predict(x)
mu1=np.array([np.mean(x[y_hat1 == i], axis=0) for i in range(3)]) 
print 'K-Means均值 = \\n', mu1
print '分类正确率为',np.mean(y_hat1==y)
gmm=GaussianMixture(n_components=3,covariance_type='full', random_state=0)
gmm.fit(x)
print 'GMM均值 = \\n', gmm.means_
y_hat2=gmm.predict(x)
print '分类正确率为',np.mean(y_hat2==y)

输出结果为

实际均值 =

[[5.006 3.418]

[5.936 2.77 ]

[6.588 2.974]]

K-Means均值 =

[[5.77358491 2.69245283]

[ 5.006 3.418 ]

[ 6.81276596 3.07446809]]

分类正确率为 0.233333333333

GMM均值 =

[[5.01494511 3.44040237]

[ 6.69225795 3.03018616]

[ 5.90652226 2.74740414]]

分类正确率为 0.533333333333怒摔键盘啊，什么破正确率呀！！！憋急啊，我看事情并不简单。机智的我们观察一下均值矩阵。K-Means给出的第一行似乎和实际的第二行很接近，第二行和实际的第一行很接近。同样，GMM给出的均值矩阵也有同样的问题，第二行和第三行似乎对调了。这不是算法的锅啊，它只管给你聚类，哪里还能保证标签和你一样啊，三个类别六种标签方式人家算法也只能随机一种好吗，所以现在我们把预测的结果的标签改一下看看实际的正确率如何。

import numpy as np
from sklearn import datasets
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
#读取数据
iris=datasets.load_iris()
x=iris.data[:,:2]
y=iris.target
mu = np.array([np.mean(x[y == i], axis=0) for i in range(3)])  
print '实际均值 = \\n', mu
#K-Means
kmeans=KMeans(n_clusters=3,init='k-means++',random_state=0)
y_hat1=kmeans.fit_predict(x)
y_hat1[y_hat1==0]=3
y_hat1[y_hat1==1]=0
y_hat1[y_hat1==3]=1
mu1=np.array([np.mean(x[y_hat1 == i], axis=0) for i in range(3)]) 
print 'K-Means均值 = \\n', mu1
print '分类正确率为',np.mean(y_hat1==y)
gmm=GaussianMixture(n_components=3,covariance_type='full', random_state=0)
gmm.fit(x)
print 'GMM均值 = \\n', gmm.means_
y_hat2=gmm.predict(x)
y_hat2[y_hat2==1]=3
y_hat2[y_hat2==2]=1
y_hat2[y_hat2==3]=2
print '分类正确率为',np.mean(y_hat2==y)

输出结果为

实际均值 =

[[5.006 3.418]

[ 5.936 2.77 ]

[ 6.588 2.974]]

K-Means均值 =

[[5.006 3.418 ]

[ 5.77358491 2.69245283]

[ 6.81276596 3.07446809]]

分类正确率为 0.82

GMM均值 =

[[5.01494511 3.44040237]