G M M 高 斯 混 合 模 型 ( 2 0 2 0 )

高斯混合模型GMM与EM算法的Python实现

GMM与EM算法的Python实现

高斯混合模型(GMM)是一种常用的聚类模型，通常我们利用最大期望算法(EM)对高斯混合模型中的参数进行估计。

1. 高斯混合模型(Gaussian Mixture models, GMM)

高斯混合模型（Gaussian Mixture Model，GMM）是一种软聚类模型。GMM也可以看作是K-means的推广，因为GMM不仅是考虑到了数据分布的均值，也考虑到了协方差。和K-means一样，我们需要提前确定簇的个数。

GMM的基本假设为数据是由几个不同的高斯分布的随机变量组合而成。如下图，我们就是用三个二维高斯分布生成的数据集。

在高斯混合模型中，我们需要估计每一个高斯分布的均值与方差。从最大似然估计的角度来说，给定某个有n个样本的数据集X，假如已知GMM 中一共有k?簇，我们就是要找到k组均值μ1,?,μk，k组方差σ1,?,σk?来最大化以下似然函数L

这里直接计算似然函数比较困难，于是我们引入隐变量(latent variable)，这里的隐变量就是每个样本属于每一簇的概率。假设W是一个n×k的矩阵，其中?Wi,j?是第?i?个样本属于第?j 簇的概率。

在已知W的情况下，我们就很容易计算似然函数LW，

将其写成

其中P(Xi?|?μj,?σj)?是样本Xi在第j个高斯分布中的概率密度函数。

以一维高斯分布为例，

2. 最大期望算法(Expectation–Maximization, EM)

有了隐变量还不够，我们还需要一个算法来找到最佳的W，从而得到GMM的模型参数。EM算法就是这样一个算法。

简单说来，EM算法分两个步骤。

第一个步骤是E（期望），用来更新隐变量WW；

第二个步骤是M（最大化），用来更新GMM中各高斯分布的参量μj,?σj。

然后重复进行以上两个步骤，直到达到迭代终止条件。

3. 具体步骤以及Python实现

完整代码在第4节。

首先，我们先引用一些我们需要用到的库和函数。

1 import numpy as np

2 import matplotlib.pyplot as plt

3 from matplotlib.patches import Ellipse

4 from scipy.stats import multivariate_normal

5 e('seaborn')?

接下来，我们生成2000条二维模拟数据，其中400个样本来自N(μ1,var1)，600个来自N(μ2,var2)，1000个样本来自N(μ3,var3)

1 # 第一簇的数据

2 num1, mu1, var1 = 400, [0.5, 0.5], [1, 3]

3 X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1)

4 # 第二簇的数据

5 num2, mu2, var2 = 600, [5.5, 2.5], [2, 2]

6 X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2)

7 # 第三簇的数据

8 num3, mu3, var3 = 1000, [1, 7], [6, 2]

9 X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3)

10 # 合并在一起

11 X = np.vstack((X1, X2, X3))

数据如下图所示：

1 plt.figure(figsize=(10, 8))

2 plt.axis([-10, 15, -5, 15])

3 plt.scatter(X1[:, 0], X1[:, 1], s=5)

4 plt.scatter(X2[:, 0], X2[:, 1], s=5)

5 plt.scatter(X3[:, 0], X3[:, 1], s=5)

6 plt.show()

3.1 变量初始化

首先要对GMM模型参数以及隐变量进行初始化。通常可以用一些固定的值或者随机值。

n_clusters?是GMM模型中聚类的个数，和K-Means一样我们需要提前确定。这里通过观察可以看出是3。（拓展阅读：如何确定GMM中聚类的个数？）

n_points?是样本点的个数。

Mu?是每个高斯分布的均值。

Var?是每个高斯分布的方差，为了过程简便，我们这里假设协方差矩阵都是对角阵。

W?是上面提到的隐变量，也就是每个样本属于每一簇的概率，在初始时，我们可以认为每个样本属于某一簇的概率都是1-3。

Pi?是每一簇的比重，可以根据W求得，在初始时，Pi?=?[1-3,?1-3,?1-3]

1 n_clusters = 3

2 n_points = len(X)

3 Mu = [[0, -1], [6, 0], [0, 9]]

4 Var = [[1, 1], [1, 1], [1, 1]]

5 Pi = [1 - n_clusters] * 3

6 W = np.ones((n_points, n_clusters)) - n_clusters

7 Pi = W.sum(axis=0) - W.sum()

3.2 E步骤

E步骤中，我们的主要目的是更新W。第i个变量属于第m簇的概率为：根据W，我们就可以更新每一簇的占比πm，

1 def update_W(X, Mu, Var, Pi):

2 n_points, n_clusters = len(X), len(Pi)

3 pdfs = np.zeros(((n_points, n_clusters)))

4 for i in range(n_clusters):

5 pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i],