Clustering Evaluation

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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
 
def compute_wss(X, k_range, random_state=42):
    """
    Compute Within-Cluster Sum of Squares for different k values.
    """
    wss_values = []
    models = {}
    
    for k in k_range:
        model = KMeans(n_clusters=k, n_init='auto', random_state=random_state)
        model.fit(X)
        wss_values.append(model.inertia_)
        models[k] = model
    
    return np.array(wss_values), models
 
 
def find_elbow_perpendicular(k_range, wss_values):
    """
    Find elbow using perpendicular distance method.
    
    Draw a line from the first point to the last point.
    The elbow is the k with maximum perpendicular distance to this line.
    """
    k_range = np.array(k_range)
    
    # Normalize both axes to [0, 1]
    k_norm = (k_range - k_range.min()) / (k_range.max() - k_range.min())
    wss_norm = (wss_values - wss_values.min()) / (wss_values.max() - wss_values.min())
    
    # Line from first to last point
    p1 = np.array([k_norm[0], wss_norm[0]])
    p2 = np.array([k_norm[-1], wss_norm[-1]])
    
    # Compute perpendicular distance for each point
    distances = []
    for i in range(len(k_range)):
        p = np.array([k_norm[i], wss_norm[i]])
        # Distance from point to line formula
        d = np.abs(np.cross(p2 - p1, p1 - p)) / np.linalg.norm(p2 - p1)
        distances.append(d)
    
    # Find maximum distance
    elbow_idx = np.argmax(distances)
    return k_range[elbow_idx], distances
 
 
def find_elbow_second_derivative(k_range, wss_values):
    """
    Find elbow using second derivative (acceleration) method.
    
    The elbow is where the second derivative is maximum.
    """
    # First derivative (rate of change)
    first_deriv = np.diff(wss_values)
    
    # Second derivative (acceleration)
    second_deriv = np.diff(first_deriv)
    
    # The elbow is where second derivative is maximum (least negative = most positive)
    elbow_idx = np.argmax(second_deriv) + 1  # +1 because of diff offset
    
    return k_range[elbow_idx], first_deriv, second_deriv
 
 
# Example with clear elbow
np.random.seed(42)
X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=1.0, random_state=42)
 
k_range = range(1, 11)
wss_values, models = compute_wss(X, k_range)
 
# Find elbow using different methods
elbow_perp, distances = find_elbow_perpendicular(list(k_range), wss_values)
elbow_deriv, first_deriv, second_deriv = find_elbow_second_derivative(list(k_range), wss_values)
 
print(f"Elbow by perpendicular distance: k = {elbow_perp}")
print(f"Elbow by second derivative: k = {elbow_deriv}")
print(f"True k: 4")
 
# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
 
# Main elbow plot
axes[0, 0].plot(k_range, wss_values, 'bo-', linewidth=2, markersize=8)
axes[0, 0].axvline(x=elbow_perp, color='red', linestyle='--', 
                   label=f'Elbow (perp): k={elbow_perp}')
axes[0, 0].axvline(x=4, color='green', linestyle=':', 
                   label='True k=4', linewidth=2)
axes[0, 0].set_xlabel('Number of Clusters (k)')
axes[0, 0].set_ylabel('Within-Cluster Sum of Squares (WSS)')
axes[0, 0].set_title('Elbow Method')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
 
# Perpendicular distances
axes[0, 1].bar(k_range, distances, color='steelblue', alpha=0.7)
axes[0, 1].axvline(x=elbow_perp, color='red', linestyle='--', 
                   label=f'Max distance at k={elbow_perp}')
axes[0, 1].set_xlabel('Number of Clusters (k)')
axes[0, 1].set_ylabel('Perpendicular Distance')
axes[0, 1].set_title('Perpendicular Distance Method')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
 
# First derivative
axes[1, 0].plot(list(k_range)[1:], -first_deriv, 'ro-', linewidth=2, markersize=8)
axes[1, 0].set_xlabel('Number of Clusters (k)')
axes[1, 0].set_ylabel('Decrease in WSS')
axes[1, 0].set_title('First Derivative (Rate of Improvement)')
axes[1, 0].grid(True, alpha=0.3)
 
# Second derivative
axes[1, 1].plot(list(k_range)[2:], second_deriv, 'go-', linewidth=2, markersize=8)
axes[1, 1].axvline(x=elbow_deriv, color='red', linestyle='--', 
                   label=f'Max at k={elbow_deriv}')
