Verwenden von BIC zum Schätzen der Anzahl von k in KMEANS

Ich versuche derzeit, den BIC für meinen Spielzeugdatensatz (ofc iris (:)) zu berechnen. Ich möchte die hier gezeigten Ergebnisse reproduzieren (Abb. 5). Dieses Papier ist auch meine Quelle für die BIC-Formeln.

Ich habe 2 Probleme damit:

• Notation:
  • $n_i$ = Anzahl der Elemente in Cluster$i$
  • $C_i$ = Mittelkoordinaten des Clusters$i$
  • $x_j$ = Datenpunkte, die dem Cluster i zugewiesen sind$i$
  • $m$ = Anzahl der Cluster

1) Die Varianz wie in Gl. (2):

Soweit ich sehe, ist es problematisch und nicht abgedeckt, dass die Varianz negativ sein kann, wenn es mehr Cluster $m$ als Elemente im Cluster gibt. Ist das richtig?

2) Ich kann meinen Code einfach nicht zum Arbeiten bringen, um den korrekten BIC zu berechnen. Hoffentlich gibt es keinen Fehler, aber es wäre sehr dankbar, wenn jemand nachsehen könnte. Die ganze Gleichung ist in Gl. (5) in der Zeitung. Ich benutze Scikit Learn für alles im Moment (um das Schlüsselwort zu rechtfertigen: 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")

Meine Ergebnisse für den BIC sehen so aus:

Was nicht annähernd dem entspricht, was ich erwartet habe und auch keinen Sinn ergibt ... Ich habe mir die Gleichungen jetzt einige Zeit angesehen und finde meinen Fehler nicht weiter):

Anscheinend haben Sie einige Fehler in Ihren Formeln, die durch den Vergleich mit folgenden Faktoren festgestellt wurden:

• https://www.cs.cmu.edu/~dpelleg/download/xmeans.pdf (das Papier enthält einige Fehler)
• 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

Hier gibt es drei Fehler in der Arbeit, in der vierten und fünften Zeile fehlt ein Faktor von d, die letzte Zeile ersetzt m für 1. Es sollte sein:

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.

Der const_term:

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

sollte sein:

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

3.

Die Varianzformel:

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

sollte ein Skalar sein:

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.

Verwenden Sie natürliche Protokolle anstelle Ihrer base10-Protokolle.

5.

Schließlich und vor allem hat der BIC, den Sie berechnen, ein umgekehrtes Vorzeichen gegenüber der regulären Definition. Sie möchten also maximieren statt minimieren

Dies ist im Grunde die Lösung von eyaler, mit ein paar Notizen. Ich habe es nur getippt, wenn jemand ein schnelles Kopieren / Einfügen wünschte:

Anmerkungen:

1. eyalers 4. Kommentar ist falsch. np.log ist bereits ein natürliches Protokoll, keine Änderung erforderlich

2. eyalers 5. Kommentar über Inverse ist korrekt. Im folgenden Code suchen Sie nach dem MAXIMUM - denken Sie daran, dass das Beispiel negative BIC-Nummern hat

Der Code lautet wie folgt (wiederum alle Kredite an 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