1、自己实现KMean算法

上一课,我们已经使用sklearn做了一次KMean聚类,但是遗憾的是,那个算法不是我们自己写的。使用别人的算法可以加快试验的进程,尽快的完成作业,但是不利于我们了解算法的原理,也不利于我们锻炼分析问题,解决问题的能力。

要知道在10多年的学习生涯中,老师们一直告诉我们分析、解决问题的能力很重要,这也是以后的绝对竞争力。人工智能能不能学好,关键在于分析、解决问题的能力,这些能力不能靠使用别人的代码提高,所以,虽然我们写教程很辛苦,但是,我们还是打算给大家讲一讲。

本章我们将自己实现KMean算法,效果如下图所示:

本例中的样本被分成了6类,每一类用不同颜色和形状的图形标记。

2、KMean算码代码详解

这是至今为止,我们写得最长的一段python代码,耐心的来看看吧,代码中大部分语句,我们都加了注释,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
# -*- coding: utf-8 -*-
 
from numpy import *
import matplotlib.pyplot as plt
import operator
import time
 
# 无限大的数
INF = 9999999.0
 
def loadDataSet(fileName, splitChar='\t'):
    """
    输入:文件名
    输出:数据集
    描述:从文件读入数据集
    """
    dataSet = []
    with open(fileName) as fr:
        for line in fr.readlines():
            curline = line.strip().split(splitChar)
            # map()是 Python 内置的高阶函数,它接收一个函数 f 和一个 list,
            # 并通过把函数 f 依次作用在 list 的每个元素上,得到一个新的 list 并返回。
            # 这里map(float, curline)的含义就是把curline中的每一个元素都转换为浮点数,没转换之前他们是字符串类型。
            # python中list就是数组
            fltline = list(map(float, curline))
            # 将样本数据放到dataSet中。
            dataSet.append(fltline)
    # 返回的数组结构[[1,2],[2,3],[1,4]......]这样的形式。
    return dataSet
 
def distEclud(vecA, vecB):
    """
    输入:向量A, 向量B
    输出:两个向量的欧式距离
    """
    # 两个向量的欧式距离
    return sqrt(sum(power(vecA - vecB, 2)))
 
def randCent(dataSet, k):
    """
    输入:数据集, 聚类个数
    输出:k个随机质心的矩阵
    """
    # 没一个数据样本的维度,本例是2
    n = shape(dataSet)[1]
    # zeros生成一个值为0的k,n维矩阵,mat是转换为矩阵
    # k目前为6,n为2,这里是随机取6个中心点
    centroids = mat(zeros((k, n)))
    for j in range(n):
        # 取各列的最小值
        minJ = min(dataSet[:, j])
        # 取各列的最大值,并用最大值减去最小值,从而得到各个维度的取值范围
        rangeJ = float(max(dataSet[:, j]) - minJ)
        # random.rand(k, 1),生成6个0到1的随机数,
        centroids[:, j] = minJ + rangeJ * random.rand(k, 1)
    return centroids
 
 
def kMeans(dataSet, k, distMeans=distEclud, createCent=randCent):
    """
    输入:数据集, 聚类个数, 距离计算函数, 生成随机质心函数
    输出:质心矩阵, 簇分配和距离矩阵
    """
    # 样本的条数
    m = shape(dataSet)[0]
 
    clusterAssment = mat(zeros((m, 2)))
    # 创建k个中心点
    centroids = createCent(dataSet, k)
    #
    clusterChanged = True
    while clusterChanged:
        clusterChanged = False
        # m个样本,
        for i in range(m): # 寻找最近的质心
            # 距离最小值
            minDist = INF
            # 距离最大值
            minIndex = -1
            # 每个点离K个中心的距离
            for j in range(k):
                # centroids[j, :]表示取第j行的数据,dataSet[i, :]表示取第i个样本
                # 计算他们之间的距离
                distJI = distMeans(centroids[j, :], dataSet[i, :])
                # 找到第i个样本,离那个中心最近
                if distJI < minDist:
                    minDist = distJI
                    minIndex = j
            # 第i和样本属于哪一个类,如果不属于minIndex这个类,那么说明该样本点的分类有所改变,
            # 那么说明还需要进行一次迭代计算
            if clusterAssment[i, 0] != minIndex:
                clusterChanged = True
            # 表示第i个样本,属于minIndex类,离这个类的中心距离是minDist**2
            clusterAssment[i, :] = minIndex, minDist**2
 
        # 更新中心的位置
        for cent in range(k):
            # nonzeros(a)返回数组a中值不为零的元素的下标
            # clusterAssment[:, 0].A == cent 表示
            ptsInClust = dataSet[nonzero(clusterAssment[:, 0].A == cent)[0]]
 
            centroids[cent, :] = mean(ptsInClust, axis=0)
 
    # 当聚类迭代结束后,返回各类的最终中心和每个样本最终属于的类别
    return centroids, clusterAssment
 
