一般讲到GMM就会讲到EM。
我不过多的介绍EM算法。这里只是举一些例子来看看真实的GMM怎么用EM算的。

一、GMM的作用

记住GMM的作用，就是聚类！

---

二、GMM有hard和soft两种

hard-GMM和soft-GMM是为了对标k-means和soft k-means。在中文互联网上搜索到的GMM，其实基本都是soft-GMM。
所以这里我先介绍一下hard-GMM。

---

三、hard-GMM

1）已知初始化簇， B 1 = { x 1 , x 2 } , B 2 = { x 3 , x 4 , x 5 , x 6 } B1 = \{x_1, x_2\}, B2 = \{x_3,x_4,x_5,x_6\} B1={x1​,x2​},B2={x3​,x4​,x5​,x6​}.

初始化两个簇权重分别为$w_1=2/6=1/3,w_2=4/6=2/3$.

由此可知，两个簇的高斯分布参数为：
${\mu }_1=(0+0.5)/2=0.25$
${\mu }_2=(1+2+3+4)/4=2.5$.
${\sigma }_1^2=0.25^2$,
${\sigma }_2^2=1.25$

概率密度
$\hat{P}(x)$
$=w_{1}N(x;{\mu }_1,{\sigma }_1^2)+w_{2}N(x;{\mu }_2,{\sigma }_2^2)$
$=1/3N(x;0.25,0.25^2)+2/3N(x;2.5,1.25)$

查表知：

$\hat{P}(x_1)=0.342$	$\hat{P}(x_2)=0.371$
$\hat{P}(x3)=0.103$	$\hat{P}(x4)=0.215$
$\hat{P}(x5)=0.215$	$\hat{P}(x6)=0.097$

$E_{in}=-\sum _{n=1}^{6}log\hat{P}(x_n)=9.748$

2）使用EM算法迭代直至收敛，需要写出每一步的聚类结果和模型参数

E-step：

计算归属度${\lambda }_{ij}$，表示数据点$i$对簇$j$的归属度。
${\lambda }_{ij}=w_{j}N(x_i;{\mu }_j,{\sigma }_j^2)$

（1）对数据点$x_1$来说，对簇1簇2的归属度分别为：
${\lambda }_{11}$
$=w_{1}N(x_1;{\mu }_1,{\sigma }_1^2)$
$=1/3N(0;0.25,0.25^2)$
$=0.323$
${\lambda }_{12}$
$=w_{2}N(x_1;{\mu }_2,{\sigma }_2^2)$
$=2/3N(0;2.5,1.25)$
$=0.019$

${\lambda }_{11}>{\lambda }_{12}$
说明该步下，数据点1归属到簇1

（2）对数据点$x_3$来说，对簇1簇2的归属度分别为：
${\lambda }_{31}$
$=w_{1}N(x_3;{\mu }_1,{\sigma }_1^2)$
$=1/3N(1;0.25,0.25^2)$
$=4/(3\sqrt{2\pi })\ast e^{-4.5)}$
$=0.006$
${\lambda }_{32}$
$=w_{2}N(x_3;{\mu }_2,{\sigma }_2^2)$
$=2/3N(1;2.5,1.25)$
$=(2/\sqrt{1.25})/(3\sqrt{2\pi })\ast e^{-0.9}$
$=0.097$
${\lambda }_{31}<{\lambda }_{32}$
说明该步下，数据点3归属到簇2

（3）同理$x_2$归属到簇1，$x_4,x_5,x_6$归属到簇2
所以最终我们还是得到 B 1 = { x 1 , x 2 } , B 2 = { x 3 , x 4 , x 5 , x 6 } B1 = \{x1, x2\}, B2 = \{x3,x4,x5,x6\} B1={x1,x2},B2={x3,x4,x5,x6}. 也就是说已经收敛，解答完毕。

但如果还没有收敛怎么办？
假设得到了 B 1 = { x 1 , x 2 , x 3 } , B 2 = { x 4 , x 5 , x 6 } B1 = \{x_1, x_2, x_3\}, B2 = \{x_4,x_5,x_6\} B1={x1​,x2​,x3​},B2={x4​,x5​,x6​}
那就再通过M步更新高斯模型参数

