python GMM二维高斯混合模型代码高斯混合模型算法步骤

转载

mob64ca14017c37 2024-01-29 17:41:38

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∼N(μj,Σj),ｋ代表着标签z(i)可以有的取值的数目。所以我们的模型实际上是当z(i)从{1,...,k}随机取值时，从ｋ个高斯分量中选一个，再根据相应的z(i)生成对应的x(i)。
z(i)是一个隐变量，就是说它是隐藏的，不可见的。
我们的似然函数可以这样写：

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

如果直接通过令偏导数为０，对其进行求解是无法得到closed form的。
这个隐变量告诉了我们，数据x(i)是从ｋ个高斯分量中的哪一个来的。如果我们一开始就知道是哪个，那么最大似然问题是很简单的。具体而言，我们就可以将似然函数展开：

L(ϕ,μ,Σ)=∑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])

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

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

叠加两部分数据：

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

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

默认迭代的次数为100次

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

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

密度估计

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

clf.fit(X_train)

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

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

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_))

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

logprob = logsumexp(lpr, axis=1)

计算结果减法完成：

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

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

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

更新协方差矩阵：

Σ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)

这里实际调用的函数为：

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

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

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

python GMM二维高斯混合模型代码高斯混合模型算法步骤_协方差矩阵