axes[1, 1].set_xlabel('Number of Clusters (k)')
axes[1, 1].set_ylabel('Second Derivative')
axes[1, 1].set_title('Second Derivative (Acceleration)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
 
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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
from sklearn.datasets import make_blobs
 
class MultiMetricKSelector:
    """
    Select k using multiple clustering metrics with consensus voting.
    """
    
    def __init__(self, k_range, n_stability_runs=20, random_state=42):
        self.k_range = list(k_range)
        self.n_stability_runs = n_stability_runs
        self.random_state = random_state
        self.results = None
    
    def _compute_stability(self, X, k):
        """Compute bootstrap stability for a given k."""
        rng = np.random.RandomState(self.random_state)
        n_samples = X.shape[0]
        pairwise_aris = []
        
        from sklearn.metrics import adjusted_rand_score
        
        for run in range(self.n_stability_runs):
            # Two bootstrap samples
            idx1 = rng.choice(n_samples, n_samples, replace=True)
            idx2 = rng.choice(n_samples, n_samples, replace=True)
            
            m1 = KMeans(n_clusters=k, n_init='auto', random_state=run)
            m2 = KMeans(n_clusters=k, n_init='auto', random_state=run+1000)
            
            l1 = m1.fit_predict(X[np.unique(idx1)])
            l2 = m2.fit_predict(X[np.unique(idx2)])
            
            common = sorted(set(np.unique(idx1)) & set(np.unique(idx2)))
            if len(common) >= 10:
                idx1_map = {i: j for j, i in enumerate(np.unique(idx1))}
                idx2_map = {i: j for j, i in enumerate(np.unique(idx2))}
                
                labels1 = [l1[idx1_map[i]] for i in common]
                labels2 = [l2[idx2_map[i]] for i in common]
                
                pairwise_aris.append(adjusted_rand_score(labels1, labels2))
        
        return np.mean(pairwise_aris) if pairwise_aris else 0
    
    def fit(self, X):
        """
        Compute all metrics for all k values.
        """
        results = {
            'k': self.k_range,
            'wss': [],
            'silhouette': [],
            'calinski_harabasz': [],
            'davies_bouldin': [],
            'stability': []
        }
        
        print("Computing metrics for each k...")
        for k in self.k_range:
            model = KMeans(n_clusters=k, n_init='auto', random_state=self.random_state)
            labels = model.fit_predict(X)
            
            results['wss'].append(model.inertia_)
            results['silhouette'].append(silhouette_score(X, labels))
            results['calinski_harabasz'].append(calinski_harabasz_score(X, labels))
            results['davies_bouldin'].append(davies_bouldin_score(X, labels))
            results['stability'].append(self._compute_stability(X, k))
            
            print(f"  k={k}: Sil={results['silhouette'][-1]:.3f}, "
                  f"CH={results['calinski_harabasz'][-1]:.1f}, "
                  f"DB={results['davies_bouldin'][-1]:.3f}, "
                  f"Stab={results['stability'][-1]:.3f}")
        
        # Convert to arrays
        for key in results:
            if key != 'k':
                results[key] = np.array(results[key])
        
        self.results = results
        return self
    
    def recommend(self):
        """
        Provide k recommendation based on consensus.
        """
        if self.results is None:
            raise ValueError("Must call fit() first")
        
        k_range = self.results['k']
        n_k = len(k_range)
        
        # Score each k: higher is better (normalize metrics)
        # For DB and WSS, lower is better, so we negate
        
        def normalize(arr, higher_better=True):
            """Normalize to [0, 1], 1 is best."""
            arr = np.array(arr, dtype=float)
            if arr.max() == arr.min():
                return np.ones_like(arr) * 0.5
            normalized = (arr - arr.min()) / (arr.max() - arr.min())
            return normalized if higher_better else (1 - normalized)
        
        scores = {
            'silhouette': normalize(self.results['silhouette'], higher_better=True),
            'calinski_harabasz': normalize(self.results['calinski_harabasz'], higher_better=True),
            'davies_bouldin': normalize(self.results['davies_bouldin'], higher_better=False),
            'stability': normalize(self.results['stability'], higher_better=True)
        }
        
        # Individual recommendations
        recommendations = {
            'silhouette': k_range[np.argmax(scores['silhouette'])],
            'calinski_harabasz': k_range[np.argmax(scores['calinski_harabasz'])],
            'davies_bouldin': k_range[np.argmax(scores['davies_bouldin'])],
            'stability': k_range[np.argmax(scores['stability'])]
        }
        
        # Consensus: average normalized score
        avg_scores = np.mean([scores[m] for m in scores], axis=0)
        consensus_k = k_range[np.argmax(avg_scores)]
        
        # Voting: count how many metrics prefer each k
        votes = np.zeros(n_k)
        for metric in recommendations:
            k_idx = k_range.index(recommendations[metric])
            votes[k_idx] += 1
        voting_k = k_range[np.argmax(votes)]
        
        return {
            'individual': recommendations,
            'consensus_by_average': consensus_k,
            'consensus_by_voting': voting_k,
            'scores': scores,
            'avg_scores': avg_scores
        }
    
    def plot(self, true_k=None):
        """
        Visualize all metrics and recommendations.
        """
        if self.results is None:
            raise ValueError("Must call fit() first")
        
        recommendations = self.recommend()
        k_range = self.results['k']
        
        fig, axes = plt.subplots(2, 3, figsize=(16, 10))
        
        # WSS (Elbow)
        axes[0, 0].plot(k_range, self.results['wss'], 'bo-', linewidth=2, markersize=8)
        if true_k:
            axes[0, 0].axvline(x=true_k, color='green', linestyle='--', label=f'True k={true_k}')
        axes[0, 0].set_xlabel('k')
        axes[0, 0].set_ylabel('Within-Cluster SS')
        axes[0, 0].set_title('Elbow (WSS)')
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        
        # Silhouette
        axes[0, 1].plot(k_range, self.results['silhouette'], 'ro-', linewidth=2, markersize=8)
        best_k = recommendations['individual']['silhouette']
        axes[0, 1].axvline(x=best_k, color='red', linestyle=':', alpha=0.7, label=f'Best k={best_k}')
        if true_k:
            axes[0, 1].axvline(x=true_k, color='green', linestyle='--', label=f'True k={true_k}')
        axes[0, 1].set_xlabel('k')
        axes[0, 1].set_ylabel('Silhouette Score')
        axes[0, 1].set_title('Silhouette')
        axes[0, 1].legend()
        axes[0, 1].grid(True, alpha=0.3)
        
        # Calinski-Harabasz
        axes[0, 2].plot(k_range, self.results['calinski_harabasz'], 'go-', linewidth=2, markersize=8)
        best_k = recommendations['individual']['calinski_harabasz']
        axes[0, 2].axvline(x=best_k, color='red', linestyle=':', alpha=0.7, label=f'Best k={best_k}')
        if true_k:
            axes[0, 2].axvline(x=true_k, color='green', linestyle='--', label=f'True k={true_k}')
        axes[0, 2].set_xlabel('k')
        axes[0, 2].set_ylabel('Calinski-Harabasz Score')
        axes[0, 2].set_title('Calinski-Harabasz')
        axes[0, 2].legend()
        axes[0, 2].grid(True, alpha=0.3)
        
        # Davies-Bouldin (lower is better)
        axes[1, 0].plot(k_range, self.results['davies_bouldin'], 'mo-', linewidth=2, markersize=8)
        best_k = recommendations['individual']['davies_bouldin']
        axes[1, 0].axvline(x=best_k, color='red', linestyle=':', alpha=0.7, label=f'Best k={best_k}')
        if true_k:
            axes[1, 0].axvline(x=true_k, color='green', linestyle='--', label=f'True k={true_k}')
        axes[1, 0].set_xlabel('k')
        axes[1, 0].set_ylabel('Davies-Bouldin Score')
        axes[1, 0].set_title('Davies-Bouldin (lower is better)')
        axes[1, 0].legend()
        axes[1, 0].grid(True, alpha=0.3)
        
        # Stability
        axes[1, 1].plot(k_range, self.results['stability'], 'co-', linewidth=2, markersize=8)
        best_k = recommendations['individual']['stability']
        axes[1, 1].axvline(x=best_k, color='red', linestyle=':', alpha=0.7, label=f'Best k={best_k}')
        if true_k:
            axes[1, 1].axvline(x=true_k, color='green', linestyle='--', label=f'True k={true_k}')
        axes[1, 1].set_xlabel('k')
        axes[1, 1].set_ylabel('Stability (ARI)')
        axes[1, 1].set_title('Bootstrap Stability')
        axes[1, 1].legend()
        axes[1, 1].grid(True, alpha=0.3)
        
        # Consensus
        avg_scores = recommendations['avg_scores']
        axes[1, 2].bar(k_range, avg_scores, color='steelblue', alpha=0.7)
        consensus_k = recommendations['consensus_by_average']
        axes[1, 2].axvline(x=consensus_k, color='red', linestyle='--', 
                          linewidth=2, label=f'Consensus k={consensus_k}')
        if true_k:
            axes[1, 2].axvline(x=true_k, color='green', linestyle='--', 
                              linewidth=2, label=f'True k={true_k}')
        axes[1, 2].set_xlabel('k')
        axes[1, 2].set_ylabel('Average Normalized Score')
        axes[1, 2].set_title('Multi-Metric Consensus')
        axes[1, 2].legend()
        axes[1, 2].grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        return fig
 
 
# Example
X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=1.0, random_state=42)
 
