计算k-means的方差度量百分比？

Question

在维基百科页面中，介绍了一种使用肘部法（elbow method）确定k-means聚类数量的方法。 Scipy内置方法提供了实现，但我不确定他们所谓的畸变度量是如何计算的。

更精确地说，如果将集群解释的方差百分比根据集群数绘制成图形，则前几个集群会添加大量信息（解释大量方差），但在某个点上，边际收益将下降，在图形中形成一个角度。

假设我有以下数据点及其关联的质心，有什么好的方法来计算此度量？

points = numpy.array([[ 0,  0],
       [ 0,  1],
       [ 0, -1],
       [ 1,  0],
       [-1,  0],
       [ 9,  9],
       [ 9, 10],
       [ 9,  8],
       [10,  9],
       [10,  8]])

kmeans(pp,2)
(array([[9, 8],
   [0, 0]]), 0.9414213562373096)

我想计算仅给定点和质心的 0.94.. 测量值。不确定是否可以使用Scipy内置方法中的任何方法，还是必须编写自己的方法。对于大量数据，如何高效地进行计算？

简而言之，我的问题（都相关）如下：

• 在给定距离矩阵和将点分配到哪个簇的映射的情况下，如何计算可用于绘制肘部图的度量标准？
• 如果使用不同的距离函数，例如余弦相似性，方法会如何改变？

编辑2：畸变程度（Distortion）

from scipy.spatial.distance import cdist
D = cdist(points, centroids, 'euclidean')
sum(numpy.min(D, axis=1))

第一组点的输出是准确的。然而，当我尝试使用另一组时：

>>> pp = numpy.array([[1,2], [2,1], [2,2], [1,3], [6,7], [6,5], [7,8], [8,8]])
>>> kmeans(pp, 2)
(array([[6, 7],
       [1, 2]]), 1.1330618877807475)
>>> centroids = numpy.array([[6,7], [1,2]])
>>> D = cdist(points, centroids, 'euclidean')
>>> sum(numpy.min(D, axis=1))
9.0644951022459797

我猜测最后一个值不匹配的原因是kmeans似乎将该值除以数据集中的总点数。

编辑1：百分比方差

到目前为止，我的代码（应添加到Denis的K-means实现中）：

centres, xtoc, dist = kmeanssample( points, 2, nsample=2,
        delta=kmdelta, maxiter=kmiter, metric=metric, verbose=0 )

print "Unique clusters: ", set(xtoc)
print ""
cluster_vars = []
for cluster in set(xtoc):
    print "Cluster: ", cluster

    truthcondition = ([x == cluster for x in xtoc])
    distances_inside_cluster = (truthcondition * dist)

    indices = [i for i,x in enumerate(truthcondition) if x == True]
    final_distances = [distances_inside_cluster[k] for k in indices]

    print final_distances
    print np.array(final_distances).var()
    cluster_vars.append(np.array(final_distances).var())
    print ""

print "Sum of variances: ", sum(cluster_vars)
print "Total Variance: ", points.var()
print "Percent: ", (100 * sum(cluster_vars) / points.var())

以下是 k=2 时的输出结果：

Unique clusters:  set([0, 1])

Cluster:  0
[1.0, 2.0, 0.0, 1.4142135623730951, 1.0]
0.427451660041

Cluster:  1
[0.0, 1.0, 1.0, 1.0, 1.0]
0.16

Sum of variances:  0.587451660041
Total Variance:  21.1475
Percent:  2.77787757437

在我的真实数据集上（我觉得不对！）：

Sum of variances:  0.0188124746402
Total Variance:  0.00313754329764
Percent:  599.592510943
Unique clusters:  set([0, 1, 2, 3])

Sum of variances:  0.0255808508714
Total Variance:  0.00313754329764
Percent:  815.314672809
Unique clusters:  set([0, 1, 2, 3, 4])

