当使用k-means聚类时，我该如何确定k值？

Question

我一直在学习k-means聚类，但有一件事不太清楚，那就是如何选择k的值。这仅仅是试错还是还有其他因素呢？

- Jason Baker

---

啊啊...那确实是关于k均值算法的核心问题。 - mjv

你能分享函数L（对数似然）的代码吗？给定一个中心点X，Y和点(x(i=1,2,3,4,...,n),y(i=1,2,3,4,..,n))，我该如何得到L？ - user653773

这是一个关于如何确定数据集中聚类数量的维基百科文章链接：http://en.wikipedia.org/wiki/Determining_the_number_of_clusters_in_a_data_set - Amro

我曾在这里使用过半打种方法（使用 R）回答了一个类似的问题：stackoverflow.com/a/15376462/1036500 - Ben

Answer 1

您可以最大化贝叶斯信息准则（BIC）：

BIC(C | X) = L(X | C) - (p / 2) * log n

L(X | C)是模型C对数据集X的对数似然，p是模型C中的参数数量，n是数据集中的点数。 请参阅Dan Pelleg和Andrew Moore在ICML 2000中的文章"X-means: extending K-means with efficient estimation of the number of clusters"。

另一种方法是从较大的k值开始，并继续删除聚类中心（减小k），直到不再减少描述长度为止。请参见Horst Bischof、Ales Leonardis和Alexander Selb在1999年发表于Pattern Analysis and Applications vol.2, p.59-72的文章"MDL principle for robust vector quantisation"。

最后，您可以从一个聚类开始，然后继续分裂聚类，直到分配给每个聚类的点具有高斯分布。在Greg Hamerly和Charles Elkan在NIPS 2003中的文章"Learning the k in k-means"中，他们展示了一些证据表明这比BIC更有效，并且BIC没有足够强烈地惩罚模型的复杂性。

- Vebjorn Ljosa

---

非常好的回答！对于X-Means，您是否知道整体BIC得分n := k*2（k个聚类，每个聚类由具有均值/方差参数的高斯模型建模）。此外，如果确定“父”BIC > “2个子”BIC，您是否会在下一次迭代中再次拆分该聚类？ - Budric

@Budric，这些可能应该是单独的问题，也许应该在stats.stackexchange.com上提问。 - Vebjorn Ljosa

Answer 2

基本上，您想要找到两个变量之间的平衡：聚类数(k)和聚类的平均方差。您希望最小化前者同时也最小化后者。当然，随着聚类数的增加，聚类的平均方差会减少(达到k=n和方差=0的平凡情况)。

像数据分析中一样，没有一种真正有效的方法适用于所有情况。最终，您必须使用自己的最佳判断力。为此，将聚类数与平均方差绘制成图表可能有所帮助(假设您已经对多个k值运行了算法)。然后，您可以使用拐点处的聚类数。

- Jan Krüger

Answer 3

是的，您可以使用肘部法找到最佳聚类数，但我发现使用脚本从肘图中找到聚类值很麻烦。您可以观察肘图并自行找到拐点，但从脚本中找到它需要很多工作。

因此，另一个选择是使用Silhouette Method来找到它。在R中，Silhouette的结果与Elbow方法完全一致。

以下是我的操作步骤。

#Dataset for Clustering
n = 150
g = 6 
set.seed(g)
d <- data.frame(x = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))), 
                y = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))))
mydata<-d
#Plot 3X2 plots
attach(mtcars)
par(mfrow=c(3,2))

#Plot the original dataset
plot(mydata$x,mydata$y,main="Original Dataset")

#Scree plot to deterine the number of clusters
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
  for (i in 2:15) {
    wss[i] <- sum(kmeans(mydata,centers=i)$withinss)
}   
plot(1:15, wss, type="b", xlab="Number of Clusters",ylab="Within groups sum of squares")

# Ward Hierarchical Clustering
d <- dist(mydata, method = "euclidean") # distance matrix
fit <- hclust(d, method="ward") 
plot(fit) # display dendogram
groups <- cutree(fit, k=5) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters 
rect.hclust(fit, k=5, border="red")

