Использование BIC для оценки количества k в KMEANS

В настоящее время я пытаюсь вычислить BIC для моего игрушечного набора данных (ofc iris (:). Я хочу воспроизвести результаты, как показано здесь (Рис. 5). Этот документ также является моим источником для формул BIC.

У меня есть 2 проблемы с этим:

• Обозначения:
  • $n_i$ я = количество элементов в кластере$i$
  • $C_i$ я = координаты центра кластера$i$
  • $x_j$ я = точки данных, назначенные кластеру$i$
  • $m$ = количество кластеров

1) Дисперсия, как определено в формуле. (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.

The 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 имеет обратный знак по сравнению с обычным определением. так что вы хотите максимизировать, а не минимизировать

Это в основном решение для глаз, с несколькими заметками ... Я просто напечатал его, если кто-то хотел быструю копию / вставку:

Примечания:

1. eyalers 4-й комментарий неверен np.log уже является естественным логом, никаких изменений не требуется

2. eyalers 5-й комментарий об обратном верен. В приведенном ниже коде вы ищете МАКСИМАЛЬНЫЙ - имейте в виду, что в примере есть отрицательные числа BIC

Код выглядит следующим образом (опять же, все заслуги перед глазами):

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