广州大学机器学习与数据挖掘实验三

实验三聚类分析

一、实验目的
本实验课程是计算机、人工智能、软件工程等专业学生的一门专业课程，通过实验，帮助学生更好地掌握数据挖掘与机器学习相关概念、技术、原理、应用等；通过实验提高学生编写实验报告、总结实验结果的能力；使学生对机器学习模型、算法等有比较深入的认识。要掌握的知识点如下：

掌握机器学习中涉及的相关概念、模型、算法；
熟悉机器学习模型训练、验证、测试的流程；
熟悉常用的数据预处理方法；
掌握聚类分析问题的表示、求解及编程。

二、基本要求

实验前，复习《数据挖掘与机器学习》课程中的有关内容。
准备好实验数据，编程完成实验内容，收集实验结果。
独立完成实验报告。

三、实验软件
推荐使用Python编程语言（允许使用numpy库，需实现详细实验步骤，不允许直接调用scikit-learn中回归、分类、聚类等高层API）。

四、实验内容：
基于IRIS鸢尾花数据集，完成关于鸢尾花的聚类分析。

1 准备数据集并认识数据
下载IRIS数据集
https://archive.ics.uci.edu/ml/datasets/iris
了解数据集各个维度特征的含义

2 探索数据并预处理数据
观察数据集各个维度特征的数值类型与分布
挑选sepal length、petal length两维特征作为聚类依据

3 求解聚类中心
编程实现k-means聚类、混合高斯聚类

4 测试和评估模型
在数据集上计算聚类的性能指标

五、学生实验报告
（1）简要介绍k-means、混合高斯聚类的原理
k-means原理：
k-means算法是常用的一种聚类算法。算法的输入为一个样本集（点集），通过该算法可以将样本进行聚类，具有相似特征的样本聚为一类。
算法思想：
假设我们要把数据分成K个类，则可分为以下步骤：
1、随机选k个点，作为聚类中心
2计算每个点分别到k个聚类中心的距离，然后将该点分给距离最近的那个中心点，这样就形成了k个簇
3、再重新计算每个簇的质心（均值）
4、重复2-3步，直到质心的位置不在发生变化或者达到设定的迭代次数

混合高斯聚类原理：
①假设观测数据y1,y2,…,yN由高斯混合模型生成，即

其中

我们用EM算法估计高斯混合模型的参数θ

②初始化模型参数：

③EM算法的E步：
计算

（有点相当于朴素贝叶斯中计算后验概率，用先验概率乘以条件概率）
这是当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据yj的影响度

④EM算法的M步：更新模型参数

重复E步和M步直到模型收敛

（2）程序清单（包含详细求解步骤）
k-means聚类算法：
①要引进的库

②导入数据集，观察数据特点

③挑选sepal length、petal length两维特征作为聚类依据，将类别’class’列的值赋给labels,并进行标签编码

④初始化聚类中心

⑤开始训练：计算每个点分别到k个聚类中心的距离，然后将该点分给距离最近的那个中心点

⑥重新计算每个簇的质心（均值）

⑦ 迭代第⑤、⑥步max_iter=1000次

⑧显示k-means算法的聚类情况，同时也显示实际数据的归类情况

⑨计算准确率（由于聚类算法只会将原始数据样本划分为K个簇，但是并不会告诉我们每个簇分别对应哪个类别，因此我们用排列组合的方法，分别计算每种情况的准确率，选最高的为最终的准确率值）

混合高斯算法：
①要引入的库

②导入数据集，观察数据特点

③挑选sepal length、petal length两维特征作为聚类依据，用data存储这两列数据。labels存储类别‘class’列的值,并进行标签编码

④实例化类GMM_EM的对象gmm

执行类的__init__函数，确定进行聚类的类数为3 即n_components=3

⑤调用对象gmm中的fit_predict函数，得到用混合高斯模型进行聚类的结果

分析fit_predict(data)函数内进行的步骤：
1’进行数据预处理，调用类内函数preprocess()

data数据集中大小为150，特征数为2

2’调用类内函数_init()初始化高斯模型的参数

3’在max_iter迭代次数为1000的情况下执行EM算法中的E、M步骤，当前后两次迭代概率的变化＜1e-6时，即可跳出迭代

E步骤的函数对应该算法

_e_step()函数内调用的guass()函数定义如下：

M步骤的函数对应算法：

⑥求解混合高斯模型聚类中心

⑦画图显示k-means算法的聚类情况及中心点，同时也显示实际数据的归类情况

⑧计算准确率，方法与K-means中的相同

（3）展示实验结果，可视化聚类结果
k-means算法：
k-means算法的聚类情况：

数据的实际归类情况：

该聚类算法得出准确率为：

混合高斯算法：（ EM算法是对初始值敏感的，修改初始化的值会发现模型性能变化很大）

