三维点云学习(3)7- 实现GMM
程序员文章站
2022-05-18 19:25:40
三维点云学习(3)7- 实现GMMgithub大神参考代码高斯混合模型的通俗理解GMM课程个人总结笔记最终效果图原图进行高斯聚类后的图代码编写流程1.输入数据集x1 x2 …xn,和K,初始化三个高斯模型的未知数2.E-step 算出后验概率 --一个点属于哪个类的概率3.M-step 算出高斯模型三个参数核心编码K-means为上一章写得KMeans.py代码,放置到统一目录下即可调用#GMM.pyclass GMM(object): # max_iter...
三维点云学习(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_clusters
self.max_iter = max_iter
self.k = n_clusters
self.mu = None #mean
self.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(3,1) #对pi进行均等分 3*1 矩阵
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]:1
self.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*2
self.prior = np.asarray([self.Nk / self.k]).reshape(3,1) #self.prior 3*1
self.fitted = True
个人心得
1.用好numpy进行矩阵运算,可以大大节省算法运算时间和提高效率
2.初始值的选取,看过大神的代码,初始值。mean选取的是K-means后的中心点,pi 初始化时 K-means后每个类在整体中所占据的比例,方差 K-means类中点的协方差矩阵,后验概率是经过前三值得计算后结果,初始中心点得选取尽可能之间分布较远,有利于聚类更快得收敛
3.理清楚矩阵之间得运算,多print shape
4.用好numpy!!!
完整代码
# 文件功能:实现 GMM 算法
import numpy as np
from numpy import *
import pylab
import random, math
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normal
import KMeans
plt.style.use('seaborn')
class GMM(object):
# max_iter 最大迭代次数为50次
def __init__(self, n_clusters, max_iter=1):
self.n_clusters = n_clusters
self.max_iter = max_iter
self.k = n_clusters
self.mu = None #mean
self.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(3,1) #对pi进行均等分 3*1 矩阵
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]:1
self.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*2
self.prior = np.asarray([self.Nk / self.k]).reshape(3,1) #self.prior 3*1
self.fitted = True
# 屏蔽结束
#计算每个点的类别
def predict(self, data):
# 屏蔽开始
result = []
if not self.fitted:
print('Unfitter. ')
return result
result = 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 X
def 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()
本文地址:https://blog.csdn.net/weixin_41281151/article/details/107440893