Scikit-Learn Clustering

Real-World Scenario

An e-commerce company wants to segment 500,000 customers for targeted marketing. They have purchase history, browsing behavior, and demographics. K-Means segments customers into behavioral groups (high-value loyalists, deal-hunters, casual browsers). DBSCAN identifies unusual shopping patterns as potential fraud. Both are clustering problems.

K-Means Clustering

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
from sklearn.datasets import make_blobs
import pandas as pd

rng = np.random.default_rng(42)

# Simulated customer feature matrix
n_customers = 2000
X = np.column_stack([
    rng.exponential(500, n_customers),        # avg_purchase_value
    rng.poisson(8, n_customers),              # num_purchases
    rng.exponential(30, n_customers),         # days_since_last_purchase
    rng.uniform(0, 1, n_customers),           # email_open_rate
])
feature_names = ["avg_value", "n_purchases", "days_inactive", "email_open_rate"]

# Scale features — K-Means is distance-based; scaling is essential
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# ─── Find optimal K using Silhouette Score ─────────────────────────────────

silhouette_scores = {}
inertias = {}

for k in range(2, 11):
    km = KMeans(n_clusters=k, random_state=42, n_init=10)
    labels = km.fit_predict(X_scaled)
    silhouette_scores[k] = silhouette_score(X_scaled, labels)
    inertias[k] = km.inertia_

print("K vs Silhouette Score:")
for k, score in silhouette_scores.items():
    bar = "█" * int(score * 30)
    print(f"  K={k}: {score:.4f} {bar}")

best_k = max(silhouette_scores, key=silhouette_scores.get)
print(f"\nBest K: {best_k} (silhouette={silhouette_scores[best_k]:.4f})")

# ─── Train final model with best K ────────────────────────────────────────

km_final = KMeans(n_clusters=best_k, random_state=42, n_init=10)
labels   = km_final.fit_predict(X_scaled)

df = pd.DataFrame(X, columns=feature_names)
df["cluster"] = labels

# Cluster profiles — interpret what each cluster represents
profile = df.groupby("cluster").agg({
    "avg_value":       "mean",
    "n_purchases":     "mean",
    "days_inactive":   "mean",
    "email_open_rate": "mean",
})
profile["size"] = df["cluster"].value_counts().sort_index()
print("\nCluster Profiles:")
print(profile.round(1))

DBSCAN — Density-Based Clustering

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_moons

# DBSCAN excels on non-spherical clusters
X, _ = make_moons(n_samples=500, noise=0.1, random_state=42)

# eps: maximum distance between two points in the same neighborhood
# min_samples: minimum points to form a dense region (core point)
db = DBSCAN(eps=0.15, min_samples=5)
labels = db.fit_predict(X)

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

print(f"Clusters found: {n_clusters}")   # should be 2
print(f"Noise points:   {n_noise}")       # outliers labeled -1

# Parameter sensitivity — run a grid to understand eps and min_samples
from sklearn.metrics import silhouette_score

best_score = -1
best_params = {}

for eps in np.arange(0.05, 0.5, 0.05):
    for min_s in [3, 5, 10]:
        db_test = DBSCAN(eps=eps, min_samples=min_s)
        lbls = db_test.fit_predict(X)
        n_c  = len(set(lbls)) - (1 if -1 in lbls else 0)
        if n_c >= 2 and n_c <= 10:
            try:
                score = silhouette_score(X, lbls)
                if score > best_score:
                    best_score  = score
                    best_params = {"eps": eps, "min_samples": min_s}
            except ValueError:
                pass

print(f"\nBest DBSCAN params: {best_params}, Silhouette: {best_score:.4f}")

Hierarchical Clustering

import numpy as np
from sklearn.cluster import AgglomerativeClustering
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import pdist

rng = np.random.default_rng(42)
X = rng.standard_normal((200, 4))

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

# Agglomerative clustering — no need to choose K upfront (use dendrogram to decide)
# Linkage methods: ward (minimizes variance), complete, average, single
for linkage_method in ["ward", "complete", "average"]:
    model = AgglomerativeClustering(n_clusters=4, linkage=linkage_method)
    labels = model.fit_predict(X_scaled)
    score  = silhouette_score(X_scaled, labels)
    print(f"linkage={linkage_method:8s}: silhouette={score:.4f}")

