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

13

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:
    • nich = Anzahl der Elemente in Clusterich
    • Cich = Mittelkoordinaten des Clustersich
    • xj = Datenpunkte, die dem Cluster i zugewiesen sindich
    • m = Anzahl der Cluster

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

ich=1nich-mj=1nichxj-Cich2

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

Kam Sen
quelle
Die Berechnung des BIC für das Clustering finden Sie hier . So macht es SPSS. Nicht unbedingt genau so wie du es zeigst.
TTNPHNS
Vielen Dank, dass Sie ttnphns. Ich habe deine Antwort schon einmal gesehen. Aber das hat keinen Bezug darauf, wie die Schritte abgeleitet werden und somit nicht das, wonach ich gesucht habe. Darüber hinaus ist diese SPSS-Ausgabe oder was auch immer die Syntax ist nicht sehr lesbar. Trotzdem danke. Aufgrund des mangelnden Interesses an diesen Fragen werde ich nach der Referenz suchen und eine andere Schätzung für die Varianz verwenden.
Kam Sen
Ich weiß, dass dies Ihre Frage nicht beantwortet (daher lasse ich es als Kommentar), aber das R-Paket mclust passt auf endliche Mischungsmodelle (eine parametrische Clustering-Methode) und optimiert automatisch die Anzahl, Form, Größe, Ausrichtung und Heterogenität von Clustern. Ich verstehe, Sie verwenden sklearn, wollten das aber nur rauswerfen.
Dreistes Gleichgewicht
1
Brash, sklearn hat auch GMM
eyaler 16.04.15
@KamSen kannst du mir bitte hier weiterhelfen? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede

Antworten:

14

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

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

eyaler
quelle
1
MK(ϕ)2
@eyaler kannst du mich bitte hier korrigieren? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede
Kannst du eine Arbeit verlinken oder in mathematischer Aufschrift schreiben?
donlan
Bitte
rnso
@ Seanny123 und eyaler siehe post stats.stackexchange.com/questions/374002/… von rnso. Diese Formel gibt etwa 9 Cluster auf Irisdaten, die 3 Cluster haben sollten
Bernardo Braga
11

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
Prabhath Nanisetty
quelle
Wenn Sie sich github.com/bobhancock/goxmeans/blob/master/doc/BIC_notes.pdf ansehen, können Sie erläutern, wie diese BIC-Formel für das MAXIMUM optimiert wird? Könnten Sie für das Minimum zeigen und erklären, was es in der verbalen Sprache tut? Es
fällt
Bitte
rnso
1
Es scheint einen Fehler in der Formel zu geben. Hat es jemand geschafft, es zu beheben?
STiGMa