def plotFeature(dataSet, centroids, clusterAssment):
    """
    绘制分类的图表
    :param dataSet:
    :param centroids:
    :param clusterAssment:
    :return:
    """
    # 中心的个数,本例是6个中心
    m = shape(centroids)[0]
    # 新建一个图表
    fig = plt.figure()
    # 每一类使用不同的标记
    scatterMarkers = ['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<']
    # 每一类使用不同的颜色
    scatterColors = ['blue', 'green', 'yellow', 'purple', 'orange', 'black', 'brown']
    # 参数111的意思是:将画布分割成1行1列,图像画在从左到右从上到下的第1块
    ax = fig.add_subplot(111)
    #
    for i in range(m):
        ptsInCurCluster = dataSet[nonzero(clusterAssment[:, 0].A == i)[0], :]
        markerStyle = scatterMarkers[i % len(scatterMarkers)]
        colorSytle = scatterColors[i % len(scatterColors)]
        ax.scatter(ptsInCurCluster[:, 0].flatten().A[0], ptsInCurCluster[:, 1].flatten().A[0], marker=markerStyle, c=colorSytle, s=90)
    # 绘制点
    ax.scatter(centroids[:, 0].flatten().A[0], centroids[:, 1].flatten().A[0], marker='+', c='red', s=300)
 
def main():
    # 加载数据,请使用绝对路径
    dataSet = loadDataSet('./data_folder/points.txt', splitChar=',')
    # mat 将list转换为矩阵
    dataSet = mat(dataSet)
    resultCentroids, clustAssing = kMeans(dataSet, 6)
    print('打印中心的位置')
    print(resultCentroids)
    print('*******************')
    plotFeature(dataSet, resultCentroids, clustAssing)
 
if __name__ == '__main__':
    # 计算该算法花费的时间
    start = time.clock()
    main()
    end = time.clock()
    print('该算法花费时间 %s' % str(end - start))
    plt.show()

下面,我们将要分析一下代码,但是,因为注释中,已经把基本原理讲清楚了,所以大家一定结合注释来看哦。

首先,从145行开始看,代码首先启动一个时钟time.clock,用来统计这个算法一共会用多少时间。

第149和150行,计算出了程序运行的时间,最终打印了出来。

第148行,是main函数,主要的算法入口。我们下面分析一下main函数中的第一个函数loadDataSet。

3、loadDataSet函数详解

代码中第11行的loadDataSet函数,是解析文件 points.txt中的数据,并将其转换为一个样本数组。

point.txt中存放的是样本数据,如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
15.55,28.65
14.9,27.55
14.45,28.35
14.15,28.8
13.75,28.05
13.35,28.45
13,29.15
13.45,27.5
13.6,26.5
12.8,27.35
12.4,27.85
......
......

这里,每一行是一个样本点数据,分别代码x,y坐标,x和y用逗号隔开。 loadDataSet函数的主要功能,就是解析这个文本文件,并转化为列表数组。loadDataSet函数代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def loadDataSet(fileName, splitChar='\t'):
    """
    输入:文件名
    输出:数据集
    描述:从文件读入数据集
    """
    dataSet = []
    with open(fileName) as fr:
        for line in fr.readlines():
            curline = line.strip().split(splitChar)
            # map()是 Python 内置的高阶函数,它接收一个函数 f 和一个 list,
            # 并通过把函数 f 依次作用在 list 的每个元素上,得到一个新的 list 并返回。
            # 这里map(float, curline)的含义就是把curline中的每一个元素都转换为浮点数,没转换之前他们是字符串类型。
            # python中list就是数组
            fltline = list(map(float, curline))
            # 将样本数据放到dataSet中。
            dataSet.append(fltline)
    # 返回的数组结构[[1,2],[2,3],[1,4]......]这样的形式。
    return dataSet

代码第1行,fileName参数表示解析哪一个文件,splitChar表示每一行的分隔符,默认为tab键。

第8行到17行是一个with语句,这里表示不管在处理文件过程中是否发生异常,都能保证 with 语句执行完毕后,能够自动关闭了打开的文件句柄。它比try/finally更简单一些。

第9行,读文件的每一行,并赋值给line变量

第10行,去掉每一行左右的空格,然后以分隔符分开

第15含,将curline中的每一个数据转换为float

