由于我刚才正在研究差距统计数据，因此我将尝试回答这个老问题。

$\newcommand{\Gap}{{\rm Gap}}$

将数据点划分为$k$使用例如k-means 的集群，您需要一些指标来描述集群的紧凑程度。这样做的一种方法是计算集群内点之间的所有成对距离并将这些距离平均。对所有集群执行此操作，将这些平均距离相加并记录日志。这是$\log(W_{k})$在Tibshirani 等人。2001 年。

如果集群是真实的并且非常紧凑，$W_{k}$会很低。但是，我们需要一些基线参考来说明“低”是什么意思。为此，我们问什么$W_{k}$我们可以从与我们拥有的数据集相似但根本没有集群的数据集中期望的值。观察到的最大差异$W_{k}$和“空”$W_{k}$我们越有信心集群的数量是正确的。

假设我们的数据集有$n$数据点在$d$方面。蒂布希拉尼等人。表明这样的参考数据集可以通过取$n$数据点，其中每个数据点的坐标是从均匀分布中得出的。每个维度中均匀分布的最小值和最大值由该维度中观察到的数据点的最小值和最大值给出。

然后使用与真实数据集相同的策略（例如，k-means with$k$簇），计算其$W_{k}$，并计算观察到的差异$\log(W_{k})$和参考$\log{W_k}$. 当然，一个随机抽取的数据集并不是可靠的参考。相反，您绘制几个数据集并对结果进行平均$\log(W_{k})$的。这是$E^{*}\{\log(W_{k})\}$和区别$E^{*}\{\log(W_{k})\} - \log(W_{k})$是差距统计量。

蒂布希拉尼等人。还表明旋转数据集并重新映射主分量轴上的数据点并使用此旋转数据集绘制参考空数据集会更好。这给出了更好的结果，因为原始变量（维度）可能相互关联。旋转后，相关性已被消除。

这是在 R 中使用模拟的示例。首先，我生成了一个示例数据集，该示例数据集具有 2 个维度和 3 个真实集群，每个集群包含 20 个观察值。

library(data.table)

options(digits= 3)
set.seed(1234)

K <- 3
ksize <- 20

M <- rbind(
    mvrnorm(n= ksize, mu= c( 5,  5), Sigma= diag(c(1,1))),
    mvrnorm(n= ksize, mu= c( 5, 10), Sigma= diag(c(1,1))),
    mvrnorm(n= ksize, mu= c(10, 10), Sigma= diag(c(1,1)))
)
data <- data.table(M, cluster= rep(LETTERS[1:K], each= ksize), observation= 1:nrow(M))

plot(data$V1, data$V2, col= rep(1:K, each= ksize), pch= 19, xlab= 'Dim 1', ylab= 'Dim 2')

# 60 observations in 2 dimensions
M
       [,1]  [,2]
 [1,]  4.87  3.79
 [2,]  5.49  5.28
 [3,]  5.44  6.08
 [4,]  4.54  2.65
       ...
[57,] 11.24  8.87
[58,]  9.83 10.88
[59,]  9.33 10.97
[60,] 10.03 12.12

这是一个方便的函数，给定观察X计算所有成对的欧几里得平方距离，并将它们平均（根据 Tibshirani等人）。

avg.dist2 <- function(X) {
    dist2 <- dist(X)^2
    avg   <- sum(dist2)/(2 * nrow(X))
    return(avg)
}

在这里，我使用 k-means 对数据集进行聚类并计算$\log(W_{k})$：

k  <- 3
km <- kmeans(M, centers= k)
data[, cluster := km$cluster]

logW <- log(sum(data[, list(d2= avg.dist2(X= cbind(V1, V2))), by= cluster]$d2))

现在我生成参考数据集。首先，我通过将其列居中（即，从每个值中减去其列的平均值）来旋转观察到的数据，并将其乘以它的特征向量（从奇异值分解中获得的特征向量）：

M_cent <- scale(M, center= TRUE, scale= FALSE)
M_rot  <- M %*% svd(M_cent)$v
plot(M_rot[,1], M_rot[,2], col= rep(1:K, each= ksize), pch= 19, xlab= 'PC 1', ylab= 'PC 2')

现在我可以绘制参考数据集，对它们进行聚类，计算$W_{k}$对于每个随机绘制的数据集 - 与真实数据相同。请注意runif(..., min= min(M_rot[,v]), max= max(M_rot[,v]))用于创建新数据点。

B <- 500
null_logW <- rep(NA, B)
for(i in 1:length(null_logW)){
    null_M <- matrix(data= NA, ncol= ncol(M_rot), nrow= nrow(M_rot))
    for(v in 1:ncol(M_rot)){
        null_M[,v] <- runif(n= nrow(M_rot), min= min(M_rot[,v]), max= max(M_rot[,v]))
    }
    null_km <- kmeans(null_M, centers= k)
    null_M  <- as.data.table(null_M)
    null_M[, cluster := null_km$cluster]
    null_logW[i] <- log(sum(null_M[, list(d2= avg.dist2(X= cbind(V1, V2))), by= cluster]$d2))
}
plot(null_M$V1, null_M$V2, col= null_M$cluster, pch= 19, xlab= 'Dim 1', ylab= 'Dim 2')

这是其中一个随机数据集的示例：

最后，平均空值$\log(W_{k})$并计算差距统计。正如 Tibshirani等人所建议的那样，我们还应该计算空分布的标准偏差，以指导选择最佳聚类数：

ElogW <- mean(null_logW) # 4.88
sdk   <- sqrt((1 + 1/B) * sum((null_logW - mean(null_logW))^2)/B) # 0.0988
gap   <- ElogW - logW  # 1.04

现在，通过替换k <- 3为 egk <- 2或 4、5 等，您可以根据间隙统计找到最佳聚类数。即，k这样$\Gap(k)\ge \Gap(k+1)-sd_{k+1}$.

信用：我发现这里的代码the-gap-statistic很容易理解我自己。