M-step

两个簇权重分别为$w_1=2/6=1/3,w_2=4/6=2/3$
两个簇的高斯分布为
${\mu }_1=(0+0.5+1)/3=0.5$
${\mu }_2=(2+3+4)/3=3$
${\sigma }_1^2=1/6$
${\sigma }_2^2=2/3$

概率密度
$\hat{P}(x)$
$=w_{1}N(x;{\mu }_1,{\sigma }_1^2)+w_{2}N(x;{\mu }_2,{\sigma }_2^2)$
$=1/3N(x;0.5,1/6)+2/3N(x;3,2/3)$

这样就可以计算新一轮的E-step，如此反复，直至数据点的归属簇不再变化

---

类比一下：hard-GMM中，E-step就像是k-means中的将各数据点归属到各个簇，M-step就像是k-means计算新的簇中心

---

四、soft-GMM

第1）问没有区别，第2）问开始重新写soft-GMM版本的

2）使用EM算法迭代直至收敛，需要写出每一步的聚类结果和模型参数

E-step：

计算归属度${\lambda }_{ij}$，表示数据点$i$对簇$j$的归属度。
${\lambda }_{ij}=w_{j}N(x_i;{\mu }_j,{\sigma }_j^2)$

（1）对数据点$x_1$来说，对簇1簇2的归属度分别为：
${\lambda }_{11}$
$=w_{1}N(x_1;{\mu }_1,{\sigma }_1^2)$
$=1/3N(0;0.25,0.25^2)$
$=0.323$
${\lambda }_{12}$
$=w_{2}N(x_1;{\mu }_2,{\sigma }_2^2)$
$=2/3N(0;2.5,1.25)$
$=0.019$

与hard-GMM不同的时，我们需要计算归一化后的归属度：
${\lambda }_{11}^{\prime }=\frac{{\lambda }_{11}}{{\lambda }_{11}+{\lambda }_{12}}=0.943$
${\lambda }_{12}^{\prime }=\frac{{\lambda }_{12}}{{\lambda }_{11}+{\lambda }_{12}}=0.057$
${\lambda }_{11}={\lambda }_{11}^{\prime }$
${\lambda }_{12}={\lambda }_{12}^{\prime }$

（2）对数据点$x_3$来说，对簇1簇2的归属度分别为：
${\lambda }_{31}$
$=w_{1}N(x_3;{\mu }_1,{\sigma }_1^2)$
$=1/3N(1;0.25,0.25^2)$
$=4/(3\sqrt{2\pi })\ast e^{-4.5)}$
$=0.006$
${\lambda }_{32}$
$=w_{2}N(x_3;{\mu }_2,{\sigma }_2^2)$
$=2/3N(1;2.5,1.25)$
$=(2/\sqrt{1.25})/(3\sqrt{2\pi })\ast e^{-0.9}$
$=0.097$

计算归一化后的归属度：
${\lambda }_{31}^{\prime }=\frac{{\lambda }_{31}}{{\lambda }_{31}+{\lambda }_{32}}=0.058$
${\lambda }_{32}^{\prime }=\frac{{\lambda }_{32}}{{\lambda }_{31}+{\lambda }_{32}}=0.942$
${\lambda }_{31}={\lambda }_{31}^{\prime }$
${\lambda }_{32}={\lambda }_{32}^{\prime }$

（3）同理计算剩下的数据点的归属度
soft-GMM如何认为算法已经收敛了呢？
当两次归属度的值变化量 < tolerance即可，比如我们这里设定 tolerance=0.02

我们这里再通过M步更新高斯模型参数

M-step