Sum of variances:  0.0588210052519
Total Variance:  0.00313754329764
Percent:  1874.74720416
Unique clusters:  set([0, 1, 2, 3, 4, 5])

Sum of variances:  0.0672406353655
Total Variance:  0.00313754329764
Percent:  2143.09824556
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6])

Sum of variances:  0.0646291452839
Total Variance:  0.00313754329764
Percent:  2059.86465055
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7])

Sum of variances:  0.0817517362176
Total Variance:  0.00313754329764
Percent:  2605.5970695
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7, 8])

Sum of variances:  0.0912820650486
Total Variance:  0.00313754329764
Percent:  2909.34837831
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Sum of variances:  0.102119601368
Total Variance:  0.00313754329764
Percent:  3254.76309585
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

Sum of variances:  0.125549475536
Total Variance:  0.00313754329764
Percent:  4001.52168834
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

Sum of variances:  0.138469402779
Total Variance:  0.00313754329764
Percent:  4413.30651542
Unique clusters:  set([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])

- Legend

Answer 1

就Kmeans而言，失真度被用作停止准则（如果两次迭代之间的变化小于某个阈值，则假定收敛）

如果您想从一组点和质心中计算它，可以执行以下操作（代码使用 pdist2 函数在MATLAB中编写，但应该很容易重写为Python / Numpy / Scipy）：

% data
X = [0 1 ; 0 -1 ; 1 0 ; -1 0 ; 9 9 ; 9 10 ; 9 8 ; 10 9 ; 10 8];

% centroids
C = [9 8 ; 0 0];

% euclidean distance from each point to each cluster centroid
D = pdist2(X, C, 'euclidean');

% find closest centroid to each point, and the corresponding distance
[distortions,idx] = min(D,[],2);

结果：

% total distortion
>> sum(distortions)
ans =
           9.4142135623731

---

编辑#1:

我有一些时间来尝试这个.. 这是应用于'费舍尔鸢尾花数据集'（4个特征，150个实例）的KMeans聚类的示例。我们迭代k = 1..10，绘制弯曲曲线，选择K = 3 作为簇数，并显示结果的散点图。

请注意，我包括了许多计算簇内方差（扭曲度）的方法，给定点和质心。 scipy.cluster.vq.kmeans函数默认返回此度量（使用欧几里得距离度量计算）。您还可以使用scipy.spatial.distance.cdist函数使用您选择的函数计算距离（前提是使用相同的距离度量获取了聚类质心：@Denis有一个解决方案），然后从中计算扭曲度。

import numpy as np
from scipy.cluster.vq import kmeans,vq
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

# load the iris dataset
fName = 'C:\\Python27\\Lib\\site-packages\\scipy\\spatial\\tests\\data\\iris.txt'
fp = open(fName)
X = np.loadtxt(fp)
fp.close()

##### cluster data into K=1..10 clusters #####
K = range(1,10)

# scipy.cluster.vq.kmeans
KM = [kmeans(X,k) for k in K]
centroids = [cent for (cent,var) in KM]   # cluster centroids
#avgWithinSS = [var for (cent,var) in KM] # mean within-cluster sum of squares

# alternative: scipy.cluster.vq.vq
#Z = [vq(X,cent) for cent in centroids]
#avgWithinSS = [sum(dist)/X.shape[0] for (cIdx,dist) in Z]

# alternative: scipy.spatial.distance.cdist
D_k = [cdist(X, cent, 'euclidean') for cent in centroids]
cIdx = [np.argmin(D,axis=1) for D in D_k]
dist = [np.min(D,axis=1) for D in D_k]
avgWithinSS = [sum(d)/X.shape[0] for d in dist]

##### plot ###
kIdx = 2

# elbow curve
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(K, avgWithinSS, 'b*-')
ax.plot(K[kIdx], avgWithinSS[kIdx], marker='o', markersize=12, 
    markeredgewidth=2, markeredgecolor='r', markerfacecolor='None')
