机器学习（二）：k-means算法（基础篇）

机器学习（二）：k-means算法（基础篇）

k均值聚类算法（k-means clustering algorithm）是一种迭代求解的聚类分析算法，其步骤是随机选取K个对象作为初始的聚类中心，然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心。聚类中心以及分配给它们的对象就代表一个聚类。每分配一个样本，聚类的聚类中心会根据聚类中现有的对象被重新计算。这个过程将不断重复直到满足某个终止条件。终止条件可以是没有（或最小数目）对象被重新分配给不同的聚类，没有（或最小数目）聚类中心再发生变化，误差平方和局部最小。

K-means算法的实现原理

计算距离的方法

每种距离的计算方式都有自己的限制，再具体使用中需要结合样本的特点和场景不同合理选择描述样本之间相似性的算法，在这里我们说一下常见的几种算法。

• 欧式距离
  上面公式是一个通用公式，当p=2时即欧式距离，欧氏距离是我们最常用的距离计算方式，也成为几何距离，表示两个点的距离差平方和，简单直观
• 曼哈顿距离
  当p=1时表示这个距离，计算的是两个点之间的实地走的距离，欧式距离时两点之间的空间直线距离，不考虑其它影响
• 切比雪夫距离
  当p趋于无穷大时，公式如下：
  
• 向量內积
  向量的內积在很多场合都在使用，向量內积是在向量空间里面，公式如下：
  
  內积的结果是一个标量，大小没有限制，上下没有界限因此不好表示两个向量的相关性，于是有人提出来使用余弦来表示相似性，它有大小限制容易表示相似性。
  

损失函数（目标函数）SSE

误差平方和（SSE）：表示样本中的每个样本点到质心的距离的平方和，最优解为：聚类结果应使 SSE 达到最小值。用通俗易懂的话讲：就是每个样本中有一个质心，当样本中的点到质心的距离平方和最小的时候，该质心为我们所求的质心。

算法实现步骤

1. List item选择合适 k 值，从数据集可以得出 k=3
2. 随机生成 k 个质点作为起始质点
3. 将数据集中的每一个点分配到每一个簇，即每一个点找到距离近的质心作为所对应的簇
4. 对每一个簇，计算簇中所有点的均值将其作为新的质心
5. 重复 cd 两个过程直到分配的所有点不再改变

K-means算法的实践

基本任务

1. 理解 k-means 算法的思想。
2. 使用 Python 实现 k-means 算法。
3. 使用 UCI 上面的 Iris 数据集进行算法测试。
4. 对参数 k 进行调整，记录结果。

扩展任务

1. 现二分 K-means 代码并进行测试（会在下一节中提到）

数据集

• 网址 http://archive.ics.uci.edu/ml/index.php
• 内容：Iris 数据集

第一步：对数据集进行分析：

1） 数据集介绍：
Iris.data 数据集主要有如下：

sl	花萼长度
sw	花萼宽度
pl	花瓣长度
pw	花瓣宽度
variety	花的品种

2）数据特征相关分析：
采用Preson相关系数对四个特征值进行两两之间的相关性分析
3）数据特征画图：
使用python对数据之间两两关系图画出来

代码如下：

# -*- coding:utf-8 -*-
import math
import warnings
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import matplotlib.patches as mpatchs

warnings.filterwarnings('ignore')   #忽略警告
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

def loadDataSet(filename):
    """
    函数说明：从文件中下载数据，并将分离除连续型变量和标签变量
    :parameter:
            data - Iris数据集
            attributes - 鸢尾花的属性
            type - 鸢尾花的类别
            sl-花萼长度 , sw-花萼宽度, pl-花瓣长度, pw-花瓣宽度
    :return:
    """
    iris_data = pd.read_csv(filename)   #打开文件
    iris_data = pd.DataFrame(data=np.array(iris_data), columns=['sl', 'sw', 'pl', 'pw', 'type'], index=range(149))   #给数据集添加列名，方便后面的操作
    attributes = iris_data[['sl', 'sw', 'pl', 'pw']]   #分离出花的属性
    iris_data['type'] = iris_data['type'].apply(lambda x: x.split('-')[1])  # 最后类别一列，感觉前面的'Iris-'有点多余即把class这一列的数据按'-'进行切分取切分后的第二个数据
    labels = iris_data['type']            #分理出花的类别
    return attributes, labels