#Silhouette analysis for determining the number of clusters
library(fpc)
asw <- numeric(20)
for (k in 2:20)
  asw[[k]] <- pam(mydata, k) $ silinfo $ avg.width
k.best <- which.max(asw)

cat("silhouette-optimal number of clusters:", k.best, "\n")
plot(pam(d, k.best))

# K-Means Cluster Analysis
fit <- kmeans(mydata,k.best)
mydata 
# get cluster means 
aggregate(mydata,by=list(fit$cluster),FUN=mean)
# append cluster assignment
mydata <- data.frame(mydata, clusterid=fit$cluster)
plot(mydata$x,mydata$y, col = fit$cluster, main="K-means Clustering results")

希望这能有所帮助！

- Udeep Shakya

---

只需为Python用户添加一个链接到Silhouette分析教程 http://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html#sphx-glr-auto-examples-cluster-plot-kmeans-silhouette-analysis-py - Chaitanya Shivade

此外，关于绘图，请参见黄砖 https://www.scikit-yb.org/en/latest/api/cluster/silhouette.html ，他们还有肘部法则。 - A_Arnold

Answer 4

也许像我这样的初学者正在寻找代码示例。有关silhouette_score的信息可在此处找到。

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

range_n_clusters = [2, 3, 4]            # clusters range you want to select
dataToFit = [[12,23],[112,46],[45,23]]  # sample data
best_clusters = 0                       # best cluster number which you will get
previous_silh_avg = 0.0

for n_clusters in range_n_clusters:
    clusterer = KMeans(n_clusters=n_clusters)
    cluster_labels = clusterer.fit_predict(dataToFit)
    silhouette_avg = silhouette_score(dataToFit, cluster_labels)
    if silhouette_avg > previous_silh_avg:
        previous_silh_avg = silhouette_avg
        best_clusters = n_clusters

# Final Kmeans for best_clusters
kmeans = KMeans(n_clusters=best_clusters, random_state=0).fit(dataToFit)

- bhargav patel

---

示例在scikit-learn版本0.24.2上无法运行。在silhouette_score（dataToFit，cluster_labels）上出现错误。 “发生异常：ValueError 标签数量为3。有效值为2到n_samples-1（含）。” - Alexander Zhukov

你可能需要查看这个链接： https://dev59.com/rlUK5IYBdhLWcg3wuB54 - bhargav patel

Answer 5

请看Greg Hamerly和Charles Elkan的论文“学习k-means中的k”，链接在这里。该论文使用高斯测试来确定正确的聚类数量。此外，作者声称这种方法比接受答案中提到的BIC更好。

- Autonomous

Answer 6

有一种叫做经验法则的东西。它说，聚类数可以通过以下公式计算得出：

k = (n/2)^0.5

其中 n 是你的样本中的总元素数。你可以在以下论文中查证此信息的真实性：

http://www.ijarcsms.com/docs/paper/volume1/issue6/V1I6-0015.pdf

还有另一种方法叫做 G-means，它适用于高斯分布或正态分布。

这种方法需要进行大量统计，但是可以完成。

以下是来源:

http://papers.nips.cc/paper/2526-learning-the-k-in-k-means.pdf

希望对你有所帮助！

- Arthur Busqueiro

Answer 7

如果您不知道要提供给k-means作为参数的聚类数k，那么有四种自动查找它的方法：

• G-means算法：该算法利用统计测试来决定是否将k-means中心点拆分成两个，从而自动发现聚类数。这种层次化的方法基于对数据子集是否遵循高斯分布（连续函数，近似精确二项分布）的假设进行统计检验，如果不是，则将其拆分为簇。它从少量中心点开始，例如仅一个簇（k=1），然后将其拆分为两个中心（k=2），再将这两个中心各自拆分（k=4），总共得到四个中心。如果G-means不接受这四个中心，则答案是上一步：在这种情况下是两个中心（k=2）。当您无法估计聚类实例数时，G-means非常有用。请注意，不恰当的“k”参数选择可能会导致错误结果。G-means的并行版本称为p-means。G-means源代码：来源1, 来源2, 来源3。

• x-means：一种新算法，它高效地搜索聚类位置和聚类数量以优化贝叶斯信息准则（BIC）或Akaike信息准则（AIC）度量。这个k-means版本不仅发现k的数量，而且还加速了k-means的执行。

