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

python

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


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

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

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

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

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

在已知,

将其写成

其中个高斯分布中的概率密度函数。

以一维高斯分布为例,

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

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

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

  • 第一个步骤是E(期望),用来更新隐变量W;
  • 第二个步骤是M(最大化),用来更新GMM中各高斯分布的参量。

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

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个样本来自

 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 是每个高斯分布的方差,为了过程简便,我们这里假设协方差矩阵都是对角阵。

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

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。第簇的概率为:

根据,

 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。 MuVar是以W的权重的样本X的均值和方差。

因为这里的数据是二维的,第个分量的均值,

第个分量的方差,

以上迭代公式写成如下函数update_Muupdate_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/

以上是 高斯混合模型GMM与EM算法的Python实现 的全部内容, 来源链接: utcz.com/z/388154.html

回到顶部