def showdatas(attributes, datinglabels):
    """
    函数说明：画出花的属性两两之间的关系图
    :parameter:
            attributes - 花的属性
            datinglabels - 花的类别
    :return:none
    """
    fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(20, 8))   #定义一个3行2列的画布
    LabelsColors = []    #建立一个颜色标签列表
    for i in datinglabels:   #遍历花的类型
        if i == 'setosa':    #setosa类型的花画成黑色的点
            LabelsColors.append('black')
        if i == 'versicolor':   #versicolor类型的花画成橙色的点
            LabelsColors.append('orange')
        if i == 'virginica':    #virginica类型的花画成红色的点
            LabelsColors.append('red')

    #在画板第一行第一列的位置绘制花萼长度和花萼宽度之间的关系
    axs[0][0].scatter(x=attributes['sl'], y=attributes['sw'], color=LabelsColors, s=15, alpha=.5) #x轴为花萼长度，y轴为花萼宽度， 点大小为15， 透明度为0.5
    axs0_title_text = axs[0][0].set_title(u'花萼长度和花萼宽度')     #设置title
    axs0_xlabel_text = axs[0][0].set_xlabel(u'花萼长度')            #设置x轴的标签
    axs0_ylabel_text = axs[0][0].set_ylabel(u'花萼宽度')            #设置y轴的标签
    plt.setp(axs0_title_text, size=9, weight='bold', color='red')
    plt.setp(axs0_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs0_ylabel_text, size=7, weight='bold', color='black')

    #在画板第一行第二列的位置绘制花萼长度和花瓣长度之间的关系
    axs[0][1].scatter(x=attributes['sl'], y=attributes['pl'], color=LabelsColors, s=15, alpha=.5) #x轴为花萼长度，y轴为花瓣长度，点的大小为15， 透明度为0.5
    axs1_title_text = axs[0][1].set_title(u'花萼长度和花瓣长度')    #设立title
    axs1_xlabel_text = axs[0][1].set_xlabel(u'花萼长度')           #设置x轴标签
    axs1_ylabel_text = axs[0][1].set_ylabel(u'花瓣长度')           #设置y轴标签
    plt.setp(axs1_title_text, size=9, weight='bold', color='red')
    plt.setp(axs1_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs1_ylabel_text, size=7, weight='bold', color='black')

    #在画板第二行第一列的位置绘制花萼长度与花瓣宽度之间的关系
    axs[1][0].scatter(x=attributes['sl'], y=attributes['pw'], color=LabelsColors, s=15, alpha=.5) #x轴为花萼长度，y轴为花瓣长度，点的大小为15， 透明度为0.5
    axs2_title_text = axs[1][0].set_title(u'花萼长度和花瓣宽度')    #设立title
    axs2_xlabel_text = axs[1][0].set_xlabel(u'花萼长度')           #设立x轴标签
    axs2_ylabel_text = axs[1][0].set_ylabel(u'花瓣宽度')          #设立y轴标签
    plt.setp(axs2_title_text, size=9, weight='bold', color='red')
    plt.setp(axs2_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs2_ylabel_text, size=7, weight='bold', color='black')

    #在画板第二行第二列的位置上绘制花萼宽度与花瓣长度之间的关系
    axs[1][1].scatter(x=attributes['sw'], y=attributes['pl'], color=LabelsColors, s=15, alpha=.5)  #x轴为花萼宽度，y轴为花瓣长度
    axs3_title_text = axs[1][1].set_title(u'花萼宽度和花瓣长度')    #设立title
    axs3_xlabel_text = axs[1][1].set_xlabel(u'花萼宽度')          #设立x轴标签
    axs3_ylabel_text = axs[1][1].set_ylabel(u'花瓣长度')          #设立y轴标签
    plt.setp(axs3_title_text, size=9, weight='bold', color='red')
    plt.setp(axs3_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs3_ylabel_text, size=7, weight='bold', color='black')

    #在画板第三行第一列的位置绘制花萼宽度与花瓣宽度之间的关系
    axs[2][0].scatter(x=attributes['sw'], y=attributes['pw'], color=LabelsColors, s=15, alpha=.5)  #x轴为花萼宽度，y轴为花瓣宽度
    axs4_title_text = axs[2][0].set_title(u'花萼宽度和花瓣宽度')   #设立title
    axs4_xlabel_text = axs[2][0].set_xlabel(u'花萼宽度')          #设立x轴坐标
    axs4_ylabel_text = axs[2][0].set_ylabel(u'花瓣宽度')          #设立y轴坐标
    plt.setp(axs4_title_text, size=9, weight='bold', color='red')
    plt.setp(axs4_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs4_ylabel_text, size=7, weight='bold', color='black')

    #在画板第三行第二列的位置绘制花瓣长度和花瓣宽度之间的关系
    axs[2][1].scatter(x=attributes['pl'], y=attributes['pw'], color=LabelsColors, s=15, alpha=.5)  #x轴花瓣长度，y轴为花瓣宽度
    axs5_title_text = axs[2][1].set_title(u'花瓣长度和花瓣宽度')
    axs5_xlabel_text = axs[2][1].set_xlabel(u'花瓣长度')
    axs5_ylabel_text = axs[2][1].set_ylabel(u'花瓣宽度')
    plt.setp(axs5_title_text, size=9, weight='bold', color='red')
    plt.setp(axs5_xlabel_text, size=7, weight='bold', color='black')
    plt.setp(axs5_ylabel_text, size=7, weight='bold', color='black')

    #设置图例
    setosa = mlines.Line2D([], [], color='black', marker='.', markersize=6, label='setosa')   #设置setosa的图例为黑色的点，大小为6
    versicolor = mlines.Line2D([], [], color='orange', marker='.', markersize=6, label='versicolor')   #设置yersicolor的图例为橙色的点，大小为6
    virginica = mlines.Line2D([], [], color='red', marker='.', markersize=6, label='virginica')  #设置virginica的图例为红色的点，大小为6

    axs[0][0].legend(handles=[setosa,versicolor,virginica])   #对每一个图形设置图例
    axs[0][1].legend(handles=[setosa,versicolor,virginica])
    axs[1][0].legend(handles=[setosa,versicolor,virginica])
    axs[1][1].legend(handles=[setosa, versicolor, virginica])
    axs[2][0].legend(handles=[setosa, versicolor, virginica])
    axs[2][1].legend(handles=[setosa, versicolor, virginica])
    #绘制图形
    plt.show()

