40

在Wikipedia 页面上，描述了一种肘部方法，用于确定 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：失真

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])
4

2 回答 2

55

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

如果要从一组点和质心计算它，可以执行以下操作（代码在 MATLAB 中使用pdist2函数，但在 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：

我有一些时间来解决这个问题。这是一个应用在“Fisher Iris 数据集”上的 KMeans 聚类示例（4 个特征，150 个实例）。我们迭代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()

弯头曲线 数字质心 PCA_比较

您可以看到一些集群实际上如何对应于可区分的数字，而另一些则不匹配单个数字。

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

于 2011-07-11T21:51:44.530 回答
6

一个简单的集群测量：
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，特别是
如何判断数据是否“聚集”到足以让聚类算法产生有意义的结果？

于 2011-07-11T11:39:28.040 回答