Pregunta:
Actualmente estoy tratando de calcular el BIC para mi conjunto de datos de juguetes (ofc iris (:)). Quiero reproducir los resultados como se muestra aquí (Fig. 5). Ese documento también es mi fuente para las fórmulas 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 clúster $ i $
- $ m $ = número de conglomerados
1) La varianza como se define en la Ec. (2):
$$ \ sum_i = \ frac {1} {n_i-m} \ sum_ {j = 1} ^ {n_i} \ Vert x_j – C_i \ Vert ^ 2 $$
Por lo que puedo ver, es problemático y no se cubre 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. Es de esperar que no haya ningún error, pero sería muy apreciado si alguien pudiera comprobarlo. La ecuación completa se puede encontrar en la ecuación. (5) en el papel. 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 se acerca a lo que esperaba y tampoco tiene sentido … miré las ecuaciones ahora durante un tiempo y no estoy obteniendo más información sobre mi error):
Respuesta:
Parece que tiene algunos errores en sus fórmulas, según se determina comparando 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 la cuarta y quinta líneas le 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 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.
Use troncos naturales, en lugar de sus troncos base10.
5.
Finalmente, y lo más importante, el BIC que está calculando tiene un signo inverso de la definición regular. por lo que busca maximizar en lugar de minimizar