Htydjtk

I'm trying to cluster some vectors with 90 features with K-means. Since this algorithm asks me the number of clusters, I want to validate my choice with some nice math.
I expect to have from 8 to 10 clusters. The features are Z-score scaled.

Elbow method and variance explained

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

From these two pictures, it seems that the number of clusters never stops :D. Strange! Where is the elbow? How can I choose K?

Bayesian information criterion

This methods comes directly from X-means and uses the BIC to choose the number of clusters. another 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')

Same problem here... What is 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()

Alleluja! Here it seems to make sense and this is what I expect. But why is this different from the others?

Just posting a summary of above comments and some more thoughts so that this question is removed from "unanswered questions".

Image_doctor's comment is right that these graphs are typical for k-means. (I am not familiar with the "Silhouette" measure though.) The in-cluster variance is expected to go down continuously with increasing k. The elbow is where the curve bends the most. (Maybe think "2nd derivative" if you want something mathematical.)

Generally, it is best to pick k using the final task. Do not use statistical measures of your cluster to make your decision but use the end-to-end performance of your system to guide your choices. Only use the statistics as a starting point.

Finding the elbow can be made more easier by computing the angles between the consecutive segments.

Replace your:

kIdx = 10-1

with:

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)

and you will see something like:

If you visualize the seg_gains, you will see something like this:

I hope you can find the tricky elbow now :)

I created a Python library that attempts to implement the Kneedle algorithim to detect the point of maximum curvature in functions like this.
It can be installed with pip install kneed.

Code and output for four different shapes of functions:

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'");

Another (simpler) solution to find the elbow point, or point of maximum curvature, to detect the number of clusters, is proposed in my article, see here.