python 混合高斯分布 混合高斯模型聚类分析

高斯混合聚类和k 均值算法（k-means）都属于原型聚类，但与k均值用原型向量来刻画聚类结构不同，高斯混合聚类采用概率模型来表达聚类原型。

一、混合模型（Mixture Model）

混合模型是一个可以用来表示在总体分布（distribution）中含有 K 个子分布的概率模型。换句话说，混合模型表示了观测数据在总体中的概率分布，它是一个由 K 个子分布组成的混合分布。混合模型不要求观测数据提供关于子分布的信息，来计算观测数据在总体分布中的概率。

二、高斯混合模型（Gaussian Mixture Model,GMM）

1. 单高斯模型

• 当样本数据 X 是一维数据（Univariate）时，高斯分布(也叫正态分布)遵从下方概率密度函数（Probability Density Function）：

$P(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^2}}exp(-\frac{(x-\mu)^2}{2\sigma^2})$

正态分布通常也记为 $N(\mu,{\sigma}^2)$，其中 $\mu$ 为数据均值（期望）， $\sigma$ 为数据标准差（Standard deviation）。

• 当样本数据 X 是多维数据（Multivariate）时，（多元）高斯分布遵从下方概率密度函数：

$P(x|\mu,\Sigma)=\frac{1}{(2\pi)^{D/2}|\Sigma|^{1/2}}exp(-\frac{1}{2}{(x-\mu)^T \Sigma^{-1}(x-\mu)})$

其中， $\mu$ 为数据均值（期望）， $\Sigma$ 为协方差矩阵，D 为数据维度。注意这里的 T 是矩阵转置符号，多维数据 $x-\mu$

2. 高斯混合模型

高斯混合模型可以看作是由 K 个单高斯模型组合而成的模型，每个 Gaussian 称为一个“Component”，这 K 个子模型是混合模型的隐变量（Hidden variable）。也就是说我们假设样本中的每个数据都由这 K 个子模型中的某一个生成。一般来说，一个混合模型可以使用任何概率分布，这里使用高斯混合模型是因为高斯分布具备很好的数学性质以及良好的计算性能。下面给出高斯混合模型的数学定义：

首先定义如下信息：

• $x_j$
• k 是混合模型中子高斯模型的编号，K是所有子模型的数量， k = 1,2,…,K
• $\alpha_k$ 是观测数据属于第 k 个子模型的概率，它其实是一个“混合系数”，和下面的$\gamma_{jk}$相比，$\alpha_k$相当于说先从k个子分布中选则一个子分布的概率。 $\alpha_k \geq 0,\sum_{i=1}^{K}{\alpha_k} = 1$
• $\gamma_{jk}$ 表示第 j 个观测数据$x_j$属于第 k 个子模型的概率
• $\phi(x|\theta_k)$ 是第 k 个子模型的高斯分布密度函数， $\theta_k = (\mu_k, \sigma_k^{2})$

那么，高斯混合模型的概率分布为：

$\begin{array}{ll} P(x) & =\sum_{k=1}^{K} \alpha_k \phi(x|\theta_k)\\ & =\sum_{k=1}^{K} \alpha_k \frac{1}{(2\pi)^{D/2}|\Sigma|_k^{1/2}}exp(-\frac{(x-\mu_k)^T \Sigma_k^{-1}(x-\mu_k)}{2})\\ & =\sum_{k=1}^{K} \alpha_k \cdot p(x|\mu_k,\Sigma_k)\\ \end{array}$

对于这个高斯混合模型（$\theta_k = ({\mu_k}, {\Sigma_k}, {\alpha_k})$）而言，参数 $\mu_k$ 表示每个子模型的期望（均值向量）、$\Sigma_k$表示协方差矩阵、${\alpha_k}$

根据上面的式子，如果我们要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：

• 首先随机地在这 K 个 Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数 $\alpha_k$
• 选中了 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

那么如何用 GMM 来做聚类呢？其实很简单，现在我们有了数据，假定它们是由 GMM 生成出来的，那么我们只要根据数据推出 GMM 的概率分布来就可以了，然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算概率密度通常被称作“密度估计(density estimation)” ，特别地，当我们在已知（或假定）了概率密度函数的形式，而要估计其中的参数的过程被称作“参数估计”。

