Clustering Algorithms

How do you choose the right number of clusters? The Elbow method and Silhouette score are two essential techniques that help you find the optimal K value for your clustering analysis.

The Challenge of Choosing K

Here's the tricky part about K-Means: you have to tell it how many groups to create BEFORE it runs. But wait — if you're using clustering because you don't know the structure of your data, how are you supposed to know how many groups exist? It's like being asked "How many colors should I use?" before you've even seen the painting!

This is one of the biggest challenges in clustering. Get the number wrong, and you'll get misleading results:

The good news? We have data-driven methods to help us find the "sweet spot" — the optimal K where clusters are meaningful and well-separated. These techniques evaluate clustering quality at different K values, so you're not just guessing. Let's explore the two most popular methods: the Elbow Method and Silhouette Score.

The Elbow Method

The Elbow method works by plotting the inertia (a measure of cluster "tightness") against different values of K. Here's the key insight: as you add more clusters, inertia ALWAYS decreases. Why? Because with more clusters, each point gets closer to its assigned center. In the extreme case, if K equals the number of data points, every point IS its own center, and inertia becomes zero!

But here's the catch: we don't want K to equal our data size — that defeats the purpose of clustering! Instead, we want to find the "elbow" point where adding more clusters stops giving us meaningful improvement. It's the sweet spot between too few clusters (high inertia) and too many (overfitting).

# ============================================
# Step 1: Import Libraries and Create Sample Data
# ============================================
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

# Generate sample data with known clusters
X, _ = make_blobs(
    n_samples=300,      # 300 data points
    centers=4,          # We know there are 4 true clusters
    cluster_std=0.60,   # Spread of points
    random_state=42     # For reproducibility
)

# ============================================
# Step 2: Test Different K Values (1 to 10)
# ============================================
print("Testing different K values...")
print("=" * 50)

# We'll store the inertia (error) for each K
inertias = []
k_range = range(1, 11)  # Test K from 1 to 10

for k in k_range:
    # Create K-Means with current K value
    kmeans = KMeans(
        n_clusters=k,       # Number of clusters to create
        random_state=42,    # Same starting point each time
        n_init=10           # Try 10 different starting positions
    )
    
    # Train the model on our data
    kmeans.fit(X)
    
    # Save the inertia (how "spread out" points are from centers)
    inertias.append(kmeans.inertia_)

# ============================================
# Step 3: Print Results to See the Pattern
# ============================================
for i, k in enumerate(k_range):
    print(f"K={k:2d}: Inertia={inertias[i]:8.2f}")

print("=" * 50)
print()

# Calculate how much inertia DECREASES when we add one more cluster
print("Inertia reduction when adding one more cluster:")
print("-" * 50)

for i in range(len(inertias) - 1):
    reduction = inertias[i] - inertias[i+1]
    percent = (reduction / inertias[i]) * 100
    
    # Mark significant drops with an arrow
    marker = " ← BIG DROP!" if percent > 30 else ""
    print(f"K={i+1} → {i+2}: Reduction = {reduction:7.2f} ({percent:5.1f}%){marker}")

# ============================================
# Step 4: Find the Elbow Point Automatically
# ============================================
print()
print("→ Look for where reduction % drops significantly!")
print("→ That's your elbow point (optimal K)")
print()

# Simple automatic detection: find biggest "drop in drops"
drops = []
for i in range(len(inertias) - 1):
    drops.append(inertias[i] - inertias[i+1])

# The elbow is where the drop suddenly becomes much smaller
# (i.e., the "second derivative" is largest)
elbow_k = 1
max_drop_change = 0

for i in range(len(drops) - 1):
    drop_change = drops[i] - drops[i+1]  # How much did the improvement slow down?
    if drop_change > max_drop_change:
        max_drop_change = drop_change
        elbow_k = i + 2  # K value at this elbow

print(f"🎯 Detected Elbow Point: K = {elbow_k}")
print(f"   (This is likely the optimal number of clusters)")

Now visualize the elbow curve to find the optimal K:

# ============================================
# Step 5: Visualize the Elbow Curve
# ============================================
plt.figure(figsize=(10, 6))

# Plot the inertia values as a line with dots
plt.plot(k_range, inertias, 'bo-', linewidth=2, markersize=8)

