Wykorzystanie BIC do oszacowania liczby k w KMEANS

13

Obecnie próbuję obliczyć BIC dla mojego zestawu danych zabawek (ofc iris (:). Chcę odtworzyć wyniki, jak pokazano tutaj (ryc. 5). Ten papier jest również moim źródłem dla formuł BIC.

Mam z tym 2 problemy:

  • Notacja:
    • ni I = liczba elementów w klastrzei
    • Ci i = współrzędne środkowe klastrai
    • xj i = punkty danych przypisane do klastrai
    • m = liczba klastrów

1) Wariancja zdefiniowana w równaniu. (2):

i=1nimj=1nixjCi2

Z tego, co widzę, jest problematyczne i nieujęte, że wariancja może być ujemna, gdy w klastrze jest więcej klastrów m niż elementów. Czy to jest poprawne?

2) Po prostu nie mogę zmusić mojego kodu do działania w celu obliczenia poprawnego BIC. Mam nadzieję, że nie ma błędu, ale byłoby bardzo mile widziane, gdyby ktoś mógł to sprawdzić. Całe równanie można znaleźć w równaniu. (5) w pracy. Używam scikit learning do wszystkiego teraz (aby uzasadnić słowo kluczowe: 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")

Moje wyniki dla BIC wyglądają tak:

Co nie jest nawet bliskie temu, czego się spodziewałem, a także nie ma sensu ... Patrzyłem teraz na równania od jakiegoś czasu i nie dostaję dalszej lokalizacji mojego błędu):

Kam Sen
źródło
Obliczenia BIC dla klastrowania można znaleźć tutaj . Tak właśnie działa SPSS. Niekoniecznie dokładnie tak samo, jak pokazujesz.
ttnphns
Dziękuję ttnphns. Widziałem już twoją odpowiedź. Ale to nie ma odniesienia do tego, jak powstają kroki, a zatem nie tego, czego szukałem. Co więcej, to wyjście SPSS lub jakakolwiek składnia nie jest bardzo czytelna. Mimo wszystko dziekuję. Z powodu braku zainteresowania tymi pytaniami poszuka referencji i użyję innego oszacowania dla wariancji.
Kam Sen
Wiem, że to nie odpowiada na twoje pytanie (dlatego zostawiam to jako komentarz), ale pakiet R mclust pasuje do modeli mieszanki skończonej (parametryczna metoda grupowania) i automatycznie optymalizuje liczbę, kształt, rozmiar, orientację i heterogeniczność klastrów. Rozumiem, że używasz sklearn, ale po prostu chciałem to tam wyrzucić.
Brash Equilibrium
1
Zuchwały, sklearn ma również GMM
eyaler
@KamSen, czy możesz mi pomóc tutaj? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede

Odpowiedzi:

14

Wygląda na to, że masz kilka błędów w formułach, które zostały określone na podstawie:

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

Tutaj są trzy błędy w dokumencie, w czwartej i piątej linii brakuje współczynnika d, ostatnia linia zastępuje m dla 1. Powinno to być:

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)

Powinien być:

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

3)

Formuła wariancji:

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

powinien być skalarem:

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

Używaj dzienników naturalnych zamiast dzienników base10.

5

Co najważniejsze, obliczany przez ciebie BIC ma odwrotny znak w stosunku do zwykłej definicji. więc chcesz zmaksymalizować zamiast zminimalizować

eyaler
źródło
1
MK(ϕ)2
@eyaler, czy możesz mnie tutaj poprawić? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede
czy możesz połączyć artykuł lub napisać to znacznikiem matematycznym?
donlan
Pl, zobacz moje powiązane pytanie tutaj: stats.stackexchange.com/questions/374002/…
rnso
@ Seanny123 i eyaler patrz post stats.stackexchange.com/questions/374002/... od rnso. Ta formuła podaje około 9 skupień danych tęczówki, które powinny mieć 3 skupienia
Bernardo Braga
11

Jest to w zasadzie rozwiązanie eyalera, z kilkoma notatkami. Napisałem to, jeśli ktoś chciałby szybko skopiować / wkleić:

Uwagi:

  1. eyalers 4. komentarz jest niepoprawny np. log jest już logiem naturalnym, nie trzeba wprowadzać żadnych zmian

  2. eyalers Piąty komentarz na temat odwrotności jest poprawny. W poniższym kodzie szukasz MAKSIMUM - pamiętaj, że przykład ma ujemne liczby BIC

Kod jest następujący (ponownie, cały kredyt dla eyalera):

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
źródło
Patrząc na github.com/bobhancock/goxmeans/blob/master/doc/BIC_notes.pdf, czy mógłbyś wyjaśnić, w jaki sposób ta formuła BIC optymalizuje się dla MAXIMUM? Czy możesz pokazać minimum i wyjaśnić, co robi w języku werbalnym? trudno jest zinterpretować formułę
305883
Pl, zobacz moje powiązane pytanie tutaj: stats.stackexchange.com/questions/374002/…
rnso
1
wydaje się, że w formule jest błąd. Czy ktoś zdołał to naprawić?
STiGMa