两个簇权重分别为$w_1=\sum _{i=1}^6{\lambda }_{i1}=0.312,w_2=\sum _{i=1}^6{\lambda }_{i2}=0.688$
两个簇的高斯分布为
${\mu }_1=\frac{\sum _{i=1}^6{\lambda }_{i1}{x}_i}{\sum _{i=1}^6{\lambda }_{i1}}=\frac{({\lambda }_{11}{x}_1+{\lambda }_{21}{x}_2+...+{\lambda }_{61}{x}_6)}{{\lambda }_{11}+{\lambda }_{21}+...+{\lambda }_{61}}=0.263$

${\mu }_2=\frac{\sum _{i=1}^6{\lambda }_{i2}{x}_i}{\sum _{i=1}^6{\lambda }_{i2}}=2.424$

${\sigma }_1^2=\frac{\sum _{i=1}^6{\lambda }_{i1}(x_i-{\mu }_1{)}^2}{\sum _{i=1}^6{\lambda }_{i1}}=0.279$

${\sigma }_2^2=\frac{\sum _{i=1}^6{\lambda }_{i2}(x_i-{\mu }_2{)}^2}{\sum _{i=1}^6{\lambda }_{i2}}=1.180$

概率密度
$\hat{P}(x)$
$=w_{1}N(x;{\mu }_1,{\sigma }_1^2)+w_{2}N(x;{\mu }_2,{\sigma }_2^2)$
$=0.312N(x;0.263,0.279)+0.688N(x;2.424,1.180)$

这样就可以计算新一轮的E-step，如此反复，直至收敛（参数是否收敛或对数似然函数收敛）
参数收敛比方说第19轮和第20轮迭代时：

iter	${\mu }_1$	${\mu }_2$	${\sigma }_1^2$	${\sigma }_1^2$
19	0.485	2.923	0.408	0.896
20	0.485	2.923	0.408	0.895

此时指定的参数已经的变化量已经 < tolerance = 0.02，可以认为已经收敛。

---

类比一下：soft-GMM中，E-step就像是k-means中的将各数据点对各个簇的归属度，M-step就像是k-means计算新的簇中心

---

五、python复现代码