# Add labels and title
plt.xlabel('Number of Clusters (K)', fontsize=12)
plt.ylabel('Inertia (Within-cluster sum of squares)', fontsize=12)
plt.title('Elbow Method for Optimal K', fontsize=14, fontweight='bold')

# Make the x-axis show each K value
plt.xticks(k_range)

# Add grid for easier reading
plt.grid(True, alpha=0.3)

# Mark the elbow point
plt.axvline(x=elbow_k, color='r', linestyle='--', label=f'Elbow at K={elbow_k}')
plt.legend()

plt.tight_layout()
plt.show()

Silhouette Score

While the Elbow method looks at how tight clusters are internally, the Silhouette Score takes a different approach: it measures both how close points are to their own cluster AND how far they are from other clusters. This gives us a more complete picture of clustering quality.

The Silhouette score measures how similar each point is to its own cluster compared to other clusters. Unlike inertia (where lower is better), higher silhouette scores are better, and the score has a meaningful, interpretable range from -1 to +1. This makes it easier to understand than inertia, which can be any positive number.

# ============================================
# Calculate Silhouette Score for Different K Values
# ============================================
from sklearn.metrics import silhouette_score

silhouette_scores = []
k_range = range(2, 11)  # NOTE: Silhouette requires at least 2 clusters (K=1 has no meaning)

print("Testing Silhouette Score for K = 2 to 10...")
print("=" * 50)

for k in k_range:
    # Create and fit K-Means
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    labels = kmeans.fit_predict(X)  # Get cluster labels for each point
    
    # Calculate silhouette score
    score = silhouette_score(X, labels)
    silhouette_scores.append(score)
    
    # Interpret the score
    if score >= 0.71:
        quality = "Excellent ⭐"
    elif score >= 0.51:
        quality = "Good ✓"
    elif score >= 0.26:
        quality = "Fair ⚠"
    else:
        quality = "Poor ✗"
    
    print(f"K={k}: Silhouette Score = {score:.3f} → {quality}")

print("=" * 50)
print(f"\\n→ Best K = {k_range[silhouette_scores.index(max(silhouette_scores))]} "
      f"(highest score: {max(silhouette_scores):.3f})")

Visualize to find the K with highest silhouette score:

# ============================================
# Visualize Silhouette Scores
# ============================================
plt.figure(figsize=(10, 6))

# Plot the scores
plt.plot(k_range, silhouette_scores, 'go-', linewidth=2, markersize=10)

# Highlight the best K
best_k = k_range[silhouette_scores.index(max(silhouette_scores))]
plt.axvline(x=best_k, color='r', linestyle='--', label=f'Best K = {best_k}')

# Labels
plt.xlabel('Number of Clusters (K)', fontsize=12)
plt.ylabel('Silhouette Score', fontsize=12)
plt.title('Silhouette Score for Optimal K\\n(Higher is Better!)', fontsize=14, fontweight='bold')
plt.xticks(k_range)
plt.ylim(0, 1)  # Score range is -1 to 1, but we expect positive values
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

Combining Both Methods

For robust cluster selection, use both methods together. They often agree, but when they differ, consider domain knowledge and the specific use case.

# Combined visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Elbow plot
axes[0].plot(range(1, 11), inertias, 'bo-', linewidth=2)
axes[0].set_xlabel('Number of Clusters (K)')
axes[0].set_ylabel('Inertia')
axes[0].set_title('Elbow Method')
axes[0].grid(True, alpha=0.3)

# Silhouette plot
axes[1].plot(range(2, 11), silhouette_scores, 'go-', linewidth=2)
axes[1].set_xlabel('Number of Clusters (K)')
axes[1].set_ylabel('Silhouette Score')
axes[1].set_title('Silhouette Analysis')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Silhouette Visualization

A silhouette plot shows the score for each individual point, helping you understand cluster quality in more detail:

# Silhouette plot for individual samples
from sklearn.metrics import silhouette_samples
import numpy as np

kmeans = KMeans(n_clusters=4, random_state=42)
labels = kmeans.fit_predict(X)

# Get silhouette values for each sample
sample_silhouette_values = silhouette_samples(X, labels)

