高斯混合模型EM算法

转自：
http://blog.csdn.net/dudubird90/article/details/49759431
k-means算法对于数据点和clusters之间的关系，是all-or-nothing的关系，这是一个hard decision，往往会导致的局部最小值，这不是理想的求解。一种常见的做法，是学习这个协方差矩阵，而不是固定它们为单位矩阵。

ＧＭＭ模型及算法流程

ＧＭＭ的全称是Gaussian Mixture Model，即高斯混合模型。

假设我们有一个训练集{x1,...,xm},在非监督学习中，这些数据都是没有标签的。我们希望可以对数据进行建模，通过定义一个联合分布p(x(i),z(i))=p(x(i)|z(i))p(z(i)) ,这里z(i)∼Multinomial(ϕ),其中ϕj≥0,∑kj=1ϕj=1,
这个参数ϕj代表着p(z(i)=j)的概率。

而xi|zi=j∼(μj,Σj),ｋ代表着标签z(i)可以有的取值的数目。所以我们的模型实际上是当z(i)从{1,...,k}随机取值时，从ｋ个高斯分量中选一个，再根据相应的z(i)生成对应的x(i)。

z(i)是一个隐变量，就是说它是隐藏的，不可见的。

我们的似然函数可以这样写：

(ϕ,μ,Σ)=∑i=1mlog p(x(i);ϕ,μ,Σ)=∑i=1mlog ∑z(i)=1kp(x(i)|z(i);μ,Σ)p(z(i);ϕ)

如果直接通过令偏导数为０，对其进行求解是无法得到closed form的。

这个隐变量告诉了我们，数据x(i)是从ｋ个高斯分量中的哪一个来的。如果我们一开始就知道是哪个，那么最大似然问题是很简单的。具体而言，我们就可以将似然函数展开：

(ϕ,μ,Σ)=∑i=1mlog p(x(i)|z(i);μ,Σ)+log p(z(i);ϕ)

随后分别对ϕ,μ,Σ求导，可以得到：

⎧⎩⎨⎪⎪⎪⎪ϕj=1m∑mi=1I{z(i)=j}μj=∑mi=1I{z(i)=j}x(i)∑mi=1I{z(i)=j}Σj=∑mi=1I{z(i)=j}x((i)−μj)(x(i)−μj)T∑mi=1I{z(i)=j}

当我们知道类别信息，事实上这就是一个高斯判别分析的参数估计问题。但是我们知道，所以我们可以分两个步骤:

一个步骤猜测这个zi的值，也就是Ｅ-step;

一个步骤利用已知的类别信息来估计高斯的参数μ,Σ，以及各个高斯分量的权重ϕj。

那么这个标签值要如何猜测呢? 我们可以在固定参数ϕ,μ,Σ的情况下，在给定数据x(i)时z(i)后验概率的来判定，再结合贝叶斯定律：

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

EM算法在k-means 聚类中也有用到，但是用的是hard decision， 而利用GMM，我们使用的其实是soft的 assignments w(i)j.
同样它也可能会收敛到局部极小值，所以需要尝试不同的初始参数。

用GMM来进行聚类，K个高斯component实际上就对应着K个cluster，我们根据数据来推算概率密度density estimation。

利用Mixtures of Gaussians模型来实现聚类算法的流程总结如下：

输入：

{x(i)}mi=1 (数据)，
k(聚类的数目)

初始化：

∀μi←1k

Σi←var({x(i)})

Repeat until convergence :

E-step:

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

w(i)j:=p(z(i)=j|x(i);ϕ,μ,Σ)

M-step:

ϕj:=1m∑i=1mw(i)j
,

μj:=∑mi=1w(i)jx(i)∑mi=1w(i)j

Σj:=∑mi=1w(i)j(x(i)−μj)(x(i)−μj)T∑mi=1w(i)j

输出：

{z(i)←argmaxj∈{1,...,k}w(i)j}mi=1

算法代码demo及关键步骤注解

