广州大学学生实验报告-聚类分析
一、实验目的
本实验课程是计算机、人工智能、软件工程等专业学生的一门专业课程，通过实验，帮助学生更好地掌握数据挖掘与机器学习相关概念、技术、原理、应用等；通过实验提高学生编写实验报告、总结实验结果的能力；使学生对机器学习模型、算法等有比较深入的认识。要掌握的知识点如下：
1.掌握机器学习中涉及的相关概念、模型、算法；
2.熟悉机器学习模型训练、验证、测试的流程；
3.熟悉常用的数据预处理方法；
4.掌握聚类分析问题的表示、求解及编程。
二、基本要求
1.实验前，复习《数据挖掘与机器学习》课程中的有关内容。
2.准备好实验数据，编程完成实验内容，收集实验结果。
3.独立完成实验报告。
三、实验软件
推荐使用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算法的原理是：通过循环，让程序不断迭代更新中心点，计算各样本点到新的中心点的距离，根据最近原则重新对样本点进行归类，直到样本点的归类不再发生变化。
混合高斯聚类的原理是：高斯混合模型是一种概率式的聚类方法，它假定所有的数据样本x由k个混合多元高斯分布组合成的混合分布生成。首先，假设样本集可分为k类且每类内符合高斯分布。然后，根据贝叶斯原理利用极大似然法同时求出决定分类的比例，再通过贝叶斯原理求出样本该分在哪个簇。
（2）程序清单（包含详细求解步骤）
import pandas as pd
import matplotlib.pyplot as plt
from numpy import *
import numpy as np
from sklearn.model_selection import train_test_split
import random
data_address='D:/XI/zhuomian/ToDo/学习/machine learning/1122/iris.data'
#查看数据分布情况
def show_data(data):
plt.axis("off")
data.value_counts('sepal_length',sort=False).plot(kind='bar',color='red')
plt.title('sepal_length')
plt.show()
data.value_counts('sepal_width',sort=False).plot(kind='bar',color='blue')
plt.title('sepal_width')
plt.show()
data.value_counts('petal_length',sort=False).plot(kind='bar',color='orange')
plt.title('petal_length')
plt.show()
data.value_counts('petal_width',sort=False).plot(kind='bar',color='green')
plt.title('petal_width')
plt.show()
data.value_counts('class',sort=False).plot(kind='bar',color='yellow')
plt.title('class')
plt.show()
#点图
def point_graph(data,label,title):
colors=label
plt.title(title)
plt.scatter(data.iloc[:,0],data.iloc[:,1],c=colors)
plt.axis("off")
plt.show()
#读取数据
def read_data():
data=pd.read_csv(data_address,header=None,sep=',',names=["sepal_length", "sepal_width", "petal_length", "petal_width","class"])#按逗号分隔读取数据
#查看数据分布情况
# show_data(data)
# 删除缺失的数据
data.replace("?",np.nan,inplace=True)
data.dropna(inplace=True)
# 将类别用数字替换
data.replace("Iris-setosa",0,inplace=True)#山鸢尾
data.replace("Iris-versicolor",1,inplace=True)#变色鸢尾
data.replace("Iris-virginica",2,inplace=True)#维吉尼亚鸢尾
#仅保留聚类特征及类别
temp_1=data.iloc[:,0]
temp_2=data.iloc[:,2]
y=data.iloc[:,-1]#所有数据label
x=pd.concat([temp_1,temp_2],axis=1)#所有数据x
train_X,_,train_y,_=train_test_split(x,y,test_size=0.3,random_state=0)#取70%数据(固定的）用于训练 random_state相同则每次运行划分是相同的
return x,y,train_X,train_y
#更新中心点
def update(c,j,x):#c:数据类别列表 j更新哪一类 x数据
temp_1=[]#分子x0
temp_2=[]#分子x1
temp_3=0#分母
temp_x=0
temp_center1=0
temp_center2=0
for i in range (len(c)):
if c[i]==j:
temp_1.append(x.iloc[i,0])
temp_2.append(x.iloc[i,1])
temp_3+=1
temp_x=i
for i in range(temp_3):
temp_center1+=temp_1[i]
temp_center2+=temp_2[i]
temp_center1/=temp_3
temp_center2/=temp_3
center=[temp_center1,temp_center2]
return center,temp_x
#更新中心点对应的类别
def update_kind(x0,x1,x2,y,center0,center1,center2):
temp_center0=center0
temp_center1=center1
temp_center2=center2
if y.values[x0]==1:
temp_center1=center0
elif y.values[x0]==2:
temp_center2=center0
if y.values[x1]==0:
temp_center0=center1
elif y.values[x1]==2:
temp_center2=center1
if y.values[x2]==0:
temp_center0=center2
elif y.values[x2]==1:
temp_center1=center2
center0=temp_center0
center1=temp_center1
center2=temp_center2
return center0,center1,center2
#欧氏距离
def distance(x0,x1,center0,center1,center2):
c=[]
for i in range(len(x0)):
d0=(x0.iloc[i]-center0[0])**2+(x1.iloc[i]-center0[1])**2
d1=(x0.iloc[i]-center1[0])**2+(x1.iloc[i]-center1[1])**2
d2=(x0.iloc[i]-center2[0])**2+(x1.iloc[i]-center2[1])**2
if d0<d1:
if d0<d2:
c.append(0)
else:
c.append(2)
else:
if d1<d2:
c.append(1)
else:
c.append(2)
return c
#损失函数
def loss_function(y,y_predict):
loss=0
for i in range(len(y)):
loss+=(y[i]-y_predict[i])**2
return loss
#正确率
def correct_rate(y,y_pre):
correct=0
for i in range (len(y)):
if y[i]==y_pre[i]:
correct+=1
correct/=len(y)
print("准确率为：",correct)
#各项性能指标
def performance(C,D):
a=0
b=0
c=0
d=0
for i in range(len(C)):
for j in range(i+1,len(C)):
if C[i]==C[j] and D[i]==D[j]:
a+=1
elif C[i]==C[j] and D[i]!=D[j]:
b+=1
elif C[i]!=C[j] and D[i]==D[j]:
c+=1
else:
d+=1
#JC:有限样本集之间的相似性 越大 样本相似性越高
jc=a/(a+b+c)
#FM
fmi=np.sqrt((a/(a+b))*(a/(a+c)))
#Rand
m=len(C)
randi=(2*(a+b))/(m*(m-1))
print("JC系数:",jc)
print("FM指数:",fmi)
print("Rand指数:",randi)
#k-means聚类
def cluster(x,y):
#在数据集中随机获取三个点作为中心
c0=random.randint(1,len(x))
c1=random.randint(1,len(x))
c2=random.randint(1,len(x))
center0=x.iloc[c0,:].tolist()
center1=x.iloc[c1,:].tolist()
center2=x.iloc[c2,:].tolist()
lossen=[]
for j in range (20):
#当前预测出来的类别
c_list=distance(x.iloc[:,0],x.iloc[:,1],center0,center1,center2)
lossen.append(loss_function(y.values.tolist(),c_list))#计算损失
# x0为0标签数据在原始数据的下标，其余同理
center0,x0=(update(c_list,0,x))
center1,x1=(update(c_list,1,x))
center2,x2=(update(c_list,2,x))
#更新预测出来的中心的实际对应的类别
center0,center1,center2=update_kind(x0,x1,x2,y,center0,center1,center2)
print(center0)
print(center1)
print(center2)
plt.axis("off")
plt.plot(lossen)
plt.show()
return center0,center1,center2
#k-means
def KMeans(x,y,train_x,train_y):
#训练中心点
center0,center1,center2=cluster(train_x,train_y)
#获得预测值
predict_y=[]
predict_y=distance(x.iloc[:,0],x.iloc[:,1],center0,center1,center2)
print("K-means:")
#计算准确度
correct_rate(y,predict_y)
#计算各项指标
performance(predict_y,y)
#绘制预测出来的散点图
point_graph(x,predict_y,"K-means_predict")
#绘制真实的散点图
point_graph(x,y,"K-means_real")
#高斯概率密度函数
def Density_function(x,Mean,Sigma):
n=3
temp_1=np.exp(-0.5*(x-Mean)*(np.linalg.inv(Sigma))*(x-Mean).T)#分子 np.inc()计算矩阵的逆
temp_2=((2*np.pi)**(n/2))*(np.linalg.det(Sigma)**0.5)# 分母 np.linalg.det计算矩阵行列式
return temp_1/temp_2
#极大似然估计
def EM(x,steps=50):
m=x.shape[0] #有m个数据
n=x.shape[1] #有n个特征
#初始化参数
Alpha=[1/3,1/3,1/3]
c0=random.randint(1,len(x))
c1=random.randint(1,len(x))
c2=random.randint(1,len(x))
Mean=[x[c0,:],x[c1,:],x[c2,:]]
Sigma = [mat([[0.1, 0], [0, 0.1]]) for x in range(3)]
Gamma = mat(zeros((m, 3)))
for i in range(steps):
#计算后验概率
for j in range(m):
Sum_Gamma=0
for k in range(3):
Gamma[j,k]=Alpha[k]*Density_function(x[j,:],Mean[k],Sigma[k])
Sum_Gamma+=Gamma[j,k]
for k in range(3):
Gamma[j,k]/=Sum_Gamma
Sum_Gamma=sum(Gamma,axis=0)
#更新
for k in range(3):
Mean[k]=mat(zeros((1,n)))
Sigma[k]=mat(zeros((n,n)))
#更新均值向量
for j in range(m):
Mean[k]+=Gamma[j,k]*x[j,:]
Mean[k]/=Sum_Gamma[0,k]
#更新协方差矩阵
for j in range(m):
Sigma[k]+=Gamma[j,k]*((x[j,:]-Mean[k]).T)*(x[j,:]-Mean[k])
Sigma[k]/=Sum_Gamma[0,k]
#更新Alpha
Alpha[k]=Sum_Gamma[0,k]/m
return Gamma
#高斯混合聚类
def GMM(x,y):
m=x.shape[0] #有m个数据
n=x.shape[1] #有n个特征
centers=zeros((m,n))
#初始化
for i in range(m):
index = int(random.uniform(0, m)) #随机生成一个在0~m的实数作为初始下标
centers[i, :] = x[index, :]
Cluster=mat(zeros((m,2)))
Gamma=EM(x)
#找到各样本对应最大概率下标
for i in range(m):
Cluster[i:]=argmax(Gamma[i,:]),amax(Gamma[i,:])
#更新分类
for i in range(m):
p=x[nonzero(Cluster[:,0].A==i)[0]] #nonzero返回非零下标
centers[i,:]=mean(p,axis=0)
print("GMM:")
# #计算准确度
# correct_rate(y,Cluster[:,0].tolist())
#计算各项指标
performance(Cluster[:,0].tolist(),y)
#可视化
colors=Cluster[:,0].tolist()
plt.title("Gauss_predict")
plt.scatter(x[:,0].tolist(),x[:,1].tolist(),c=colors)
plt.axis("off")
plt.show()
colors=y
plt.title("Gauss_real")
plt.scatter(x[:,0].tolist(),x[:,1].tolist(),c=colors)
plt.axis("off")
plt.show()
return centers,Cluster
if __name__ =='__main__':
#拆分训练集
x,y,train_x,train_y=read_data()
#初始化类别全为0
label_y=np.zeros(len(train_x))
# 查看数据点分布
point_graph(train_x,label_y,"inits")
# k-means
KMeans(x,y,train_x,train_y)
# GMM
GMM(mat(x),y)
（3）展示实验结果，可视化聚类结果
数据集各个维度特征的数值类型与分布如下（仅查看作为特征的两个维度）
原始数据聚类图：
K-means算法预测数据聚类结果：
高斯混合聚类结果：
各项指标：
（4）讨论实验结果，分析k-means聚类数量与聚类指标的关系
从结果看出这两种算法都能将样本较好的分为三类。
关于聚类数量与聚类指标：随着聚类数目增多，每一个类别中数量越来越少，距离越来越近，聚类的效果不一定会增强。
508
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
