[Kaggle] *Clustering and Gaussian Mixture Models (GMMs)
Ref: [Scikit-learn] 2.1 Gaussian mixture models & EM
The PCA representation of the data
It's a low-dimensional representation of the full images - in this case, 256-pixel images are squished down into 3 dimensions. 【什么技术?高大上的样子】
If you are curious about this compression trick, you can use the variables mu and E to explore it.
漫谈 Clustering 系列 <-- 转自
- 漫谈 Clustering (1): k-means
- 漫谈 Clustering (2): k-medoids
- 漫谈 Clustering (番外篇): Vector Quantization
- 漫谈 Clustering (3): Gaussian Mixture Model
- 漫谈 Clustering (番外篇): Expectation Maximization
- 漫谈 Clustering (4): Spectral Clustering
- 漫谈 Clustering (番外篇): Dimensionality Reduction
- 漫谈 Clustering (5): Hierarchical Clustering
- 漫谈 Clustering (6): Evaluation & Visualization
- 漫谈 Clustering (番外篇): Deciding the Number of Clusterings
漫谈 Clustering (3): Gaussian Mixture Model
- k-means 的结果是每个数据点被 assign 到其中某一个 cluster 了,
- GMM 则给出这些数据点被 assign 到每个 cluster 的概率,又称作 soft assignment .
相对于KNN的优势:
得出一个概率有很多好处,因为它的信息量比简单的一个结果要多,比如,我可以把这个概率转换为一个 score ,表示算法对自己得出的这个结果的把握。
也许我可以对同一个任务,用多个方法得到结果,最后选取“把握”最大的那个结果;
另一个很常见的方法是在诸如疾病诊断之类的场所,机器对于那些很容易分辨的情况(患病或者不患病的概率很高)可以自动区分,而对于那种很难分辨的情况,比如,49% 的概率患病,51% 的概率正常,如果仅仅简单地使用 50% 的阈值将患者诊断为“正常”的话,风险是非常大的,因此,在机器对自己的结果把握很小的情况下,会“拒绝发表评论”,而把这个任务留给有经验的医生去解决。
归纳偏置:【对事物的归纳总结】
在机器学习中,一个学习算法也会有一个前提假设,这里被称作“归纳偏执 (bias)”(bias 这个英文词在机器学习和统计里还有其他许多的意思)。
例如线性回归,目的是要找一个函数尽可能好地拟合给定的数据点,它的归纳偏执就是“满足要求的函数必须是线性函数”。一个没有归纳偏执的学习算法从某种意义上来说毫无用处,就像一个完全没有归纳能力的人一样,在第一次看到鱼的时候有人告诉他那是鱼,下次看到另一条鱼了,他并不知道那也是鱼,因为两条鱼总有一些地方不一样的,或者就算是同一条鱼,在河里不同的地方看到,或者只是看到的时间不一样,也会被他认为是不同的,因为他无法归纳,无法提取主要矛盾、忽略次要因素,只好要求所有的条件都完全一样──然而哲学家已经告诉过我们了:世界上不会有任何样东西是完全一样的,所以这个人即使是有无比强悍的记忆力,也绝学不到任何一点知识。
过拟合:
这个问题在机器学习中称作“过拟合 (Overfitting)”,例如前面的回归的问题,如果去掉“线性函数”这个归纳偏执,因为对于 N 个点,我们总是可以构造一个 N-1 次多项式函数,让它完美地穿过所有的这 N 个点,或者如果我用任何大于 N-1 次的多项式函数的话,我甚至可以构造出无穷多个满足条件的函数出来。如果假定特定领域里的问题所给定的数据个数总是有个上限的话,我可以取一个足够大的 N ,从而得到一个(或者无穷多个)“超级函数”,能够 fit 这个领域内所有的问题。然而这个(或者这无穷多个)“超级函数”有用吗?只要我们注意到学习的目的(通常)不是解释现有的事物,而是从中归纳出知识,并能应用到新的事物上,结果就显而易见了。
平衡点:
没有归纳偏执或者归纳偏执太宽泛会导致 Overfitting ,然而另一个极端──限制过大的归纳偏执也是有问题的:如果数据本身并不是线性的,强行用线性函数去做回归通常并不能得到好结果。难点正在于在这之间寻找一个平衡点。
不过人在这里相对于(现在的)机器来说有一个很大的优势:人通常不会孤立地用某一个独立的系统和模型去处理问题,一个人每天都会从各个来源获取大量的信息,并且通过各种手段进行整合处理,归纳所得的所有知识最终得以统一地存储起来,并能有机地组合起来去解决特定的问题。这里的“有机”这个词很有意思,搞理论的人总能提出各种各样的模型,并且这些模型都有严格的理论基础保证能达到期望的目的,然而绝大多数模型都会有那么一些“参数”(例如 K-means 中的 k ),通常没有理论来说明参数取哪个值更好,而模型实际的效果却通常和参数是否取到最优值有很大的关系,我觉得,在这里“有机”不妨看作是所有模型的参数已经自动地取到了最优值。另外,虽然进展不大,但是人们也一直都期望在计算机领域也建立起一个统一的知识系统(例如语意网就是这样一个尝试)。
开始GMM:
按照我们前面的讨论,作为一个流行的算法,GMM 肯定有它自己的一个相当体面的归纳偏执了。其实它的假设非常简单,顾名思义,Gaussian Mixture Model ,就是假设数据服从 Mixture Gaussian Distribution ,换句话说,数据可以看作是从数个 Gaussian Distribution 中生成出来的。实际上,我们在 K-means 和 K-medoids 两篇文章中用到的那个例子就是由三个 Gaussian 分布从随机选取出来的。实际上,从中心极限定理可以看出,Gaussian 分布(也叫做正态 (Normal) 分布)这个假设其实是比较合理的,除此之外,Gaussian 分布在计算上也有一些很好的性质,所以,虽然我们可以用不同的分布来随意地构造 XX Mixture Model ,但是还是 GMM 最为流行。另外,Mixture Model 本身其实也是可以变得任意复杂的,通过增加 Model 的个数,我们可以任意地逼近任何连续的概率密分布。
每个 GMM 由
个 Gaussian 分布组成,每个 Gaussian 称为一个“Component”,这些 Component 线性加成在一起就组成了 GMM 的概率密度函数:

那么如何用 GMM 来做 clustering 呢?其实很简单,现在我们有了数据,假定它们是由 GMM 生成出来的,那么我们只要根据数据推出 GMM 的概率分布来就可以了,然后 GMM 的
个 Component 实际上就对应了
个 cluster 了。根据数据来推算概率密度通常被称作 density estimation ,特别地,当我们在已知(或假定)了概率密度函数的形式,而要估计其中的参数的过程被称作“参数估计”。
现在假设我们有
个数据点,并假设它们服从某个分布(记作
),现在要确定里面的一些参数的值,例如,在 GMM 中,我们就需要确定
、
和
这些参数。 我们的想法是:
- 找到这样一组参数,它所确定的概率分布生成这些给定的数据点的概率最大,而这个概率实际上就等于
,我们把这个乘积称作似然函数 (Likelihood Function)。 - 通常单个点的概率都很小,许多很小的数字相乘起来在计算机里很容易造成浮点数下溢,因此我们通常会对其取对数,把乘积变成加和
,得到 log-likelihood function 。 - 接下来我们只要将这个函数最大化(通常的做法是求导并令导数等于零,然后解方程),亦即找到这样一组参数值,它让似然函数取得最大值,我们就认为这是最合适的参数,这样就完成了参数估计的过程。
下面让我们来看一看 GMM 的 log-likelihood function :