第17行,将解析出来的2个数字列表放到dataSet中

第19行,最后返回dataSet。

4、randCent函数详解,随机生成中心点

回到kMean函数的69行,首先是初始化中心的位置,调用了createCent函数,这个函数等于randCent,它的原型如下:

def randCent(dataSet, k):

randCent生成随机的中心点,参数dataSet表示样本数据,k表示需要产生几个中心点。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def randCent(dataSet, k):
    """
    输入:数据集, 聚类个数
    输出:k个随机质心的矩阵
    """
    # 没一个数据样本的维度,本例是2
    n = shape(dataSet)[1]
    # zeros生成一个值为0的k,n维矩阵,mat是转换为矩阵
    # k目前为6,n为2,这里是随机取6个中心点
    centroids = mat(zeros((k, n)))
    for j in range(n):
        # 取各列的最小值
        minJ = min(dataSet[:, j])
        # 取各列的最大值,并用最大值减去最小值,从而得到各个维度的取值范围
        rangeJ = float(max(dataSet[:, j]) - minJ)
        # random.rand(k, 1),生成6个0到1的随机数,
        centroids[:, j] = minJ + rangeJ * random.rand(k, 1)
    return centroids

第10行,生成6行2列的0矩阵,用于存放聚类中心。第10行生成的centroids打印的结果是:

5、distEclud函数详解,两个样本之间的距离

回到第140行,其调用了kMeans函数,传入了数据集,并准备将数据聚为6类。其返回值是聚类中心和每个样本属于哪一个类。

kMeans函数中大部分代码注释都很详细,请仔细看代码,直到84行distMeans函数,计算2个点的距离。这个函数被distEculd实现。

distEclud函数用来计算两个坐标(或者叫向量)之间的距离。这里使用的欧式距离,就是A点到B点的距离。

distEclud的代码如下:

1
2
3
4
5
6
7
def distEclud(vecA, vecB):
    """
    输入:向量A, 向量B
    输出:两个向量的欧式距离
    """
    # 两个向量的欧式距离
    return sqrt(sum(power(vecA - vecB, 2)))

第7行,power(vecA – vecB, 2)表示计算两个向量差的平方。

sum表示将平方差加起来。

sqrt表示开根号求距离,这段代码就是计算出的距离,如下公式:

6、plotFeature函数

plotFeature函数用来绘制分类的结果,函数原型如下:

def plotFeature(dataSet, centroids, clusterAssment)

参数dataSet: 表示数据集。

参数centroids:每个类的中心。

参数clusterAssment:每个样本属于哪一类。

plotFeature的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def plotFeature(dataSet, centroids, clusterAssment):
    """
    绘制分类的图表
    :param dataSet:
    :param centroids:
    :param clusterAssment:
    :return:
    """
    # 中心的个数,本例是6个中心
    m = shape(centroids)[0]
    # 新建一个图表
    fig = plt.figure()
    # 每一类使用不同的标记
    scatterMarkers = ['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<']
    # 每一类使用不同的颜色
    scatterColors = ['blue', 'green', 'yellow', 'purple', 'orange', 'black', 'brown']
    # 参数111的意思是:将画布分割成1行1列,图像画在从左到右从上到下的第1块
    ax = fig.add_subplot(111)
    #
    for i in range(m):
        ptsInCurCluster = dataSet[nonzero(clusterAssment[:, 0].A == i)[0], :]
        markerStyle = scatterMarkers[i % len(scatterMarkers)]
        colorSytle = scatterColors[i % len(scatterColors)]
        ax.scatter(ptsInCurCluster[:, 0].flatten().A[0], ptsInCurCluster[:, 1].flatten().A[0], marker=markerStyle, c=colorSytle, s=90)
    # 绘制点
    ax.scatter(centroids[:, 0].flatten().A[0], centroids[:, 1].flatten().A[0], marker='+', c='red', s=300)

第10行, 取centroids的有多少行,有多少行就代表有多少个中心点。

第12行,调用figure新建一个图表,在这个图表上可以绘制点、线等数据。

第14行,marker这个单词表示标记的意思,也就是每一类使用什么形状的标记绘制在图中。scatterMarkers中每一个字母表示一种形状,关于形状的具体情况,请参考matplot的定义。举一个例来说,s表示方块的意思。

第16行,表示每一个分类的颜色,可以直接使用英文的颜色来表示,这个matplot是支持的。

第18行,add_subplot的参数是一个三位数,前2位表示将一个图分为几行几列,第三个数,表示在第几个位置绘图。这里的111,就是绘制一个图的意思。

第20到24行,将分类点用scatter散点图绘制出来。