Mixture-of-Gaussian Clustering详解
第二十次写博客,本人数学基础不是太好,如果有幸能得到读者指正,感激不尽,希望能借此机会向大家学习。这一篇文章是原型聚类中介绍的第四个算法,主要是谈一谈“高斯混合分布”(Mixture-of-Gaussian Distribution)。其他有关于原型聚类算法的讨论可以移步到该类算法的导航页《原型聚类算法综述(原型聚类算法开篇)》。
高斯混合分布(Mixture-of-Gaussian Distribution)与样本生成
同样是基于原型的聚类方法,混合高斯聚类(Mixture-of-Gaussian)不同于K-Means和LVQ中使用向量作为聚类原型,他使用概率模型(高斯混合模型)来表示原型,与一元高斯分布类似,多元高斯分布可以表示为
其中, 是 维均值向量, 是 的协方差矩阵,这两个参数决定了整个高斯分布的概率分布,为了体现出不同模型参数对样本点概率的影响,将该概率密度函数记为 。由此引申出如下高斯混合分布 (假设由 个多元高斯模型混合)的定义式
上式中 和 代表第 个多元高斯混合成分的模型参数, 为该混合成分的“混合系数”(mixture coefficient),且 。高斯混合分布中样本的生成过程,首先根据先验概率 选择对应的高斯混合成分,然后根据该混合成分的模型参数 和 进行采样,进而生成该样本点。
混合模型参数的确定
令 表示样本点 所属的混合成分,那么样本 由第 个多元高斯混合成分生成的概率可以表示为如下后验概率,并且由贝叶斯定理可得
其中, 是第 个混合系数,同时代表了第 个混合成分被选中的先验概率, 简写为 。当高斯混合分布的各模型参数已知后,可以通过下式决定该样本点所属的簇
式中
是第
个样本的簇标记。从上式中可以看出,簇划分是依据原型的后验概率
进行的,因此,问题的关键是如何求得原型参数
。
给定数据集
,原型参数应使得每个样本的高斯混合概率密度之和
最大,因此取对数似然后,可以得到如下优化目标
下面分别将优化目标对各模型参数求偏导,来求得模型参数的定义式
(1) 对
求偏导并置零,得到
(2) 对 求偏导并置零,得到
(3) 由于参数 要满足约束 ,因此在对 求偏导时要引入拉格朗日乘子
将上式对 求导并置零可以得到
注意到上面各式中存在一个未知量 ,且 的计算需要用到上述原型参数,因此与K-Means类似,根据EM算法对 和 进行交替迭代更新,具体的算法描述如下图所示
代码实现及对比
下面是我自己实现的代码,建议不要用load_Data2()中加载的数据集,可以用load_Data1()自己加载一个比较小的数据集,这样可以更快的看出算法效果,下面我还会对sklearn中的实现进行比较。
代码细节
"""
@author: Ἥλιος
@CSDN:https://blog.csdn.net/qq_40793975/article/details/82382583
"""
print(__doc__)
import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
# 加载数据集(从文件中)
def load_Data1(filename):
data_set = []
with open(filename) as fi:
for line in fi.readlines():
cur_line = line.strip().split('\t')
flt_line = []
for i in cur_line:
flt_line.append(float(i))
data_set.append(flt_line)
data_mat = np.mat(data_set) # 转化为矩阵形式
return data_mat
# 加载数据集(自建数据集)
def load_Data2(n_samples=1500):
# 带噪声的圆形数据
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5, noise=.05)
# 带噪声的月牙形数据
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
# 随机分布数据
no_structure = np.random.rand(n_samples, 2), np.ones((1, n_samples), dtype=np.int32).tolist()[0]
# 各向异性分布数据(Anisotropicly distributed data)
random_state = 170
X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
aniso = (X_aniso, y)
# 不同方差的气泡形数据(blobs with varied variances)
varied = datasets.make_blobs(n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state)
# 相同方差的气泡形数据
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)
# 合并数据集
data_sets = [noisy_circles, noisy_moons, no_structure, aniso, varied, blobs]
cluster_nums = [2, 2, 3, 3, 3, 3]
data_mats = []
for i in range(data_sets.__len__()):
X, y = data_sets[i]
X = StandardScaler().fit_transform(X) # 对数据集进行标准化处理
X_mat = np.mat(X)
y_mat = np.mat(y)
data_mats.append((X_mat, y_mat))
# # 展示数据集
# plt.figure(figsize=(2.5, 14))
# plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
# for i in range(data_sets.__len__()):
# X, y = data_sets[i]
# X = StandardScaler().fit_transform(X) # 对数据集进行标准化处理
# plt.subplot(len(data_sets), 1, i+1)
# if i == 0:
# plt.title("Self-built Data Set", size=18)
# color = np.array(y)
# plt.scatter(X[:, 0], X[:, 1], c=color, s=10)
# plt.xlim(-2.5, 2.5)
# plt.ylim(-2.5, 2.5)
#
# plt.show()
return data_mats, cluster_nums
# 高斯模型概率计算(P(x_j|Z_i))
def calculate_GM(x_j, mu_i, Sigma_i):
n = np.shape(x_j)[1]
p = (np.exp((-1/2)*(x_j-mu_i)*Sigma_i.I*(x_j-mu_i).T)[0, 0] /
((np.power(2*np.pi, n/2))*np.sqrt(np.linalg.det(Sigma_i))))
return p
# 高斯混合聚类
def mixture_Gauss(data_mat, k, maxIter=100):
# 变量声明
m, n = np.shape(data_mat)
cluster_assment = np.mat(np.zeros((m, 1)), dtype=np.int32) # 存储每个样本点的簇隶属、簇隶属度和距离所属簇质心的长度
gamma = np.mat(np.random.rand(m, k), dtype=np.float32) # 存储每个样本点的后验概率
gamma = gamma / np.sum(gamma, axis=1) # 初始化后验概率
gamma_temp = np.mat(np.zeros((m, k)), dtype=np.float32)
alpha = np.mat(np.random.rand(k, 1), dtype=np.float32) # 存储混合系数
mu = [] # 存储均值向量
Sigma = [] # 存储协方差矩阵
pts_in_cluster = [] # 存储每个簇中的样本点数
# 初始化
alpha /= np.sum(alpha) # 初始化混合系数
cluster_assment = np.argmax(gamma, axis=1) # 根据后验概率对数据集进行划分
for i in range(k): # 初始化均值向量
mu.append(np.sum(data_mat[np.nonzero(cluster_assment[:, 0] == i)[0], :], axis=0) /
np.sum(cluster_assment[:, 0] == i))
pts_in_cluster.append(np.sum(cluster_assment[:, 0] == i))
for i in range(k): # 初始化协方差矩阵
sigma_temp = np.mat(np.zeros((n, n)), dtype=np.float32)
for j in range(m):
sigma_temp += gamma[j, i]*(data_mat[j, :]-mu[i]).T*(data_mat[j, :]-mu[i])
Sigma.append(sigma_temp / np.sum(gamma[:, i]))
# cluster_changed = True # 后验概率发生变化
# 迭代更新上述参数
for _ in range(maxIter): # 迭代停止条件可以选择最大迭代次数
# while cluster_changed: # 迭代停止条件可以选择后验概率不发生变化
# cluster_changed = False
print(_)
for j in range(m): # E步:更新每个样本点的后验概率
for i in range(k):
gamma_temp[j, i] = alpha[i, :] * calculate_GM(data_mat[j, :], mu[i], Sigma[i])
gamma_temp = gamma_temp / np.sum(gamma_temp, axis=1)
gamma = gamma_temp.copy() # 迭代停止条件(最大迭代次数)
cluster_assment = np.argmax(gamma, axis=1) # 划分数据集(由于动态显示和计算SSE)
SSE_temp = 0
for i in range(k): # 计算SSE
mean = np.sum(data_mat[np.nonzero(cluster_assment[:, 0] == i)[0], :], axis=0) /\
np.sum(cluster_assment == i)
SSE_temp += np.sum(np.power(data_mat[np.nonzero(cluster_assment[:, 0] == i)[0], :] - mean, 2))
print(SSE_temp)
# if (gamma != gamma_temp).any(): # 迭代停止条件(后验概率不发生变化)
# cluster_changed = True
# gamma = gamma_temp.copy()
for i in range(k): # M步:更新每个混合成分
mu_i = np.sum(np.multiply(gamma[:, i], data_mat), axis=0) / np.sum(gamma[:, i])
Sigma_i = np.mat(np.zeros((n, n)), dtype=np.float32)
for j in range(m):
Sigma_i = Sigma_i + gamma[j, i]*(data_mat[j, :]-mu_i).T*(data_mat[j, :]-mu_i)
Sigma_i = Sigma_i / np.sum(gamma[:, i])
alpha_i = np.sum(gamma[:, i]) / m
alpha[i, :] = alpha_i
mu[i] = mu_i
Sigma[i] = Sigma_i
# 动态显示(!不要再Pycharm中运行,会变成幻灯片)
plt.scatter(data_mat[:, 0].T.A[0], data_mat[:, 1].T.A[0], c=cluster_assment[:, 0].T.A[0])
plt.show()
plt.pause(1)
plt.clf()
cluster_assment = np.argmax(gamma, axis=1) # 划分数据集
# # 静态显示
# plt.scatter(data_mat[:, 0].T.A[0], data_mat[:, 1].T.A[0], c=cluster_assment[:, 0].T.A[0])
# plt.show()
return cluster_assment
# GMM
# data_mats, cluster_nums = load_Data2()
# plt.figure(figsize=(2.5, 14))
# plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
# for i in range(len(data_mats)):
# data_mat = data_mats[i][0] # 获取自建数据集
# k = cluster_nums[i] # 获取自建数据集的簇标记
# t0 = time.time() # 计算运行时间
# cluster_assment = mixture_Gauss(data_mat, k) # 自己实现的高斯混合聚类
# # gmm = GaussianMixture(n_components=k, covariance_type='full', init_params="kmeans") # sklearn中自带的高斯混合聚类
# # gmm.fit(data_mat) # sklearn中自带的高斯混合聚类
# # y_pred = gmm.predict(data_mat) # sklearn中自带的高斯混合聚类
# t1 = time.time()
#
# plt.subplot(len(data_mats), 1, i + 1)
# if i == 0:
# plt.title("Mixture of Gaussian(Self-programming Implementation)", size=10)
# # plt.title("Mixture of Gaussian(Sklearn Implementation)", size=10)
# im = plt.scatter(data_mat[:, 0].T.A[0], data_mat[:, 1].T.A[0], c=cluster_assment.T.A[0], s=10, alpha=0.5)
# # im = plt.scatter(data_mat[:, 0].T.A[0], data_mat[:, 1].T.A[0], c=y_pred, s=10, alpha=0.5)
# plt.xlim(-2.5, 2.5)
# plt.ylim(-2.5, 2.5)
# plt.text(.99, .01, ('%.2fs' % (t1 - t0)).lstrip('0'), transform=plt.gca().transAxes, size=15,
# horizontalalignment='right')
# plt.colorbar(im)
# plt.show()
data_mat = load_Data1("C:\\Users\\Administrator\\Desktop\\testSet.txt")
cluster_assment = mixture_Gauss(data_mat, 4)
算法对比
上图从左至右依次是:原始数据集及标签、自己实现的GMM以及Sklearn中实现的GMM,可以清楚地看出,自己实现的算法在数据集比较大时,时间复杂度是很高的,而且迭代次数设置的比较低时算法效果比较差….QAQ,由于每次参数的迭代更新都要对数据集进行遍历,算法效率比较低,另外参数的初始化也影响收敛速度,在Sklearn中初始化方式是一个用户可选的参数,可以设置为“kmeans”或者“random”,当然前者的收敛速度更快(如上图所示)。因此,在运行该代码时,建议在大小为100左右的数据集上查看算法效果(如下图所示),或者以K-Meanas作为前驱算法,之后再运行GMM。
上图所示是在大小为80的数据集上运行的效果,代码还对聚类结果的SSE进行计算并可视化在命令行中,因此可以将迭代停止条件修改为类似“SSE变化很小就停止迭代更新”,如果执行过程中有警告“RuntimeWarning: invalid value encountered in true_divide”,这是由于以零作为分母报的错,可以无视他,也可以添加如下代码“np.seterr(divide=’ignore’, invalid=’ignore’)”。