好的聚类:类内凝聚度高,类间分离度高。
本文介绍两种聚类评估方法,轮廓系数(Silhouette Coefficient)以及标准化互信息(NMI),并且用Python实现。
导航
- 效果评估综述
- 轮廓系数
- 互信息
- 参考文章
效果评估综述
这里直接贴上 聚类算法初探(七)聚类分析的效果评测
它摘自于中国科学院计算技术研究所周昭涛的硕士论文《文本聚类分析效果评价及文本表示研究》的第三章。建议先看看原文,可以对聚类评估有一个很好的了解。
综合来说,我们希望最终的聚类结果是:同一个簇内的点是紧密的,而不同簇之间的距离是较远的;同时,它也要与我们人工的判断相一致。
类内距 和 类间距 ;而后者是需要对于数据进行 标注 ,判断聚类的 准确性 。
轮廓系数
类别未知 时,可以选择轮廓系数作为聚类性能的评估指标。从 百度百科 摘取的轮廓系数的介绍:
轮廓系数是聚类效果好坏的一种评价方式。它结合内聚度和分离度两种因素。可以用来在相同原始数据的基础上用来评价不同算法、或者算法不同运行方式对聚类结果所产生的影响。
的值是需要事先给定的。如果我们不太清楚数据分几个簇比较好,我们可以尝试用不同的 值来计算轮廓系数,从而选择 最优 的一个
- 对于一个簇中的一个点
- ( ) = (
- ( ) = (
- 那么 的轮廓系数为
- 将所有点的轮廓系数求平均,就是该聚类结果总的轮廓系数。
( ) 可以理解为一点与簇内其他点不相似性的平均值,即类内的 凝聚度 ; ( ) 可以理解为一点与最近簇的不相似性的平均值,即类间的 分离度 。
是介于[-1, 1]
之间的。 ( ) 越小, ( ) 越大,说明聚类效果越好,轮廓系数 就越接近1
。
那么,如果用Python来实现呢?给出它的伪代码:
输入:数据的个数m,k的取值范围
输出:每个k所对应的轮廓系数
在取值范围内的每一个k:
对于每一个簇:
对于每一个点i:
计算i与簇内其他点距离之和的平均值,记为a(i)
计算i与其他簇每个点距离之和的平均值,选择最小的那一个,记为b(i)
s(i) = (b(i)-a(i)) / max(a(i), b(i))
加总每一个点的轮廓系数,记为 sum
轮廓系数 s = sum / m
接下来用代码实现,K-Means 具体细节可以参考 利用K-Means聚类算法对未标注数据分组。
因此,运行K-Means的代码我就不贴啦,我就直接附上轮廓系数的代码。代码中的一些细节在上面的博客中得到了解释,有什么疑问可以先看看那篇博文哦。
首先,我的数据集为testSet.txt
,运行 K-Means 的代码为biKMeans(data, k)
, 为传入的数据集矩阵, 是簇的个数。这个函数返回两个值,一个是中心点,另一个是clusterAssment
,它一共有 行(
if __name__ == '__main__':
data = np.mat(loadDataSet('testSet.txt'))
m = np.shape(data)[0] # 一共有m行数据
for k in range(2, 10): # 簇的个数取值在2到9之间
clusterAssment = biKMeans(data, k)[1] # 进行二分类,返回保存簇索引的矩阵
s_sum = 0 # 所有簇的s值
cluster_s = 0 # 一个簇所有点的的s值
for cent in range(k): # 对于每一个簇
category = np.nonzero(clusterAssment[:, 0] == cent)[0] # 得到簇索引为cent的值的位置,形式类似为[1, 4 ,6]
clusterNum = len(category) # 该簇中点的个数
s = 0
for index, lineNum in enumerate(category): # 对于簇中的每一个点,index为索引,lineNum为点所在的行数
# 计算该点到簇内其他点的距离之和的平均值
innerSum = 0
for i in range(clusterNum):
if i == index:
continue # 若为当前该点,则跳出本次循环
dis = distCal(data[category[i]], data[lineNum]) # 若二者为不同点,计算二者之间的距离
innerSum += dis # 将之保存到内部距离
a = innerSum / clusterNum
# 计算该点到其他簇所有点距离之和的最小平均值
minDis = np.inf # 设定初始最小值为无穷大
for other_cent in range(k): # 对于每一个簇
if other_cent != cent: # 如果和上面给定的簇不一样
other_category = np.nonzero(clusterAssment[:, 0] == other_cent)[0] # 得到簇里面点对应的行数
other_clusterNum = len(other_category) # 该簇中点的个数
other_sum = 0
for other_lineNum in other_category: # 对于簇中的每一个点
other_dis = distCal(data[other_lineNum], data[lineNum])
other_sum += other_dis
other_sum = other_sum / other_clusterNum # 求平均
if other_sum < minDis:
minDis = other_sum # 如果一个点距离另外一个簇所有点的距离小于当前最小值,则更新
b = minDis
s += (b-a) / max(a, b) # 每一个点的轮廓系数
cluster_s += s # 每一个簇的s值
s_sum = cluster_s / m # 取平均
print("当前k的值为:%d" % k)
print("轮廓系数为:%s" % str(s_sum))
print('***' * 20)
取值范围在[2, 10)
,最终打印的结果如下:
可以看到,当
不应该 用来评估 不同聚类算法 之间的优劣,比如Kmeans聚类结果与DBSCAN聚类结果之间的比较。
画图的代码是用sklearn做的,请看原文博主:sklearn之聚类评估指标—轮廓系数。
还是贴出代码吧!sklearn有现成的函数可以调用,就不用写轮廓系数的代码啦!很方便的!
import numpy as np
import matplotlib.pyplot as mp
import sklearn.cluster as sc
import sklearn.metrics as sm
# 读取数据,绘制图像
x = np.loadtxt('testSet.txt', unpack=False, dtype='f8', delimiter='\t')
print(x.shape)
# 基于Kmeans完成聚类
model = sc.KMeans(n_clusters=4)
model.fit(x) # 完成聚类
pred_y = model.predict(x) # 预测点在哪个聚类中
print(pred_y) # 输出每个样本的聚类标签
# 打印轮廓系数
print(sm.silhouette_score(x, pred_y, sample_size=len(x), metric='euclidean'))
# 获取聚类中心
centers = model.cluster_centers_
print(centers)
# 绘制分类边界线
l, r = x[:, 0].min() - 1, x[:, 0].max() + 1
b, t = x[:, 1].min() - 1, x[:, 1].max() + 1
n = 500
grid_x, grid_y = np.meshgrid(np.linspace(l, r, n), np.linspace(b, t, n))
bg_x = np.column_stack((grid_x.ravel(), grid_y.ravel()))
bg_y = model.predict(bg_x)
grid_z = bg_y.reshape(grid_x.shape)
# 画图显示样本数据
mp.figure('Kmeans', facecolor='lightgray')
mp.title('Kmeans', fontsize=16)
mp.xlabel('X', fontsize=14)
mp.ylabel('Y', fontsize=14)
mp.tick_params(labelsize=10)
mp.pcolormesh(grid_x, grid_y, grid_z, cmap='gray')
mp.scatter(x[:, 0], x[:, 1], s=80, c=pred_y, cmap='brg', label='Samples')
mp.scatter(centers[:, 0], centers[:, 1], s=300, color='red', marker='+', label='cluster center')
mp.legend()
mp.show()
互信息
熵 的知识,建议看这位博主写的 信息论基础 。
这里仍然参考了百度百科 互信息 。
互信息 (Mutual Information) 是信息论里一种有用的信息度量,它可以看成是一个随机变量中包含的关于另一个随机变量的信息量,或者说是一个随机变量由于已知另一个随机变量而减少的不肯定性。
时,由于通道往往会有噪音和干扰,使得信宿往往接受到的是变形的 。 收到后推测由 发出的概率,称为后验概率;而信源发出 的概率,我们称为先验概率。定义 的后验概率与先验概率比值的对数为 对
我们看如何计算,首先是熵的计算公式:
由熵的连锁规则:
所以:
是互信息。
可以这样理解: 是 变量的不确定性, 是知道了 之后, 的不确定性, 就是由于知道 而使得 不确定性减少 的量(即 透露了多少 的信息?),这部分就是 互信息 。
这个就是互信息的计算公式。那么如何运用到聚类中呢?举一个例子,这里参考了 标准化互信息NMI计算步骤及其Python实现 。
是经过某种算法聚类后,得到的结果;
、 都将数据聚成了两类:,,由
首先计算联合概率分布:
再分别计算各自的概率分布:
,
,
要对计算的值进行归一化,公式为:
其中,
是经过某种算法聚类后,得到的结果;
import math
import numpy as np
from sklearn import metrics
def NMI(A, B):
# 样本点数
total = len(A) # 数据及的个数
A_ids = set(A) # 创建一个无序不重复的集合,保存A簇的个数
B_ids = set(B) # 保存B簇的个数
# 互信息计算
MI = 0
eps = 1.4e-45 # 防止出现log 0的情况
for idA in A_ids: # 对于A中的每一个簇
for idB in B_ids: # 对于B中的每一个簇
idAOccur = np.where(A == idA) # np.where(condition),则输出满足条件 (即非0) 元素的坐标
idBOccur = np.where(B == idB) # 一个簇的元素所对应的位置
idABOccur = np.intersect1d(idAOccur, idBOccur) # 找到A/B位置的交集
px = 1.0 * len(idAOccur[0])/total
py = 1.0 * len(idBOccur[0])/total
pxy = 1.0 * len(idABOccur)/total
MI = MI + pxy * math.log(pxy/(px * py) + eps, 2)
# 标准化互信息
Hx = 0
for idA in A_ids:
idAOccurCount = 1.0 * len(np.where(A == idA)[0])
Hx = Hx - (idAOccurCount/total)*math.log(idAOccurCount/total+eps, 2)
Hy = 0
for idB in B_ids:
idBOccurCount = 1.0 * len(np.where(B == idB)[0])
Hy = Hy - (idBOccurCount/total)*math.log(idBOccurCount/total+eps, 2)
MIhat = 2.0 * MI / (Hx + Hy)
return MIhat
if __name__ == '__main__':
A = np.array([1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3])
B = np.array([1,2,1,1,1,1,1,2,2,2,2,3,1,1,3,3,3])
print(NMI(A, B))
print(metrics.normalized_mutual_info_score(A, B)) # 直接调用函数计算的互信息
结果如下: