使用 BIC 估計 KMEANS 中的 k 個數

我目前正在嘗試為我的玩具數據集（ofc iris (: ) 計算 BIC。我想重現結果，如圖所示（圖 5）。那篇論文也是我 BIC 公式的來源。

我有兩個問題：

• 符號：
  • $ n_i $ = 簇中的元素數 $ i $
  • $ C_i $ = 集群的中心坐標 $ i $
  • $ x_j $ = 分配給集群的數據點 $ i $
  • $ m $ = 聚類數

1. 方程式中定義的方差。(2):

$$ \sum_i = \frac{1}{n_i-m}\sum_{j=1}^{n_i}\Vert x_j - C_i \Vert^2 $$

據我所知，這是有問題的，並且沒有涵蓋當有更多集群時方差可能為負 $ m $ 比集群中的元素。這個對嗎？

2. 我只是無法讓我的代碼工作來計算正確的 BIC。希望沒有錯誤，但如果有人可以檢查，將不勝感激。整個方程可以在方程中找到。(5) 在論文中。我現在正在使用 scikit learn 進行所有操作（以證明關鍵字 :P 的合理性）。

from sklearn import cluster
from scipy.spatial import distance
import sklearn.datasets
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import numpy as np

def compute_bic(kmeans,X):
   """
Computes the BIC metric for a given clusters

Parameters:
-----------------------------------------
kmeans: List of clustering object from scikit learn

X : multidimension np array of data points

Returns:
-----------------------------------------
BIC value
"""
   # assign centers and labels
   centers = [kmeans.cluster_centers_]
   labels  = kmeans.labels_
   #number of clusters
   m = kmeans.n_clusters
   # size of the clusters
   n = np.bincount(labels)
   #size of data set
   N, d = X.shape

   #compute variance for all clusters beforehand
   cl_var = [(1.0 / (n[i] - m)) * sum(distance.cdist(X[np.where(labels == i)], [centers[0][i]], 'euclidean')**2)  for i in xrange(m)]

   const_term = 0.5 * m * np.log10(N)

   BIC = np.sum([n[i] * np.log10(n[i]) -
          n[i] * np.log10(N) -
        ((n[i] * d) / 2) * np.log10(2*np.pi) -
         (n[i] / 2) * np.log10(cl_var[i]) -
        ((n[i] - m) / 2) for i in xrange(m)]) - const_term

   return(BIC)



# IRIS DATA
iris = sklearn.datasets.load_iris()
X = iris.data[:, :4]  # extract only the features
#Xs = StandardScaler().fit_transform(X)
Y = iris.target

ks = range(1,10)

# run 9 times kmeans and save each result in the KMeans object
KMeans = [cluster.KMeans(n_clusters = i, init="k-means++").fit(X) for i in ks]

# now run for each cluster the BIC computation
BIC = [compute_bic(kmeansi,X) for kmeansi in KMeans]

plt.plot(ks,BIC,'r-o')
plt.title("iris data (cluster vs BIC)")
plt.xlabel("# clusters")
plt.ylabel("# BIC")

我的 BIC 結果如下所示：

這甚至不符合我的預期，也沒有任何意義……我現在看了一段時間的方程，並沒有進一步定位我的錯誤）：

您的公式中似乎有一些錯誤，通過比較來確定：

• https://www.cs.cmu.edu/~dpelleg/download/xmeans.pdf（論文有一些錯誤）
• https://github.com/bobhancock/goxmeans/blob/master/km.go
• https://github.com/mynameisfiber/pyxmeans/blob/master/pyxmeans/xmeans.py
• https://github.com/bobhancock/goxmeans/blob/master/doc/BIC_notes.pdf

1.

np.sum([n[i] * np.log(n[i]) -
              n[i] * np.log(N) -
            ((n[i] * d) / 2) * np.log(2*np.pi) -
             (n[i] / 2) * np.log(cl_var[i]) -
            ((n[i] - m) / 2) for i in range(m)]) - const_term

這裡論文中出現了三個錯誤，第四行和第五行缺少一個因子 d，最後一行將 m 替換為 1。應該是：

np.sum([n[i] * np.log(n[i]) -
              n[i] * np.log(N) -
            ((n[i] * d) / 2) * np.log(2*np.pi*cl_var) -
            ((n[i] - 1) * d/ 2) for i in range(m)]) - const_term

2.

const_term：

const_term = 0.5 * m * np.log(N)

應該：

const_term = 0.5 * m * np.log(N) * (d+1)

3.

方差公式：

cl_var = [(1.0 / (n[i] - m)) * sum(distance.cdist(p[np.where(label_ == i)], [centers[0][i]], 'euclidean')**2)  for i in range(m)]

應該是一個標量：

cl_var = (1.0 / (N - m) / d) * sum([sum(distance.cdist(p[np.where(labels == i)], [centers[0][i]], 'euclidean')**2) for i in range(m)])

4.

使用自然日誌，而不是 base10 日誌。

5.

最後，也是最重要的一點，您正在計算的 BIC 具有與常規定義相反的符號。所以你正在尋找最大化而不是最小化

引用自：https://stats.stackexchange.com/questions/90769