def autoNorm(attributes):
    """
    函数说明： 对数据进行归一化
    :parameter
            attributes - 特征矩阵
    :return:  nonormAttributes - 归一化后的矩阵
    """
    attributes = attributes.values  #将DataFrame类型转变为array类型
    minVal = attributes.min()      #找出数据中的最小值
    maxVal = attributes.max()      #找出数据中的最大值
    range = maxVal - minVal        #数据范围
    normAttributes = np.zeros(np.shape(attributes))  #初始化归一化数据
    m = attributes.shape[0]     #获取数据的行数
    normAttributes = attributes - np.tile(minVal, (m, 1))  #创建一个全是最小值得数组
    normAttributes = normAttributes / np.tile(range, (m, 1))  #创建一个全是范围值得数组
    return normAttributes   #返回归一化后的数据

def colcPer(attributes):
    """
    函数说明：计算每个特征之间的Person相关系数
    :parameter:
        attributes - 特征值
    :return: none
    """
    attributes = attributes.astype('float32')  #将attributes数据类型更改成float32
    colcPerAttributes = np.corrcoef(attributes, rowvar=0)   #以每一列计算相关系数
    colcPerAttributes = pd.DataFrame(data=colcPerAttributes, columns=['sl', 'sw', 'pl', 'pw'], index=['sl', 'sw', 'pl', 'pw'])  # 给数据集添加列名，方便观察
    print(colcPerAttributes)   #输出相关系数矩阵

if __name__ == '__main__':
    df = "iris.data"      #文件路径
    attributes, labels = loadDataSet(df)   #得到特征矩阵和标签变量
    showdatas(attributes, labels)      #输出相关的散点图
    normAttributes = autoNorm(attributes)  #进行归一化处理
    colcPer(normAttributes)    #计算相关系数

计算出来的相关系数矩阵记录如下：

数据特征图如下：

第二步：使用算法进行聚类

1）算法步骤：

1. 选择合适 k 值，从数据集可以得出 k=3
2. 随机生成 k 个质点作为起始质点
3. 将数据集中的每一个点分配到每一个簇，即每一个点找到距离近的质心作为所对应的簇
4. 对每一个簇，计算簇中所有点的均值将其作为新的质心
5. 重复 cd 两个过程直到分配的所有点不再改变