由于在对数函数里面又有加和,我们没法直接用求导解方程的办法直接求得最大值。【这就是为什么使用迭代逼近的原因之一】
为了解决这个问题,我们采取之前从 GMM 中随机选点的办法:分成两步,实际上也就类似于 K-means 的两步。
- 估计数据由每个 Component 生成的概率(并不是每个 Component 被选中的概率):对于每个数据
来说,它由第
个 Component 生成的概率为:
由于式子里的
和
也是需要我们估计的值,我们采用迭代法,在计算
的时候我们假定
和
均已知,我们将取上一次迭代所得的值(或者初始值)。 - 估计每个 Component 的参数:现在我们假设上一步中得到的
就是正确的“数据
由 Component
生成的概率”,亦可以当做该 Component 在生成这个数据上所做的贡献,或者说,我们可以看作
这个值其中有
这部分是由 Component
所生成的。集中考虑所有的数据点,现在实际上可以看作 Component 生成了
这些点。由于每个 Component 都是一个标准的 Gaussian 分布,可以很容易分布求出最大似然所对应的参数值:
其中
,并且
也顺理成章地可以估计为
。 - 重复迭代前面两步,直到似然函数的值收敛为止。
当然,上面给出的只是比较“直观”的解释,想看严格的推到过程的话,可以参考 Pattern Recognition and Machine Learning 这本书的第九章。有了实际的步骤,再实现起来就很简单了。Matlab 代码如下:

相对于之前 K-means 给出的结果,这里的结果更好一些,左下角的比较稀疏的那个 cluster 有一些点跑得比较远了。
当然,因为这个问题原本就是完全有 Mixture Gaussian Distribution 生成的数据,GMM (如果能求得全局最优解的话)显然是可以对这个问题做到的最好的建模。
【体现了GMM的优势】
K-means init + GMM
另外,从上面的分析中我们可以看到 GMM 和 K-means 的迭代求解法其实非常相似(都可以追溯到 EM 算法,下一次会详细介绍),因此也有和 K-means 同样的问题 ── 并不能保证总是能取到全局最优,如果运气比较差,取到不好的初始值,就有可能得到很差的结果。
对于 K-means 的情况,我们通常是重复一定次数然后取最好的结果,
不过 GMM 每一次迭代的计算量比 K-means 要大许多,一个更流行的做法是先用 K-means (已经重复并取最优值了)得到一个粗略的结果,然后将其作为初值(只要将 K-means 所得的 centroids 传入 gmm 函数即可),再用 GMM 进行细致迭代。
Gaussian Mixture Model 的迭代求解方法可以算是 EM 算法最典型的应用,而最开始说的 K-means 其实也可以看作是 Gaussian Mixture Model 的一个变种(固定所有的
,并令
即可)。
【方差无限小,某种意义上,是抹掉了“概率”这个所含信息的优势】
然而 EM 实际上是一种通用的算法(或者说是框架),可以用来解决很多类似的问题,我们最后将以一个中文分词的例子来说明这一点。
- Expectation 一步 - 对应于求关于后验概率的期望亦即
; - Maximization 一步 - 则对应于接下来的正常的最大似然的方法估计相关的参数亦即
、
和
。

EM为什么会work?
从偏理论一点的角度来简略地证明一下 EM 这个迭代的过程每一步都会对结果有所改进,除非已经达到了一个(至少是局部的)最优解。
现在我们的讨论将不局限于 GMM ,并使用一些稍微紧凑一点的符号。用
表示所有的 sample ,同时用
表示所有对应的隐含变量。我们的问题是要通过最大似然的方法估计出
中的参数
。
在这里我们假设这个问题很困难,但是要优化
却很容易。这就是 EM 算法能够解决的那一类问题。【加了个Z就完全不一样了,看上去】
现在我们引入一个关于隐含变量的分布
。注意到对于任意的
,我们都可以对 Log-likelihood Function 作如下分解:

其中
是分布
和
之间的 Kullback-Leibler divergence 。【分布之间的相似性度量 or KL-divergence,俗称KL距离,常用来衡量两个概率分布的距离】
由于 Kullback-Leibler divergence 是非负的,并且只有当两个分布完全相同的时候才会取到 0 。因此我们可以得到关系
,亦即
是
的一个下界。【其实就是对公式的理解】
换个角度看“E”:
现在考虑 EM 的迭代过程,记上一次迭代得出的参数为
,现在我们要选取
以令
最大,
由于
并不依赖于
,因此
的上限(在
固定的时候)是一个定值,它取到这个最大值的条件就是 Kullback-Leibler divergence 为零,亦即
等于后验概率
。
把它带入到
的表达式中可以得到

其中 const 是常量,而
则正是我们之前所得到的“同时包含了 sample 和隐含变量的 Log-likelihood function 关于后验概率的期望”,因此这个对应到 EM 中的“E”一步。
换个角度看“M”:
在接下来的“M”一步中,我们讲固定住分布
,再选取合适的
以将
最大化,这次其上界
也依赖于变量
,并会随着
的增大而增大(因为我们有前面的不等式成立)。
一般情况下
增大的量会比
要多一些,这时候 Kullback-Leibler divergence 在新的参数
下又不为零了,因此我们可以进入下一轮迭代,重新回到“E”一步去求新的
;
另一方面,如果这里 Kullback-Leibler divergence 在新的参数下还是等于 0 ,那么说明我们已经达到了一个(至少是局部的)最优解,迭代过程可以结束了。
小总结,以及扩展方案
上面的推导中我们可以看到每一次迭代 E 和 M 两个步骤都是在对解进行改进,因此迭代的过程中得到的 likelihood 会逐渐逼近(至少是局部的)最优值。
另外,在 M 一步中除了用最大似然之外,也可以引入先验使用 Maximum a Posteriori (MAP) 的方法来做。
还有一些很困难的问题,甚至在迭代的过程中 M 一步也不能直接求出最大值,这里我们通过把要求放宽──并不一定要求出最大值,只要能够得到比旧的值更好的结果即可,这样的做法通常称作 Generalized EM (GEM)。
【貌似,徐的变分推断也讲这个】
用 EM 来做中文分词
当然,一开始我们就说了,EM 是一个通用的算法,并不只是用来求解 GMM 。下面我们就来举一个用 EM 来做中文分词的例子。
这个例子来源于论文 Self-supervised Chinese Word Segmentation 。我在上次 MSTC 搜索引擎系列小课堂讲文本处理的时候提到过。这里为了把注意力集中到 EM 上,略去一些细节的东西,简单地描述一下基本的模型。
现在我们有一个字符序列
,并希望得到一个模型
(依赖于参数
)能用来将其词的序列
。按照生成模型的方式来考虑,将
看成是由
生成的序列的话,我们自然希望找到那个生成
的可能性最大的
,亦即通过最大似然的方式来估计参数
。
然而我们不知道似然函数
应该如何去优化,因此我们引入 latent variable
,如果我们知道
的话,似然值很容易求得:

其中
的值是从模型
中已知的。
但是现在我们不知道
的值,因此我们转而取
的后验概率的期望:【没有准确值,就用后验方法去估计,估计值就是期望值】

