三维点云学习（3）7- 实现GMM

github大神参考代码
高斯混合模型的通俗理解
GMM课程个人总结笔记

最终效果图

原图

进行高斯聚类后的图

代码编写流程

1.输入数据集x1 x2 …xn,和K，初始化三个高斯模型的未知数
2.E-step 算出后验概率 --一个点属于哪个类的概率
3.M-step 算出高斯模型三个参数

核心编码

K-means为上一章写得KMeans.py代码,放置到统一目录下即可调用

#GMM.py
class GMM(object):# max_iter  最大迭代次数为50次def __init__(self, n_clusters, max_iter=1):self.n_clusters = n_clustersself.max_iter = max_iterself.k = n_clustersself.mu = None            #meanself.cov = None           #协方差矩阵self.prior = None         #先验概率 -> 权重self.posteriori = None   #后验概率       # k*N矩阵

聚类流程
1.初始化
类的个数：k=3；
协方差矩阵 cov为三个 2*2的单位矩阵；
后验概率 pi 为 1/k （全概率/类的个数）；
mean 为K-Means选取的center点

2.E-step 算出后验概率 --一个点属于哪个类的概率

3.M-step 算出高斯模型三个参数

4.step2 step3 往复迭代

    def fit(self, data):# 作业3# 屏蔽开始# step1 初始化 Mu pi cov##init Mu ，使用K-means中心点k_means = KMeans.K_Means(n_clusters=self.k)k_means.fit(data)self.mu = np.asarray(k_means.centers_)   # 将mean的初始值为 k-means 的中心点   3*2 矩阵self.cov = np.asarray([eye(2,2)]*self.k)                      #初始化的cov为 3*2*2 的单位矩阵self.prior = np.asarray([1/self.k ]*self.k).reshape(1,3)                          #对pi进行均等分  1*3 矩阵self.posteriori = np.zeros((self.k,len(data)))                                    #初始化 后验概率为 K*N的 矩阵for _ in range(self.max_iter):   #迭代#step2 E-step 算出后验概率posteriori --一个点属于哪个类的概率for k in range(self.k):self.posteriori[k] = multivariate_normal.pdf(x=data,mean=self.mu[k],cov=self.cov[k])    #提取每个点的概率密度分布self.posteriori = np.dot(diag(self.prior.ravel()),self.posteriori)                                  #diag 将一维数组元素放在对角线上，方便进行对应的数据乘法运算 变为3*3对角矩阵 3*3 * 3*N = 3*N,ravel 将矩阵里所有元素变为列表##归一化self.posteriori /= np.sum(self.posteriori,axis=0)   #后验概率,3*N矩阵#step3 M-step 使用MLE 算出高斯模型三个参数  mu:mean cov：协方差  prior：先验概率self.Nk = np.sum(self.posteriori,axis=1)self.mu = np.asarray([np.dot(self.posteriori[k],data) / self.Nk[k] for k in range(self.k)])        #self.posteriori[k]: 3*2  data:n*2  self.Nk[k]:1self.cov = np.asarray([np.dot((data-self.mu[k]).T, np.dot(np.diag(self.posteriori[k].ravel()),data-self.mu[k])) / self.Nk[k] for k in range(self.k)])   #sel.cov : 3*2*2self.prior = np.asarray([self.Nk / self.k])               #self.prior  1*3self.fitted = True

个人心得

1.用好numpy进行矩阵运算，可以大大节省算法运算时间和提高效率
2.初始值的选取，看过大神的代码，初始值。mean选取的是K-means后的中心点，pi 初始化时 K-means后每个类在整体中所占据的比例，方差 K-means类中点的协方差矩阵，后验概率是经过前三值得计算后结果，初始中心点得选取尽可能之间分布较远，有利于聚类更快得收敛
3.理清楚矩阵之间得运算，多print shape
4.用好numpy！！！

完整代码

python实现K-means

