icaoberg / Statistics - Clustering

Created Sun, 17 May 2026 00:00:00 +0000 Modified Sun, 17 May 2026 17:14:32 -0400

Jupyter Notebook python

statistics clustering k-means dbscan hierarchical gmm scikit-learn

Statistics - Clustering¶

Clustering is an unsupervised learning technique that groups similar observations together without using predefined labels. It is one of the most widely used exploratory tools in data science, biology, marketing, image analysis, and many other fields.

This notebook covers:

Definition and background
When clustering is useful
A real-life example: customer segmentation
Common clustering algorithms — K-Means, Hierarchical, DBSCAN, Gaussian Mixture Models
Choosing the right algorithm
Evaluating clustering quality

Definition¶

According to the NIST/SEMATECH e-Handbook of Statistical Methods:

"Cluster analysis is a method of multivariate analysis in which observations are classified into groups (clusters) on the basis of a set of measured characteristics. The goal is to form clusters such that observations within each cluster are more similar to each other than they are to observations in other clusters."

Source: NIST/SEMATECH e-Handbook of Statistical Methods, Section 7.5 — Cluster Analysis

Key concepts¶

Intra-cluster similarity: observations within a cluster should be as similar as possible.
Inter-cluster dissimilarity: clusters should be as different from each other as possible.
Clustering is unsupervised — no labeled output is used during fitting.
The "right" number of clusters is often unknown and must be inferred from the data.

When Is Clustering Useful?¶

Scenario	Example
Customer segmentation	Group shoppers by purchase behavior for targeted marketing
Document/topic grouping	Cluster news articles by topic without predefined categories
Anomaly detection	Points far from every cluster may be outliers or fraud
Image segmentation	Group pixels by color/texture to identify regions in an image
Gene expression analysis	Find co-regulated gene groups in transcriptomics data
Spatial analysis	Identify geographic hotspots of disease, crime, or sales

Clustering is less useful when:

The data has no natural grouping structure.
The number of clusters must be exact and is unknown (some algorithms require it).
You have labeled data — supervised classification will typically outperform clustering.

Setup¶

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_blobs, make_moons
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn.metrics import silhouette_score, davies_bouldin_score
from scipy.cluster.hierarchy import dendrogram, linkage

sns.set_theme(style='whitegrid')
plt.rcParams['figure.dpi'] = 120
rng = np.random.default_rng(42)

print('Libraries loaded.')

Libraries loaded.

Real-Life Example: Customer Segmentation¶

Imagine a retail company that wants to understand its customer base better. They have collected two measurements for each customer:

Annual income (thousands USD)
Annual spending score (0–100, higher = spends more)

The company wants to segment customers into groups to tailor marketing campaigns. We will simulate this dataset and apply multiple clustering algorithms, comparing their behavior.

In [2]:

# Simulate customer data: 5 natural segments
centers = [
    [25, 80],   # low income, high spender
    [55, 55],   # mid income, mid spender
    [85, 85],   # high income, high spender
    [85, 20],   # high income, low spender
    [30, 20],   # low income, low spender
]
X, y_true = make_blobs(n_samples=300, centers=centers, cluster_std=6.0, random_state=42)
X[:, 0] = np.clip(X[:, 0], 15, 135)   # income range
X[:, 1] = np.clip(X[:, 1], 1, 99)      # spending score range

df = pd.DataFrame(X, columns=['Annual Income (k$)', 'Spending Score (1-100)'])

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(df['Annual Income (k$)'], df['Spending Score (1-100)'],
           alpha=0.7, edgecolors='white', s=60)
ax.set_xlabel('Annual Income (k$)')
ax.set_ylabel('Spending Score (1-100)')
ax.set_title('Customer Data — No Labels')
plt.tight_layout()
plt.show()

print(f'Dataset shape: {df.shape}')

No description has been provided for this image

Dataset shape: (300, 2)

Standardize Features¶

Distance-based clustering algorithms (K-Means, hierarchical, DBSCAN) are sensitive to feature scale. We standardize to zero mean and unit variance.

In [3]:

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
print('Features standardized.')

Features standardized.

Algorithm 1 — K-Means¶

K-Means partitions data into k clusters by iteratively assigning each point to its nearest centroid and recomputing centroids. It minimizes within-cluster sum of squares (WCSS).

Best for: spherical, similarly-sized clusters with a known k.

Choosing k with the Elbow Method¶

Plot WCSS against k and look for the "elbow" — the point where additional clusters stop providing meaningful improvement.

In [4]:

wcss = []
k_range = range(1, 11)
for k in k_range:
    km = KMeans(n_clusters=k, random_state=42, n_init='auto')
    km.fit(X_scaled)
    wcss.append(km.inertia_)

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(k_range, wcss, marker='o', color='steelblue')
ax.axvline(5, linestyle='--', color='red', alpha=0.6, label='k=5 (elbow)')
ax.set_xlabel('Number of Clusters (k)')
ax.set_ylabel('WCSS (Inertia)')
ax.set_title('Elbow Method for Optimal k')
ax.legend()
plt.tight_layout()
plt.show()