GMM与K-means比较

• 相同点：都是可用于聚类的算法；都需要指定K值。
• 不同点：GMM可以给出一个样本属于某类的概率是多少。

三、极大似然估计

上面只是给出了 GMM 的概率密度，但是它的参数还是未知的，那么如何该如何求概率密度的参数呢？

先来聊聊似然函数。现在假设我们有 N 个数据点，并假设它们服从某个分布（记作 p(x) ），要确定里面的一些参数的值。现在这里的分布就是高斯混合分布，我们就需要确定 $\alpha_k,\mu_k$ 和$\Sigma_k$ 这些参数。 我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于 $\prod_{i=1}^N p(x_i)$ ，我们把这个乘积称作似然函数 (Likelihood Function)。

通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和 $\sum_{i=1}^N \log p(x_i)$，得到对数似然函数（log-likelihood function） 。接下来我们只要将这个函数最大化（通常的做法是求导并令导数等于零，然后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，我们就认为这是最合适的参数，这样就完成了参数估计的过程。下面让我们来看一看 GMM 的对数似然函数的具体形式 ：

$l(\Theta)=ln\prod_{j=1}^{N}p(x_{i};\Theta)=\sum_{j=1}^{N}ln(\sum_{t=1}^{K}\alpha_{t} \cdot p(x_{j};\mu_{t},\Sigma_{t}))$

由于在对数函数里面又有加和，我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们要用到 EM 算法。

四、EM算法

EM算法全程是 Expectation Maximization，即期望最大化算法，专门用来迭代求解极大似然估计。我们就是要找到最佳的模型参数，使得期望最大，“期望最大化算法”名字由此而来。

EM算法每次迭代包含两个步骤：

• E-step：求期望 $E(\gamma_{jk} | X, \theta)$
• M-step：求极大，计算新一轮迭代的模型参数。

这里不具体介绍一般性的 EM 算法（通过 Jensen 不等式得出似然函数的下界 Lower bound，通过极大化下界做到极大化似然函数），只介绍怎么在高斯混合模型里应用从来推算出模型参数。通过 EM 迭代更新高斯混合模型参数的方法如下：

假设：我们有样本数据 $x_{1}, x_{2}, ...,x_{N}$

1. 首先初始化参数

2. E-step：依据当前参数，计算每个数据 j 来自子模型 k 的可能性

首先定义$\gamma_{jk}$ 表示估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率 $\alpha_k$ ）：对于给定的样本 $x_j$ 来说，$z_j$是生成这个样本$x_j$的子高斯分布，那么$z_j$的先验概率分布（大p）$P(z_j=k)$ 就是前面定义的$\alpha_k$表示属于第k个子分布的概率。根据贝叶斯定理，$z_j$的后验分布为：

$\begin{array}{ll} \gamma_{jk} & =PM(z_j=k|x_j)\\ & =\frac{P(z_j=k) \cdot PM(x_j|z_j=k)}{PM(x_j)}\\ & = \frac{\alpha_{k} \cdot p(x_{j}|\mu_k,\Sigma_k)}{\sum_{t=1}^{K}{\alpha_{k} \cdot p(x_{j}|\mu_t,\Sigma_t)}}, j = 1,2,...,N; t = 1,2,...,K\\ \end{array}$

也就是说，PM 表示给定样本数据$x_j$的情况下，它是由第k个Component生成的后验概率。（小p）$p$就是高斯混合概率分布。

此外如果要做聚类的话，设每个$x_j$样本的簇标记为$\lambda_j$，则：

$\lambda_j=argmax \: \gamma_{jk},k \in\{1,2,...,K\}$

3. M-step：计算新一轮迭代的模型参数，并更新