# 文件功能：实现 GMM 算法import numpy as np
from numpy import *
import pylab
import random, mathimport matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normal
import KMeansplt.style.use('seaborn')class GMM(object):# max_iter  最大迭代次数为50次def __init__(self, n_clusters, max_iter=1):self.n_clusters = n_clustersself.max_iter = max_iterself.k = n_clustersself.mu = None            #meanself.cov = None           #协方差矩阵self.prior = None         #先验概率 -> 权重self.posteriori = None   #后验概率       # k*N矩阵def fit(self, data):# 作业3# 屏蔽开始# step1 初始化 Mu pi cov##init Mu ，使用K-means中心点k_means = KMeans.K_Means(n_clusters=self.k)k_means.fit(data)self.mu = np.asarray(k_means.centers_)   # 将mean的初始值为 k-means 的中心点   3*2 矩阵self.cov = np.asarray([eye(2,2)]*self.k)                      #初始化的cov为 3*2*2 的单位矩阵self.prior = np.asarray([1/self.k ]*self.k).reshape(1,3)                          #对pi进行均等分  1*3 矩阵self.posteriori = np.zeros((self.k,len(data)))                                    #初始化 后验概率为 K*N的 矩阵for _ in range(self.max_iter):   #迭代#step2 E-step 算出后验概率posteriori --一个点属于哪个类的概率for k in range(self.k):self.posteriori[k] = multivariate_normal.pdf(x=data,mean=self.mu[k],cov=self.cov[k])    #提取每个点的概率密度分布self.posteriori = np.dot(diag(self.prior.ravel()),self.posteriori)                                  #diag 将一维数组元素放在对角线上，方便进行对应的数据乘法运算 变为3*3对角矩阵 3*3 * 3*N = 3*N,ravel 将矩阵里所有元素变为列表##归一化self.posteriori /= np.sum(self.posteriori,axis=0)   #后验概率,3*N矩阵#step3 M-step 使用MLE 算出高斯模型三个参数  mu:mean cov：协方差  prior：先验概率self.Nk = np.sum(self.posteriori,axis=1)self.mu = np.asarray([np.dot(self.posteriori[k],data) / self.Nk[k] for k in range(self.k)])        #self.posteriori[k]: 3*2  data:n*2  self.Nk[k]:1self.cov = np.asarray([np.dot((data-self.mu[k]).T, np.dot(np.diag(self.posteriori[k].ravel()),data-self.mu[k])) / self.Nk[k] for k in range(self.k)])   #sel.cov : 3*2*2self.prior = np.asarray([self.Nk / self.k])               #self.prior  1*3self.fitted = True# 屏蔽结束#计算每个点的类别def predict(self, data):# 屏蔽开始result = []if not self.fitted:print('Unfitter. ')return resultresult = np.argmax(self.posteriori,axis=0)    #比较每个点的后验概率return result# 屏蔽结束# 生成仿真数据
def generate_X(true_Mu, true_Var):# 第一簇的数据num1, mu1, var1 = 400, true_Mu[0], true_Var[0]X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1)# 第二簇的数据num2, mu2, var2 = 600, true_Mu[1], true_Var[1]X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2)# 第三簇的数据num3, mu3, var3 = 1000, true_Mu[2], true_Var[2]X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3)# 合并在一起X = np.vstack((X1, X2, X3))# 显示数据plt.figure(figsize=(10, 8))plt.axis([-10, 15, -5, 15])plt.scatter(X1[:, 0], X1[:, 1], s=5)plt.scatter(X2[:, 0], X2[:, 1], s=5)plt.scatter(X3[:, 0], X3[:, 1], s=5)plt.show()return Xdef points_show(point,color):x = []y = []point = np.asarray(point)for i in range(len(point)):x.append(point[i][0])y.append(point[i][1])plt.scatter(x, y,color=color)#plt.show()if __name__ == '__main__':# 生成数据true_Mu = [[0.5, 0.5], [5.5, 2.5], [1, 7]]true_Var = [[1, 3], [2, 2], [6, 2]]X = generate_X(true_Mu, true_Var)gmm = GMM(n_clusters=3)gmm.fit(X)  # 确定中心点的位置cat = gmm.predict(X)  # 提取点的索引K = 3# visualize:color = ['red', 'blue', 'green', 'cyan', 'magenta']labels = [f'Cluster{k:02d}' for k in range(K)]cluster = [[] for i in range(K)]  # 用于分类所有数据点for i in range(len(X)):if cat[i] == 0:cluster[0].append(X[i])elif cat[i] == 1:cluster[1].append(X[i])elif cat[i] == 2:cluster[2].append(X[i])points_show(cluster[0], color="red")points_show(cluster[1], color="blue")points_show(cluster[2], color="yellow")plt.show()# for k in range(K):#     plt.scatter(X[cat == k][:, 0], X[cat == k][:, 1], c=color[k], label=labels[k])## mu = gmm.mu# plt.scatter(mu[:, 0], mu[:, 1], s=300, c='grey', marker='P', label='Centroids')## plt.xlabel('X')# plt.ylabel('Y')# plt.legend()     #显示图例# plt.title('GMM Testcase')# plt.show()