算法-聚类-K均值与模糊K均值：原理+python代码

这篇文章是根据作业修改后得到的，个人感觉写的比较详细了。但还有许多不足，希望大家评论指出。

K均值聚类与模糊K均值

1. 算法原理及流程

相关名词解释如表1。
表1-相关名词解释

1.1 K均值算法原理

K均值是基于原型的、划分的聚类计数。K均值使用质心定义原型，质心是一组点的均值。
1）首先，选择K个初始质心（通常随机初始化），其中K是用户预期的簇个数。
2）每个点指派到最近的质心，而指派到一个质心的点集合为一个簇。这里的距离量度有多种，通常采用欧几里得距离：

SSE为误差的平方和（Sum of Squared Error, SSE），为第个簇，为中的点，为第个簇的均值。
3）然后，根据指派到簇的点，更新每个簇的质心。

4）重复指派和更新步骤，直到簇不发生变化，或等价为质心不发生变化。

1.2 K均值算法伪代码流程

1.3 模糊K均值算法原理

在大部分情况下，数据集中的对象并不能划分为明显分离的簇，且指派一个对象到一个特定的簇也有一定的任意性。此时采用模糊的概念，对每个对象和每个簇赋予一个权值，指明该对象属于该簇的程度。这样可以更好的对该点进行划分。

1.4 模糊K均值算法伪代码流程

2. 算法相关问题

2.1 质心定义问题

2.2 聚类效果评估方法

（1）轮廓系数
轮廓系数（Silhouette Coefficient），是聚类效果好坏的一种评价方式。它结合了凝聚度和分离度两种因素。可以用来在相同原始数据的基础上用来评价不同算法、或者算法不同运行方式对聚类结果所产生的影响。

轮廓系数的值处于-1~1之间，值越大，表示聚类效果越好。具体计算方法如下：

（2）手肘法
手肘法利用SSE曲线判别最佳簇数K，SSE公式在上面已给出。通过对不同簇K进行统计，可以得到手肘样的曲线，在肘处的SSE较小，且随着K的增加，SSE逐渐趋于平稳。
（3）相似度矩阵
相似度矩阵按照簇标号调整数据对象距离矩阵的行列次序，并将数据对象间的距离转换成相似度。明显分离的簇会在会显示很强的块对角模式，可以较形象的反应聚类效果。但是对于大型数据集，该方法开销较大，复杂度为O（n2）。

注：不同距离量度下的聚类效果本应采用相似度矩阵等宏观检测，但是由于时间关系仅采用了SSE与轮廓系数来一定程度上反应聚类效果，SSE与轮廓系数中的距离量度均统一为欧式距离。

3. K均值仿真测试

本文代码在python3.7环境下运行，数据集自定义在.txt文件中。

3.1 不同簇数量K下的仿真

初始质心的选取采用随机数生成的方式，计算输入数据集x与y维的最小值与最大值，然后计算两者的范围。则每个质心的x与y维值为最小值加上范围之间的某个值。距离度量统一采用欧式距离，分类结果的评价采用手肘法与轮廓系数结合的方法进行评价。从簇中可以看到存在较小簇时，不容易让其单独分出来（簇K=5），个人感觉K=4合K=6时的分类效果较好，K=4时能够将较大的簇单独分类处理，K=6时能够将小的簇分类出来，但是有较大的簇被分为两类。

从图中曲线可以较为明显看出较高的轮廓系数出现在K=4时，第2个峰值出现在K=6时。而对于SSE曲线则呈不断衰减的趋势，到K=6以后衰减速率较慢。因此综合来看，簇K=6时的分类效果较好。

3.2 不同距离量度下的仿真

除了欧式距离，这里再采用三种其他距离作为对比。

其中对于曼哈顿距离要对质心的选取要采用中位数，而不是取均值。这里分析了四种不同的距离量度，初始质心的选取与3.1一致，簇K按照数据分布分为5簇数据。对分类好的数据分别拿出来对比了一下，并记录了每种距离量度下的SSE与轮廓系数，然后进行了对比。

从图3中可以看出前三种距离量度分类较好，但是对于余弦距离来说，分类结果比较不理想。但是再看下SSE与轮廓系数曲线。

图中1，2，3，4分别对应欧几里得距离、曼哈顿距离、切比雪夫距离和余弦距离。从图中可以看出曼哈顿距离轮廓系数最小，且SSE最高，说明该方法对这种数据的距离效果不佳。而欧几里得与切比雪夫的轮廓系数与SSE相差无几，对于余弦距离的数据就有点夸张了，轮廓系数有0.729，而SSE有0.003，但从肉眼观测来看该分类并不是很好。因为余弦距离是两向量之间的角度余弦值，故这种数据的分类其实并不适合。