2）算法的相关函数实现：

• loadDataSet()： 读入数据，得到四个特征变量
• distEclud(): 欧式距离公式
• randCent(): 生成随机 k 个质心
• Kmeans(): 主要算法，得到局部最优质心，和分类结果
• chooseK(): 画出肘部图
• writeTxt(): 写入文件

3）代码实现：

# -*- coding:utf-8 -*-
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

warnings.filterwarnings('ignore')   #忽略警告
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

def loadDataSet(filename):
    """
    函数说明：从文件中下载数据，并将分离除连续型变量和标签变量
    :parameter:
            data - Iris数据集
            attributes - 鸢尾花的属性
            type - 鸢尾花的类别
            sl-花萼长度 , sw-花萼宽度, pl-花瓣长度, pw-花瓣宽度
    :return:
    """
    iris_data = pd.read_csv(filename)   #打开文件
    iris_data = pd.DataFrame(data=np.array(iris_data), columns=['sl', 'sw', 'pl', 'pw', 'type'], index=range(149))   #给数据集添加列名，方便后面的操作
    attributes = iris_data[['sl', 'sw', 'pl', 'pw']]   #分离出花的属性
    iris_data['type'] = iris_data['type'].apply(lambda x: x.split('-')[1])  # 最后类别一列，感觉前面的'Iris-'有点多余即把class这一列的数据按'-'进行切分取切分后的第二个数据
    labels = iris_data['type']     #分理出花的类别
    attriLabels = []      #建立一个标签列表
    for label in labels:        #为了更方便操作，将三中不同的类型分别设为1，2，3
        if label == 'setosa':    #如果类别为setosa的话，设为1
            attriLabels.append(1)
        elif label == 'versicolor':  #如果是versicolor的时候设为2
            attriLabels.append(2)
        elif label == 'virginica':  #如果是virginica的时候设为3
            attriLabels.append(3)
    attributes = attributes.values
    return attributes, attriLabels

def randCent(dataSet, k):
    """
    函数说明：随机初始化k个质心（质心满足数据边界之内）
    :param dataSet:数据集
    :param k:  质点个数
    :return:centroids 返回初始化得到的k个质心向量
    """
    m, n = dataSet.shape  # 得到行列数目
    centroids = np.zeros((k, n))  # 定义一个k行n列的全为0的质点集合
    for j in range(n):
        #得到该列数据的最小值
        minJ = min(dataSet[:, j])
        #得到该列数据的范围
        rangeJ = float(max(dataSet[:, j]) - minJ)
        #k个质心向量的第j维数据值随机为位于范围内的某一个值
        centroids[:, j] = np.mat(minJ + rangeJ * np.random.rand(k))
    return centroids


def distEclud(vecA,vecB):
    """
    函数说明：欧式距离公式
    :param vecA:样本点
    :param vecB：质心
    :return:距离
    """
    return np.sqrt(np.sum((vecA - vecB) ** 2))

def KMeans(dataSet, k, distMeas=distEclud, createCent=randCent):
    """
    函数说明：k均值算法
    :param dataSet: 特征集
    :param k: k个值
    :return:centroids, clusterAssment
    """
    m = np.shape(dataSet)[0]     #得到行的数目
    clusterAssment = np.mat(np.zeros((m, 2))) #初始化一个（m，2）矩阵
    clusterChange = True      #定义聚类是否发生变化
    centroids = createCent(dataSet, k)        #创建初始化的k个质心向量
    while clusterChange:          #只要聚类结果一直发生变化，就一直执行聚类算法，直到所有数据点聚类结果不发生变化
        clusterChange = False  #聚类结果定义为False
        for i in range(m):     #遍历数据集中的每一个样本
            minDist = np.inf       #初始化最小距离为100000
            minIndex = -1    #最小距离索引定为-1
            for j in range(k): #循环k个类的质心
                distance = distMeas(centroids[j, :], dataSet[i, :]) #计算数据点到质心的欧式距离
                if distance < minDist:  #如果当前距离少于最小距离
                    minDist = distance#当前距离定为最小距离，最小距离的索引定为j
                    minIndex = j
            if clusterAssment[i, 0] != minIndex:  #当前聚类结果中第i个样本的结果发生变化
                clusterChange = True    #把clusterChange定义为Ture，代表发生了变化
            clusterAssment[i, :] = minIndex, minDist**2     #更新当前新变化的聚类结果和错误平方
        for j in range(k):                #遍历每一个质心
            #因此首先先比较clusterAssment[:,0].A==cent的真假，如果为真则记录了他所在的行，因此在用切片进行取值。
            # print(clusterAssment[:, 0].A == j)
            if (clusterAssment[:, 0].A == j).all() == False:   #再chooseK防止报错
                continue
            pointsInCluster = dataSet[np.nonzero(clusterAssment[:, 0].A == j)[0]]
            # 计算这些数据的均值（axis=0：求列的均值），作为该类质心向量
            centroids[j, :] = np.mean(pointsInCluster, axis=0)
    # #返回k个聚类，聚类结果和误差
    return centroids, clusterAssment

