看了很多博文,包括《统计学习知识》和西瓜书上对GMM算法的推导,总有些重要的步骤被略去(比如从公式一推到公式二,书上直接给出结果,却没有具体步骤),导致理解整个算法非常困难。后来幸运地发现一篇博文,使用了对我而言易于理解的语言,重要把整个推导过程疏通成功,最后在纸上手推了一遍,真是酣畅淋漓!算法实现很简单,结构跟K-均值形似,参数的推导过程不用体现在代码上,直接根据推导出来的公式计算就行(所以说,仅写代码还容易些呢,但是不懂原理,背公式有什么意思呢)。以下是一元(即数据集为m*1)的聚类算法,扩展成多元,则要改变miu, delta2扩展成矩阵和协方差矩阵。以后有需求了再写吧,核心思想是一样的。供参考
另感谢这篇博文的博主! https://www.cnblogs.com/mindpuzzle/archive/2013/04/24/3036447.html
import numpy as np
class GMM():
'''
一元高斯混合聚类算法
'''
def randTheta(self, dataSet, k):
'''
初始化模型参数
:param dataSet: 数据集,m*1, np.array
:param k: 指定的聚类数 k>=2
:return: alpha(k,), miu(k,), delta2(k,)
'''
alpha = np.ones(k) / k # 初始化每个alpha为1/k
miu = dataSet.min() + (dataSet.max() - dataSet.min()) * np.random.rand(k, 1) # 从[最小值,最大值)中随机取k个数
delta2 = np.ones(k) / 10.0 # 初始化每个delta2为0.1,delta2意为delta的平方
return alpha, miu.ravel(), delta2
def calcLambda(self, X, alpha, miu, delta2):
'''
计算数据X在每个高斯分布下的lambda(根据推导,此lambda就是后验概率)
:param X: 一元数据(即单个实数,X属于R)
:param alpha: 模型参数,相当于各高斯分布的权重
:param miu: 模型参数,高斯分布的参数min
:param delta2: 模型参数,高斯分布的参数delta2
:return: lambda(k,)
'''
k = len(alpha)
lamb = np.zeros(k) # 定义好lamb,此处针对一个样本数据
temp = np.sum(alpha / (np.sqrt(2 * np.pi) * np.sqrt(delta2)) * np.exp(-(X - miu) ** 2 / (2 * delta2)))
for i in range(k):
lamb[i] = alpha[i] / (np.sqrt(2 * np.pi) * np.sqrt(delta2[i])) * np.exp(-(X - miu[i]) ** 2 / (2 * delta2[i])) / temp # 根据公式计算
return lamb
def gMM(self, dataSet, k):
'''
一元高斯混合分布算法
:param dataSet: 数据集,m*1, np.array
:param k: 指定的聚类数 k>=2
:return: alpha(k,), miu(k,), delta2(k,), 样本的信息矩阵(m*2),存储所属高斯分布的索引和后验概率。
'''
m = dataSet.shape[0]
alpha, miu, delta2 = self.randTheta(dataSet, k) # 初始化参数
lamb = np.zeros((m, k)) # 定义好lamb,此处针对所有样本数据
clusterAss = np.zeros((m, 2)) # 定义最终需要输出的样本信息:所属高斯分布的索引和后验概率
clusterChanged = True # 循环控制开关
while clusterChanged:
clusterChanged = False
for i in range(m):
lamb[i] = self.calcLambda(dataSet[i], alpha, miu, delta2)
maxLamb = lamb[i].max() # 后验概率最大的高斯分布,即为样本所属分布
if int(clusterAss[i][0]) != list(lamb[i]).index(maxLamb): # 如果样本所属分布有变动,那么继续循环
clusterChanged = True
clusterAss[i, :] = list(lamb[i]).index(maxLamb), maxLamb # 更新样本信息
alpha = np.sum(lamb, axis=0) / m # 根据公式更新alpha, 利用array数组,计算非常方便
miu = np.sum((lamb * dataSet), axis=0) / np.sum(lamb, axis=0) # 同上
delta2 = np.sum((lamb * np.power(dataSet - miu, 2)), axis=0) / np.sum(lamb, axis=0) # 同上
return alpha, miu, delta2, clusterAss
# 以下是测试数据,观察分类结果,可以发现数值相近的值被很好的分在一起
dataSet = np.random.rand(100, 1) # 创建100*1个[0,1)的浮点数
gmm = GMM()
alpha, miu, delta, clusterAss = gmm.gMM(dataSet, 4)
print(dataSet[np.nonzero(clusterAss[:, 0] == 0)[0], :].ravel())
print(dataSet[np.nonzero(clusterAss[:, 0] == 1)[0], :].ravel())
print(dataSet[np.nonzero(clusterAss[:, 0] == 2)[0], :].ravel())
print(dataSet[np.nonzero(clusterAss[:, 0] == 3)[0], :].ravel())