# Average silhouette score
avg_score = silhouette_score(X, labels)
print(f"Average Silhouette Score: {avg_score:.3f}")

from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=200, centers=5, random_state=42)
kmeans = KMeans(n_clusters=5, random_state=42)
kmeans.fit(X)

print(f"Inertia: {kmeans.inertia_:.2f}")

import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

inertias = []
for k in range(1, 11):
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(X)
    inertias.append(kmeans.inertia_)

plt.plot(range(1, 11), inertias, 'bo-')
plt.xlabel('K')
plt.ylabel('Inertia')
plt.title('Elbow Method')
plt.show()
# Elbow appears at K=4

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

kmeans = KMeans(n_clusters=4, random_state=42)
labels = kmeans.fit_predict(X)

score = silhouette_score(X, labels)
print(f"Silhouette Score: {score:.3f}")

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

best_k = 2
best_score = -1

for k in range(2, 9):
    kmeans = KMeans(n_clusters=k, random_state=42)
    labels = kmeans.fit_predict(X)
    score = silhouette_score(X, labels)
    print(f"K={k}: Score={score:.3f}")
    if score > best_score:
        best_score = score
        best_k = k

print(f"\nBest K: {best_k} with score {best_score:.3f}")

Hierarchical clustering builds a tree of clusters, letting you explore data at multiple granularity levels. Unlike K-Means, you do not need to specify the number of clusters upfront, giving you more flexibility in exploratory analysis.

How Hierarchical Clustering Works

Agglomerative (bottom-up) hierarchical clustering starts with each data point as its own cluster, then progressively merges the closest clusters until only one remains. This creates a tree structure called a dendrogram that shows how clusters are related.

Imagine you have 100 customers. Initially, each customer is their own "cluster" (100 clusters). The algorithm finds the two most similar customers and merges them (99 clusters). It repeats this process, always merging the two closest clusters, until all customers form one giant cluster. The dendrogram records this entire merge history, letting you "rewind" to any stage and extract the cluster structure you need.

Linkage Methods

The linkage method determines how to measure distance between clusters. Different linkage methods produce different clustering results:

Implementing Hierarchical Clustering

Scikit-learn provides AgglomerativeClustering, and scipy offers tools for dendrograms:

# ============================================
# Hierarchical (Agglomerative) Clustering
# ============================================

# Import required libraries
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs
import numpy as np

print("Setting up Hierarchical Clustering...")
print("=" * 50)

# Generate sample data with 4 natural groups
X, _ = make_blobs(
    n_samples=150,        # 150 data points
    centers=4,            # 4 true clusters
    cluster_std=0.60,     # Moderate spread
    random_state=42       # Reproducibility
)

print(f"Data shape: {X.shape}")  # (150, 2)
print(f"Starting with {len(X)} individual clusters (each point is its own cluster)")
print()

# ============================================
# STEP 1: Apply Agglomerative Clustering
# ============================================

# Create the clustering model
agg_clustering = AgglomerativeClustering(
    n_clusters=4,         # Cut the tree to get 4 final clusters
    linkage='ward'        # Ward minimizes within-cluster variance (good default)
)

print("Applying Agglomerative Clustering...")
print("Process: Start with 150 clusters → Merge closest pairs → Repeat until 4 remain")
print()

# Fit and predict cluster labels
# Under the hood: Builds full dendrogram, then cuts at the right height
labels = agg_clustering.fit_predict(X)

print("✓ Clustering complete!")
print()

# ============================================
# STEP 2: Analyze Results
# ============================================

print(f"Cluster labels: {np.unique(labels)}")  # [0 1 2 3]
print()

# Count points per cluster
for cluster_id in range(4):
    count = np.sum(labels == cluster_id)
    percentage = (count / len(X)) * 100
    print(f"Cluster {cluster_id}: {count} points ({percentage:.1f}%)")

print()
print(f"Total merges performed: {len(X) - 4}")  # Started with 150, ended with 4
print(f"Number of clusters: {agg_clustering.n_clusters_}")

Creating Dendrograms

Dendrograms visualize the hierarchical structure and help you decide where to cut for a specific number of clusters:

# ============================================
# Creating a Dendrogram Visualization
# ============================================

# Import visualization tools
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage

print("Building dendrogram...")
print()