然后将这个期望针对
最大化即完成了 EM 的一次迭代。
具体的做法通常是先把一个初始文本(比如所有的
的集合)按照 N-gram 分割(N-gram 在 讲 K-medoids 的那篇文章中介绍过)为
,形成一个最初的辞典,而模型
的参数
实际上就描述了各个 N-gram 的概率
,初始值可以直接取为频率值。
有了辞典之后对于任意的
,我们可以根据辞典枚举出所有可能的分割
,而每个分割的后验概率
就是其中的单词的概率乘积。其他的就是标准的 EM 算法的迭代步骤了。
实际上在实际产品中我们会使用一些带了更多启发式规则的分词方法(比如 MMSEG),而这个方法更大的用处在于从一堆文本中“学习”出一个词典来(也就是
),并且这个过程是全自动的,主要有两个优点:
- 不需要人参与手工分词、标记等。
- 能自动从文本中学习,换句话说,你给它一些某个领域的专业文章,它能从中学习到这个领域的专业词汇来。
不管怎样,以上两点看起来都非常诱人。不过理论离实际往往还是有不少差距的。我不知道实际分词系统中有没有在用这样的方法来做训练的。之前我曾经用 Wikipedia (繁体中文)上抓下来的一些数据跑过一下小规模的试验,结果还可以,但是也并不如想像中的完美。因为当时也并没有弄明白 EM 是怎么一回事,加上这个过程本身计算负荷还是非常大的,也就没有深入下去。也许以后有时间可以再尝试一下。
总之,EM 是一种应用广泛的算法,虽然讲了点理论也讲了例子,但是一没有贴代码二没有贴图片,似乎实在是不像我的作风,不过精力和篇幅有限,也只能到此为止了。
【原作者很敬业】
Gaussian mixtures
Now we will fit a Gaussian mixture model to the data.
We don't need to write our own code here for one; instead we can use the convenient scikit-learn Gaussian Mixture models.
Goto: [Scikit-learn] 2.1 Gaussian mixture models & EM
以及如何恢复256维的问题。
通过GMM预测数字,采用了PCA的降维的思想。
# Make division default to floating-point, saving confusion from __future__ import division from __future__ import print_function # Necessary libraries import numpy as np import matplotlib.pyplot as plt from random import random from scipy.io import loadmat import seaborn as sns from collections import OrderedDict as odict # Put the graphs where we can see them %matplotlib inline sns.set(style="ticks") # easier debugging display np.set_printoptions(edgeitems=5, precision=3, suppress=False) from pprint import pprint
数据:usps_gmm_3d.mat
data = loadmat('./usps_gmm_3d.mat')
xtrain3d = data['xtrain3d'] # 1500*3
ytrain = data['ytrain'] # 1500*1
xtest3d = data['xtest3d']
ytest = data['ytest']
pca_mu = data['mu']
pca_e = data['E']
data.keys()
del(data)
理解:
一个样本的数据本来是28*28维度。这里通过PCA降维后,变成了三维的数据。过程原理详见PCA。
添加了index策略,容易提取特定数据
x2test = (ytest==2).ravel()
x3test = (ytest==3).ravel()
x5test = (ytest==5).ravel()
x2train = (ytrain==2).ravel()
x3train = (ytrain==3).ravel()
x5train = (ytrain==5).ravel()
#x2test
#output: array([False, False, False, ..., True, True, False], dtype=bool)
# 可提取特定的数据
xtest3d[x2test,:]
3D: 一个样本是三维,所以可用3D space表示:
# enable interactive plots
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# 这里的星号是什么个意思?
ax.scatter(*xtrain3d[x2train, :].T, alpha=0.3, label='2', c='red');
ax.scatter(*xtrain3d[x3train, :].T, alpha=0.3, label='3', c='green');
ax.scatter(*xtrain3d[x5train, :].T, alpha=0.3, label='5', c='blue');
ax.legend();
生成的这个图是可以用鼠标点击翻转的:)

2D: 只表示两维的数据:
# disable interactive plots
%matplotlib inline
plt.scatter(*xtrain3d[x2train, 0:2].T, alpha=0.3, label='2');
ax.legend();

