Usando BIC para estimar el número de k en KMEANS

Actualmente estoy tratando de calcular el BIC para mi conjunto de datos de juguete (ofc iris (:). Quiero reproducir los resultados como se muestra aquí (Fig. 5). Ese documento también es mi fuente para las fórmulas de BIC.

Tengo 2 problemas con esto:

• Notación:
  • $n_i$ = número de elementos en el clúster$i$
  • $C_i$ = coordenadas centrales del grupo$i$
  • $x_j$ = puntos de datos asignados al grupo$i$
  • $m$ = número de grupos

1) La varianza como se define en la ecuación. (2):

Hasta donde puedo ver, es problemático y no está cubierto que la varianza puede ser negativa cuando hay más grupos $m$ que elementos en el grupo. ¿Es esto correcto?

2) Simplemente no puedo hacer que mi código funcione para calcular el BIC correcto. Esperemos que no haya ningún error, pero sería muy apreciado si alguien pudiera verificarlo. Toda la ecuación se puede encontrar en la ecuación. (5) en el periódico. Estoy usando scikit learn para todo en este momento (para justificar la palabra clave: 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")

Mis resultados para el BIC se ven así:

Lo cual ni siquiera está cerca de lo que esperaba y tampoco tiene sentido ... Miré las ecuaciones ahora por un tiempo y no encuentro más mi error):

Parece que tiene algunos errores en sus fórmulas, según lo determinado en comparación con:

• https://www.cs.cmu.edu/~dpelleg/download/xmeans.pdf (hay algunos errores en el documento)
• 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

Aquí hay tres errores en el documento, a las líneas cuarta y quinta les falta un factor de d, la última línea sustituye m por 1. Debería ser:

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)

El const_term:

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

debiera ser:

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

3)

La fórmula de la varianza:

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

debería ser un escalar:

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)

Utilice registros naturales, en lugar de sus registros de base10.

5)

Finalmente, y lo más importante, el BIC que está computando tiene un signo inverso de la definición regular. entonces busca maximizar en lugar de minimizar

Esta es básicamente la solución de eyaler, con algunas notas. Simplemente la escribí si alguien quería copiar / pegar rápidamente:

Notas:

1. el 4º comentario de eyalers es incorrecto np.log ya es un registro natural, no se necesita ningún cambio

2. El quinto comentario de eyalers sobre el inverso es correcto. En el siguiente código, está buscando el MÁXIMO: tenga en cuenta que el ejemplo tiene números BIC negativos

El código es el siguiente (de nuevo, todo crédito para 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