Internal Metrics: Silhouette and Calinski-Harabasz

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import silhouette_samples, silhouette_score
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
 
# Generate sample data
X, y_true = make_blobs(n_samples=500, centers=4, cluster_std=0.8, random_state=42)
 
# Perform clustering
kmeans = KMeans(n_clusters=4, random_state=42, n_init='auto')
cluster_labels = kmeans.fit_predict(X)
 
# Compute silhouette scores
silhouette_avg = silhouette_score(X, cluster_labels)
sample_silhouette_values = silhouette_samples(X, cluster_labels)
 
print(f"Overall Silhouette Score: {silhouette_avg:.4f}")
 
# Per-cluster analysis
n_clusters = len(np.unique(cluster_labels))
for i in range(n_clusters):
    cluster_silhouette = sample_silhouette_values[cluster_labels == i]
    cluster_size = len(cluster_silhouette)
    cluster_mean = np.mean(cluster_silhouette)
    negative_count = np.sum(cluster_silhouette < 0)
    
    print(f"
Cluster {i}:")
    print(f"  Size: {cluster_size}")
    print(f"  Mean silhouette: {cluster_mean:.4f}")
    print(f"  Min: {np.min(cluster_silhouette):.4f}, Max: {np.max(cluster_silhouette):.4f}")
    print(f"  Potentially misassigned (s < 0): {negative_count} ({100*negative_count/cluster_size:.1f}%)")
 
# Create silhouette plot
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
y_lower = 10
 
for i in range(n_clusters):
    cluster_silhouette = sample_silhouette_values[cluster_labels == i]
    cluster_silhouette.sort()
    
    y_upper = y_lower + len(cluster_silhouette)
    color = plt.cm.viridis(float(i) / n_clusters)
    
    ax.fill_betweenx(np.arange(y_lower, y_upper), 0, cluster_silhouette,
                     facecolor=color, edgecolor=color, alpha=0.7)
    ax.text(-0.05, y_lower + 0.5 * len(cluster_silhouette), str(i))
    y_lower = y_upper + 10
 
ax.axvline(x=silhouette_avg, color="red", linestyle="--", label=f"Mean: {silhouette_avg:.3f}")
ax.set_xlabel("Silhouette Coefficient")
ax.set_ylabel("Cluster")
ax.set_title("Silhouette Plot for Cluster Analysis")
ax.legend()
plt.tight_layout()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

import numpy as np
from sklearn.metrics import silhouette_score
from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
 
def simplified_silhouette(X, labels, centroids):
    """
    Simplified silhouette using distances to centroids instead of 
    full pairwise distances. Complexity: O(n * k) instead of O(n²).
    
    Parameters:
    -----------
    X : ndarray of shape (n_samples, n_features)
    labels : ndarray of shape (n_samples,) - cluster assignments
    centroids : ndarray of shape (n_clusters, n_features)
    
    Returns:
    --------
    float : simplified silhouette score
    """
    n_samples = X.shape[0]
    n_clusters = len(centroids)
    
    # Distance from each point to each centroid: O(n * k)
    distances_to_centroids = cdist(X, centroids, metric='euclidean')
    
    silhouette_scores = np.zeros(n_samples)
    
    for i in range(n_samples):
        cluster_i = labels[i]
        
        # a(i): distance to own centroid
        a_i = distances_to_centroids[i, cluster_i]
        
        # b(i): minimum distance to other centroids
        other_distances = distances_to_centroids[i, np.arange(n_clusters) != cluster_i]
        b_i = np.min(other_distances)
        
        # Silhouette score
        if max(a_i, b_i) == 0:
            silhouette_scores[i] = 0
        else:
            silhouette_scores[i] = (b_i - a_i) / max(a_i, b_i)
    
    return np.mean(silhouette_scores)
 
 
def sampled_silhouette(X, labels, sample_size=1000, n_bootstrap=10, random_state=42):
    """
    Estimate silhouette score using stratified sampling with bootstrap confidence intervals.
    
    Parameters:
    -----------
    X : ndarray of shape (n_samples, n_features)
    labels : ndarray of shape (n_samples,)
    sample_size : int, number of samples to use
    n_bootstrap : int, number of bootstrap iterations
    random_state : int
    
    Returns:
    --------
    tuple : (mean estimate, standard error)
    """
    rng = np.random.RandomState(random_state)
    n_samples = X.shape[0]
    unique_labels = np.unique(labels)
    
    bootstrap_scores = []
    
    for b in range(n_bootstrap):
        # Stratified sampling
        sample_indices = []
        for label in unique_labels:
            label_indices = np.where(labels == label)[0]
            # Sample proportionally from each cluster
            n_from_cluster = max(1, int(sample_size * len(label_indices) / n_samples))
            sampled = rng.choice(label_indices, size=min(n_from_cluster, len(label_indices)), replace=False)
            sample_indices.extend(sampled)
        
        sample_indices = np.array(sample_indices)
        X_sample = X[sample_indices]
        labels_sample = labels[sample_indices]
        
        score = silhouette_score(X_sample, labels_sample)
        bootstrap_scores.append(score)
    
    return np.mean(bootstrap_scores), np.std(bootstrap_scores)
 
 
# Example usage
from sklearn.datasets import make_blobs
 
# Generate larger dataset
X, y = make_blobs(n_samples=10000, centers=5, random_state=42)
kmeans = KMeans(n_clusters=5, random_state=42, n_init='auto')
labels = kmeans.fit_predict(X)
 
# Exact silhouette (for comparison)
exact_score = silhouette_score(X, labels)
print(f"Exact Silhouette Score: {exact_score:.4f}")
 
# Simplified silhouette
simplified_score = simplified_silhouette(X, labels, kmeans.cluster_centers_)
print(f"Simplified Silhouette Score: {simplified_score:.4f}")
 
# Sampled silhouette with CI
mean_est, std_err = sampled_silhouette(X, labels, sample_size=1000, n_bootstrap=20)
print(f"Sampled Silhouette Score: {mean_est:.4f} ± {1.96*std_err:.4f} (95% CI)")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import calinski_harabasz_score
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
 
def compute_variance_decomposition(X, labels):
    """
    Compute the within-cluster (W), between-cluster (B), 
    and total (T) sum of squares.
    """
    n = X.shape[0]
    global_centroid = np.mean(X, axis=0)
    unique_labels = np.unique(labels)
    
    # Total sum of squares
    T = np.sum((X - global_centroid) ** 2)
    
    # Within-cluster sum of squares
    W = 0
    for label in unique_labels:
        cluster_points = X[labels == label]
        cluster_centroid = np.mean(cluster_points, axis=0)
        W += np.sum((cluster_points - cluster_centroid) ** 2)
    
    # Between-cluster sum of squares
    B = 0
    for label in unique_labels:
        cluster_points = X[labels == label]
        cluster_centroid = np.mean(cluster_points, axis=0)
        B += len(cluster_points) * np.sum((cluster_centroid - global_centroid) ** 2)
    
    return T, W, B
 
 
def manual_calinski_harabasz(X, labels):
    """
    Manual implementation of CH index for understanding.
    """
    n = X.shape[0]
    k = len(np.unique(labels))
    
    T, W, B = compute_variance_decomposition(X, labels)
    
    # CH = (B / (k-1)) / (W / (n-k))
    ch = (B / (k - 1)) / (W / (n - k))
    
    return ch, T, W, B
 
 
# Generate data with 4 true clusters
X, y_true = make_blobs(n_samples=500, centers=4, cluster_std=1.0, random_state=42)
 
# Compute CH for different values of k
k_range = range(2, 11)
ch_scores = []
variance_data = []
 
for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init='auto')
    labels = kmeans.fit_predict(X)
    
    # Using sklearn
    ch_sklearn = calinski_harabasz_score(X, labels)
    
    # Manual computation for verification
    ch_manual, T, W, B = manual_calinski_harabasz(X, labels)
    
    ch_scores.append(ch_sklearn)
    variance_data.append({'k': k, 'T': T, 'W': W, 'B': B, 'CH': ch_sklearn})
    
    print(f"k={k}: CH={ch_sklearn:.2f}, W={W:.2f}, B={B:.2f}, W+B={W+B:.2f}, T={T:.2f}")
 
# Plot CH vs k
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
 
# CH score plot
axes[0].plot(k_range, ch_scores, 'bo-', linewidth=2, markersize=8)
axes[0].axvline(x=4, color='red', linestyle='--', label='True k=4')
optimal_k = k_range[np.argmax(ch_scores)]
axes[0].scatter([optimal_k], [max(ch_scores)], color='green', s=200, zorder=5, 
                label=f'Optimal k={optimal_k}')
axes[0].set_xlabel('Number of Clusters (k)')
axes[0].set_ylabel('Calinski-Harabasz Score')
axes[0].set_title('CH Score vs Number of Clusters')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
 
# Variance decomposition plot
W_values = [v['W'] for v in variance_data]
B_values = [v['B'] for v in variance_data]
 
axes[1].bar(k_range, W_values, label='Within (W)', alpha=0.7)
axes[1].bar(k_range, B_values, bottom=W_values, label='Between (B)', alpha=0.7)
axes[1].set_xlabel('Number of Clusters (k)')
axes[1].set_ylabel('Sum of Squares')
axes[1].set_title('Variance Decomposition: W + B = T')
axes[1].legend()
 
plt.tight_layout()
plt.show()