聚类分析实验

广州大学学生实验报告-聚类分析

一、实验目的

本实验课程是计算机、人工智能、软件工程等专业学生的一门专业课程,通过实验,帮助学生更好地掌握数据挖掘与机器学习相关概念、技术、原理、应用等;通过实验提高学生编写实验报告、总结实验结果的能力;使学生对机器学习模型、算法等有比较深入的认识。要掌握的知识点如下:
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聚类数量与聚类指标的关系
从结果看出这两种算法都能将样本较好的分为三类。
关于聚类数量与聚类指标:随着聚类数目增多,每一个类别中数量越来越少,距离越来越近,聚类的效果不一定会增强。

  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值