机器算法验证 - 使用 BIC 估计 KMEANS 中的 k 个数 - 吾爱随笔录

使用 BIC 估计 KMEANS 中的 k 个数

机器算法验证 Python k-均值 scikit-学习比克

2022-02-09 01:47:58

我目前正在尝试为我的玩具数据集（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}} ‖ x_{j} - C_{i} ‖^{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 结果如下所示：

这甚至不符合我的预期，也没有任何意义......我现在看了一段时间的方程式，并没有进一步定位我的错误）：

3个回答

您的公式中似乎有一些错误，通过比较来确定：

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

const_term：

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

应该：

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

方差公式：

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

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

最后，也是最重要的一点，您正在计算的 BIC 具有与常规定义相反的符号。所以你正在寻找最大化而不是最小化

这基本上是 eyaler 的解决方案，带有一些注释。如果有人想要快速复制/粘贴，我只是将其输入：

笔记：

eyalers 第 4 条评论不正确 np.log 已经是自然日志，无需更改
eyalers 关于逆的第 5 条评论是正确的。在下面的代码中，您正在寻找 MAXIMUM - 请记住，该示例具有负 BIC 数字

代码如下（再次感谢 eyaler）：

from sklearn import cluster
from scipy.spatial import distance
import sklearn.datasets
from sklearn.preprocessing import StandardScaler
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 - m) / d) * sum([sum(distance.cdist(X[np.where(labels == i)], [centers[0][i]], 
             'euclidean')**2) for i in range(m)])

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

    BIC = 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

    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]

print BIC

在我的环境中，Prabhath 的回答不起作用，因为np.where()子句无法散列，X无法确定要引用的记录计算cl_var。让我修复该错误并重新发布代码，因为我没有足够的声誉来添加评论。

# Almost all credits to elayer and Prabhath
from sklearn import cluster

from scipy.spatial import distance
import sklearn.datasets
from sklearn.preprocessing import StandardScaler
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
    #fixed
    cl_var = (1.0 / (N - m) / d) * sum([sum(distance.cdist(X.iloc[labels == i], [centers[0][i]], 
             'euclidean')**2) for i in range(m)])

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

    BIC = 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

    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]

print BIC

其它你可能感兴趣的问题

上一篇如何检测两极分化的用户意见（高星级和低星级）下一篇为什么叫几何分布和超几何分布？