GMM与EM算法的Python实现

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

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

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

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

在高斯混合模型中，我们需要估计每一个高斯分布的均值与方差。从最大似然估计的角度来说，给定某个有 $n$

这里直接计算似然函数比较困难，于是我们引入隐变量(latent variable)，这里的隐变量就是每个样本属于每一簇的概率。假设 $W$

在已知 $W$

将其写成

其中 $P (X_{i} | μ_{j}, σ_{j})$

以一维高斯分布为例，

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

有了隐变量还不够，我们还需要一个算法来找到最佳的 $W$

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

第一个步骤是E（期望），用来更新隐变量 $W$
第二个步骤是M（最大化），用来更新GMM中各高斯分布的参量 $μ_{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 plt.style.use(‘seaborn‘)

接下来，我们生成2000条二维模拟数据，其中400个样本来自 $N (μ_{1}, {var}_{1})$

 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,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[:,X2[:,s=5) 5 plt.scatter(X3[:,X3[:,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],[0,9]] 4 Var = [[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$

根据 $W$

 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],np.diag(Var[i]))  6     W = pdfs / pdfs.sum(axis=1).reshape(-1,1)  7     return W  8 
 9 
10 def update_Pi(W): 11     Pi = W.sum(axis=0) / W.sum() 12     return Pi

以下是计算对数似然函数的logLH以及用来可视化数据的plot_clusters。

 1 def logLH(X,Pi,Var):  2     n_points,np.diag(Var[i]))  6     return np.mean(np.log(pdfs.sum(axis=1)))  7 
 8 
 9 def plot_clusters(X,Mu_true=None,Var_true=None): 10     colors = [‘b‘,‘g‘,‘r‘] 11     n_clusters = len(Mu) 12     plt.figure(figsize=(10,8)) 13     plt.axis([-10,15]) 14     plt.scatter(X[:,X[:,s=5) 15     ax = plt.gca() 16     for i in range(n_clusters): 17         plot_args = {‘fc‘: ‘None‘,‘lw‘: 2,‘edgecolor‘: colors[i],‘ls‘: ‘:‘} 18         ellipse = Ellipse(Mu[i],3 * Var[i][0],3 * Var[i][1],**plot_args) 19  ax.add_patch(ellipse) 20     if (Mu_true is not None) & (Var_true is not None): 21         for i in range(n_clusters): 22             plot_args = {‘fc‘: ‘None‘,‘alpha‘: 0.5} 23             ellipse = Ellipse(Mu_true[i],3 * Var_true[i][0],3 * Var_true[i][1],**plot_args) 24  ax.add_patch(ellipse) 25     plt.show()

3.2 M步骤

M步骤中，我们需要根据上面一步得到的W来更新均值Mu和方差Var。 Mu和Var是以W的权重的样本X的均值和方差。

因为这里的数据是二维的，第 $m$

第 $m$

以上迭代公式写成如下函数update_Mu和update_Var。

 1 def update_Mu(X,W):  2     n_clusters = W.shape[1]  3     Mu = np.zeros((n_clusters,2))  4     for i in range(n_clusters):  5         Mu[i] = np.average(X,axis=0,weights=W[:,i])  6     return Mu  7 
 8 def update_Var(X,W):  9     n_clusters = W.shape[1] 10     Var = np.zeros((n_clusters,2)) 11     for i in range(n_clusters): 12         Var[i] = np.average((X - Mu[i]) ** 2,i]) 13     return Var

3.3 迭代求解

下面我们进行迭代求解。

图中实现是真实的高斯分布，虚线是我们估计出的高斯分布。可以看出，经过5次迭代之后，两者几乎完全重合。

1 loglh = [] 2 for i in range(5): 3  plot_clusters(X,[mu1,mu3],[var1,var2,var3]) 4  loglh.append(logLH(X,Var)) 5     W = update_W(X,Pi) 6     Pi = update_Pi(W) 7     Mu = update_Mu(X,W) 8     print(‘log-likehood:%.3f‘%loglh[-1]) 9     Var = update_Var(X,W)

1 log-likehood:-8.054

1 log-likehood:-4.731

1 log-likehood:-4.729

1 log-likehood:-4.728

1 log-likehood:-4.728

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 plt.style.use(‘seaborn‘)  6 
 7 # 生成数据
 8 def generate_X(true_Mu,true_Var):  9     # 第一簇的数据
10     num1,var1 = 400,true_Mu[0],true_Var[0] 11     X1 = np.random.multivariate_normal(mu1,num1) 12     # 第二簇的数据
13     num2,true_Mu[1],true_Var[1] 14     X2 = np.random.multivariate_normal(mu2,num2) 15     # 第三簇的数据
16     num3,true_Mu[2],true_Var[2] 17     X3 = np.random.multivariate_normal(mu3,num3) 18     # 合并在一起
19     X = np.vstack((X1,X3)) 20     # 显示数据
21     plt.figure(figsize=(10,8)) 22     plt.axis([-10,15]) 23     plt.scatter(X1[:,s=5) 24     plt.scatter(X2[:,s=5) 25     plt.scatter(X3[:,s=5) 26  plt.show() 27     return X 28 
29 
30 # 更新W
31 def update_W(X,Pi): 32     n_points,len(Pi) 33     pdfs = np.zeros(((n_points,n_clusters))) 34     for i in range(n_clusters): 35         pdfs[:,np.diag(Var[i])) 36     W = pdfs / pdfs.sum(axis=1).reshape(-1,1) 37     return W 38 
39 
40 # 更新pi
41 def update_Pi(W): 42     Pi = W.sum(axis=0) / W.sum() 43     return Pi 44 
45 
46 # 计算log似然函数
47 def logLH(X,Var): 48     n_points,len(Pi) 49     pdfs = np.zeros(((n_points,n_clusters))) 50     for i in range(n_clusters): 51         pdfs[:,np.diag(Var[i])) 52     return np.mean(np.log(pdfs.sum(axis=1))) 53 
54 
55 # 画出聚类图像
56 def plot_clusters(X,Var_true=None): 57     colors = [‘b‘,‘r‘] 58     n_clusters = len(Mu) 59     plt.figure(figsize=(10,8)) 60     plt.axis([-10,15]) 61     plt.scatter(X[:,s=5) 62     ax = plt.gca()for i in range(n_clusters): 63         plot_args ={‘fc‘:‘None‘,‘lw‘:2,‘ls‘:‘:‘} 64         ellipse =Ellipse(Mu[i],3*Var[i][0],3*Var[i][1],**plot_args) 65         ax.add_patch(ellipse)if(Mu_trueisnotNone)&(Var_trueisnotNone):for i in range(n_clusters): 66             plot_args ={‘fc‘:‘None‘,‘alpha‘:0.5} 67             ellipse =Ellipse(Mu_true[i],3*Var_true[i][0],3*Var_true[i][1],**plot_args) 68  ax.add_patch(ellipse) 69     plt.show()# 更新Mudef update_Mu(X,W):
70     n_clusters = W.shape[1]Mu= np.zeros((n_clusters,2))for i in range(n_clusters):Mu[i]= np.average(X,weights=W[:,i])returnMu# 更新Vardef update_Var(X,W):
71     n_clusters = W.shape[1]Var= np.zeros((n_clusters,2))for i in range(n_clusters):Var[i]= np.average((X -Mu[i])**2,i])returnVarif __name__ ==‘__main__‘:# 生成数据
72     true_Mu =[[0.5,7]] 73     true_Var =[[1,3],2],2]] 74     X = generate_X(true_Mu,true_Var)# 初始化
75     n_clusters =3
76     n_points = len(X)Mu=[[0,9]]Var=[[1,1]]Pi=[1/ n_clusters]*3
77     W = np.ones((n_points,n_clusters))/ n_clusters 78     Pi= W.sum(axis=0)/ W.sum()# 迭代
79     loglh =[]for i in range(5): 80  plot_clusters(X,true_Mu,true_Var) 81  loglh.append(logLH(X,Var)) 82         W = update_W(X,Pi)Pi= update_Pi(W)Mu= update_Mu(X,W)print(‘log-likehood:%.3f‘%loglh[-1])Var= update_Var(X,W)

本教程基于Python 3.6。

原创者：u_u | 修改校对：SofaSofa TeamM | 转自： http://sofasofa.io/tutorials/gmm_em/

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