3.3 不同初始质心选取条件下的仿真

原来是采用范围内取随机数的方法进行指派初始质心，这样的效果有时较好有时较差，并不是特别稳定。

层次聚类是另一种指派质心的方式，通过层次聚类，划分k个层次，计算出每个簇对应的质心作为K-Means算法的初始质心。这种方法可以很好地解决初始质心指派不合理的问题。但是也有局限性。

K-Means++方法中初始质心的选取也是一种较新的方式，第一个质心是随机选择的，接下来的质心基于样本点与最近质心的距离，距离越大越可能被选为下一个质心，直到选择完k个质心。该方法有效地解决了关于初始质心的选取问题，目前已经成为了一种硬聚类算法的标准。但是该方法无法解决离群点问题。

这里仅采用了K-Means++中初始质心选取的方法作为对比，在K=5条件下，采用欧几里得距离度量，并统计轮廓系数，K-Means++普遍比随机质心的轮廓系数高。

4. 模糊K均值仿真测试

在原来函数体的框架下新建了新的函数，改变了质心更新公式、SSE公式，并添加了每个点对不同簇的权值计算公式。在与上面仿真类似的情况下，对模糊K均值进行了同样的仿真。

4.1 不同簇数量K下的仿真

这里并没有一一列举不同簇数量K下的仿真图，其效果与K均值类似。这里将K均值与模糊K均值在K=5的情况下进行了仿真对比，其实效果相差不大。但是模糊K均值较多出现将下面的那个大簇分为两类，而较少识别出较小的簇，这样来看这个对比其实效果不大。再看一下SSE曲线与轮廓曲线，SSE曲线与K均值的类似，也是在K=6时开始趋向平稳，SSE较低的原因是因为每个点与质心的距离乘了权值的缘故。再看轮廓系数曲线，与K均值在K=6时明显不同，且呈不断降低的趋势，整体数值相差不大。K取4或5轮廓系数会更好。

4.2不同距离量度下的仿真

从图中可以看出余弦度量还是与K均值一样，轮廓系数极高，SSE极低，但是分类效果极差，因为它本身不适合这种数据的聚类。而欧几里得的轮廓系数是与其他两个相比较高的距离量度，且SSE与曼哈顿相同，故在模糊K均值中采用欧几里得距离更好。

5. 小结

K均值简单并可以用于各种数据类型。然而，K均值并不适合所有的数据类型。它不能处理非球形簇、不同尺寸和不同密度的簇，尽管指定足够大的簇个数K时它通常可以发现纯子簇。对包含离群点的数据进行聚类时，K均值也有问题。在这种情况下，离群点检测与删除大有帮助。最后，K均值仅限于具有质心概念的数据。

本次主要对比了K均值与模糊K均值在不同簇数K的情况下、不同距离量度下和不同初始质心选取的情况下的分类效果，采用了SSE和轮廓系数评价聚类效果。对于不同的数据类型应选取适合的距离量度，再根据数据大致分类缺点簇K的大小，初始质心选取可以采用与K-Means++中的方法，这样可以获得较稳定质量较高的聚类效果。

参考文献

[1] 数据库https://scikit-learn.org/stable/datasets/index.html#iris-plants-database
[2] Pang-Ning Tan, Michael Steinbach, Vipin Kumar, 数据挖掘导论,2011
[3] 各种距离度量https://blog.csdn.net/lj6052317/article/details/78770383

代码

