美文网首页生信算法生信工具生信算法
利用Kmeans聚类分析两类问题

利用Kmeans聚类分析两类问题

作者: apricoter | 来源:发表于2019-04-04 17:47 被阅读13次

    聚类分析是一种无监督的学习方法,根据一定条件将相对同质的样本归到一个类总(俗话说人以类聚,物以群分)
    正式一点的:聚类是对点集进行考察并按照某种距离测度将他们聚成多个“簇”的过程。聚类的目标是使得同一簇内的点之间的距离较短,而不同簇中点之间的距离较大。

    两种方法对比:

    在K-means聚类中,是预先规定出要产生多少个类别的数量,再根据类别数量自动聚成相应的类。对K-means而言,首先是随机产生于类别数相同的初始点,然后判断每个点与初始点的距离,每个点选择最近的一个初始点,作为其类别。

    当类别产生后,在计算各个类别的中心点,然后计算每个点到中心点的距离,并根据距离再次选择类别。当新类别产生后,再次根据中心点重复选择类别的过程,直到中心点的变化不再明显。最终根据中心点产生的类别,就是聚类的结果。正如图中所示,一组对象中需要生成三个类别,各个类别之间都自然聚焦在一起。

    在层次聚类中,不需要规定出类别的数量,最终聚类的数量可以根据人为要求进行划分。对层次聚类,首先每个对象都是单独的类别,通过比较两两之间距离,首先把距离最小的两个对象聚成一类。接着把距离次小的聚成一类,然后就是不断重复按距离最小的原则,不断聚成一类的过程,直到所有对象都被聚成一类。

    在层次聚类中,可以以一张树状图来表示聚类的过程,如果要讲对象分类的话,就可以从根节点触发,按照树状图的分叉情况,划分出不同的类别来。在图中,把一组对象分成了三个类别,可见这三个类别就是构成了树状图最开始的三个分支。

    k-means

    首先,随机选择K个对象,并且所选择的每个对象都代表一个组的初始均值或初始的组中心值;对剩余的每个对象,根据其与各个组初始均值的距离,将它们分配给最近的(最相似)小组;然后,重新计算每个小组新的均值;这个过程不断重复,直到所有的对象在K组分布中都找到离自己最近的组。

    优点:容易实现。

    缺点:可能收敛到局部最小值,在大规模数据集上收敛较慢;
    先指定k,同时对异常值很敏感。

    聚类技术在数据分析和数据化运营中的主要用途表现在:既可以直接作为模型对观察对象进行群体划分,为业务方的精细化运营提供具体的细分依据和相应的运营方案建议,又可在数据处理阶段用作数据探索的工具,包括发现离群点、孤立点,数据降维的手段和方法,通过聚类发现数据间的深层次的关系等。

    一、已知聚类簇数的iris数据集

    # 导入第三方包
    import pandas as pd
    import numpy as np  
    import matplotlib.pyplot as plt
    from sklearn.cluster import KMeans
    from sklearn import metrics
    # 读取iris数据集
    iris = pd.read_csv(r'F:\iris.csv')
    # 查看数据集的前几行
    iris.head()
    

    不存在量纲上的差异,无需做标准化处理

    聚类簇数为3,

    # 提取出用于建模的数据集X
    X = iris.drop(labels = 'Species', axis = 1)
    # 构建Kmeans模型
    kmeans = KMeans(n_clusters = 3)
    kmeans.fit(X)
    # 聚类结果标签
    X['cluster'] = kmeans.labels_
    # 各类频数统计
    X.cluster.value_counts()
    

    各簇样本量分别为62,50,38

    对比建模前后差异

    # 导入第三方模块
    import seaborn as sns
    # 中文和负号的正常显示
    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
    plt.rcParams['axes.unicode_minus'] = False
    import matplotlib.pyplot as plt
    # 设置绘图风格
    plt.style.use('ggplot')
    # 三个簇的簇中心
    centers = kmeans.cluster_centers_
    # 绘制聚类效果的散点图
    sns.lmplot(x = 'Petal_Length', y = 'Petal_Width', hue = 'cluster', markers = ['^','s','o'], 
               data = X, fit_reg = False, scatter_kws = {'alpha':0.8}, legend_out = False)
    plt.scatter(centers[:,2], centers[:,3], marker = '*', color = 'black', s = 130)
    plt.xlabel('花瓣长度')
    plt.ylabel('花瓣宽度')
    # 图形显示
    plt.show()
    

    以上为聚类效果的散点图,五角星为每个簇的簇中心

    # 增加一个辅助列,将不同的花种映射到0,1,2三种值,目的方便后面图形的对比
    iris['Species_map'] = iris.Species.map({'virginica':0,'setosa':1,'versicolor':2})
    # 绘制原始数据三个类别的散点图
    sns.lmplot(x = 'Petal_Length', y = 'Petal_Width', hue = 'Species_map', data = iris, markers = ['^','s','o'],
               fit_reg = False, scatter_kws = {'alpha':0.8}, legend_out = False)
    plt.xlabel('花瓣长度')
    plt.ylabel('花瓣宽度')
    # 图形显示
    plt.show()
    

    以上为原始数据的散点图,与聚类图对比,标记为1的与原始数据吻合,0和2存在一些错误分割,但还是比较一致

    对比样本差异使用雷达图,导入pygal模块

    # 导入第三方模块
    import pygal
    # 调用Radar这个类,并设置雷达图的填充,及数据范围
    radar_chart = pygal.Radar(fill = True)
    # 添加雷达图各顶点的名称
    radar_chart.x_labels = ['花萼长度','花萼宽度','花瓣长度','花瓣宽度']
    
    # 绘制三个雷达图区域,代表三个簇中心的指标值
    radar_chart.add('C1', centers[0])
    radar_chart.add('C2', centers[1])
    radar_chart.add('C3', centers[2])
    # 保存图像
    radar_chart.render_to_file('radar_chart.svg')
    

    雷达图无法通过plt.show展示,通过浏览器打开svg文件


    二、未知聚类簇数的NBA球员数据集

    # 读取球员数据
    players = pd.read_csv(r'F:\players.csv')
    players.head()
    
    # 绘制得分与命中率的散点图
    sns.lmplot(x = '得分', y = '命中率', data = players, 
               fit_reg = False, scatter_kws = {'alpha':0.8, 'color': 'steelblue'})
    plt.show()
    
    from sklearn import preprocessing
    # 数据标准化处理
    X = preprocessing.minmax_scale(players[['得分','罚球命中率','命中率','三分命中率']])
    # 将数组转换为数据框
    X = pd.DataFrame(X, columns=['得分','罚球命中率','命中率','三分命中率'])
    X.head()
    

    重点在于选择最佳k值

    1.拐点法

    # 构造自定义函数,用于绘制不同k值和对应总的簇内离差平方和的折线图
    def k_SSE(X, clusters):
        # 选择连续的K种不同的值
        K = range(1,clusters+1)
        # 构建空列表用于存储总的簇内离差平方和
        TSSE = []
        for k in K:
            # 用于存储各个簇内离差平方和
            SSE = []
            kmeans = KMeans(n_clusters=k)
            kmeans.fit(X)
            # 返回簇标签
            labels = kmeans.labels_
            # 返回簇中心
            centers = kmeans.cluster_centers_
            # 计算各簇样本的离差平方和,并保存到列表中
            for label in set(labels):
                SSE.append(np.sum((X.loc[labels == label,]-centers[label,:])**2))
            # 计算总的簇内离差平方和 
            TSSE.append(np.sum(SSE))
    
        # 中文和负号的正常显示
        plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
        plt.rcParams['axes.unicode_minus'] = False
        # 设置绘图风格
        plt.style.use('ggplot')
        # 绘制K的个数与GSSE的关系
        plt.plot(K, TSSE, 'b*-')
        plt.xlabel('簇的个数')
        plt.ylabel('簇内离差平方和之和')
        # 显示图形
        plt.show()
        
    # 使用拐点法选择最佳的K值
    k_SSE(X, 15)
    

    当k在4附近,折线斜率的变动不是很大,故k为3,或4或5

    2,轮廓系数法

    # 构造自定义函数,用于绘制不同k值和对应轮廓系数的折线图
    def k_silhouette(X, clusters):
        K = range(2,clusters+1)
        # 构建空列表,用于存储个中簇数下的轮廓系数
        S = []
        for k in K:
            kmeans = KMeans(n_clusters=k)
            kmeans.fit(X)
            labels = kmeans.labels_
            # 调用字模块metrics中的silhouette_score函数,计算轮廓系数
            S.append(metrics.silhouette_score(X, labels, metric='euclidean'))
    
        # 中文和负号的正常显示
        plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
        plt.rcParams['axes.unicode_minus'] = False
        # 设置绘图风格
        plt.style.use('ggplot')    
        # 绘制K的个数与轮廓系数的关系
        plt.plot(K, S, 'b*-')
        plt.xlabel('簇的个数')
        plt.ylabel('轮廓系数')
        # 显示图形
        plt.show()
        
    # 使用轮廓系数选择最佳的K值
    k_silhouette(X, 15)
    

    k=2时轮廓系数最大

    3,间隔统计量法

    # 自定义函数,计算簇内任意两样本之间的欧氏距离
    def short_pair_wise_D(each_cluster):
        mu = each_cluster.mean(axis = 0)
        Dk = sum(sum((each_cluster - mu)**2)) * 2.0 * each_cluster.shape[0]
        return Dk
    
    # 计算簇内的Wk值
    def compute_Wk(data, classfication_result):
        Wk = 0
        label_set = set(classfication_result)
        for label in label_set:
            each_cluster = data[classfication_result == label, :]
            Wk = Wk + short_pair_wise_D(each_cluster)/(2.0*each_cluster.shape[0])
        return Wk
    
    # 计算GAP统计量 
    def gap_statistic(X, B=10, K=range(1,11), N_init = 10):
        # 将输入数据集转换为数组
        X = np.array(X)
        # 生成B组参照数据
        shape = X.shape
        tops = X.max(axis=0)
        bots = X.min(axis=0)
        dists = np.matrix(np.diag(tops-bots))
        rands = np.random.random_sample(size=(B,shape[0],shape[1]))
        for i in range(B):
            rands[i,:,:] = rands[i,:,:]*dists+bots
        
        # 自定义0元素的数组,用于存储gaps、Wks和Wkbs
        gaps = np.zeros(len(K))
        Wks = np.zeros(len(K))
        Wkbs = np.zeros((len(K),B))
        # 循环不同的k值,
        for idxk, k in enumerate(K):
            k_means =  KMeans(n_clusters=k)
            k_means.fit(X)
            classfication_result = k_means.labels_
            # 将所有簇内的Wk存储起来
            Wks[idxk] = compute_Wk(X,classfication_result)
            
            # 通过循环,计算每一个参照数据集下的各簇Wk值
            for i in range(B):
                Xb = rands[i,:,:]
                k_means.fit(Xb)
                classfication_result_b = k_means.labels_
                Wkbs[idxk,i] = compute_Wk(Xb,classfication_result_b)
    
        # 计算gaps、sd_ks、sk和gapDiff
        gaps = (np.log(Wkbs)).mean(axis = 1) - np.log(Wks)        
        sd_ks = np.std(np.log(Wkbs), axis=1)
        sk = sd_ks*np.sqrt(1+1.0/B)
        # 用于判别最佳k的标准,当gapDiff首次为正时,对应的k即为目标值
        gapDiff = gaps[:-1] - gaps[1:] + sk[1:]
        
        # 中文和负号的正常显示
        plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
        plt.rcParams['axes.unicode_minus'] = False
        # 设置绘图风格
        plt.style.use('ggplot')
        # 绘制gapDiff的条形图
        plt.bar(np.arange(len(gapDiff))+1, gapDiff, color = 'steelblue')
        plt.xlabel('簇的个数')
        plt.ylabel('k的选择标准')
        plt.show()
        
    # 使用间隙统计量选择最佳的K值
    gap_statistic(X, B = 20, K=range(1, 16))
    

    纵坐标首次为正时k=3

    综合考虑以上3种,选择k=3

    基于k值进行聚类

    # 将球员数据集聚为3类
    kmeans = KMeans(n_clusters = 3)
    kmeans.fit(X)
    # 将聚类结果标签插入到数据集players中
    players['cluster'] = kmeans.labels_
    # 构建空列表,用于存储三个簇的簇中心
    centers = []
    for i in players.cluster.unique():
        centers.append(players.ix[players.cluster == i,['得分','罚球命中率','命中率','三分命中率']].mean())
    # 将列表转换为数组,便于后面的索引取数
    centers = np.array(centers)
    centers
    
    # 绘制散点图
    sns.lmplot(x = '得分', y = '命中率', hue = 'cluster', data = players, markers = ['^','s','o'],
               fit_reg = False, scatter_kws = {'alpha':0.8}, legend = False)
    # 添加簇中心
    plt.scatter(centers[:,0], centers[:,2], c='k', marker = '*', s = 180)
    plt.xlabel('得分')
    plt.ylabel('命中率')
    # 图形显示
    plt.show()
    

    需要注意的是,由于对原数据做了标准化处理,簇中心不能直接使用cluster_centers_得到,返回的是原数据标准化后的中心,需要通过For循环重新找到原始数据下的簇中心,即五角星



    可以得到高得分高命中率型诸如此类

    再看四个指标上的差异,由于四个维度上量纲不一致,需要使用标准化后的中心点绘制雷达图

    # 雷达图
    # 调用模型计算出来的簇中心
    centers_std = kmeans.cluster_centers_
    # 设置填充型雷达图
    radar_chart = pygal.Radar(fill = True)
    # 添加雷达图各顶点的名称
    radar_chart.x_labels = ['得分','罚球命中率','命中率','三分命中率']
    
    # 绘制雷达图代表三个簇中心的指标值
    radar_chart.add('C1', centers_std[0])
    radar_chart.add('C2', centers_std[1])
    radar_chart.add('C3', centers_std[2])
    # 保存图像
    radar_chart.render_to_file('radar_chart.svg')
    

    C2、C3得分没有差异,但命中率C2比C3高很多诸如此类结论

    相关文章

      网友评论

        本文标题:利用Kmeans聚类分析两类问题

        本文链接:https://www.haomeiwen.com/subject/vdqdvqtx.html