plt.grid(True)
plt.xlabel('Number of clusters')
plt.ylabel('Average within-cluster sum of squares')
plt.title('Elbow for KMeans clustering')

# scatter plot
fig = plt.figure()
ax = fig.add_subplot(111)
#ax.scatter(X[:,2],X[:,1], s=30, c=cIdx[k])
clr = ['b','g','r','c','m','y','k']
for i in range(K[kIdx]):
    ind = (cIdx[kIdx]==i)
    ax.scatter(X[ind,2],X[ind,1], s=30, c=clr[i], label='Cluster %d'%i)
plt.xlabel('Petal Length')
plt.ylabel('Sepal Width')
plt.title('Iris Dataset, KMeans clustering with K=%d' % K[kIdx])
plt.legend()

plt.show()

---

编辑#2:

针对评论，我在下面提供了另一个完整的例子，使用NIST手写数字数据集：它包含1797个数字图像，从0到9，每个图像大小为8×8像素。我稍微修改了上面的实验重复实验：应用主成分分析将维度从64降至2：

import numpy as np
from scipy.cluster.vq import kmeans
from scipy.spatial.distance import cdist,pdist
from sklearn import datasets
from sklearn.decomposition import RandomizedPCA
from matplotlib import pyplot as plt
from matplotlib import cm

##### data #####
# load digits dataset
data = datasets.load_digits()
t = data['target']

# perform PCA dimensionality reduction
pca = RandomizedPCA(n_components=2).fit(data['data'])
X = pca.transform(data['data'])

##### cluster data into K=1..20 clusters #####
K_MAX = 20
KK = range(1,K_MAX+1)

KM = [kmeans(X,k) for k in KK]
centroids = [cent for (cent,var) in KM]
D_k = [cdist(X, cent, 'euclidean') for cent in centroids]
cIdx = [np.argmin(D,axis=1) for D in D_k]
dist = [np.min(D,axis=1) for D in D_k]

tot_withinss = [sum(d**2) for d in dist]  # Total within-cluster sum of squares
totss = sum(pdist(X)**2)/X.shape[0]       # The total sum of squares
betweenss = totss - tot_withinss          # The between-cluster sum of squares

##### plots #####
kIdx = 9        # K=10
clr = cm.spectral( np.linspace(0,1,10) ).tolist()
mrk = 'os^p<dvh8>+x.'

# elbow curve
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(KK, betweenss/totss*100, 'b*-')
ax.plot(KK[kIdx], betweenss[kIdx]/totss*100, marker='o', markersize=12, 
    markeredgewidth=2, markeredgecolor='r', markerfacecolor='None')
ax.set_ylim((0,100))
plt.grid(True)
plt.xlabel('Number of clusters')
plt.ylabel('Percentage of variance explained (%)')
plt.title('Elbow for KMeans clustering')

# show centroids for K=10 clusters
plt.figure()
for i in range(kIdx+1):
    img = pca.inverse_transform(centroids[kIdx][i]).reshape(8,8)
    ax = plt.subplot(3,4,i+1)
    ax.set_xticks([])
    ax.set_yticks([])
    plt.imshow(img, cmap=cm.gray)
    plt.title( 'Cluster %d' % i )

# compare K=10 clustering vs. actual digits (PCA projections)
fig = plt.figure()
ax = fig.add_subplot(121)
for i in range(10):
    ind = (t==i)
    ax.scatter(X[ind,0],X[ind,1], s=35, c=clr[i], marker=mrk[i], label='%d'%i)
plt.legend()
plt.title('Actual Digits')
ax = fig.add_subplot(122)
for i in range(kIdx+1):
    ind = (cIdx[kIdx]==i)
    ax.scatter(X[ind,0],X[ind,1], s=35, c=clr[i], marker=mrk[i], label='C%d'%i)
plt.legend()
plt.title('K=%d clusters'%KK[kIdx])

plt.show()