Train:
这里使用的是新的api: mixture.GMM拟合数据:
from sklearn import mixture
from collections import OrderedDict as odict
models = odict()
for n_components in (1, 2, 4, 6, 8, 10):
gmm = mixture.GMM(
n_components=n_components,
covariance_type='full')
gmm.fit(xtrain3d[x2train, :])
models[n_components] = gmm
# let's look at the 6-component mixture as an example:
gmm = models[6]
print('cluster weights: {}'.format()) # 混合比例
print('cluster centres: {}'.format(gmm.means_))
print('cluster covars : {}'.format(gmm.covars_))
cluster weights:
[ 0.14066026 0.09201028 0.18819883 0.08630041 0.22760965 0.26522056]
cluster centres:
[[ 525.34206155 229.82762931 -84.56957136] [ 668.83752127 -211.13219553 -556.70353001] [ 708.82292215 97.37435102 -162.24295728] [ -67.54755363 -498.65315683 205.45762409] [ 184.48129824 -225.04402981 -259.68106114] [ 468.51537602 -449.72214396 16.11527165]]
cluster covars :
[[[ 28844.36131231 2248.0120006 -16091.78553269] [ 2248.0120006 15739.98515107 -9287.48344138] [ -16091.78553269 -9287.48344138 52983.91662063]] [[ 20536.87170667 31705.03194212 -13878.5693464 ] [ 31705.03194212 75164.4199765 -14529.40953692] [ -13878.5693464 -14529.40953692 100473.7693028 ]] [[ 15275.85814349 2454.43190215 -21078.0571328 ] [ 2454.43190215 13564.86221594 2581.56705871] [ -21078.0571328 2581.56705871 56779.47115667]] [[ 16530.12269771 4728.04364736 18134.65810034] [ 4728.04364736 27730.24884845 6822.27940972] [ 18134.65810034 6822.27940972 40211.56343292]] [[ 45326.27458989 15412.885695 -7565.9600019 ] [ 15412.885695 106608.10060722 385.32567872] [ -7565.9600019 385.32567872 70784.44410959]] [[ 17985.11679955 656.5057687 -18187.7813835 ] [ 656.5057687 45790.83090656 3976.8033297 ] [ -18187.7813835 3976.8033297 60795.35810013]]]
Prediction:
# let's look at the predict_proba method on 10 test samples.
# The result is a 6-by-10 array
gmm.predict_proba(xtest3d[x2test, :][0:10])
# 得到结果类似softmax
gmm.predict_proba(xtest3d[x2test, :][0:10]).sum(1)
# 通过这个验证可得,Pr之和当然为1
如果要看所有6个假设的综合结果:
for n, gmm in models.items():
print('test subset log likelihood with {} components: {}'.format(n, gmm.score(xtest3d[x2test, :][0:10])))
The score method calculates the average log likelihood of a given sample.
test subset log likelihood with 1 components:
[-20.15743676 -20.89899758 -20.57470118 -21.38488846 -21.80538554 -22.85135431 -19.99027883 -20.69164439 -20.39108309 -20.57871817]
test subset log likelihood with 2 components:
[-20.50944149 -20.85671343 -19.47743735 -20.51780262 -21.49312112 -22.32285379 -20.51512882 -21.23519172 -20.82624165 -19.56157603]
test subset log likelihood with 4 components:
[-20.82218924 -20.81328986 -19.62220702 -21.3270596 -21.65611915 -22.36649503 -20.73613886 -20.92524886 -21.02064744 -19.510683 ]
test subset log likelihood with 6 components:
[-21.28179157 -21.62172376 -19.28294185 -20.64373138 -21.64069071 -22.21199095 -20.10286151 -21.00646018 -20.31653127 -19.08305399]
test subset log likelihood with 8 components:
[-21.21808576 -21.51171094 -19.12618852 -20.5857113 -21.07915049 -22.12327529 -20.11120626 -21.26127906 -20.59065857 -19.07384354]
test subset log likelihood with 10 components:
[-20.99695722 -21.77845886 -19.05018406 -20.40976237 -21.21723584 -21.81931593 -20.04737598 -21.10808762 -20.5069136 -18.94941105]
还原数字图像:

浙公网安备 33010602011771号