我一直在使用 k-means 在 R 中对我的数据进行聚类，但我希望能够使用贝西信息准则 (BIC) 和 AIC 评估我的聚类的拟合与模型复杂性。目前我在 R 中使用的代码是：
KClData <- kmeans(Data, centers=2, nstart= 100)
但我希望能够提取 BIC 和对数似然。任何帮助将不胜感激！
问问题
21659 次
5 回答
17
对于其他登陆这里的人，Sherry Towers 在http://sherrytowers.com/2013/10/24/k-means-clustering/提出了一种方法，该方法使用来自stats::kmeans. 我引用：
可以使用以下函数计算 AIC：
kmeansAIC = function(fit){ m = ncol(fit$centers) n = length(fit$cluster) k = nrow(fit$centers) D = fit$tot.withinss return(D + 2*m*k) }
从 的帮助中stats::AIC还可以看到，BIC 的计算方法与 AIC 类似。获取 BIC 的一种简单方法是return()将上述函数中的 替换为：
return(data.frame(AIC = D + 2*m*k,
BIC = D + log(n)*m*k))
因此，您将按如下方式使用它：
fit <- kmeans(x = data,centers = 6)
kmeansAIC(fit)
于 2014-08-28T20:09:09.103 回答
7
要计算 BIC，只需将.5*k*d*log(n)（其中k是平均值的数量，d是数据集中向量的长度，以及n数据点的数量）添加到标准 k-means 误差函数。
标准的 k-means 惩罚是\sum_n (m_k(n)-x_n)^2，其中m_k(n)是与第 n 个数据点相关的平均值。这个惩罚可以解释为对数概率，所以 BIC 是完全有效的。
BIC 只是在与 成比例的 k 均值误差上添加了一个额外的惩罚项k。
于 2013-04-07T17:40:20.917 回答
4
只是为了添加 user1149913 所说的内容（我没有足够的声誉来发表评论），因为您在 R 中使用了 kmeans 函数，\sum_n (m_k(n)-x_n)^2因此已经为您计算了KClData$tot.withinss.
于 2014-07-31T16:13:58.310 回答
2
我们可以为对象定义一个对数似然函数，而不是重新实现AICor ；这将被包中的函数使用。BICkmeansBICstats
logLik.kmeans <- function(object) structure(
-object$tot.withinss/2,
df = nrow(object$centers)*ncol(object$centers),
nobs = length(object$cluster)
)
然后使用它，BIC正常调用。例如：
example(kmeans, local=FALSE)
BIC(cl)
# [1] 26.22842084
此方法将在stackoverflow软件包的下一个版本中提供。
于 2015-10-18T19:22:19.430 回答
0
一个函数qualityCriterion::longitudinalData可以计算 k-means 集群的 BIC 和 AIC。它首先计算每个集群中心在每个个体上的可能性，然后再与集群大小的权重合并。正态密度函数的 sd 基于 RSS。
虽然原始代码给出了 -BIC 而不是 BIC，但我采用了 BIC 的代码：
qualityCriterion <- function (traj, clusters)
{
if (nrow(traj) != length(clusters)) {
stop("[qualityCriterion] the cluster and the number of trajectory should be the same.")
}
clusters <- as.integer(clusters)
nbIndiv <- nrow(traj)
nbTime <- ncol(traj)
nbClusters <- length(unique(clusters))
# Cluster frequency
preProba <- as.numeric(table(clusters))
preProba <- preProba / sum(preProba)
# Centers as cluster means
moy <- matrix(, nbClusters, nbTime)
for (i in 1:nbClusters) {
moy[i, ] <- apply(traj[as.numeric(clusters) == i, , drop = FALSE], 2, meanNA)
}
# sd of residuals
ecart <- sqrt(mean(as.numeric(traj - moy[clusters, ])^2, na.rm=TRUE))
# likelihood
vraisIndivXcluster <- matrix(, nbIndiv, nbClusters)
for (i in 1:nbClusters) {
vraisIndivXcluster[, i] <- preProba[i] * apply(dnorm(t(traj), moy[i, ], ecart),
2, prod, na.rm = TRUE
)
}
vraisIndivXcluster <- apply(vraisIndivXcluster, 1, sum)
logVraisemblance <- sum(log(vraisIndivXcluster))
nbParam <- nbClusters * nbTime + 1 # cluster centers and sd
#BIC <- -2 * logVraisemblance + nbParam * log(nbIndiv) # BIC for time series
BIC2 <- -2 * logVraisemblance + nbParam * log(nbIndiv * nbTime) # BIC for independent columns
AIC <- 2 * nbParam - 2 * logVraisemblance
#AICc <- AIC + (2 * nbParam * (nbParam + 1)) / (nbIndiv - nbParam - 1) # AICc for time series
#AICc2 <- AIC + (2 * nbParam * (nbParam + 1)) / (nbIndiv * nbTime - nbParam - 1) # AICc for independent columns
return(list(criters = c(BIC2 = BIC2, AIC = AIC)))
}
100 个个体的示例，每个个体有 2 个数据点，形成 3 个集群：
set.seed(1)
dat <- matrix(rnorm(100 * 2), nrow = 100, ncol = 2) # data of 100 individuals
dat[34:66,] <- dat[34:66,] + 4
dat[67:100,] <- dat[67:100,] + 8
plot(dat[,1], dat[,2]) # 3 cluster centers at (0,0), (4,4), (8,8)
正如预期的那样，k = 3 个集群时的最小 BIC：
# k-means with k = 2:5
clusters_2 <- kmeans(dat, centers = 2)
clusters_3 <- kmeans(dat, centers = 3)
clusters_4 <- kmeans(dat, centers = 4)
clusters_5 <- kmeans(dat, centers = 5)
BIC <- c(qualityCriterion(dat, rep(1, nrow(dat)))$criters["BIC2"],
qualityCriterion(dat, clusters_2$cluster)$criters["BIC2"],
qualityCriterion(dat, clusters_3$cluster)$criters["BIC2"],
qualityCriterion(dat, clusters_4$cluster)$criters["BIC2"],
qualityCriterion(dat, clusters_5$cluster)$criters["BIC2"]
)
plot(1:5, BIC, xlab = "k", ylab = "BIC")
lines(1:5, BIC)
于 2021-09-21T12:27:03.647 回答