• 在线K-means或流式K-means：它允许仅通过扫描整个数据一次来执行K-means，并自动找到最优的k值。Spark实现了该算法。

• MeanShift算法：它是一种非参数聚类技术，不需要先验知识，也不限制簇的形状。均值漂移聚类旨在发现光滑样本密度中的“blob”。它是一种基于质心的算法，通过更新候选中心点为给定区域内点的平均值来工作。然后，在后处理阶段过滤这些候选中心点以消除近似重复项，形成最终的质心集。来源：来源1，来源2，来源3

- user9562553

Answer 8

我很惊讶没有人提到这篇出色的文章： http://www.ee.columbia.edu/~dpwe/papers/PhamDN05-kmeans.pdf

在尝试了几个其他建议之后，我最终在阅读这篇博客时找到了这篇文章： https://datasciencelab.wordpress.com/2014/01/21/selection-of-k-in-k-means-clustering-reloaded/

之后，我使用Scala实现了这个算法，对于我的用例来说，这个实现提供了真正优秀的结果。以下是代码：

import breeze.linalg.DenseVector
import Kmeans.{Features, _}
import nak.cluster.{Kmeans => NakKmeans}

import scala.collection.immutable.IndexedSeq
import scala.collection.mutable.ListBuffer

/*
https://datasciencelab.wordpress.com/2014/01/21/selection-of-k-in-k-means-clustering-reloaded/
 */
class Kmeans(features: Features) {
  def fkAlphaDispersionCentroids(k: Int, dispersionOfKMinus1: Double = 0d, alphaOfKMinus1: Double = 1d): (Double, Double, Double, Features) = {
    if (1 == k || 0d == dispersionOfKMinus1) (1d, 1d, 1d, Vector.empty)
    else {
      val featureDimensions = features.headOption.map(_.size).getOrElse(1)
      val (dispersion, centroids: Features) = new NakKmeans[DenseVector[Double]](features).run(k)
      val alpha =
        if (2 == k) 1d - 3d / (4d * featureDimensions)
        else alphaOfKMinus1 + (1d - alphaOfKMinus1) / 6d
      val fk = dispersion / (alpha * dispersionOfKMinus1)
      (fk, alpha, dispersion, centroids)
    }
  }

  def fks(maxK: Int = maxK): List[(Double, Double, Double, Features)] = {
    val fadcs = ListBuffer[(Double, Double, Double, Features)](fkAlphaDispersionCentroids(1))
    var k = 2
    while (k <= maxK) {
      val (fk, alpha, dispersion, features) = fadcs(k - 2)
      fadcs += fkAlphaDispersionCentroids(k, dispersion, alpha)
      k += 1
    }
    fadcs.toList
  }

  def detK: (Double, Features) = {
    val vals = fks().minBy(_._1)
    (vals._3, vals._4)
  }
}

object Kmeans {
  val maxK = 10
  type Features = IndexedSeq[DenseVector[Double]]
}

- eirirlar

---

使用Scala 2.11.7实现，使用Breeze 0.12和Nak 1.3。 - eirirlar

嗨@eirirlar，我正在尝试使用Python实现相同的代码-但我无法按照网站上的代码。请查看我的帖子：http://stackoverflow.com/questions/36729826/python-k-means-clustering - piccolo

@ImranRashid 抱歉，我只测试了二维，并且我不是 Python 专家。 - eirirlar

Answer 9

首先构建数据的最小生成树。 删除K-1个最昂贵的边将树分成K个簇，
因此您可以构建MST一次，查看各种K的簇间距离/度量标准， 并取曲线的拐点。

这仅适用于单链接聚类，但对于此方法来说，速度快且易于操作。此外，MST可视化效果良好。
例如，请参见stats.stackexchange聚类可视化软件下的MST图。

- denis

Answer 10

如果你使用MATLAB，任何2013b及以后的版本，你都可以利用函数evalclusters来找出给定数据集的最佳k值。

这个函数让你从三个聚类算法中选择-kmeans、linkage和gmdistribution。

它还让你从四个聚类评估标准中选择-CalinskiHarabasz、DaviesBouldin、gap和silhouette。

- Kristada673