from numpy import *
import numpy as np
import matplotlib.pyplot as plt
import operator
import random
from scipy.spatial.distance import pdist
from pylab import *
# 正常显示中文
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
INF = 9999999.0
def loadDataSet(fileName):'''加载数据集:param fileName::return:'''# 初始化一个空列表dataSet = []# 读取文件fr = open(fileName)# 循环遍历文件所有行for line in fr.readlines():# 切割每一行的数据curLine = line.strip(',').split(',')fltLine = list(map(float, curLine))    # 映射所有的元素为 float（浮点数）类型dataSet.append(fltLine)return dataSet
def distEclud(vecA, vecB):"""欧式距离输入：向量A, 向量B输出：两个向量的欧式距离"""a = (vecA[0, 0] - vecB[0, 0]) ** 2b = (vecA[0, 1] - vecB[0, 1]) ** 2return sqrt(a+b)
def distManhattan(vecA, vecB):'''曼哈顿距离输入：向量A, 向量B输出：两个向量的曼哈顿距离'''return np.sum(np.abs(vecA-vecB))
def distCosine(vecA, vecB):'''余弦距离输入：向量A, 向量B输出：两个向量的余弦距离'''return pdist(np.vstack([vecA, vecB]), 'cosine')
def distChebyshev(vecA, vecB):'''切比雪夫距离输入：向量A, 向量B输出：两个向量的切比雪夫距离'''return np.max(np.abs(vecA-vecB))
def randCent(dataSet, k):"""生成随机质心输入：数据集, 聚类个数输出：k个随机质心的矩阵"""n = dataSet.shape[1] #每个数据的维度centroids = mat(zeros((k, n))) # 生成k*n维数据for j in range(n):minJ = min(dataSet[:, j]) # 第j列最小值rangeJ = float(max(dataSet[:, j]) - minJ) # 第j列最大值与最小值的差值centroids[:, j] = minJ + rangeJ * np.random.rand(k, 1) # 最小值加上差值的（0，1）之间的倍数return centroids
def w_update(dataSet, centroids, clusterAssment, distMeans=distEclud):"""更新权重值输入：数据集, 质心输出：更新后的权值"""for j in range(dataSet.shape[0]):dist_all = 0for cen in range(len(centroids)):dist = distMeans(dataSet[j, :], centroids[cen, :])dist_all += distfor cen in range(len(centroids)):dist_self = distMeans(dataSet[j, :], centroids[cen, :])clusterAssment[j, cen+2] = dist_self/dist_allreturn clusterAssment
def kMeans(dataSet, k, distMeans=distEclud, createCent=randCent):"""输入：数据集, 聚类个数, 距离计算函数, 生成随机质心函数输出：质心矩阵, 簇分配和距离矩阵"""m = dataSet.shape[0]clusterAssment = mat(zeros((m, 2))) # 初始化聚类矩阵centroids = createCent(dataSet, k) # 生成随机质心clusterChanged = True #启动while clusterChanged:clusterChanged = Falsefor i in range(m): # 寻找最近的质心minDist = INFminIndex = -1for j in range(k):distJI = distMeans(centroids[j, :], dataSet[i, 0:2]) #距离方法,计算该点与质心的距离if distJI < minDist:   # 得到与质心距离最少的距离下标minDist = distJIminIndex = jif clusterAssment[i, 0] != minIndex: # 只要有一个点的簇发生变化，就继续进行clusterChanged = TrueclusterAssment[i, :] = minIndex, minDist #存储该点的簇类，该点到簇心的距离for cent in range(k): # 更新质心的位置ptsInClust = dataSet[nonzero(clusterAssment[:, 0].A == cent)[0]]centroids[cent, :] = mean(ptsInClust, axis=0) #计算均值return centroids, clusterAssment
def fuzzy_kmeans(dataSet, k, distMeans=distEclud, createCent=randCent):"""输入：数据集, 聚类个数, 距离计算函数, 生成随机质心函数输出：质心矩阵, 簇分配和距离矩阵"""m = dataSet.shape[0]clusterAssment = mat(zeros((m, k+2)))  # 初始化聚类矩阵centroids = createCent(dataSet, k)  # 生成随机质心clusterChanged = True  # 启动clusterAssment[:, 2:k+2] = 1/kwhile clusterChanged:clusterChanged = Falsefor i in range(m):  # 寻找最近的质心minDist = INFminIndex = -1for j in range(k):distJI = (clusterAssment[i, j+2]**2)*distMeans(centroids[j, :], dataSet[i, 0:2])  # 距离方法,计算该点与质心的距离if distJI < minDist:  # 得到与质心距离最少的距离下标minDist = distJIminIndex = jif clusterAssment[i, 0] != minIndex:  # 只要有一个点的簇发生变化，就继续进行clusterChanged = TrueclusterAssment[i, 0:2] = minIndex, minDist  # 存储该点的簇类，该点到簇心的距离# 更新权值clusterAssment = w_update(dataSet[i, :], centroids, clusterAssment)c = mat(zeros((3, k)))for cent in range(dataSet.shape[0]): #更新质心的位置c[2, int(clusterAssment[cent, 0])] += clusterAssment[cent, int(clusterAssment[cent, 0])+2]c[0, int(clusterAssment[cent, 0])] += clusterAssment[cent, int(clusterAssment[cent, 0])+2]*dataSet[cent, 0]c[1, int(clusterAssment[cent, 0])] += clusterAssment[cent, int(clusterAssment[cent, 0])+2]*dataSet[cent, 1]for cent in range(k):cc0 = c[0, cent]/c[2, cent]cc1 = c[1, cent]/c[2, cent]centroids[cent, :] = cc0, cc1  # 计算均值return centroids, clusterAssment
def plotFeature(dataSet, centroids, clusterAssment):m = shape(centroids)[0]fig = plt.figure()scatterMarkers = ['v', '^']scatterColors = ['blue', 'green', 'black', 'purple', 'orange', 'black', 'yellow']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)] # 选择颜色# flatten为返回一维数组ax.scatter(ptsInCurCluster[:, 0].flatten().A[0], ptsInCurCluster[:, 1].flatten().A[0], marker=markerStyle, c=colorSytle, s=30) # 添加每类数据点ax.scatter(centroids[:, 0].flatten().A[0], centroids[:, 1].flatten().A[0], marker='h', c='red', s=200) # 添加质心的位置
def silhouette(dataset, clustAssing, K, distMeans=distEclud):'''计算每个点的轮廓系数:param dataset:所有数据集:param clustAssing:输入为每个簇的标签和其到簇心的距离:return:返回平均轮廓系数，和数组:k:簇数量'''mean_s = 0clustAssing_new = mat(zeros((dataset.shape[0], 2)))for i in range(dataset.shape[0]):a = 0 ; n1 = 0buffer = []for k in range(K):buffer.append([])for j in range(dataset.shape[0]):if i == j :continueelif clustAssing[i, 0] == clustAssing[j, 0]: #计算同簇内的距离n1 += 1distJI = distMeans(dataset[j, 0:2], dataset[i, 0:2])  # 距离方法,计算该点与质心的距离a += distJIelif clustAssing[i, 0] != clustAssing[j, 0]: #计算不同簇内的距离distJI = distMeans(dataset[j, 0:2], dataset[i, 0:2])buffer[int(clustAssing[j, 0])].append(distJI)b_min = INFfor p in range(K):if p == clustAssing[i, 0]:continueelse:b_buff = 0for q in range(len(buffer[p])):b_buff += buffer[p][q]try:b_buff = np.round(b_buff / len(buffer[p]), 5)except Exception as e:print(e)if b_buff <= b_min:b_min = b_buffaa = np.round(a/n1, 5)s = (b_min - aa)/max(aa, b_min)clustAssing_new[i, :] = clustAssing[i, 0], smean_s += smean_s = mean_s/dataset.shape[0]return mean_s, clustAssing
def SSE(K, clustAssing):value = []SSE_mean = 0for k in range(K):value.append([])for i in range(clustAssing.shape[0]):value[int(clustAssing[i, 0])].append(clustAssing[i, 1])for k in range(K):SSE_mean += np.mean(value[k])SSE_mean /= Kreturn SSE_mean
def main():dataset = loadDataSet('Dataset.txt')dataSet = mat(dataset)resultCentroids, clustAssing = fuzzy_kmeans(dataSet, 5)#  resultCentroids, clustAssing = kMeans(dataSet, 5)plotFeature(dataSet, resultCentroids, clustAssing)mean_s, clustAssing_new = silhouette(dataSet, clustAssing, 5)mean_sse = SSE(5, clustAssing)print(mean_s, mean_sse)plt.show()#   plt.xlabel('四种不同距离量度')  # 将坐标系x轴命名为x1#   plt.ylabel('轮廓系数')  # 将坐标系y轴命名为y1#   plt.plot([1,2,3,4],[0.557,0.579,0.497,0.725])#   plt.show()#    performence = [[], []]#    for i in range(3, 11):#        resultCentroids, clustAssing = fuzzy_kmeans(dataSet, i)# #      print(i, resultCentroids)#        mean_s, clustAssing_new = silhouette(dataSet, clustAssing, i)#         # mean_s = SSE(i, clustAssing)#        # print('簇K%-2d的SSE为 %.5f'%(i, mean_s))#        print('簇K %-2d 的轮廓系数为 %.5f'%(i, mean_s))#        performence[0].append(i)#        performence[1].append(mean_s)#    plt.plot(performence[0], performence[1])#    plt.show()
if __name__ == '__main__':main()

如果觉得写得不错，请点赞、关注、评论。非常感谢~~