# ============================================
# STEP 1: Calculate Linkage Matrix
# ============================================
# The linkage matrix records the merge history:
# - Which clusters were merged at each step
# - The distance between them at merge time
# - How many original points each new cluster contains

linkage_matrix = linkage(
    X,                # Data to cluster
    method='ward'     # Ward linkage (minimizes variance)
)

print(f"Linkage matrix shape: {linkage_matrix.shape}")  # (n_samples-1, 4)
print(f"Number of merges recorded: {len(linkage_matrix)}")
print()
print("Each row represents one merge:")
print("[cluster1_id, cluster2_id, merge_distance, num_points_in_new_cluster]")
print()
print("First 3 merges:")
print(linkage_matrix[:3])
print()

# ============================================
# STEP 2: Plot Dendrogram
# ============================================

plt.figure(figsize=(12, 6))

# Create the dendrogram plot
# truncate_mode='level' limits display to top 5 levels (for readability)
dendrogram(
    linkage_matrix,         # The merge history
    truncate_mode='level',  # Show only top levels
    p=5                     # Show 5 levels (prevents clutter with large datasets)
)

plt.xlabel('Sample Index (or Cluster Size)', fontsize=12)
plt.ylabel('Distance (Height of Merge)', fontsize=12)
plt.title('Hierarchical Clustering Dendrogram\n(Higher merges = more different clusters)', 
          fontsize=14, fontweight='bold')
plt.axhline(y=10, color='r', linestyle='--', label='Potential cut point')
plt.legend()
plt.grid(True, alpha=0.3, axis='y')
plt.show()

print("✓ Dendrogram created!")
print()
print("How to read it:")
print("- X-axis: Data points or clusters")
print("- Y-axis: Distance at which clusters merge")
print("- Vertical lines: Clusters being merged")
print("- Height of U-shape: How different the clusters are")
print("- Cut horizontally to get desired number of clusters")

You can cut the dendrogram at different heights to get different numbers of clusters:

# Cut dendrogram at specific distance threshold
from scipy.cluster.hierarchy import fcluster

# Get clusters by cutting at distance threshold
distance_threshold = 10
clusters = fcluster(linkage_matrix, t=distance_threshold, criterion='distance')

print(f"Number of clusters at threshold {distance_threshold}: {len(np.unique(clusters))}")

# Or cut to get specific number of clusters
n_clusters = 4
clusters = fcluster(linkage_matrix, t=n_clusters, criterion='maxclust')
print(f"Cluster labels: {np.unique(clusters)}")

Visualizing Cluster Results

# Visualize clustering results
plt.figure(figsize=(12, 5))

# Without predefined n_clusters - explore structure
plt.subplot(1, 2, 1)
agg = AgglomerativeClustering(n_clusters=4, linkage='ward')
labels = agg.fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', alpha=0.6)
plt.title('Agglomerative Clustering (4 clusters)')

# Compare with different n_clusters
plt.subplot(1, 2, 2)
agg = AgglomerativeClustering(n_clusters=6, linkage='ward')
labels = agg.fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', alpha=0.6)
plt.title('Agglomerative Clustering (6 clusters)')

plt.tight_layout()
plt.show()

When to Use Hierarchical Clustering

from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=150, centers=3, random_state=42)

agg = AgglomerativeClustering(n_clusters=3, linkage='ward')
labels = agg.fit_predict(X)

print(f"Labels: {labels[:10]}")

import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=50, centers=3, random_state=42)

linkage_matrix = linkage(X, method='ward')

plt.figure(figsize=(10, 5))
dendrogram(linkage_matrix)
plt.xlabel('Sample Index')
plt.ylabel('Distance')
plt.title('Dendrogram')
plt.show()

import matplotlib.pyplot as plt
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=150, centers=4, random_state=42)
linkages = ['ward', 'complete', 'average']

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for ax, link in zip(axes, linkages):
    agg = AgglomerativeClustering(n_clusters=4, linkage=link)
    labels = agg.fit_predict(X)
    ax.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis')
    ax.set_title(f'Linkage: {link}')

plt.tight_layout()
plt.show()

from scipy.cluster.hierarchy import linkage, fcluster
from sklearn.datasets import make_blobs
import numpy as np