固定一个分类较好的随机情况：

此时混合高斯算法聚类情况：

数据实际分类情况：

此时该聚类算法得出准确率为：

（4）讨论实验结果，分析k-means聚类数量与聚类指标的关系

由于在根据鸢尾花数据集进行编码时，我就先默认聚类数量为3，代码中的许多运算都是固定死了针对3这个聚类数量，代码很不灵活，所以该小节的关系分析没法做出。
在网上查找了一下针对k-means聚类数量与聚类指标的关系问题，但找到的信息无法理解分析。针对k-means聚类数目的确定，得到的相关信息是：k-means分类数目k值很难估计，不确定分成多少类才最合适。

（5）源代码
k-means

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
from sklearn.preprocessing import LabelEncoderiris_data=pd.read_csv("Iris/iris.data",header=None,names=['sepal length','sepal width','petal length','petal width','class'])print(iris_data.info())#发现有iris-setosa   iris-versicolor   iris-virginica三种
print(iris_data['class'].value_counts())#对labels进行标签编码
labels=iris_data['class'].values
label_encoder=LabelEncoder()
labels=label_encoder.fit_transform(labels)
# print(labels)#挑选sepal length、petal length两维特征作为聚类依据
x_axis=iris_data['sepal length']  #series(150,)
y_axis=iris_data['petal length']print(x_axis.shape)
print(y_axis.shape)#随机三个index值，在150条数据集中随机选中三个种类的开始点的标号
indexList=random.sample(range(0,150),3)
print(indexList)#随机初始的中心点
x_center1=x_axis[indexList[0]]
y_center1=y_axis[indexList[0]]
x_center2=x_axis[indexList[1]]
y_center2=y_axis[indexList[1]]
x_center3=x_axis[indexList[2]]
y_center3=y_axis[indexList[2]]print(x_center1)
print(x_axis[0])#---------------------开始训练  训练100次-------------------------for i in range(100):# 用来装分属于三类的数据的index值belong1 = []belong2 = []belong3 = []#计算每条分数据别到3个聚类中心的距离for j in range(150):belong=0          #belong用来记录该条数据属于的类别dis_1=pow((x_axis[j]-x_center1),2)+pow((y_axis[j]-y_center1),2)dis_2=pow((x_axis[j]-x_center2),2)+pow((y_axis[j]-y_center2),2)dis_3=pow((x_axis[j]-x_center3),2)+pow((y_axis[j]-y_center3),2)#比较离三类中心点哪个更近，将该条数据归于距离更近的中心点所属类if dis_2<dis_1:belong=2if dis_3<dis_2:belong=3else:belong=1if dis_3<dis_1:belong=3# print(belong)if belong==1:belong1.append(j)elif belong==2:belong2.append(j)else:belong3.append(j)#进行center点的位置更新for k in range(len(belong1)):x_center1+=x_axis[belong1[k]]y_center1+=y_axis[belong1[k]]for k in range(len(belong2)):x_center2 += x_axis[belong2[k]]y_center2 += y_axis[belong2[k]]for k in range(len(belong3)):x_center3 += x_axis[belong3[k]]y_center3 += y_axis[belong3[k]]x_center1=x_center1/(1+len(belong1))x_center2=x_center2/(1+len(belong2))x_center3=x_center3/(1+len(belong3))y_center1 = y_center1 / (1 + len(belong1))y_center2 = y_center2 / (1 + len(belong2))y_center3 = y_center3 / (1 + len(belong3))#y_pred用来装k-means聚类算法对每条数据所归到的类
#注意这里的类值1、2、3没有实际意义，跟实际数据labels中用标签编码所得的值0、1、2值没有对应关系
#仅仅是为了区分类别
y_pred=np.array(np.zeros(150))for i in range(len(belong1)):y_pred[belong1[i]]=1for i in range(len(belong2)):y_pred[belong2[i]]=2for i in range(len(belong3)):y_pred[belong3[i]]=3#k-means求出的聚类中心
x_center=[x_center1,x_center2,x_center3]
y_center=[y_center1,y_center2,y_center3]
#数据集的实际中心
x_ac_center=[x_axis[0:50].mean(),x_axis[50:100].mean(),x_axis[100:150].mean()]
y_ac_center=[y_axis[0:50].mean(),y_axis[50:100].mean(),y_axis[100:150].mean()]
#画图
#聚类算法的归类情况
plt.scatter(x_axis,y_axis,c=y_pred)
plt.scatter(x_center,y_center,c='r',marker='x')
plt.show()#实际数据的归类情况
plt.scatter(x_axis,y_axis,c=labels)
plt.scatter(x_ac_center,y_ac_center,c='r',marker='x')
plt.show()#计算准确率
#分别计算三种类别组合  0 1 2    1 0 2   0 2 1    1 2 0    2 1 0   2 0 1
y_pred_1=np.array(np.zeros(150))
y_pred_2=np.array(np.zeros(150))
y_pred_3=np.array(np.zeros(150))
y_pred_4=np.array(np.zeros(150))
y_pred_5=np.array(np.zeros(150))
y_pred_6=np.array(np.zeros(150))
for i in range(150):if y_pred[i]==1:y_pred_1[i]=0y_pred_2[i] = 1y_pred_3[i] = 0y_pred_4[i] = 1y_pred_5[i] = 2y_pred_6[i] = 2if y_pred[i]==2:y_pred_1[i] = 1y_pred_2[i] = 0y_pred_3[i] = 2y_pred_4[i] = 2y_pred_5[i] = 1y_pred_6[i] = 0if y_pred[i]==3:y_pred_1[i] = 2y_pred_2[i] = 2y_pred_3[i] = 1y_pred_4[i] = 0y_pred_5[i] = 0y_pred_6[i] = 1def correct_rate(lei_list):correct_num = 0for i in range(150):if (lei_list[i] == labels[i]):correct_num += 1rate = correct_num / 150return raterate1=correct_rate(y_pred_1)
rate2=correct_rate(y_pred_2)
rate3=correct_rate(y_pred_3)
rate4=correct_rate(y_pred_4)
rate5=correct_rate(y_pred_5)
rate6=correct_rate(y_pred_6)#比较
rate=[rate1,rate2,rate3,rate4,rate5,rate6]
max_rate=0
for i in range(6):if rate[i]>max_rate:max_rate=rate[i]print('准确率为：',max_rate)

