Utilisation de BIC pour estimer le nombre de k dans KMEANS

13

J'essaie actuellement de calculer le BIC pour mon jeu de données de jouets (ofc iris (:). Je veux reproduire les résultats comme indiqué ici (Fig. 5). Ce papier est également ma source pour les formules BIC.

J'ai 2 problèmes avec ceci:

  • Notation:
    • ni = nombre d'éléments dans le clusteri
    • Ci = coordonnées centrales du clusteri
    • xj = points de données affectés au clusteri
    • m = nombre de clusters

1) La variance telle que définie dans l'équation. (2):

i=1nimj=1nixjCi2

Pour autant que je puisse voir, il est problématique et non couvert que la variance peut être négative lorsqu'il y a plus de clusters m que d'éléments dans le cluster. Est-ce correct?

2) Je ne peux tout simplement pas faire fonctionner mon code pour calculer le BIC correct. J'espère qu'il n'y a pas d'erreur, mais il serait très apprécié que quelqu'un puisse vérifier. L'équation entière peut être trouvée dans l'équation. (5) dans le document. J'utilise scikit learn pour tout en ce moment (pour justifier le mot-clé: 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")

Mes résultats pour le BIC ressemblent à ceci:

Ce qui n'est même pas proche de ce à quoi je m'attendais et n'a également aucun sens ... J'ai regardé les équations maintenant depuis un certain temps et je n'arrive pas à trouver plus loin mon erreur):

Kam Sen
la source
Vous pouvez trouver le calcul du BIC pour le clustering ici . C'est ainsi que SPSS le fait. Pas nécessairement exactement de la même manière que vous le montrez.
ttnphns
Merci ttnphns. J'ai déjà vu ta réponse. Mais cela n'a aucune référence sur la façon dont les étapes sont dérivées et donc pas ce que je cherchais. De plus, cette sortie SPSS ou quelle que soit la syntaxe n'est pas très lisible. Merci quand même. En raison du manque d'intérêt pour ces questions, je vais chercher la référence et utiliser une autre estimation pour la variance.
Kam Sen
Je sais que cela ne répond pas à votre question (donc je le laisse comme un commentaire), mais le package R mclust s'adapte aux modèles de mélange fini (une méthode de clustering paramétrique) et optimise automatiquement le nombre, la forme, la taille, l'orientation et l'hétérogénéité des clusters. Je comprends que vous utilisez sklearn, mais je voulais juste jeter cela là-bas.
Brash Equilibrium
1
Brash, sklearn a également GMM
eyaler
@KamSen pouvez-vous m'aider ici? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede

Réponses:

14

Il semble que vous ayez quelques erreurs dans vos formules, déterminées en les comparant à:

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

Ici, il y a trois erreurs dans le papier, il manque un facteur de d aux quatrième et cinquième lignes, la dernière ligne remplace m par 1. Il devrait être:

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.

Le const_term:

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

devrait être:

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

3.

La formule de variance:

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

devrait être un scalaire:

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.

Utilisez des journaux naturels, au lieu de vos journaux base10.

5.

Enfin, et surtout, le BIC que vous calculez a un signe inverse de la définition régulière. vous cherchez donc à maximiser au lieu de minimiser

eyaler
la source
1
MK(ϕ)2
@eyaler pouvez-vous me corriger ici? : - stats.stackexchange.com/questions/342258/…
Pranay Wankhede
pouvez-vous lier un article ou l'écrire en balisage mathématique?
donlan
Veuillez voir ma question connexe ici: stats.stackexchange.com/questions/374002/…
rnso
@ Seanny123 et eyaler, veuillez consulter le post stats.stackexchange.com/questions/374002/… de rnso. Cette formule donne environ 9 clusters sur des données d'iris qui devraient avoir 3 clusters
Bernardo Braga
11

C'est essentiellement la solution d'Eyaler, avec quelques notes .. Je viens de le taper si quelqu'un voulait un copier / coller rapide:

Remarques:

  1. eyalers 4ème commentaire est incorrect np.log est déjà un journal naturel, aucun changement nécessaire

  2. eyalers 5ème commentaire sur l'inverse est correct. Dans le code ci-dessous, vous recherchez le MAXIMUM - gardez à l'esprit que l'exemple a des nombres BIC négatifs

Le code est le suivant (encore une fois, tout crédit à 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
la source
En regardant github.com/bobhancock/goxmeans/blob/master/doc/BIC_notes.pdf pourriez-vous expliquer comment cette formule BIC optimise pour le MAXIMUM? Pourriez-vous montrer le minimum et expliquer ce qu'il fait en langage verbal? trouver difficile d'interpréter la formule
user305883
Veuillez voir ma question connexe ici: stats.stackexchange.com/questions/374002/…
rnso
1
il semble y avoir un bug dans la formule. Quelqu'un a-t-il réussi à le réparer?
STiGMa