"""
import numpy as np
import timedef CreateData(mu1, sigma1, alpha1, mu2, sigma2, alpha2, length):"""生成数据集输入两个高斯分布的真实参数:param mu1::param sigma1::param alpha1::param mu2::param sigma2::param alpha2::param length::return:        通过两个高斯分布生成的数据集"""# 生成第一个分模型的数据First = np.random.normal(loc=mu1, scale=sigma1, size=int(length * alpha1))# 生成第二个分模型的数据Second = np.random.normal(loc=mu2, scale=sigma2, size=int(length * alpha2))# 混合两个数据ObservedData = np.concatenate((First, Second), axis=0)  # 行拼接，生成的还是一列# 打乱顺序(打不打乱其实对算法没有影响)np.random.shuffle(ObservedData)return ObservedDatadef Cal_Gaussian(Data, mu, sigma):"""根据高斯密度函数计算:param Data:数据集:param mu:  当前的高斯分布的参数:param sigma::return:    GaussianP：概率"""GaussianP = 1 / (sigma * np.sqrt(2 * np.pi)) \* np.exp(-1 * np.square(Data - mu) / (2 * np.square(sigma)))return GaussianPdef E_step(Data, mu1, sigma1, alpha1, mu2, sigma2, alpha2):"""E步:param Data:    使用的是ndarray格式:param mu1::param sigma1::param alpha1::param mu2::param sigma2::param alpha2::return:        Gamma1； Gamma2：响应度"""# 计算样本对每个高斯混合模型的响应度（归属度，gamma）# 因此得到的数列结果中，包含该分模型对每个样本的计算响应度公式的分子项print(f" Cal_Gaussian(Data, mu1, sigma1) = { Cal_Gaussian(Data, mu1, sigma1)}")Gamma1 = alpha1 * Cal_Gaussian(Data, mu1, sigma1)  # 是一个列表，保存了每个样本对高斯模型1的响应度print(f" Cal_Gaussian(Data, mu2, sigma2) = { Cal_Gaussian(Data, mu2, sigma2)}")Gamma2 = alpha2 * Cal_Gaussian(Data, mu2, sigma2)# 响应度归一化Gamma1 = Gamma1 / (Gamma1 + Gamma2)# Gamma2 = Gamma2 / (Gamma1 + Gamma2)   # 注意原来是这样写的，但是Gamma1已经在上一行变化了，不能再用来归一化Gamma2 = 1 - Gamma1return Gamma1, Gamma2def M_step(Data, mu1_old, mu2_old, Gamma1, Gamma2):"""M步:param Data::param mu1_old: 当前高斯分布的参数:param mu2_old::param Gamma1:  对高斯分布1的响应度:param Gamma2:  对高斯分布2的响应度:return:    新的高斯分布的参数"""# 计算新的参数mumu1_new = np.dot(Gamma1, Data) / np.sum(Gamma1)mu2_new = np.dot(Gamma2, Data) / np.sum(Gamma2)# 计算新的参数sigma# sigma1_new = np.sqrt(np.dot(Gamma1, np.square(Data - mu1_old)) / np.sum(Gamma1))sigma1_new = np.sqrt(np.dot(Gamma1, np.square(Data - mu1_new)) / np.sum(Gamma1))# sigma2_new = np.sqrt(np.dot(Gamma2, np.square(Data - mu2_old)) / np.sum(Gamma2))sigma2_new = np.sqrt(np.dot(Gamma2, np.square(Data - mu2_new)) / np.sum(Gamma2))# 计算新的参数alpham = len(Data)alpha1_new = np.sum(Gamma1) / malpha2_new = np.sum(Gamma2) / mreturn mu1_new, sigma1_new, alpha1_new, mu2_new, sigma2_new, alpha2_newdef EM(Data, mu1, sigma1, alpha1, mu2, sigma2, alpha2, itertime):"""EM算法:param Data::param mu1::param sigma1::param alpha1::param mu2::param sigma2::param alpha2::param itertime:    迭代次数:return:            模型中分模型的参数"""# 迭代for i in range(itertime):print("************************************************")print(f"\n\n第{i}次迭代")# 计算响应度，E步Gamma1, Gamma2 = E_step(Data, mu1, sigma1, alpha1, mu2, sigma2, alpha2)print(f"对1类的隶属度 = {Gamma1}, \n对2类的隶属度 = {Gamma2}")# 更新参数，M步mu1, sigma1, alpha1, mu2, sigma2, alpha2 \= M_step(Data, mu1, mu2, Gamma1, Gamma2)print(f"mu1={mu1}, sigma1={sigma1}, alpha1={alpha1},\n" f"mu2={mu2}, sigma2={sigma2}, alpha2={alpha2}")return mu1, sigma1, alpha1, mu2, sigma2, alpha2if __name__ == '__main__':# 生成数据Data = np.array([0,0.5,1,2,3,4])# 初始化数据init_mu1 = 0.25init_sigma1 = 0.25init_theta1, init_alpha1 = [init_mu1, init_sigma1], 1/3init_mu2 = 2.5init_sigma2 = np.sqrt(1.25)init_theta2, init_alpha2 = [init_mu2, init_sigma2], 2/3itertime = 20tol = 2e-3# 训练模型start = time.time()print("start training")model_mu1, model_sigma1, model_alpha1, model_mu2, model_sigma2, model_alpha2 = \EM(Data, init_mu1, init_sigma1, init_alpha1, init_mu2, init_sigma2, init_alpha2, itertime)print('end training')end = time.time()print('training time: ', end - start)# print('真实参数：mu1 ', mu1, 'sigma1 ', sigma1, 'alpha1 ', alpha1, 'mu2 ', mu2, 'sigma2 ', sigma2, 'alpha2 ', alpha2)print('模型参数：mu1 ', model_mu1, 'sigma1 ', model_sigma1, 'alpha1 ', model_alpha1, 'mu2 ', model_mu2, 'sigma2 ',model_sigma2, 'alpha2 ', model_alpha2)

---

参考：

[1]《Dempster A P . Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society, 1977, 39.》
[2]《高斯混合模型(GMM)》
[3]《EM算法python复现 - 高斯混合模型》