估计每个 Component 的参数：现在我们假设上一步中得到的 $\gamma(j, k)$ 就是正确的“数据 $x_j$ 由 Component k 生成的概率”，亦可以当做该 Component 在生成 $x_j$ 这个数据上所做的贡献。或者说，我们可以看作 $x_j$ 这个值其中有 $\gamma(j, k)x_j$ 这部分是由 Component k 所生成的。集中考虑所有的数据点，现在实际上可以看作 Component 生成了 $\gamma(1, k)x_1, \ldots, \gamma(N, k)x_N$ 这些点。由于每个 Component 都是一个标准的 Gaussian 分布，可以很容易求出最大似然所对应的参数($\mu,\Sigma,\alpha$)的值（分别对似然函数求导，见下面的详细推导部分）：

$\mu_{k} = \frac{\sum_{j}^{N}{(\gamma_{jk}}x_{j})}{\sum_{j}^{N}{\gamma_{jk}}}, k=1,2,...,K$

$\Sigma_{k} = \frac{\sum_{j}^{N}{\gamma_{jk}}(x_{j}-\mu_{k})(x_{j}-\mu_{k})^{T}}{\sum_{j}^{N}{\gamma_{jk}}}, k = 1,2,...,K$

$\alpha_{k} = \frac{\sum_{j=1}^{N}{\gamma_{jk}}}{N}, k=1,2,...,K$

3.1 $\mu_k$

首先，对于给定的样本集D,GMM的对数似然函数$l(D)$为：
$\begin{array}{ll} l(D) & =l(\alpha,\mu,\Sigma)\\ & =\sum_{j=1}^{N}ln(\sum_{t=1}^{K}\alpha_{t} \cdot p(x_j|\mu_{t},\Sigma_{t}))\\ \end{array}$

具体看看如何求$\mu_k$。首先，若参数{$(\alpha_k,\mu_k,\Sigma_k)|1\leq k \leq K$}能使上式最大化，则对$\mu_k$求导可得：$\frac{\partial l(D)}{\partial \mu_k}=0$。则：

$\begin{array}{ll} \frac{\partial l(D)}{\partial \mu_k} & = \frac{\partial l(\alpha,\mu,\Sigma)}{\partial \mu_k} \\ & =\frac{\partial}{\partial \mu_k} \sum_{j=1}^{N}ln(\sum_{t=1}^{K}\alpha_{t} \cdot p(x_j|\mu_{t},\Sigma_{t}))\\ & =\sum_{j=1}^{N} \frac{\alpha_k \cdot \partial p(x_j|\mu_k,\Sigma_k)/\partial \mu_k}{\sum_{t=1}^{K}\alpha_{t} \cdot p(x_{j}|\mu_{t},\Sigma_{t})}\\ & =\sum_{j=1}^{N} \frac{\alpha_k \cdot \frac{1}{(2\pi)^{n/2} |\Sigma_t|^{1/2} } \cdot e^{-(1/2)(x_j-\mu_k)^T \Sigma_k^{-1}(x_j-\mu_k)}}{\sum_{t=1}^{K}\alpha_{t} \cdot p(x_{j}|\mu_{t},\Sigma_{t})} \cdot \frac{\partial }{\partial \mu_k}(-1/2)(x_j-\mu_k)^T \Sigma_t^{-1}(x_j-\mu_k)\\ & =\sum_{j=1}^{N} \frac{\alpha_k \cdot p(x_j|\mu_k,\Sigma_k)}{\sum_{t=1}^{K}\alpha_{t} \cdot p(x_{j}|\mu_{t},\Sigma_{t})} (x_j-\mu_k)\\ & =0\\ \end{array}$

注意：

• 第三行等号右边，分子因为除$\mu_t$之外都不含$\mu_t$，所以原来的求和只剩一项。
• 一个矩阵乘A以一个向量x，对向量求导等于矩阵的逆$A^T$，也就是$\frac{\partial Ax}{\partial x}=A^T$。

又因为 $\gamma_{jk}=\frac{\alpha_{k} \cdot p(x_{j}|\mu_k,\Sigma_k)}{\sum_{t=1}^{K}{\alpha_{k} \cdot p(x_{j}|\mu_t,\Sigma_t)}}$，所以: $\sum_{j=1}^{N} \gamma_{jk}(x_j-\mu_k)=0$，里面的式子乘法分配率一下，再移动一项到等号右边得：$\sum_{j=1}^{N} \gamma_{jk} \cdot x_j=\sum_{j=1}^{N} \gamma_{jk} \cdot \mu_k$，于是最终得到：