In [5]:

km5 = KMeans(n_clusters=5, random_state=42, n_init='auto')
labels_kmeans = km5.fit_predict(X_scaled)

fig, ax = plt.subplots(figsize=(7, 5))
scatter = ax.scatter(X[:, 0], X[:, 1], c=labels_kmeans,
                     cmap='tab10', alpha=0.8, edgecolors='white', s=60)
centers_orig = scaler.inverse_transform(km5.cluster_centers_)
ax.scatter(centers_orig[:, 0], centers_orig[:, 1],
           c='black', marker='X', s=200, zorder=5, label='Centroids')
ax.set_xlabel('Annual Income (k$)')
ax.set_ylabel('Spending Score (1-100)')
ax.set_title('K-Means Clustering (k=5)')
ax.legend()
plt.tight_layout()
plt.show()

sil = silhouette_score(X_scaled, labels_kmeans)
print(f'Silhouette Score (K-Means, k=5): {sil:.3f}')

Silhouette Score (K-Means, k=5): 0.749

Algorithm 2 — Hierarchical (Agglomerative) Clustering¶

Hierarchical clustering builds a tree of clusters (dendrogram) by successively merging the two closest clusters (agglomerative) or splitting clusters (divisive). No need to specify k in advance — you cut the dendrogram at a desired level.

Best for: exploring the natural hierarchy in data; varying cluster sizes; small-to-medium datasets.

In [6]:

Z = linkage(X_scaled, method='ward')

fig, ax = plt.subplots(figsize=(10, 4))
dendrogram(Z, ax=ax, truncate_mode='level', p=5,
           color_threshold=0.7 * max(Z[:, 2]))
ax.set_title('Dendrogram — Ward Linkage (truncated)')
ax.set_xlabel('Sample index (or cluster size)')
ax.set_ylabel('Distance')
ax.axhline(y=11, linestyle='--', color='red', alpha=0.6, label='Cut → 5 clusters')
ax.legend()
plt.tight_layout()
plt.show()

In [7]:

agg = AgglomerativeClustering(n_clusters=5, linkage='ward')
labels_agg = agg.fit_predict(X_scaled)

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(X[:, 0], X[:, 1], c=labels_agg,
           cmap='tab10', alpha=0.8, edgecolors='white', s=60)
ax.set_xlabel('Annual Income (k$)')
ax.set_ylabel('Spending Score (1-100)')
ax.set_title('Agglomerative Clustering (k=5, Ward)')
plt.tight_layout()
plt.show()

sil = silhouette_score(X_scaled, labels_agg)
print(f'Silhouette Score (Agglomerative, k=5): {sil:.3f}')

Silhouette Score (Agglomerative, k=5): 0.749

Algorithm 3 — DBSCAN¶

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) groups points that are closely packed together and marks outliers (low-density regions) as noise. It does not require specifying k and can find arbitrarily-shaped clusters.

Key parameters:

eps: maximum distance between two points to be considered neighbors.
min_samples: minimum points in a neighborhood to form a core point.

Best for: irregular cluster shapes, datasets with noise/outliers, unknown number of clusters.

To demonstrate DBSCAN's advantage, we use a crescent moon dataset where K-Means fails.

In [8]:

# DBSCAN shines on non-convex shapes
X_moons, _ = make_moons(n_samples=300, noise=0.08, random_state=42)
X_moons_scaled = StandardScaler().fit_transform(X_moons)

km_moons = KMeans(n_clusters=2, random_state=42, n_init='auto').fit_predict(X_moons_scaled)
db_moons = DBSCAN(eps=0.2, min_samples=5).fit_predict(X_moons_scaled)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].scatter(X_moons[:, 0], X_moons[:, 1], c=km_moons,
                cmap='tab10', alpha=0.8, edgecolors='white', s=50)
axes[0].set_title('K-Means on Moons (fails — assumes spherical clusters)')

axes[1].scatter(X_moons[:, 0], X_moons[:, 1], c=db_moons,
                cmap='tab10', alpha=0.8, edgecolors='white', s=50)
axes[1].set_title('DBSCAN on Moons (succeeds — density-based)')

for ax in axes:
    ax.set_xlabel('x'); ax.set_ylabel('y')

plt.tight_layout()
plt.show()

n_noise = (db_moons == -1).sum()
print(f'DBSCAN: {len(set(db_moons)) - (1 if -1 in db_moons else 0)} clusters, {n_noise} noise points')

DBSCAN: 7 clusters, 10 noise points

In [9]:

# Apply DBSCAN to customer data
db_cust = DBSCAN(eps=0.45, min_samples=8).fit_predict(X_scaled)
n_clusters_db = len(set(db_cust)) - (1 if -1 in db_cust else 0)
n_noise_db = (db_cust == -1).sum()

fig, ax = plt.subplots(figsize=(7, 5))
scatter = ax.scatter(X[:, 0], X[:, 1], c=db_cust,
                     cmap='tab10', alpha=0.8, edgecolors='white', s=60)
