数据选取和数据情况
本次聚类实验仍然选取鸢尾花数据集（http://archive.ics.uci.edu/ml/datasets/Iris）
数据包含5列，分别是花萼长度、花萼宽度、花瓣长度、花瓣宽度、鸢尾花种类。
鸢尾花属种类包含三种：iris-setosa, iris-versicolour, iris-virginica。
鸢尾花的种类不参与聚类，最后用于计算准确率作为聚类效果的评价
每一类分别是50条，共150条数据；每一类在四个属性的分布情况如下图所示
聚类前拿前三个特征可视化（四个特征没法画），可以看到区分度还是很明显的
利用高斯混合模型聚类
模型原理
高斯混合模型（GMM）是假设数据是由N个高斯分布采样所得到的，并不知道每个数据是来自哪个高斯分布，均有概率来自不同的高斯分布，数据来自每个模型的概率即为隐变量。这是含有隐变量的概率模型，故采用EM算法求解，流程如下：
- 先随机生成N个高斯分布，然后判断每个数据属于各个高斯分布的概率，将其归入概率最大的高斯分布中
- 根据上一步中的数据重新估计各个高斯分布的参数
- 循环以上2步，直至似然函数值达到设定的阈值，即不再明显变化
- 对每组数据求属于这N个高斯分布的概率，取概率最大的为该组数据最终聚类的类别
由于初始化高斯分布是随机生成的，可能会导致最后效果不好或者迭代很慢，考虑使用k-means聚类作为GMM的初始分布，可以加速收敛速度和提高聚类效果。
GMM实现
k-means生成模型初始参数
利用k-means聚为N类，根据N类数据估计这N个高斯分布的参数作为初始高斯分布
def get_init_param(self, data):
#k-means聚类来生成初始高斯分布
KMEANS = KMeans(n_clusters=self.K).fit(data)
clusters = defaultdict(list)
for index, label in enumerate(KMEANS.labels_):
clusters[label].append(index)
mu = []
alpha = []
sigma = []
for indexs in clusters.values():
partial_data = data[indexs]
#本次高斯分布的均值向量（共D个特征)
mu.append(partial_data.mean(axis=0))
#本次高斯分布的权重
alpha.append(len(indexs) / self.N)
#本次高斯分布的协方差（D个特征间）
sigma.append(np.cov(partial_data.T))
self.mu = np.array(mu)
self.alpha = np.array(alpha)
self.sigma = np.array(sigma)
return
EM算法迭代训练
利用上面的k-means得到的参数初始化GMM模型参数，接着进行EM算法迭代训练，基本思想是根据目前的各类高斯分布参数重新计算每组数据所属各类的概率，取最大概率的类为所属类；进而根据各类数据重新估计高斯分布参数并更新，循环以上步骤，直至似然函数收敛，训练结束。
def fit(self, data):
# 迭代训练
self.init_param(data)
step = 0
gamma_all_arr = None
while step < self.maxstep:
step += 1
old_alpha = copy.copy(self.alpha)
# E步
gamma_all = []
# 依次求每个样本对K个分模型的响应度
for j in range(self.N):
gamma_j = []
for k in range(self.K):
gamma_j.append(self.alpha[k]*
self._phi(data[j],self.mu[k],self.sigma[k]))
s = sum(gamma_j)
gamma_j = [item/s for item in gamma_j]
gamma_all.append(gamma_j)
gamma_all_arr = np.array(gamma_all)
# M步
for k in range(self.K):
gamma_k = gamma_all_arr[:, k]
SUM = np.sum(gamma_k)
# 更新权重
self.alpha[k] = SUM / self.N
# 更新均值向量, 1*D
new_mu = sum([gamma * y for gamma, y in
zip(gamma_k, data)]) / SUM
self.mu[k] = new_mu
# 更新协方差阵, N*D
delta_ = data - new_mu
# D*D
self.sigma[k] = sum([gamma*(np.outer(np.transpose([delta]), delta)) for gamma, delta in zip(gamma_k, delta_)]) / SUM
alpha_delta = self.alpha - old_alpha
if np.linalg.norm(alpha_delta, 1) < self.epsilon:
break
self.gamma_all_final = gamma_all_arr
return
模型预测
经过EM算法迭代训练后，得到各个高斯分布的参数，计算每组数据求属于这N个高斯分布的概率，取概率最大的为该组数据最终聚类的类别
def predict(self):
#应用模型对数据聚类
pre_label=[]
for j in range(self.N):
max_ind = np.argmax(self.gamma_all_final[j])
pre_label.append(max_ind)
return pre_label
模型聚类效果
以上过程详细代码见附件GMM_cluster.py
根据聚类后的数据，依然用鸢尾花数据集的前三个特征可视化
利用k-means初始化来定义子模型的概率及生成高斯分布的参数，进而利用EM算法迭代训练的GMM聚类算法最终用原鸢尾花数据集的标签计算得到96.67%的准确率，说明效果很好。
绘制混淆矩阵查看各类数据利用GMM模型聚类的情况：
参考资料
利用Python实现高斯混合模型（GMM）
聚类——高斯混合模型及EM算法以及Python实现
本文全部代码可以关注下方博主公众号获取，后台回复：GMM
5万+
