高斯混合模型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

产生原点处的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

叠加两部分数据：

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

训练数据是一个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

默认迭代的次数为100次

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

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

密度估计

随后要做的事情，是生成模型。 
调用方式为：

clf.fit(X_train)
 
 
1
 
 
1

里面的关键步骤，就看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

第一个函数就是计算后验概率的函数， 
其中作者对 

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

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

logprob = logsumexp(lpr, axis=1)
 
 
1
 
 
1

计算结果减法完成：

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

函数的返回值中， 
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

前面已经说了： 
ϕ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

注意每一行两个分量的权重求和为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

更新协方差矩阵： 

Σ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

这里实际调用的函数为：

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

具体可见参考文献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