def chooseK(dataSet, i):
    """
    函数说明：肘部图的绘画
    :param dataSet: 数据集
    :param i:k从1开始迭代到k次
    :return:none
    """
    list = []  #定义存放距离的列表，即y轴
    x= []      #定义存放1，2，3，4，5，即x轴
    for j in range(1, i):
        #得到聚点， 聚类结果，和误差
        cenList, clusterAssment = KMeans(dataSet, j)
        #计算每个点的误差平方和，并加入到列表中
        list.append(sum(clusterAssment[:, 1]))
        #添加x轴
        x.append(j)
    #将list转变为列表
    a = np.array(list)
    list = a.reshape(-1)
    #以x为画图
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(x, list)
    plt.show()


def writeTxt(cenList1, clusterAssment, type=1):
    """
    函数说明：将质心和分类结果写入文件的操作。
    :param cenList1:质心
    :param clusterAssment:分类结果
    :param type: 1为二分k-means， 2为k-means
    :return: none
    """
    #将clusterAssment变成列表
    clusterAssment = clusterAssment.tolist()
    #获得cenList的长度
    n = len(cenList1)
    #获得clusterAssment的长度
    m = len(clusterAssment)
    #判断是否为k-means
    if type == 1:
        file = open('consequeue01.txt', mode='w')   #如果是二分k-means则打开01文件
        file.write("二分k-means的聚类后的质心")
    else:
        file = open('consequeue02.txt', mode='w')   #如果是k-means则打开02文件
        file.write("k-means的聚类后的质心")
    for j in range(n):
        file.write("\n第%d个质心为：" % (j+1))  #输入质心
        file.write(str(cenList1[j]))           #输入质心
    file.write("\n聚类结果:\n")
    for i in range(m):
        file.write('第%d个属性被归类为：' % (i+1))     #输入类别
        file.write(str(int(clusterAssment[i][0])))
        file.write("       ")
        file.write("距离为：")
        file.write(str(clusterAssment[i][1]))    #输入距离
        file.write("\n")
    file.close()


#k-means算法
if __name__ == '__main__':
     filename = "iris.data"      #文件路径
     attributes, labels = loadDataSet(filename)   #得到数据
     centList, clusterAssment = KMeans(attributes, 3)   #k=3时得到质心和分类结果
     chooseK(attributes, 7)     #画图
     writeTxt(centList, clusterAssment, 2)   #写入文件

我们看看代码运行后的效果图：

很棒，聚类结果和三个质心都写进了记事本里面。

第四步：对K-means参数值的优化

k-means存在着一定的不足：
1, K 值的取值，簇的数目是用户预先定义的参数，如何选择正确的 k 是这个算法需要改进的地方
2, K-均值算法收敛到局部最小值，而不是全局最小值。原因是由于随机初始化质心导致收敛效果比较差。

于是对于k值优化，我采用了肘部图对k进行选择。原理：让 k 从 1 取值到 7，计算
其 SSE，利用绘图很明显的看出当 k 取一定值的时候，SSE 函数下降
最快，而后下降最小。（函数：chooseK()）:

缺陷：可能由于数据原因，当 k 取值为 4 的时候，分类完后数据只有 3 个类别，导致了 k 取 3，4，5 时候 SSE 并不一定会递减。这个问题改了很多种方法，都不能够解决。

PS：下一节会实现k-means的进阶任务：二分k-means。