大费苦心写了编辑了这么多公式，还没有进入正题。本文一大重点是要研读scikit包中gaussian mixture 部分的代码，已进一步熟悉算法和这个工具包。

官网给出的密度估计的源码可见：

http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_pdf.html#example-mixture-plot-gmm-pdf-py

算法就是随机产生两个不同中心的高斯分布的数据点，再用mixture.GMM进行拟合。

关键部分代码拆分如下。

1. 准备数据

产生以20,20 为中心的满足标准正态分布的300个样本点。

shifted_gaussian = np.random.randn(n_samples, 2) + np.array([20, 20])

1
1
[/code]

产生原点处的stretch后高斯分布的300个样本点。通过乘以一个矩阵来完成。

C = np.array([[0., -0.7], [3.5, .7]])
stretched_gaussian = np.dot(np.random.randn(n_samples, 2), C)

1
2
1
2
[/code]

叠加两部分数据：

X_train = np.vstack([shifted_gaussian, stretched_gaussian])

1
1
[/code]

训练数据是一个600*2的矩阵。

2. 初始化高斯混合模型

clf = mixture.GMM(n_components=2, covariance_type=’full’)

第二个参数这样设置，实际得到的协方差的参数 covars_是

(n_components, n_features, n_features)

这样允许了不同的高斯分量有不同的协方差矩阵。

初始化就做了这些事情：

def __init__(self, n_components=1, covariance_type='diag',
random_state=None, thresh=1e-2, min_covar=1e-3,
n_iter=100, n_init=1, params='wmc', init_params='wmc'):
self.n_components = n_components
self.covariance_type = covariance_type
self.thresh = thresh
self.min_covar = min_covar
self.random_state = random_state
self.n_iter = n_iter
self.n_init = n_init
self.params = params
self.init_params = init_params

1
2
3
4
5
6
7
8
9
10
11
12
1
2
3
4
5
6
7
8
9
10
11
12
[/code]

默认迭代的次数为100次

self.weights_ = np.ones(self.n_components) / self.n_components

1
1
[/code]

每个分量一开始赋值的权重是相同的，这里的例子中有两个分量，所以得到的是一个行向量[0.5,0.5].

密度估计

随后要做的事情，是生成模型。

调用方式为：

clf.fit(X_train)

1
1
[/code]

里面的关键步骤，就看Expection Step和Maximization Step。

Expection Step:

curr_log_likelihood, responsibilities = self.score_samples(X)
log_likelihood.append(curr_log_likelihood.sum())

# Check for convergence.
if i > 0 and abs(log_likelihood[-1] - log_likelihood[-2]) < \
self.thresh:
self.converged_ = True
break

1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
[/code]

第一个函数就是计算后验概率的函数，

其中作者对

p(z(i)=j|x(i);ϕ,μ,Σ)=p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

表示成exp(log(.))的形式，于是把除法就转换成了减法。

下面这一步是计算分子。相当于计算log(p(x(i)|z(i)=j;μ,Σ))+log(p(z(i)=j;ϕ))

lpr = (log_multivariate_normal_density(X, self.means_, self.covars_,
self.covariance_type)
+ np.log(self.weights_))

1
2
3
1
2
3
[/code]

上面得到的矩阵为600*2，一列对应着一个高斯分量。下面这一步是计算分母，先用exp得到原值，再求和，随后再取log相当于计算log∑kl=1p(x(i)|z(i)=j;μ,Σ)p(z(i)=j;ϕ)

万分精简的只用了一个函数就完成了：

logprob = logsumexp(lpr, axis=1)

1
1
[/code]

计算结果减法完成：

responsibilities = np.exp(lpr - logprob[:, np.newaxis])

1
1
[/code]

函数的返回值中，

responsibilities : array_like, shape (n_samples, n_components)

Posterior probabilities of each mixture component for each

observation

而curr_log_likelihood则是分母的值，作者通过记录这个值所有600个样本的和在每一次迭代中的变化来与阈值self.thresh比较来判断是否达到了收敛。

Maximization step