X, _ = make_blobs(n_samples=100, centers=4, random_state=42)

linkage_matrix = linkage(X, method='ward')
clusters = fcluster(linkage_matrix, t=15, criterion='distance')

n_clusters = len(np.unique(clusters))
print(f"Clusters at distance 15: {n_clusters}")

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) finds clusters of arbitrary shapes and automatically detects outliers. It is ideal when your data has non-spherical clusters or contains noise that should not belong to any cluster.

How DBSCAN Works

Unlike K-Means which uses distance to centroids, DBSCAN defines clusters as dense regions of points separated by sparse regions. It groups together points that are closely packed and marks points in low-density regions as outliers.

The algorithm has a simple logic: If you have enough neighbors nearby, you're in a cluster. It doesn't care about cluster shape (spherical, crescent, serpentine) - only about local density. This makes DBSCAN perfect for real-world data where clusters aren't always nice, round blobs.

DBSCAN Parameters

DBSCAN has two main parameters that control cluster formation:

DBSCAN(eps=0.5, min_samples=5)  # Good default for scaled data

Implementing DBSCAN

# ============================================
# DBSCAN - Density-Based Clustering
# ============================================

# Import required libraries
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
import numpy as np
import matplotlib.pyplot as plt

print("Setting up DBSCAN for non-spherical clusters...")
print("=" * 50)

# ============================================
# STEP 1: Generate Non-Spherical Data
# ============================================
# K-Means struggles with moon/crescent shapes
# DBSCAN excels at these!

X, _ = make_moons(
    n_samples=300,        # 300 data points
    noise=0.05,           # Small amount of noise (outliers)
    random_state=42       # Reproducibility
)

print(f"Data shape: {X.shape}")  # (300, 2)
print(f"Data has 2 crescent-shaped clusters (like two moons)")
print()

# ============================================
# STEP 2: Apply DBSCAN Algorithm
# ============================================

# Create DBSCAN model with two key parameters
dbscan = DBSCAN(
    eps=0.2,              # Maximum distance between neighbors (neighborhood radius)
                          # "2 neighbors within 0.2 units are in same neighborhood"
    
    min_samples=5         # Minimum points needed to form a dense region (core point)
                          # "Need at least 5 points nearby to be a 'crowd'"
)

print("DBSCAN Parameters:")
print(f"  eps (neighborhood radius): {dbscan.eps}")
print(f"  min_samples (crowd size): {dbscan.min_samples}")
print()

print("Running DBSCAN...")
print("Process: For each point, count neighbors within eps")
print("         If >= min_samples, mark as core point")
print("         Connect core points to form clusters")
print("         Label remaining points as border or noise")
print()

# Fit and predict cluster labels
# Label -1 means NOISE (outlier)
# Labels 0, 1, 2, ... are cluster IDs
labels = dbscan.fit_predict(X)

print("✓ DBSCAN complete!")
print()

# ============================================
# STEP 3: Analyze Results
# ============================================

# Calculate number of clusters (excluding noise)
# Noise is labeled as -1, so we subtract 1 if it exists
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)

# Count noise points (labeled as -1)
n_noise = list(labels).count(-1)

# Get unique cluster labels
unique_labels = set(labels)
print(f"Unique labels: {sorted(unique_labels)}")
print(f"Clusters found: {n_clusters}")
print(f"Noise points: {n_noise} ({n_noise/len(X)*100:.1f}% of data)")
print()

# Show cluster sizes
for label in sorted(unique_labels):
    if label == -1:
        continue  # Skip noise, already counted
    count = list(labels).count(label)
    percentage = (count / len(X)) * 100
    print(f"Cluster {label}: {count} points ({percentage:.1f}%)")

print()
print(f"Core samples (dense points): {len(dbscan.core_sample_indices_)}")
print(f"Components (connected regions): {dbscan.components_.shape}")

Visualize how DBSCAN handles non-spherical clusters:

# Visualize DBSCAN results
plt.figure(figsize=(12, 5))

# DBSCAN on moon data
plt.subplot(1, 2, 1)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', alpha=0.6)
plt.title(f'DBSCAN: {n_clusters} clusters, {n_noise} noise points')

