高斯混合模型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组均值μ₁,⋯,μ_k，k组方差σ₁,⋯,σ_k 来最大化如下似然函数L

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

在已知W的状况下，咱们就很容易计算似然函数L_W，

将其写成

其中P(X_i | μ_j, σ_j) 是样本X_i在第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 plt.style.use('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], 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, Mu, Var):  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 return np.mean(np.log(pdfs.sum(axis=1)))  7  8  9 def plot_clusters(X, Mu, Var, 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, -5, 15]) 14 plt.scatter(X[:, 0], X[:, 1], 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', 'lw': 2, 'edgecolor': colors[i], '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簇的第k个份量的均值，

 

第m簇的第k个份量的方差，

 

 

以上迭代公式写成以下函数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, Mu, 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, axis=0, weights=W[:, i]) 13 return Var

 

3.3 迭代求解

下面咱们进行迭代求解。

图中实现是真实的高斯分布，虚线是咱们估计出的高斯分布。能够看出，通过5次迭代以后，二者几乎彻底重合。

1 loglh = [] 2 for i in range(5): 3  plot_clusters(X, Mu, Var, [mu1, mu2, mu3], [var1, var2, var3]) 4  loglh.append(logLH(X, Pi, Mu, Var)) 5 W = update_W(X, Mu, Var, Pi) 6 Pi = update_Pi(W) 7 Mu = update_Mu(X, W) 8 print('log-likehood:%.3f'%loglh[-1]) 9 Var = update_Var(X, Mu, 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, mu1, var1 = 400, true_Mu[0], true_Var[0] 11 X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1) 12 # 第二簇的数据 13 num2, mu2, var2 = 600, true_Mu[1], true_Var[1] 14 X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2) 15 # 第三簇的数据 16 num3, mu3, var3 = 1000, true_Mu[2], true_Var[2] 17 X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3) 18 # 合并在一块儿 19 X = np.vstack((X1, X2, X3)) 20 # 显示数据 21 plt.figure(figsize=(10, 8)) 22 plt.axis([-10, 15, -5, 15]) 23 plt.scatter(X1[:, 0], X1[:, 1], s=5) 24 plt.scatter(X2[:, 0], X2[:, 1], s=5) 25 plt.scatter(X3[:, 0], X3[:, 1], s=5) 26  plt.show() 27 return X 28 29 30 # 更新W 31 def update_W(X, Mu, Var, Pi): 32 n_points, n_clusters = len(X), len(Pi) 33 pdfs = np.zeros(((n_points, n_clusters))) 34 for i in range(n_clusters): 35 pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i], 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, Pi, Mu, Var): 48 n_points, n_clusters = len(X), len(Pi) 49 pdfs = np.zeros(((n_points, n_clusters))) 50 for i in range(n_clusters): 51 pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i], np.diag(Var[i])) 52 return np.mean(np.log(pdfs.sum(axis=1))) 53 54 55 # 画出聚类图像 56 def plot_clusters(X, Mu, Var, Mu_true=None, Var_true=None): 57 colors = ['b','g','r'] 58 n_clusters = len(Mu) 59 plt.figure(figsize=(10,8)) 60 plt.axis([-10,15,-5,15]) 61 plt.scatter(X[:,0], X[:,1], s=5) 62 ax = plt.gca()for i in range(n_clusters): 63 plot_args ={'fc':'None','lw':2,'edgecolor': colors[i],'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','lw':2,'edgecolor': colors[i],'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, axis=0, weights=W[:, i])returnMu# 更新Vardef update_Var(X,Mu, 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, axis=0, weights=W[:, i])returnVarif __name__ =='__main__':# 生成数据 72 true_Mu =[[0.5,0.5],[5.5,2.5],[1,7]] 73 true_Var =[[1,3],[2,2],[6,2]] 74 X = generate_X(true_Mu, true_Var)# 初始化 75 n_clusters =3 76 n_points = len(X)Mu=[[0,-1],[6,0],[0,9]]Var=[[1,1],[1,1],[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,Mu,Var, true_Mu, true_Var) 81  loglh.append(logLH(X,Pi,Mu,Var)) 82 W = update_W(X,Mu,Var,Pi)Pi= update_Pi(W)Mu= update_Mu(X, W)print('log-likehood:%.3f'%loglh[-1])Var= update_Var(X,Mu, W)

 

本教程基于Python 3.6。

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