def _do_mstep(self, X, responsibilities, params, min_covar=0):
""" Perform the Mstep of the EM algorithm and return the class weihgts.
"""
weights = responsibilities.sum(axis=0)
weighted_X_sum = np.dot(responsibilities.T, X)
inverse_weights = 1.0 / (weights[:, np.newaxis] + 10 * EPS)

if 'w' in params:
self.weights_ = (weights / (weights.sum() + 10 * EPS) + EPS)
if 'm' in params:
self.means_ = weighted_X_sum * inverse_weights
if 'c' in params:
covar_mstep_func = _covar_mstep_funcs[self.covariance_type]
self.covars_ = covar_mstep_func(
self, X, responsibilities, weighted_X_sum, inverse_weights,
min_covar)
return weights

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
[/code]

前面已经说了：

ϕj代表着p(z(i)=j)的概率，之类对应着weights_,
每个分量的权重。

计算的方法是

ϕj:=1m∑i=1mw(i)j
,

也等于600个样本的分量1的后验概率和/（600个样本的分量1的后验概率和+600个样本的分量2的后验概率和）。

weights = responsibilities.sum(axis=0)
if 'w' in params:
self.weights_ = (weights / (weights.sum() + 10 * EPS) + EPS)

1
2
3
1
2
3
[/code]

注意每一行两个分量的权重求和为1,600个全部加起来为600，正好等于左右两列归总值的和。所以这么计算也是对的。

更新均值：

μj:=∑mi=1w(i)jx(i)∑mi=1w(i)j

%得到2*2的矩阵，（1,1）代表分量1第一维，（1,2）代表分量1第二维
weighted_X_sum = np.dot(responsibilities.T, X)
weights = 1.0 / (weights[:, np.newaxis] + 10 * EPS)
if 'm' in params:
self.means_ = weighted_X_sum * inverse_weights % 同样是2*2的矩阵

1
2
3
4
5
1
2
3
4
5
[/code]

更新协方差矩阵：

Σj:=∑mi=1w(i)j(x(i)−μj)(x(i)−μj)T∑mi=1w(i)j

if 'c' in params:
covar_mstep_func = _covar_mstep_funcs[self.covariance_type]
self.covars_ = covar_mstep_func(
self, X, responsibilities, weighted_X_sum, inverse_weights,
min_covar)

1
2
3
4
5
1
2
3
4
5
[/code]

这里实际调用的函数为：

def _covar_mstep_full(gmm, X, responsibilities, weighted_X_sum, norm,
min_covar):
"""Performing the covariance M step for full cases"""
# Eq. 12 from K. Murphy, "Fitting a Conditional Linear Gaussian
# Distribution"
n_features = X.shape[1]
cv = np.empty((gmm.n_components, n_features, n_features))
for c in range(gmm.n_components):
post = responsibilities[:, c]
# Underflow Errors in doing post * X.T are  not important
np.seterr(under='ignore')
avg_cv = np.dot(post * X.T, X) / (post.sum() + 10 * EPS)
mu = gmm.means_[c][np.newaxis]
cv[c] = (avg_cv - np.dot(mu.T, mu) + min_covar * np.eye(n_features))
return cv

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
[/code]

具体可见参考文献2的公式12. 其实我看参考文献2脱离上下文不怎么容易懂，基本上应该说这是协方差矩阵的另一种写法。https://en.wikipedia.org/wiki/Covariance_matrix

Σ=E[(X−EX)(X−EX)T]=E(XXT)−μμT

所以上面的公式就进一步写成：

Σj:=∑mi=1w(i)jx(i)x(i)T∑mi=1w(i)j−μjμTj

下一次再来梳理关键步骤的推导吧，最后的结果是好看的。



参考文献：

１. andrew ng cs229 handouts

http://cs229.stanford.edu/notes/cs229-notes7b.pdf

2. K. Murphy, “Fitting a Conditional Linear Gaussian Distribution”

http://www.cs.ubc.ca/~murphyk/Papers/learncg.pdf