$\mu_k=\frac{\sum_{j=1}^{N} \gamma_{jk} \cdot x_j}{\sum_{j=1}^{N} \gamma_{jk} }$

这个式子也可以理解为各混合成分的均值可以通过加权平均来计算，样本权重是每个样本属于该成分的后验概率$\gamma$。

3.2 $\Sigma_{k}$

类似的，可以得到$\Sigma_{k} = \frac{\sum_{j}^{N}{\gamma_{jk}}(x_{j}-\mu_{k})(x_{j}-\mu_{k})^{T}}{\sum_{j}^{N}{\gamma_{jk}}}, k = 1,2,...,K$

3.3 $\alpha_k$

除了要最大化$l(D)$之外，还需满足$\alpha \geq 0, \sum_{k=1}^{K} \alpha_k=1$，考虑拉格朗日乘子法：$\psi(D,\lambda) = l(D)+\lambda(\sum_{k=1}^{K} \alpha_k-1)$，再对 $\alpha_k$

$\sum_{j=1}^{N}\frac{p(x_j|\mu_k,\Sigma_k)}{\sum_{t}^{K}\alpha_t \cdotp(x_j|\mu_t,\Sigma_t)} + \lambda=0$

等式两边同时乘以$\alpha_k$得：$\sum_{j=1}^{N}\gamma_{jk}+\alpha_k\lambda=0$，然后对所有混合成分（子分布）求和得：$\sum_{k=1}^{K}\sum_{j=1}^{N}\gamma_{jk}+\sum_{k=1}^{K}\alpha_k\lambda=0=>N+\lambda=0$，于是得到$\lambda=-N$。再代入前面的这个式子可得：

$\alpha_k=\frac{\sum_{j=1}^{N}\gamma_{jk}}{N}$

这个式子也可以理解为每个高斯成分的混合系数由（所有）样本属于该成分的后验概率确定。

至此，我们就得到了高斯混合模型的三个参数。

这样我们可以根据上一步得到的$\gamma_{jk}$（E步），在这一步中更新模型的三个参数的值（M步。

4. 重复计算 E-step 和 M-step 直至收敛

若达到最大迭代次数，或者似然函数$l(D)$增长很少甚至不再增长，则得到最终的高斯混合模型。否则重复2，3两步。得到模型之后可以根据高斯混合分布对簇进行划分，也就是根据$argmax \lambda_{jk}$

需要注意的是，EM 算法具备收敛性，但并不保证找到全局最大值，有可能找到局部最大值。解决方法是初始化几次不同的参数进行迭代，取结果最好的那次。

5. GMM 实践

首先是生成随机的4个簇。然后使用 sklearn 内置的 GaussianMixture 模块进行聚类。

Python 代码实现

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from sklearn import mixture

# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，簇中心在[-1,-1], [0,0],[1,1], [2,2]， 簇方差分别为[0.4, 0.2, 0.2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.2, 0.2, 0.2], 
                  random_state =9)
plt.scatter(X[:, 0], X[:, 1], marker='o')
print(X.shape)
plt.show()

# concatenate the two datasets into the final training set
X_train = X #np.vstack([shifted_gaussian, stretched_gaussian])

# fit a Gaussian Mixture Model with two components
clf = mixture.GaussianMixture(n_components=4, covariance_type='full')
clf.fit(X_train)

Z = -clf.predict(X)
score= metrics.calinski_harabaz_score(X, y_pred) 

plt.scatter(X_train[:, 0], X_train[:, 1],c=Z)
plt.text(.99, .01, ('score: %.2f' % (score)),
                 transform=plt.gca().transAxes, size=10,
                 horizontalalignment='right')

plt.title('Negative log-likelihood predicted by a GMM')
plt.axis('tight')
plt.show()

看一下效果：

---

参考：

[0]. 高斯混合模型

[1]. 高斯混合模型（GMM）

[2]. 漫谈 Clustering (3): Gaussian Mixture Model

[3]. EM及高斯混合模型

[4]. 如何通俗理解EM算法

[5]. 《机器学习（周志华）》，第9.4.3节：高斯混合聚类

---

THE END.