你可以看到，有些聚类实际上对应于可区分的数字，而其他聚类则不匹配任何一个数字。

注意：scikit-learn 中包含 K-means 的实现（以及许多其他聚类算法和各种聚类度量）。这里还有另一个类似的例子。

- Amro

---

+1 谢谢您的解释。根据您所提到的，我现在唯一需要确认的是这个畸变值是否用于评估 k 的价值。在这篇文章中：http://stats.stackexchange.com/questions/9850/how-to-plot-data-output-of-clustering 作者直接使用了畸变值，但我并不真正理解他为什么这样做。您对此有什么想法吗？ - Legend

是的，在最小化簇内平方和（这里称为失真度）和最小化簇数之间存在权衡。换句话说，我们希望模型能够很好地拟合数据（小失真度），但同时，我们也希望模型尽可能简单（不要过于复杂，有太多的簇）。肘部法则是一种简单的启发式方法，可以在两者之间取得平衡。这个答案也很好地解释了它：https://dev59.com/J3I-5IYBdhLWcg3wj5JG#1793572 - Amro

Amro，不错。然而Iris数据集很小，从中推断可能会有问题。在scikits.learn的1797 x 64数字数据上运行kmeans算法，应该可以得到10个明显分离的簇 :) 当k = 7 .. 13时，我得到的平均距离点 - 簇中心为27.7 26.2 25.3 26.2 24.6 24.5 24.1。拐点在10吗？ - denis

@Denis：我增加了一个手写数字数据集的示例。 - Amro

Amro，又见面了，+1（添加到scikits.learn示例中？）这清楚地表明k = 10并不是一个很好的拐点，自动拐点有时会出现问题。 - denis

@Denis：肘部法是一种启发式方法，远非完美。还有其他方法，如AIC/BIC...此外，您必须记住Kmeans是一种无监督学习技术，这意味着它不知道数据的实际类别。相反，它试图从数据本身自然地发现聚类。因此，如果两个数字在特征空间中看起来相似，则可能会像您在上面的示例中看到的那样被分组在一起。此外，通过使用PCA，我们为了减少维数而失去了一些信息...正如您现在可能已经发现的那样，聚类是一项困难的任务 :) - Amro

Answer 2

一个简单的聚类度量方法： 1）从每个点到它最近的簇中心绘制“日光”射线， 2）查看所有射线的长度——距离（点，中心，度量=...）。

对于metric="sqeuclidean"和1个簇，平均长度平方是总方差X.var()；对于2个簇，它会减少...直到N个簇，长度都为0。"解释的方差百分比"是100%减去这个平均值。

此代码在is-it-possible-to-specify-your-own-distance-function-using-scikits-learn-k-means下可用。

def distancestocentres( X, centres, metric="euclidean", p=2 ):
    """ all distances X -> nearest centre, any metric
            euclidean2 (~ withinss) is more sensitive to outliers,
            cityblock (manhattan, L1) less sensitive
    """
    D = cdist( X, centres, metric=metric, p=p )  # |X| x |centres|
    return D.min(axis=1)  # all the distances

像任何长列表一样，这些距离可以以不同的方式查看：np.mean()，np.histogram() ...绘图和可视化并不容易。
另请参见stats.stackexchange.com/questions/tagged/clustering，特别是
如何确定数据是否“足够聚类”，使聚类算法产生有意义的结果？

- denis

---

+1 感谢您的时间和解释！我尝试编写了您在帖子中解释的内容，并将其添加到了我的问题末尾。如果您有空，能否请看一下？ - Legend

当然，足够好。真正的问题是，对于您的实际数据，这如何随着k的变化而变化--请给出数字？如果k = 5和6接近，请继续前进。 - denis

我猜测我的函数出了问题。我已经在EDIT 1下面的问题中发布了观察到的值。我得到的百分比超过了100％，甚至达到了千分之几。我现在确定我的实现是错误的。 - Legend