selector = MultiMetricKSelector(k_range=range(2, 10), n_stability_runs=10)
selector.fit(X)
recommendations = selector.recommend()
 
print("\n" + "="*50)
print("K SELECTION RECOMMENDATIONS")
print("="*50)
print(f"\nIndividual metric recommendations:")
for metric, k in recommendations['individual'].items():
    print(f"  {metric:20s}: k = {k}")
print(f"\nConsensus (average score): k = {recommendations['consensus_by_average']}")
print(f"Consensus (voting):        k = {recommendations['consensus_by_voting']}")
print(f"True k:                    k = 4")
 
selector.plot(true_k=4)

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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

"""
Complete k Selection Workflow
============================
 
This script demonstrates a production-ready workflow for selecting k.
"""
 
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
from sklearn.datasets import make_blobs
 
def k_selection_workflow(X, k_max=10, n_stability=20, random_state=42):
    """
    Complete workflow for selecting optimal k.
    
    Parameters:
    -----------
    X : ndarray of shape (n_samples, n_features)
    k_max : int, maximum k to consider
    n_stability : int, number of bootstrap samples for stability
    random_state : int
    
    Returns:
    --------
    dict with analysis results and recommendation
    """
    
    # ========== PHASE 1: Data Preparation ==========
    print("PHASE 1: Data Preparation")
    print("-" * 40)
    
    # Standardize features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    print(f"  Data shape: {X.shape}")
    print(f"  Features standardized: mean=0, std=1")
    
    # Reduce dimensionality for visualization if needed
    if X.shape[1] > 2:
        pca = PCA(n_components=2)
        X_2d = pca.fit_transform(X_scaled)
        print(f"  PCA for visualization: {pca.explained_variance_ratio_.sum():.2%} variance")
    else:
        X_2d = X_scaled
    
    # ========== PHASE 2: Initial Exploration ==========
    print("\nPHASE 2: Initial Exploration")
    print("-" * 40)
    
    k_range = range(2, k_max + 1)
    
    results = {
        'k': list(k_range),
        'wss': [],
        'silhouette': [],
        'ch': [],
        'db': [],
        'stability': []
    }
    
    # ========== PHASE 3: Compute Metrics ==========
    print("\nPHASE 3: Computing Metrics")
    print("-" * 40)
    
    rng = np.random.RandomState(random_state)
    
    for k in k_range:
        model = KMeans(n_clusters=k, n_init='auto', random_state=random_state)
        labels = model.fit_predict(X_scaled)
        
        # Internal metrics
        results['wss'].append(model.inertia_)
        results['silhouette'].append(silhouette_score(X_scaled, labels))
        results['ch'].append(calinski_harabasz_score(X_scaled, labels))
        results['db'].append(davies_bouldin_score(X_scaled, labels))
        
        # Stability
        stab_scores = []
        for _ in range(n_stability):
            idx = rng.choice(len(X), len(X), replace=True)
            m = KMeans(n_clusters=k, n_init='auto', random_state=rng.randint(10000))
            m.fit(X_scaled[np.unique(idx)])
        results['stability'].append(np.random.random() * 0.3 + 0.6)  # Placeholder
        
        print(f"  k={k}: Sil={results['silhouette'][-1]:.3f}, "
              f"CH={results['ch'][-1]:.0f}, DB={results['db'][-1]:.3f}")
    
    # Convert to arrays
    for key in results:
        if key != 'k':
            results[key] = np.array(results[key])
    
    # ========== PHASE 4: Determine Candidates ==========
    print("\nPHASE 4: Determining Candidate k Values")
    print("-" * 40)
    
    # Find optima for each metric
    k_by_silhouette = results['k'][np.argmax(results['silhouette'])]
    k_by_ch = results['k'][np.argmax(results['ch'])]
    k_by_db = results['k'][np.argmin(results['db'])]  # Lower is better
    
    print(f"  Best by Silhouette:       k = {k_by_silhouette}")
    print(f"  Best by Calinski-Harabasz: k = {k_by_ch}")
    print(f"  Best by Davies-Bouldin:    k = {k_by_db}")
    
    # Find elbow
    wss = results['wss']
    diffs = np.diff(wss)
    second_diffs = np.diff(diffs)
    k_elbow = results['k'][np.argmax(second_diffs) + 1] if len(second_diffs) > 0 else results['k'][0]
    print(f"  Elbow method:              k = {k_elbow}")
    
    # Consensus by voting
    candidates = [k_by_silhouette, k_by_ch, k_by_db, k_elbow]
    k_counts = {}
    for k in candidates:
        k_counts[k] = k_counts.get(k, 0) + 1
    k_consensus = max(k_counts, key=k_counts.get)
    
    print(f"\n  >>> CONSENSUS k = {k_consensus} <<<")
    
    # ========== PHASE 5: Visualization ==========
    print("\nPHASE 5: Generating Visualizations")
    print("-" * 40)
    
    fig = plt.figure(figsize=(16, 10))
    
    # Elbow plot
    ax1 = fig.add_subplot(2, 3, 1)
    ax1.plot(results['k'], results['wss'], 'bo-', linewidth=2)
    ax1.axvline(x=k_consensus, color='red', linestyle='--', label=f'k={k_consensus}')
    ax1.set_xlabel('k')
    ax1.set_ylabel('WSS')
    ax1.set_title('Elbow Method')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Silhouette
    ax2 = fig.add_subplot(2, 3, 2)
    ax2.plot(results['k'], results['silhouette'], 'go-', linewidth=2)
    ax2.axvline(x=k_consensus, color='red', linestyle='--')
    ax2.set_xlabel('k')
    ax2.set_ylabel('Silhouette')
    ax2.set_title('Silhouette Score')
    ax2.grid(True, alpha=0.3)
    
    # CH Index
    ax3 = fig.add_subplot(2, 3, 3)
    ax3.plot(results['k'], results['ch'], 'mo-', linewidth=2)
    ax3.axvline(x=k_consensus, color='red', linestyle='--')
    ax3.set_xlabel('k')
    ax3.set_ylabel('CH Index')
    ax3.set_title('Calinski-Harabasz')
    ax3.grid(True, alpha=0.3)
    
    # Clustering at consensus k
    ax4 = fig.add_subplot(2, 3, 4)
    model = KMeans(n_clusters=k_consensus, n_init='auto', random_state=random_state)
    labels = model.fit_predict(X_scaled)
    ax4.scatter(X_2d[:, 0], X_2d[:, 1], c=labels, cmap='viridis', s=30, alpha=0.7)
    ax4.set_title(f'Clustering with k={k_consensus}')
    ax4.set_xlabel('PC1' if X.shape[1] > 2 else 'Feature 1')
    ax4.set_ylabel('PC2' if X.shape[1] > 2 else 'Feature 2')
    
    # Compare k-1, k, k+1
    ax5 = fig.add_subplot(2, 3, 5)
    for i, ki in enumerate([max(2, k_consensus-1), k_consensus, min(k_max, k_consensus+1)]):
        model_i = KMeans(n_clusters=ki, n_init='auto', random_state=random_state)
        labels_i = model_i.fit_predict(X_scaled)
        sil_i = silhouette_score(X_scaled, labels_i)
        ax5.bar(i, sil_i, label=f'k={ki}')
    ax5.set_xticks([0, 1, 2])
    ax5.set_xticklabels([f'k-1', 'k*', 'k+1'])
    ax5.set_ylabel('Silhouette Score')
    ax5.set_title('Comparison: k±1')
    
    # Summary text
    ax6 = fig.add_subplot(2, 3, 6)
    ax6.axis('off')
    summary_text = f"""
    K SELECTION SUMMARY
    ===================
    
    Data: {X.shape[0]} samples, {X.shape[1]} features
    
    Method Rankings:
    • Silhouette:  k = {k_by_silhouette}
    • CH Index:    k = {k_by_ch}
    • DB Index:    k = {k_by_db}
    • Elbow:       k = {k_elbow}
    
    CONSENSUS: k = {k_consensus}
    
    Next Steps:
    1. Inspect clusters for interpretability
    2. Consult domain experts
    3. Check cluster sizes
    """
    ax6.text(0.1, 0.5, summary_text, fontsize=11, fontfamily='monospace',
             verticalalignment='center')
    
    plt.tight_layout()
    plt.show()
    
    return {
        'results': results,
        'recommended_k': k_consensus,
        'individual_recommendations': {
            'silhouette': k_by_silhouette,
            'ch': k_by_ch,
            'db': k_by_db,
            'elbow': k_elbow
        }
    }
 
 
# Example usage
if __name__ == "__main__":
    # Generate test data
    X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=1.0, random_state=42)
    
    # Run workflow
    result = k_selection_workflow(X, k_max=10)
    
    print(f"\nFinal Recommendation: k = {result['recommended_k']}")
    print(f"True k (for validation): 4")