# Compare with K-Means (struggles with moons)
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=2, random_state=42)
km_labels = kmeans.fit_predict(X)

plt.subplot(1, 2, 2)
plt.scatter(X[:, 0], X[:, 1], c=km_labels, cmap='viridis', alpha=0.6)
plt.title('K-Means (struggles with non-spherical)')

plt.tight_layout()
plt.show()

Finding Optimal eps

The k-distance plot helps find a good eps value. Calculate the distance to the k-th nearest neighbor for each point and look for the "elbow" where distances sharply increase:

# ============================================
# k-Distance Plot for Optimal eps Selection
# ============================================

from sklearn.neighbors import NearestNeighbors
import numpy as np
import matplotlib.pyplot as plt

print("Finding optimal eps using k-distance plot...")
print("=" * 50)

# ============================================
# STEP 1: Calculate k-Nearest Neighbor Distances
# ============================================

# Choose k = min_samples (or min_samples - 1)
# This is the same k we'll use in DBSCAN's min_samples
k = 4

print(f"Calculating distance to {k}th nearest neighbor for each point...")

# Create a nearest neighbors model
neighbors = NearestNeighbors(n_neighbors=k)
neighbors.fit(X)  # Fit on our data

# Calculate distances to k nearest neighbors for each point
# distances shape: (n_samples, k) - distances to each of k neighbors
# indices shape: (n_samples, k) - indices of k nearest neighbors
distances, indices = neighbors.kneighbors(X)

print(f"Distances shape: {distances.shape}")  # (300, 4)
print()

# ============================================
# STEP 2: Extract and Sort k-th Neighbor Distance
# ============================================

# Get distance to the k-th (farthest) neighbor for each point
# [:, k-1] gets the last column (4th neighbor distance)
k_distances = distances[:, k-1]

print(f"Sample k-distances (first 10 points):")
print(k_distances[:10])
print()

# Sort distances in ascending order for visualization
k_distances_sorted = np.sort(k_distances)

print(f"Min k-distance: {k_distances_sorted[0]:.4f}")
print(f"Max k-distance: {k_distances_sorted[-1]:.4f}")
print(f"Median k-distance: {np.median(k_distances_sorted):.4f}")
print()

# ============================================
# STEP 3: Plot k-Distance Graph
# ============================================

plt.figure(figsize=(10, 6))

# Plot sorted k-distances
plt.plot(k_distances_sorted, linewidth=2, color='steelblue')

# Add suggested eps line (visual helper)
# Rule of thumb: eps at the "knee" or elbow of the curve
suggested_eps = 0.2  # This would be determined visually or algorithmically
plt.axhline(y=suggested_eps, color='red', linestyle='--', 
            linewidth=2, label=f'Suggested eps = {suggested_eps}')

