K-signifie un comportement incohérent choisissant K avec la méthode du coude, le BIC, la variance expliquée et la silhouette

23

J'essaie de regrouper certains vecteurs avec 90 fonctionnalités avec K-means. Étant donné que cet algorithme me demande le nombre de clusters, je veux valider mon choix avec de belles mathématiques. Je m'attends à avoir de 8 à 10 grappes. Les fonctionnalités sont à l'échelle Z-score.

Explication de la méthode et de la variance du coude

from scipy.spatial.distance import cdist, pdist
from sklearn.cluster import KMeans

K = range(1,50)
KM = [KMeans(n_clusters=k).fit(dt_trans) for k in K]
centroids = [k.cluster_centers_ for k in KM]

D_k = [cdist(dt_trans, cent, 'euclidean') for cent in centroids]
cIdx = [np.argmin(D,axis=1) for D in D_k]
dist = [np.min(D,axis=1) for D in D_k]
avgWithinSS = [sum(d)/dt_trans.shape[0] for d in dist]

# Total with-in sum of square
wcss = [sum(d**2) for d in dist]
tss = sum(pdist(dt_trans)**2)/dt_trans.shape[0]
bss = tss-wcss

kIdx = 10-1

# elbow curve
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(K, avgWithinSS, 'b*-')
ax.plot(K[kIdx], avgWithinSS[kIdx], marker='o', markersize=12, 
markeredgewidth=2, markeredgecolor='r', markerfacecolor='None')
plt.grid(True)
plt.xlabel('Number of clusters')
plt.ylabel('Average within-cluster sum of squares')
plt.title('Elbow for KMeans clustering')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(K, bss/tss*100, 'b*-')
plt.grid(True)
plt.xlabel('Number of clusters')
plt.ylabel('Percentage of variance explained')
plt.title('Elbow for KMeans clustering')

Méthode du coude Variance

D'après ces deux images, il semble que le nombre de clusters ne s'arrête jamais: D. Étrange! Où est le coude? Comment choisir K?

Critère d'information bayésien

Cette méthode provient directement de X-means et utilise le BIC pour choisir le nombre de clusters. un autre ref

    from sklearn.metrics import euclidean_distances
from sklearn.cluster import KMeans

def bic(clusters, centroids):
    num_points = sum(len(cluster) for cluster in clusters)
    num_dims = clusters[0][0].shape[0]
    log_likelihood = _loglikelihood(num_points, num_dims, clusters, centroids)
    num_params = _free_params(len(clusters), num_dims)
    return log_likelihood - num_params / 2.0 * np.log(num_points)


def _free_params(num_clusters, num_dims):
    return num_clusters * (num_dims + 1)


def _loglikelihood(num_points, num_dims, clusters, centroids):
    ll = 0
    for cluster in clusters:
        fRn = len(cluster)
        t1 = fRn * np.log(fRn)
        t2 = fRn * np.log(num_points)
        variance = _cluster_variance(num_points, clusters, centroids) or np.nextafter(0, 1)
        t3 = ((fRn * num_dims) / 2.0) * np.log((2.0 * np.pi) * variance)
        t4 = (fRn - 1.0) / 2.0
        ll += t1 - t2 - t3 - t4
    return ll

def _cluster_variance(num_points, clusters, centroids):
    s = 0
    denom = float(num_points - len(centroids))
    for cluster, centroid in zip(clusters, centroids):
        distances = euclidean_distances(cluster, centroid)
        s += (distances*distances).sum()
    return s / denom

from scipy.spatial import distance
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)



sns.set_style("ticks")
sns.set_palette(sns.color_palette("Blues_r"))
bics = []
for n_clusters in range(2,50):
    kmeans = KMeans(n_clusters=n_clusters)
    kmeans.fit(dt_trans)

    labels = kmeans.labels_
    centroids = kmeans.cluster_centers_

    clusters = {}
    for i,d in enumerate(kmeans.labels_):
        if d not in clusters:
            clusters[d] = []
        clusters[d].append(dt_trans[i])

    bics.append(compute_bic(kmeans,dt_trans))#-bic(clusters.values(), centroids))

plt.plot(bics)
plt.ylabel("BIC score")
plt.xlabel("k")
plt.title("BIC scoring for K-means cell's behaviour")
sns.despine()
#plt.savefig('figures/K-means-BIC.pdf', format='pdf', dpi=330,bbox_inches='tight')

entrez la description de l'image ici

Même problème ici ... Qu'est-ce que K?

Silhouette

    from sklearn.metrics import silhouette_score

s = []
for n_clusters in range(2,30):
    kmeans = KMeans(n_clusters=n_clusters)
    kmeans.fit(dt_trans)

    labels = kmeans.labels_
    centroids = kmeans.cluster_centers_

    s.append(silhouette_score(dt_trans, labels, metric='euclidean'))

plt.plot(s)
plt.ylabel("Silouette")
plt.xlabel("k")
plt.title("Silouette for K-means cell's behaviour")
sns.despine()

entrez la description de l'image ici

Alleluja! Ici, cela semble logique et c'est ce que j'attends. Mais pourquoi est-ce différent des autres?

marcodena
la source
1
Pour répondre à votre question sur le genou dans le cas de la variance, il semble que ce soit autour de 6 ou 7, vous pouvez l'imaginer comme le point de rupture entre deux segments approximatifs linéaires de la courbe. La forme du graphique n'est pas inhabituelle, le% de variance approchera souvent asymptotiquement 100%. Je mettrais k dans votre graphique BIC un peu plus bas, environ 5.
image_doctor
mais je devrais avoir (plus ou moins) les mêmes résultats dans toutes les méthodes, non?
marcodena
Je ne pense pas en savoir assez pour le dire. Je doute beaucoup que les trois méthodes soient mathématiquement équivalentes avec toutes les données, sinon elles n'existeraient pas en tant que techniques distinctes, donc les résultats comparatifs dépendent des données. Deux des méthodes donnent un nombre de grappes proches, la troisième est plus élevée mais pas énormément. Avez-vous des informations a priori sur le nombre réel de clusters?
image_doctor
Je ne suis pas sûr à 100% mais je m'attends à avoir de 8 à 10 grappes
marcodena
2
Vous êtes déjà dans le trou noir de "Curse of Dimensionality". Nothings fonctionne avant une réduction de dimensionnalité.
Kasra Manshaei

Réponses:

7

Il suffit de publier un résumé des commentaires ci-dessus et quelques réflexions supplémentaires pour que cette question soit supprimée des "questions sans réponse".

Le commentaire de Image_doctor est juste que ces graphiques sont typiques pour k-means. (Je ne connais pas la mesure "Silhouette" cependant.) La variance dans le cluster devrait baisser continuellement avec l'augmentation de k. Le coude est l'endroit où la courbe se plie le plus. (Pensez peut-être à "dérivée 2e" si vous voulez quelque chose de mathématique.)

En règle générale, il est préférable de choisir k en utilisant la tâche finale. N'utilisez pas de mesures statistiques de votre cluster pour prendre votre décision, mais utilisez les performances de bout en bout de votre système pour guider vos choix. N'utilisez les statistiques que comme point de départ.

Joachim Wagner
la source
5

Il est plus facile de trouver le coude en calculant les angles entre les segments consécutifs.

Remplacez votre:

kIdx = 10-1

avec:

seg_threshold = 0.95 #Set this to your desired target

#The angle between three points
def segments_gain(p1, v, p2):
    vp1 = np.linalg.norm(p1 - v)
    vp2 = np.linalg.norm(p2 - v)
    p1p2 = np.linalg.norm(p1 - p2)
    return np.arccos((vp1**2 + vp2**2 - p1p2**2) / (2 * vp1 * vp2)) / np.pi

#Normalize the data
criterion = np.array(avgWithinSS)
criterion = (criterion - criterion.min()) / (criterion.max() - criterion.min())

#Compute the angles
seg_gains = np.array([0, ] + [segments_gain(*
        [np.array([K[j], criterion[j]]) for j in range(i-1, i+2)]
    ) for i in range(len(K) - 2)] + [np.nan, ])

#Get the first index satisfying the threshold
kIdx = np.argmax(seg_gains > seg_threshold)

et vous verrez quelque chose comme: entrez la description de l'image ici

Si vous visualisez les seg_gains, vous verrez quelque chose comme ceci: entrez la description de l'image ici

J'espère que vous pouvez trouver le coude délicat maintenant :)

Sahloul
la source
3

J'ai créé une bibliothèque Python qui tente d'implémenter l' algorithme Kneedle pour détecter le point de courbure maximale dans des fonctions comme celle-ci. Il peut être installé avec pip install kneed.

Code et sortie pour quatre formes différentes de fonctions:

from kneed.data_generator import DataGenerator
from kneed.knee_locator import KneeLocator

import numpy as np

import matplotlib.pyplot as plt

# sample x and y
x = np.arange(0,10)
y_convex_inc = np.array([1,2,3,4,5,10,15,20,40,100])
y_convex_dec = y_convex_inc[::-1]
y_concave_dec = 100 - y_convex_inc
y_concave_inc = 100 - y_convex_dec

# find the knee points
kn = KneeLocator(x, y_convex_inc, curve='convex', direction='increasing')
knee_yconvinc = kn.knee

kn = KneeLocator(x, y_convex_dec, curve='convex', direction='decreasing')
knee_yconvdec = kn.knee

kn = KneeLocator(x, y_concave_inc, curve='concave', direction='increasing')
knee_yconcinc = kn.knee

kn = KneeLocator(x, y_concave_dec, curve='concave', direction='decreasing')
knee_yconcdec = kn.knee

# plot
f, axes = plt.subplots(2, 2, figsize=(10,10));
yconvinc = axes[0][0]
yconvdec = axes[0][1]
yconcinc = axes[1][0]
yconcdec = axes[1][1]

yconvinc.plot(x, y_convex_inc)
yconvinc.vlines(x=knee_yconvinc, ymin=0, ymax=100, linestyle='--')
yconvinc.set_title("curve='convex', direction='increasing'")

yconvdec.plot(x, y_convex_dec)
yconvdec.vlines(x=knee_yconvdec, ymin=0, ymax=100, linestyle='--')
yconvdec.set_title("curve='convex', direction='decreasing'")

yconcinc.plot(x, y_concave_inc)
yconcinc.vlines(x=knee_yconcinc, ymin=0, ymax=100, linestyle='--')
yconcinc.set_title("curve='concave', direction='increasing'")

yconcdec.plot(x, y_concave_dec)
yconcdec.vlines(x=knee_yconcdec, ymin=0, ymax=100, linestyle='--')
yconcdec.set_title("curve='concave', direction='decreasing'");

entrez la description de l'image ici

Kevin
la source