# Reading a dendrogram to choose K
# The linkage matrix encodes the merge history
Z = linkage(X_scaled, method="ward")
# Cut the dendrogram at a threshold that gives meaningful clusters
# Rule of thumb: cut where the gap between merge distances is largest
print("\nTop 10 merge distances (look for a large gap):")
print(np.sort(Z[:, 2])[-10:].round(3))

Evaluating Cluster Quality

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import (
    silhouette_score,
    calinski_harabasz_score,
    davies_bouldin_score,
)
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler

X, y_true = make_blobs(n_samples=500, centers=4, cluster_std=1.5, random_state=42)

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

km = KMeans(n_clusters=4, random_state=42, n_init=10)
labels = km.fit_predict(X_scaled)

# Internal cluster metrics — don't require ground truth labels

# Silhouette: how similar a point is to its cluster vs the nearest other cluster
# Range: [-1, 1], higher is better
sil = silhouette_score(X_scaled, labels)
print(f"Silhouette score:          {sil:.4f}   (higher = better, max 1.0)")

# Calinski-Harabász: ratio of between-cluster to within-cluster dispersion
# Higher is better; no upper bound
ch = calinski_harabasz_score(X_scaled, labels)
print(f"Calinski-Harabász score:   {ch:.1f}  (higher = better)")

# Davies-Bouldin: average similarity of each cluster to its most similar cluster
# Lower is better, minimum is 0
db = davies_bouldin_score(X_scaled, labels)
print(f"Davies-Bouldin index:      {db:.4f}  (lower = better, min 0)")

# Per-sample silhouette to find poorly assigned points
from sklearn.metrics import silhouette_samples
sample_scores = silhouette_samples(X_scaled, labels)
for cluster_id in range(4):
    mask = labels == cluster_id
    cluster_sil = sample_scores[mask].mean()
    n_negative  = (sample_scores[mask] < 0).sum()
    print(f"Cluster {cluster_id}: avg_sil={cluster_sil:.3f}, "
          f"negative={n_negative}/{mask.sum()} (poorly placed points)")

Real-World: Customer Segmentation Pipeline

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

rng = np.random.default_rng(42)
n = 5000

# Simulate customer features
customers = pd.DataFrame({
    "customer_id":     range(n),
    "total_spend":     rng.lognormal(6.5, 1.0, n),
    "num_orders":      rng.poisson(6, n) + 1,
    "avg_order_value": rng.lognormal(4.5, 0.8, n),
    "days_since_last": rng.exponential(45, n),
    "email_open_rate": rng.beta(2, 5, n),
    "return_rate":     rng.beta(1, 10, n),
})

features = ["total_spend", "num_orders", "avg_order_value",
            "days_since_last", "email_open_rate", "return_rate"]

X = customers[features].values
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Choose K via silhouette
best_k, best_sil = 2, -1
for k in range(2, 9):
    sil = silhouette_score(X_scaled, KMeans(k, n_init=5, random_state=42).fit_predict(X_scaled))
    if sil > best_sil:
        best_k, best_sil = k, sil

km = KMeans(n_clusters=best_k, n_init=10, random_state=42)
customers["segment"] = km.fit_predict(X_scaled)

# Name segments based on profiles
profile = customers.groupby("segment")[features].mean()
print(profile.round(1))

# Assign business-friendly names by highest spending segments
spend_rank = profile["total_spend"].rank(ascending=False).astype(int)
segment_names = {
    spend_rank[spend_rank == 1].index[0]: "High-Value Champions",
    spend_rank[spend_rank == 2].index[0]: "Loyal Customers",
    spend_rank[spend_rank == 3].index[0]: "Potential Loyalists",
    spend_rank[spend_rank == 4].index[0]: "At-Risk Customers",
}
customers["segment_name"] = customers["segment"].map(segment_names).fillna("Other")
print("\nSegment distribution:")
print(customers["segment_name"].value_counts())