Skills › Scientific Agent Skills › scikit-learn › clustering_analysis.py
clustering_analysis.py
357 lines · skills/scikit-learn/scripts/clustering_analysis.py at 5f2bc03
1"""
2Clustering analysis example with multiple algorithms, evaluation, and visualization.
3"""
5import numpy as np
6import matplotlib.pyplot as plt
7from sklearn.preprocessing import StandardScaler
8from sklearn.decomposition import PCA
9from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
10from sklearn.mixture import GaussianMixture
11from sklearn.metrics import (
12 silhouette_score, calinski_harabasz_score, davies_bouldin_score
13)
14from sklearn.utils.validation import check_array
17def preprocess_for_clustering(X, scale=True, pca_components=None):
18 """
19 Preprocess data for clustering.
21 Parameters:
22 -----------
23 X : array-like
24 Feature matrix
25 scale : bool
26 Whether to standardize features
27 pca_components : int or None
28 Number of PCA components (None to skip PCA)
30 Returns:
31 --------
32 array
33 Preprocessed data
34 """
35 X_processed = check_array(X, ensure_min_samples=2, copy=True)
37 if scale:
38 scaler = StandardScaler()
39 X_processed = scaler.fit_transform(X_processed)
41 if pca_components is not None:
42 pca = PCA(n_components=pca_components)
43 X_processed = pca.fit_transform(X_processed)
44 print(f"PCA: Explained variance ratio = {pca.explained_variance_ratio_.sum():.3f}")
46 return X_processed
49def find_optimal_k_kmeans(X, k_range=range(2, 11)):
50 """
51 Explore candidate K values with inertia and silhouette; no truth guarantee.
53 Parameters:
54 -----------
55 X : array-like
56 Feature matrix (should be scaled)
57 k_range : range
58 Range of K values to test
60 Returns:
61 --------
62 dict
63 Dictionary with inertia and silhouette scores for each K
64 """
65 X = check_array(X, ensure_min_samples=3)
66 k_range = list(k_range)
67 if not k_range or any(not isinstance(k, (int, np.integer)) or
68 not 2 <= k < len(X) for k in k_range):
69 raise ValueError("k_range must contain integers with 2 <= k < n_samples")
70 inertias = []
71 silhouette_scores = []
73 for k in k_range:
74 kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
75 labels = kmeans.fit_predict(X)
77 inertias.append(kmeans.inertia_)
78 # Duplicate rows can collapse requested clusters; an undefined metric
79 # must not be ranked as evidence for a particular K.
80 score = clustering_metrics(X, labels)['silhouette']
81 silhouette_scores.append(np.nan if score is None else score)
83 if not np.isfinite(silhouette_scores).any():
84 raise ValueError("No candidate produced 2 <= distinct labels < n_samples")
86 # Plot results
87 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
89 # Elbow plot
90 ax1.plot(k_range, inertias, 'bo-')
91 ax1.set_xlabel('Number of clusters (K)')
92 ax1.set_ylabel('Inertia')
93 ax1.set_title('Elbow Method')
94 ax1.grid(True)
96 # Silhouette plot
97 ax2.plot(k_range, silhouette_scores, 'ro-')
98 ax2.set_xlabel('Number of clusters (K)')
99 ax2.set_ylabel('Silhouette Score')
100 ax2.set_title('Silhouette Analysis')
101 ax2.grid(True)
103 plt.tight_layout()
104 plt.savefig('clustering_optimization.png', dpi=300, bbox_inches='tight')
105 print("Saved: clustering_optimization.png")
106 plt.close()
108 # Find best K based on silhouette score
109 best_k = k_range[int(np.nanargmax(silhouette_scores))]
110 print(f"\nCandidate K with highest silhouette (exploratory): {best_k}")
112 return {
113 'k_values': list(k_range),
114 'inertias': inertias,
115 'silhouette_scores': silhouette_scores,
116 'best_k': best_k
117 }
120def clustering_metrics(X, labels):
121 """Score an explicit set of rows, preserving undefined metrics as None."""
122 n_labels = len(np.unique(labels))
123 if not 2 <= n_labels < len(labels):
124 return dict(silhouette=None, calinski_harabasz=None, davies_bouldin=None,
125 metric_reason="Requires 2 <= distinct labels < scored samples")
126 return dict(silhouette=float(silhouette_score(X, labels)),
127 calinski_harabasz=float(calinski_harabasz_score(X, labels)),
128 davies_bouldin=float(davies_bouldin_score(X, labels)),
129 metric_reason=None)
132def compare_clustering_algorithms(X, n_clusters=3):
133 """Compare exploratory Euclidean clusterings on finite numeric features.
135 DBSCAN scores exclude noise and therefore describe a different population;
136 use n_scored/coverage and metric_reason when interpreting its scores.
137 Internal scores alone do not establish scientific clusters.
138 """
139 X = check_array(X, ensure_min_samples=2)
140 if not isinstance(n_clusters, (int, np.integer)) or not 1 <= n_clusters <= len(X):
141 raise ValueError("n_clusters must be an integer between 1 and n_samples")
142 algorithms = {
143 'K-Means': KMeans(n_clusters=n_clusters, random_state=42, n_init=10),
144 'Agglomerative': AgglomerativeClustering(n_clusters=n_clusters, linkage='ward'),
145 'Gaussian Mixture': GaussianMixture(n_components=n_clusters, random_state=42),
146 'DBSCAN': DBSCAN(eps=0.5, min_samples=5),
147 }
148 results = {}
149 for name, algorithm in algorithms.items():
150 labels = algorithm.fit_predict(X)
151 mask = labels != -1
152 n_scored = int(mask.sum())
153 result = dict(labels=labels, n_clusters=len(np.unique(labels[mask])),
154 n_noise=int((~mask).sum()), n_scored=n_scored,
155 coverage=n_scored / len(X),
156 **clustering_metrics(X[mask], labels[mask]))
157 results[name] = result
158 print(f"\n{name}: {result['n_clusters']} clusters; "
159 f"{result['n_noise']} noise; coverage={result['coverage']:.1%}")
160 for metric in ('silhouette', 'calinski_harabasz', 'davies_bouldin'):
161 value = result[metric]
162 print(f" {metric}: {value:.4f}" if value is not None else f" {metric}: undefined")
163 if result['metric_reason']:
164 print(f" {result['metric_reason']}")
165 return results
168def visualize_clusters(X, results, true_labels=None):
169 """
170 Visualize clustering results using PCA for 2D projection.
172 Parameters:
173 -----------
174 X : array-like
175 Feature matrix
176 results : dict
177 Dictionary with clustering results
178 true_labels : array-like or None
179 True labels (if available) for comparison
180 """
181 X = check_array(X, ensure_min_samples=2, ensure_min_features=2)
182 if not results and true_labels is None:
183 raise ValueError("At least one clustering result or true_labels is required")
184 # Reduce to 2D using PCA (display only; scores use the clustering space).
185 pca = PCA(n_components=2)
186 X_2d = pca.fit_transform(X)
188 # Determine number of subplots
189 n_plots = len(results)
190 if true_labels is not None:
191 n_plots += 1
193 n_cols = min(3, n_plots)
194 n_rows = (n_plots + n_cols - 1) // n_cols
196 fig, axes = plt.subplots(n_rows, n_cols, figsize=(5*n_cols, 4*n_rows))
197 if n_plots == 1:
198 axes = np.array([axes])
199 axes = axes.flatten()
201 plot_idx = 0
203 # Plot true labels if available
204 if true_labels is not None:
205 ax = axes[plot_idx]
206 scatter = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=true_labels, cmap='viridis', alpha=0.6)
207 ax.set_title('True Labels')
208 ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
209 ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
210 plt.colorbar(scatter, ax=ax)
211 plot_idx += 1
213 # Plot clustering results
214 for name, result in results.items():
215 ax = axes[plot_idx]
216 labels = result['labels']
218 scatter = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=labels, cmap='viridis', alpha=0.6)
220 # Highlight noise points for DBSCAN
221 if name == 'DBSCAN' and -1 in labels:
222 noise_mask = labels == -1
223 ax.scatter(X_2d[noise_mask, 0], X_2d[noise_mask, 1],
224 c='red', marker='x', s=100, label='Noise', alpha=0.8)
225 ax.legend()
227 title = f"{name} (K={result['n_clusters']})"
228 if result.get('silhouette') is not None:
229 title += f"\nSilhouette: {result['silhouette']:.3f}"
230 title += f"\nCoverage: {result.get('coverage', 1.0):.0%}"
231 ax.set_title(title)
232 ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
233 ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
234 plt.colorbar(scatter, ax=ax)
236 plot_idx += 1
238 # Hide unused subplots
239 for idx in range(plot_idx, len(axes)):
240 axes[idx].axis('off')
242 plt.tight_layout()
243 plt.savefig('clustering_results.png', dpi=300, bbox_inches='tight')
244 print("\nSaved: clustering_results.png")
245 plt.close()
248def complete_clustering_analysis(X, true_labels=None, scale=True,
249 find_k=True, k_range=range(2, 11), n_clusters=3):
250 """
251 Complete clustering analysis workflow.
253 Parameters:
254 -----------
255 X : array-like
256 Feature matrix
257 true_labels : array-like or None
258 True labels (for comparison only, not used in clustering)
259 scale : bool
260 Whether to scale features
261 find_k : bool
262 Whether to search for optimal K
263 k_range : range
264 Range of K values to test
265 n_clusters : int
266 Number of clusters to use in comparison
268 Returns:
269 --------
270 dict
271 Dictionary with all analysis results
272 """
273 print("="*60)
274 print("Clustering Analysis")
275 print("="*60)
276 print(f"Data shape: {X.shape}")
278 # Preprocess data
279 X_processed = preprocess_for_clustering(X, scale=scale)
281 # Find optimal K if requested
282 optimization_results = None
283 if find_k:
284 print("\n" + "="*60)
285 print("Finding Optimal Number of Clusters")
286 print("="*60)
287 optimization_results = find_optimal_k_kmeans(X_processed, k_range=k_range)
289 # Use recommended K
290 if optimization_results:
291 n_clusters = optimization_results['best_k']
293 # Compare clustering algorithms
294 comparison_results = compare_clustering_algorithms(X_processed, n_clusters=n_clusters)
296 # Visualize results
297 print("\n" + "="*60)
298 print("Visualizing Results")
299 print("="*60)
300 visualize_clusters(X_processed, comparison_results, true_labels=true_labels)
302 return {
303 'X_processed': X_processed,
304 'optimization': optimization_results,
305 'comparison': comparison_results
306 }
309# Example usage
310if __name__ == "__main__":
311 from sklearn.datasets import load_iris, make_blobs
313 print("="*60)
314 print("Example 1: Iris Dataset")
315 print("="*60)
317 # Load Iris dataset
318 iris = load_iris()
319 X_iris = iris.data
320 y_iris = iris.target
322 results_iris = complete_clustering_analysis(
323 X_iris,
324 true_labels=y_iris,
325 scale=True,
326 find_k=True,
327 k_range=range(2, 8),
328 n_clusters=3
329 )
331 print("\n" + "="*60)
332 print("Example 2: Synthetic Dataset with Noise")
333 print("="*60)
335 # Create synthetic dataset
336 X_synth, y_synth = make_blobs(
337 n_samples=500, n_features=2, centers=4,
338 cluster_std=0.5, random_state=42
339 )
341 # Add noise points
342 noise = np.random.default_rng(42).normal(size=(50, 2)) * 3
343 X_synth = np.vstack([X_synth, noise])
344 y_synth_with_noise = np.concatenate([y_synth, np.full(50, -1)])
346 results_synth = complete_clustering_analysis(
347 X_synth,
348 true_labels=y_synth_with_noise,
349 scale=True,
350 find_k=True,
351 k_range=range(2, 8),
352 n_clusters=4
353 )
355 print("\n" + "="*60)
356 print("Analysis Complete!")
357 print("="*60)