ax.set_xlabel('Annual Income (k$)')
ax.set_ylabel('Spending Score (1-100)')
ax.set_title(f'DBSCAN — Customer Data ({n_clusters_db} clusters, {n_noise_db} noise points)')
plt.tight_layout()
plt.show()

if n_clusters_db > 1:
    mask = db_cust != -1
    sil = silhouette_score(X_scaled[mask], db_cust[mask])
    print(f'Silhouette Score (DBSCAN, no noise points): {sil:.3f}')

Silhouette Score (DBSCAN, no noise points): 0.749

Algorithm 4 — Gaussian Mixture Models (GMM)¶

GMM assumes the data is generated from a mixture of k Gaussian distributions. Unlike K-Means (hard assignment), GMM produces soft assignments — each point has a probability of belonging to each cluster. Fitted via the Expectation-Maximization (EM) algorithm.

Best for: elliptical clusters; when you need probabilistic cluster membership; clusters with different sizes/shapes.

In [10]:

gmm = GaussianMixture(n_components=5, random_state=42)
labels_gmm = gmm.fit_predict(X_scaled)
proba = gmm.predict_proba(X_scaled).max(axis=1)   # max probability per point

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

axes[0].scatter(X[:, 0], X[:, 1], c=labels_gmm,
                cmap='tab10', alpha=0.8, edgecolors='white', s=60)
axes[0].set_title('GMM — Hard Assignments')
axes[0].set_xlabel('Annual Income (k$)')
axes[0].set_ylabel('Spending Score (1-100)')

sc = axes[1].scatter(X[:, 0], X[:, 1], c=proba,
                     cmap='RdYlGn', alpha=0.9, edgecolors='white', s=60,
                     vmin=0.5, vmax=1.0)
plt.colorbar(sc, ax=axes[1], label='Max cluster probability')
axes[1].set_title('GMM — Assignment Confidence')
axes[1].set_xlabel('Annual Income (k$)')
axes[1].set_ylabel('Spending Score (1-100)')

plt.tight_layout()
plt.show()

sil = silhouette_score(X_scaled, labels_gmm)
print(f'Silhouette Score (GMM, k=5): {sil:.3f}')
print(f'BIC (lower is better): {gmm.bic(X_scaled):.1f}')

Silhouette Score (GMM, k=5): 0.749
BIC (lower is better): 930.9

Evaluating Clustering Quality¶

Because clustering is unsupervised, evaluation metrics do not use ground-truth labels. Two common internal metrics:

Metric	Range	Interpretation
Silhouette Score	−1 to +1	Higher → denser, better-separated clusters
Davies–Bouldin Index	0 to ∞	Lower → better-separated clusters

In [11]:

results = []
for name, labels in [('K-Means', labels_kmeans),
                      ('Agglomerative', labels_agg),
                      ('GMM', labels_gmm)]:
    sil = silhouette_score(X_scaled, labels)
    db_idx = davies_bouldin_score(X_scaled, labels)
    results.append({'Algorithm': name, 'Silhouette ↑': round(sil, 3),
                    'Davies-Bouldin ↓': round(db_idx, 3)})

# DBSCAN: exclude noise
mask = db_cust != -1
if mask.sum() > 1 and len(set(db_cust[mask])) > 1:
    sil = silhouette_score(X_scaled[mask], db_cust[mask])
    db_idx = davies_bouldin_score(X_scaled[mask], db_cust[mask])
    results.append({'Algorithm': 'DBSCAN (excl. noise)', 'Silhouette ↑': round(sil, 3),
                    'Davies-Bouldin ↓': round(db_idx, 3)})

pd.DataFrame(results).set_index('Algorithm')

Out[11]:

	Silhouette ↑	Davies-Bouldin ↓
Algorithm
K-Means	0.749	0.346
Agglomerative	0.749	0.351
GMM	0.749	0.346
DBSCAN (excl. noise)	0.749	0.346

Choosing the Right Clustering Algorithm¶

Algorithm	Requires k?	Handles noise?	Cluster shapes	Scalability
K-Means	Yes	No	Spherical	Very fast (O(n))
Agglomerative	No (cut dendrogram)	No	Any (with right linkage)	Slow (O(n²))
DBSCAN	No	Yes	Arbitrary	Fast with spatial index
GMM	Yes	Soft (via low probability)	Elliptical	Moderate

Rule of thumb:

Start with K-Means — fast, simple, widely understood.
Use DBSCAN if you expect noise or non-spherical clusters.
Use GMM if you need probabilistic membership or elliptical clusters.
Use hierarchical if you want to explore cluster structure at multiple granularities.

Summary¶

We demonstrated four clustering algorithms on a customer segmentation task:

Step	What we did
Data generation	Simulated 300 customers with income + spending score
K-Means	Elbow method identified k=5; fast and clean segmentation
Agglomerative	Dendrogram revealed hierarchical structure
DBSCAN	Identified non-convex shapes; handled noise automatically
GMM	Provided soft probabilistic assignments
Evaluation	Compared via Silhouette and Davies–Bouldin scores