plt.xlabel('Points (sorted by k-distance)', fontsize=12)
plt.ylabel(f'{k}th Nearest Neighbor Distance', fontsize=12)
plt.title(f'k-Distance Plot (k={k})\nLook for the elbow point!', 
          fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.legend(fontsize=11)
plt.tight_layout()
plt.show()

print("✓ k-distance plot created!")
print()
print("How to interpret:")
print("- Flat region (left): Points with close neighbors (core/border points)")
print("- Sharp rise (right): Points with distant neighbors (noise points)")
print("- Elbow point: Optimal eps value (where curve bends sharply)")
print()
print(f"Suggested eps for this data: {suggested_eps}")

# The elbow point suggests optimal eps

Handling Noise Points

DBSCAN labels noise points as -1. You can filter them out or analyze them separately:

# Identify noise points
noise_mask = labels == -1
noise_points = X[noise_mask]
clustered_points = X[~noise_mask]

print(f"Total points: {len(X)}")
print(f"Clustered points: {len(clustered_points)}")
print(f"Noise points: {len(noise_points)}")

# Visualize with noise highlighted
plt.figure(figsize=(10, 6))
plt.scatter(clustered_points[:, 0], clustered_points[:, 1], 
            c=labels[~noise_mask], cmap='viridis', alpha=0.6, label='Clustered')
plt.scatter(noise_points[:, 0], noise_points[:, 1], 
            c='red', marker='x', s=50, label='Noise')
plt.legend()
plt.title('DBSCAN with Noise Points Highlighted')
plt.show()

Algorithm Comparison Guide

Decision Flowchart

Follow this visual decision tree to choose the best clustering algorithm for your data. Start at the top and answer each question to find your optimal choice!

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons

X, _ = make_moons(n_samples=200, noise=0.05, random_state=42)

dbscan = DBSCAN(eps=0.3, min_samples=5)
labels = dbscan.fit_predict(X)

print(f"Unique labels: {set(labels)}")

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

dbscan = DBSCAN(eps=0.5, min_samples=5)
labels = dbscan.fit_predict(X)

n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = list(labels).count(-1)

print(f"Clusters: {n_clusters}")
print(f"Noise points: {n_noise}")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import NearestNeighbors
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

neighbors = NearestNeighbors(n_neighbors=5)
neighbors.fit(X)
distances, _ = neighbors.kneighbors(X)

distances = np.sort(distances[:, 4])

plt.figure(figsize=(10, 6))
plt.plot(distances)
plt.xlabel('Points')
plt.ylabel('5th NN Distance')
plt.title('k-Distance Plot')
plt.show()

import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons

X, _ = make_moons(n_samples=300, noise=0.1, random_state=42)
dbscan = DBSCAN(eps=0.2, min_samples=5)
labels = dbscan.fit_predict(X)

noise_mask = labels == -1

plt.figure(figsize=(10, 6))
plt.scatter(X[~noise_mask, 0], X[~noise_mask, 1], 
            c=labels[~noise_mask], cmap='viridis')
plt.scatter(X[noise_mask, 0], X[noise_mask, 1], 
            c='red', marker='x', s=50, label='Noise')
plt.legend()
plt.show()

Evaluating clustering results is challenging because there are no ground truth labels. Learn how to assess cluster quality using internal metrics, visualizations, and domain knowledge to ensure your clusters are meaningful and actionable.

Internal vs External Evaluation

Clustering evaluation falls into two categories depending on whether you have ground truth labels:

Internal Evaluation Metrics

These metrics evaluate cluster quality based only on the data and cluster assignments:

# Import evaluation metrics
from sklearn.metrics import silhouette_score, calinski_harabasz_score
from sklearn.metrics import davies_bouldin_score
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

# Generate sample data
X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

# Fit K-Means
kmeans = KMeans(n_clusters=4, random_state=42)
labels = kmeans.fit_predict(X)

# Calculate internal metrics
silhouette = silhouette_score(X, labels)
calinski = calinski_harabasz_score(X, labels)
davies_bouldin = davies_bouldin_score(X, labels)

print(f"Silhouette Score: {silhouette:.3f}")  # Higher is better (-1 to 1)
print(f"Calinski-Harabasz: {calinski:.2f}")   # Higher is better
print(f"Davies-Bouldin: {davies_bouldin:.3f}")  # Lower is better

External Evaluation Metrics

When you have ground truth labels (for validation or benchmarking), use these metrics:

# External metrics (when true labels are available)
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score
from sklearn.metrics import homogeneity_score, completeness_score, v_measure_score

# Generate data with known labels
X, y_true = make_blobs(n_samples=300, centers=4, random_state=42)

# Cluster the data
kmeans = KMeans(n_clusters=4, random_state=42)
y_pred = kmeans.fit_predict(X)

# Calculate external metrics
ari = adjusted_rand_score(y_true, y_pred)
nmi = normalized_mutual_info_score(y_true, y_pred)
homogeneity = homogeneity_score(y_true, y_pred)
completeness = completeness_score(y_true, y_pred)
v_measure = v_measure_score(y_true, y_pred)

print(f"Adjusted Rand Index: {ari:.3f}")
print(f"Normalized MI: {nmi:.3f}")
print(f"Homogeneity: {homogeneity:.3f}")
print(f"Completeness: {completeness:.3f}")
print(f"V-Measure: {v_measure:.3f}")

Comparing Multiple Algorithms

Evaluate different clustering algorithms on the same data to find the best approach:

# Compare clustering algorithms
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler

# Prepare data
X, _ = make_blobs(n_samples=300, centers=4, random_state=42)
X_scaled = StandardScaler().fit_transform(X)

# Test different algorithms
results = {}

# K-Means
km = KMeans(n_clusters=4, random_state=42)
km_labels = km.fit_predict(X_scaled)
results['K-Means'] = silhouette_score(X_scaled, km_labels)

# Agglomerative
agg = AgglomerativeClustering(n_clusters=4)
agg_labels = agg.fit_predict(X_scaled)
results['Agglomerative'] = silhouette_score(X_scaled, agg_labels)

# DBSCAN
db = DBSCAN(eps=0.5, min_samples=5)
db_labels = db.fit_predict(X_scaled)
if len(set(db_labels)) > 1:  # Need at least 2 clusters for silhouette
    results['DBSCAN'] = silhouette_score(X_scaled, db_labels)

for algo, score in results.items():
    print(f"{algo}: {score:.3f}")

Cluster Interpretation

After clustering, analyze what makes each cluster unique to derive actionable insights:

# Analyze cluster characteristics
import pandas as pd
import numpy as np

# Create sample customer data
np.random.seed(42)
customers = pd.DataFrame({
    'age': np.random.randint(18, 70, 100),
    'income': np.random.randint(20000, 150000, 100),
    'spending_score': np.random.randint(1, 100, 100)
})

# Scale and cluster
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(customers)

kmeans = KMeans(n_clusters=3, random_state=42)
customers['cluster'] = kmeans.fit_predict(X_scaled)

# Analyze cluster profiles
cluster_summary = customers.groupby('cluster').agg({
    'age': ['mean', 'std'],
    'income': ['mean', 'std'],
    'spending_score': ['mean', 'std']
}).round(2)

print("Cluster Profiles:")
print(cluster_summary)

# Visualize cluster profiles with radar/bar charts
import matplotlib.pyplot as plt

cluster_means = customers.groupby('cluster')[['age', 'income', 'spending_score']].mean()

# Normalize for comparison
cluster_normalized = (cluster_means - cluster_means.min()) / (cluster_means.max() - cluster_means.min())

# Bar chart comparison
cluster_normalized.T.plot(kind='bar', figsize=(10, 6))
plt.title('Cluster Profile Comparison (Normalized)')
plt.xlabel('Features')
plt.ylabel('Normalized Value')
plt.legend(title='Cluster', labels=['Young Budget', 'High Spenders', 'Conservative'])
plt.tight_layout()
plt.show()

Best Practices for Cluster Analysis

from sklearn.metrics import silhouette_score, calinski_harabasz_score
from sklearn.metrics import davies_bouldin_score
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)
kmeans = KMeans(n_clusters=4, random_state=42)
labels = kmeans.fit_predict(X)