混合高斯聚类：

from scipy.stats import multivariate_normal
from sklearn import preprocessing
from sklearn.metrics import accuracy_score
import numpy as np
import pandas as pd
from sklearn.preprocessing import LabelEncoder
import matplotlib.pyplot as pltclass GMM_EM():def __init__(self, n_components, max_iter=1000, error=1e-6):self.n_components = n_components  # 混合模型由几个gauss模型组成self.max_iter = max_iter  # 最大迭代次数self.error = error  # 收敛误差self.samples = 0   #样本个数self.features = 0 #存储特征个数self.alpha = []  # 存储模型权重self.mu = []  # 存储均值self.sigma = []  # 存储标准差def _init(self, data):  # 初始化参数np.random.seed(4)self.mu = np.array(np.random.rand(self.n_components, self.features))#sigma中为每个分模型都初始化了一个二维随机变量的协方差矩阵self.sigma = np.array([np.eye(self.features) / self.features] * self.n_components)self.alpha = np.array([1.0 / self.n_components] * self.n_components)print(self.alpha.shape, self.mu.shape, self.sigma.shape)print(self.alpha,self.mu,self.sigma)def gauss(self, Y, mu, sigma):  # 直接调用多元正态分布的概率密度函数，计算高斯函数值return multivariate_normal.pdf(Y,mean=mu, cov=sigma )def preprocess(self, data):  # 数据预处理self.samples = data.shape[0]   #定义数据集大小self.features = data.shape[1]   #定义数据集的特征数pre = preprocessing.MinMaxScaler()  #进行了特征归一化return pre.fit_transform(data)def fit_predict(self, data):  # 拟合数据data = self.preprocess(data)   #进行数据预处理self._init(data)     #初始化模型参数weighted_probs = np.zeros((self.samples, self.n_components))print(weighted_probs.shape)  #用来存放进行E步后算出来的当前模型参数下每个观测数据来自第k个分模型的概率 shape(150,3)for i in range(self.max_iter):prev_weighted_probs = weighted_probs#e步weighted_probs = self._e_step(data)#change 当后验概率没有什么改变了的时候，即收敛的时候，停止迭代change = np.linalg.norm(weighted_probs - prev_weighted_probs)if change < self.error:break#m步self._m_step(data, weighted_probs)#比较每个观测数据来自三个分模型的概率，返回概率最大的分模型所代表列的列号return weighted_probs.argmax(axis=1)def _e_step(self, data):  # E步probs = np.zeros((self.samples, self.n_components))   #shape(150,3)for i in range(self.n_components):#调用类定义的gauss函数，计算不同高斯模型下数据集对应高斯函数值probs[:, i] = self.gauss(data, self.mu[i, :], self.sigma[i, :, :])weighted_probs = np.zeros(probs.shape)for i in range(self.n_components):weighted_probs[:, i] = self.alpha[i] * probs[:, i]for i in range(self.samples):#后验概率某类概率再除以三类的概率和weighted_probs[i, :] /= np.sum(weighted_probs[i, :])return weighted_probsdef _m_step(self, data, weighted_probs):  # M步  进行mu,sigma,alpha的数值更新for i in range(self.n_components):#求出每一列的概率和 即每行数据属于某一具体类的概率和sum_probs_i = np.sum(weighted_probs[:, i])#axis=0计算每列的sumself.mu[i, :] = np.sum(np.multiply(data, np.mat(weighted_probs[:, i]).T), axis=0) / sum_probs_iself.sigma[i, :, :] = (data - self.mu[i, :]).T * np.multiply((data - self.mu[i, :]),np.mat(weighted_probs[:, i]).T) / sum_probs_i#行数  shape[0]self.alpha[i] = sum_probs_i / data.shape[0]iris_data=pd.read_csv("Iris/iris.data",header=None,names=['sepal length','sepal width','petal length','petal width','class'])#发现有iris-setosa   iris-versicolor   iris-virginica三种
print(iris_data['class'].value_counts())
labels=iris_data['class'].values#对labels进行标签编码
label_encoder=LabelEncoder()
labels=label_encoder.fit_transform(labels)
# print(labels)#挑选sepal length、petal length两维特征作为聚类依据
x_axis=iris_data['sepal length']  #series(150,)
y_axis=iris_data['petal length']data=np.array(pd.concat([x_axis,y_axis],axis=1))gmm = GMM_EM(3)
pre_label = gmm.fit_predict(data)print(pre_label)
print(labels)#混合高斯算法求出的聚类中心
num_0,num_1,num_2=[0,0,0]
xsum_0,xsum_1,xsum_2=[0,0,0]
ysum_0,ysum_1,ysum_2=[0,0,0]for i in range(len(pre_label)):if pre_label[i]==0:num_0+=1xsum_0+=x_axis[i]ysum_0 += y_axis[i]elif pre_label[i]==1:num_1+=1xsum_1+=x_axis[i]ysum_1 += y_axis[i]else:num_2+=1xsum_2+=x_axis[i]ysum_2 += y_axis[i]x_center_0=xsum_0/num_0
y_center_0=ysum_0/num_0
x_center_1=xsum_1/num_1
y_center_1=ysum_1/num_1
x_center_2=xsum_2/num_2
y_center_2=ysum_2/num_2x_center=[x_center_0,x_center_1,x_center_2]
y_center=[y_center_0,y_center_1,y_center_2]
#数据集的实际中心
x_ac_center=[x_axis[0:50].mean(),x_axis[50:100].mean(),x_axis[100:150].mean()]
y_ac_center=[y_axis[0:50].mean(),y_axis[50:100].mean(),y_axis[100:150].mean()]#画图
#画混合高斯聚类图
plt.scatter(x_axis,y_axis,c=pre_label)
plt.scatter(x_center,y_center,c='r',marker='x')
plt.xlabel('sepal length')
plt.ylabel('petal length')
plt.show()#实际数据图
plt.scatter(x_axis,y_axis,c=labels)
plt.scatter(x_ac_center,y_ac_center,c='r',marker='x')
plt.xlabel('sepal length')
plt.ylabel('petal length')
plt.show()# EM算法是对初始值敏感的，修改初始化的值会发现模型性能变化很大#计算准确率
#分别计算三种类别组合  0 1 2    1 0 2   0 2 1    1 2 0    2 1 0   2 0 1
y_pred_1=np.array(np.zeros(150))
y_pred_2=np.array(np.zeros(150))
y_pred_3=np.array(np.zeros(150))
y_pred_4=np.array(np.zeros(150))
y_pred_5=np.array(np.zeros(150))
y_pred_6=np.array(np.zeros(150))
for i in range(150):if pre_label[i]==0:y_pred_1[i]=0y_pred_2[i] = 1y_pred_3[i] = 0y_pred_4[i] = 1y_pred_5[i] = 2y_pred_6[i] = 2if pre_label[i]==1:y_pred_1[i] = 1y_pred_2[i] = 0y_pred_3[i] = 2y_pred_4[i] = 2y_pred_5[i] = 1y_pred_6[i] = 0if pre_label[i]==2:y_pred_1[i] = 2y_pred_2[i] = 2y_pred_3[i] = 1y_pred_4[i] = 0y_pred_5[i] = 0y_pred_6[i] = 1def correct_rate(lei_list):correct_num = 0for i in range(150):if (lei_list[i] == labels[i]):correct_num += 1rate = correct_num / 150return raterate1=correct_rate(y_pred_1)
rate2=correct_rate(y_pred_2)
rate3=correct_rate(y_pred_3)
rate4=correct_rate(y_pred_4)
rate5=correct_rate(y_pred_5)
rate6=correct_rate(y_pred_6)#比较
rate=[rate1,rate2,rate3,rate4,rate5,rate6]
max_rate=0
for i in range(6):if rate[i]>max_rate:max_rate=rate[i]print('准确率为：',max_rate)