print(f"Silhouette: {silhouette_score(X, labels):.3f}")
print(f"CH Index: {calinski_harabasz_score(X, labels):.2f}")
print(f"DB Index: {davies_bouldin_score(X, labels):.3f}")

from sklearn.metrics import adjusted_rand_score
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, y_true = make_blobs(n_samples=300, centers=4, random_state=42)

kmeans = KMeans(n_clusters=4, random_state=42)
y_pred = kmeans.fit_predict(X)

ari = adjusted_rand_score(y_true, y_pred)
print(f"Adjusted Rand Index: {ari:.3f}")

from sklearn.cluster import KMeans, AgglomerativeClustering
from sklearn.metrics import silhouette_score
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=4, random_state=42)

km = KMeans(n_clusters=4, random_state=42)
km_labels = km.fit_predict(X)
km_score = silhouette_score(X, km_labels)

agg = AgglomerativeClustering(n_clusters=4)
agg_labels = agg.fit_predict(X)
agg_score = silhouette_score(X, agg_labels)

print(f"K-Means: {km_score:.3f}")
print(f"Agglomerative: {agg_score:.3f}")
print(f"Winner: {'K-Means' if km_score > agg_score else 'Agglomerative'}")

import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
df = pd.DataFrame({
    'age': np.random.randint(20, 60, 100),
    'income': np.random.randint(30000, 120000, 100)
})

scaler = StandardScaler()
X_scaled = scaler.fit_transform(df)

kmeans = KMeans(n_clusters=3, random_state=42)
df['cluster'] = kmeans.fit_predict(X_scaled)

print(df.groupby('cluster').mean().round(2))