[Data Science]

Statistics Fundamentals Chapter 7: Unsupervised Learning

March 4, 2026·13 min read·

#Statistics #Unsupervised Learning #PCA #Clustering #Python

Chapters 1–6 covered supervised learning, where labeled answers are provided. Chapter 7 turns to unsupervised learning — discovering structure and patterns in data without any labels. The three main pillars are dimensionality reduction, clustering, and density estimation.

This post is based on Chapter 7 of Practical Statistics for Data Scientists.

7.1 Principal Component Analysis (PCA)#

PCA compresses high-dimensional data into a lower-dimensional representation while minimizing information loss.

Term	Description
Principal Component (PC)	A new axis defined in the direction of greatest variance — a linear combination of the original variables
Loading	The coefficient that describes how much each original variable contributes to a principal component — used to interpret the PC
Eigenvalue ( $\lambda_i$ )	The amount of variance explained by each principal component
Eigenvector	The eigenvector of the covariance matrix — defines the direction of each principal component
Scree Plot	A graph of eigenvalues in descending order — the number of components is chosen at the elbow point

Expressing a Principal Component#

$\text{PC}_1 = a_{11}x_1 + a_{12}x_2 + \cdots + a_{1p}x_p$

The coefficients $a_{1j}$ are elements of the eigenvector, representing how much variable $x_j$ contributes to the first principal component.

Cumulative Explained Variance#

$\text{Cumulative Explained Variance} = \frac{\sum_{i=1}^{k} \lambda_i}{\sum_{j=1}^{p} \lambda_j}$

Standardization is mandatory before running PCA. If variables are on different scales, those with larger magnitudes will dominate the principal components.

$z_j = \frac{x_j - \bar{x}_j}{s_j}$

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_breast_cancer

matplotlib.rcParams['font.family'] = 'AppleGothic'
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(42)

# Breast cancer dataset (30 features → reduce with PCA)
data = load_breast_cancer()
X, y = data.data, data.target

# Standardize (required before PCA)
scaler = StandardScaler()
X_sc = scaler.fit_transform(X)

# Run PCA
pca = PCA()
pca.fit(X_sc)

explained_var = pca.explained_variance_ratio_
cumulative_var = np.cumsum(explained_var)

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

# 1. Scree Plot
axes[0].bar(range(1, 11), explained_var[:10] * 100, color='steelblue', alpha=0.8, label='Individual explained variance')
axes[0].plot(range(1, 11), cumulative_var[:10] * 100, 'ro-', markersize=5, label='Cumulative explained variance')
axes[0].axhline(80, color='gray', linestyle='--', alpha=0.7, label='80% threshold')
axes[0].set_xlabel('Principal Component')
axes[0].set_ylabel('Explained Variance (%)')
axes[0].set_title('Scree Plot')
axes[0].legend()
axes[0].grid(alpha=0.3)

n_components_80 = np.argmax(cumulative_var >= 0.80) + 1
print(f"Components needed to explain 80% variance: {n_components_80} (out of {X.shape[1]} total)")

# 2. 2D visualization (PC1 vs PC2)
pca_2d = PCA(n_components=2)
X_pca = pca_2d.fit_transform(X_sc)

colors = ['red' if label == 0 else 'blue' for label in y]
axes[1].scatter(X_pca[:, 0], X_pca[:, 1], c=colors, alpha=0.6, s=20)
axes[1].set_xlabel(f'PC1 ({explained_var[0]*100:.1f}% explained)')
axes[1].set_ylabel(f'PC2 ({explained_var[1]*100:.1f}% explained)')
axes[1].set_title('PCA 2D Visualization (red=malignant, blue=benign)')
axes[1].grid(alpha=0.3)

# 3. Loading heatmap (top 10 features × PC1–3)
loadings = pd.DataFrame(
    pca.components_[:3].T,
    columns=['PC1', 'PC2', 'PC3'],
    index=data.feature_names
)
top_features = loadings['PC1'].abs().nlargest(10).index

im = axes[2].imshow(loadings.loc[top_features].values, cmap='coolwarm',
                    aspect='auto', vmin=-1, vmax=1)
axes[2].set_xticks([0, 1, 2])
axes[2].set_xticklabels(['PC1', 'PC2', 'PC3'])
axes[2].set_yticks(range(len(top_features)))
axes[2].set_yticklabels(top_features, fontsize=8)
axes[2].set_title('Loading Heatmap\n(Top 10 features by PC1 contribution)')
plt.colorbar(im, ax=axes[2])

plt.tight_layout()
plt.show()

# Top 5 features contributing to PC1
print("\nTop 5 features by PC1 loading (largest contribution to variance direction):")
print(loadings['PC1'].abs().nlargest(5).round(4))

Practical applications: Face recognition (Eigenfaces), gene expression analysis, visualization of high-dimensional text embeddings, and resolving multicollinearity. Choose the number of components based on 80–95% cumulative explained variance or the elbow point of the scree plot.

7.2 K-Means Clustering#

The most widely used clustering algorithm — partitions data into $K$ clusters.

Term	Description
Cluster	A group of similar data points
Centroid	The mean position of all points in a cluster
Inertia	Sum of squared distances from each point to its cluster center — lower is better

K-Means Algorithm#

Step 1 — Randomly initialize $K$ centroids

Step 2 — Assign each point to the nearest centroid: $\text{cluster}(x_i) = \arg\min_k \|x_i - \mu_k\|^2$

Step 3 — Recompute each centroid: $\mu_k = \frac{1}{|C_k|} \sum_{x_i \in C_k} x_i$

Step 4 — Repeat steps 2–3 until centroids no longer change

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

matplotlib.rcParams['font.family'] = 'AppleGothic'
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(42)

# Sample data
X, y_true = make_blobs(n_samples=400, centers=4, cluster_std=0.8, random_state=42)
X = StandardScaler().fit_transform(X)

# Find optimal K: Elbow method + Silhouette score
k_range = range(2, 11)
inertias, silhouettes = [], []

for k in k_range:
    km = KMeans(n_clusters=k, random_state=42, n_init=10)
    labels = km.fit_predict(X)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X, labels))

fig, axes = plt.subplots(1, 3, figsize=(16, 4))

# 1. Elbow method (inertia)
axes[0].plot(k_range, inertias, 'bo-', markersize=6, linewidth=2)
axes[0].set_xlabel('Number of Clusters K')
axes[0].set_ylabel('Inertia')
axes[0].set_title('Elbow Method — Inertia Decrease')
axes[0].grid(alpha=0.3)

# 2. Silhouette score
best_k = k_range[np.argmax(silhouettes)]
axes[1].plot(k_range, silhouettes, 'ro-', markersize=6, linewidth=2)
axes[1].axvline(best_k, color='gray', linestyle='--', label=f'Optimal K={best_k}')
axes[1].set_xlabel('Number of Clusters K')
axes[1].set_ylabel('Silhouette Score')
axes[1].set_title('Selecting Optimal K via Silhouette Score')
axes[1].legend()
axes[1].grid(alpha=0.3)

# 3. Clustering result with optimal K
km_best = KMeans(n_clusters=best_k, random_state=42, n_init=10)
labels = km_best.fit_predict(X)
centers = km_best.cluster_centers_

scatter = axes[2].scatter(X[:, 0], X[:, 1], c=labels, cmap='tab10', alpha=0.6, s=20)
axes[2].scatter(centers[:, 0], centers[:, 1], c='black', marker='X',
                s=200, linewidths=1.5, label='Centroids')
axes[2].set_title(f'K-Means Result (K={best_k})')
axes[2].legend()
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Optimal K: {best_k}")
print(f"Silhouette score: {max(silhouettes):.4f}  (closer to 1 = better cluster separation)")
print(f"Inertia: {km_best.inertia_:.2f}")

# K-Means limitation: sensitivity to initialization
print("\n=== K-Means initialization sensitivity (n_init=1, 5 runs) ===")
for _ in range(5):
    km_tmp = KMeans(n_clusters=4, n_init=1, random_state=None)
    km_tmp.fit(X)
    print(f"  Inertia: {km_tmp.inertia_:.2f}")
print("→ KMeans(n_init=10) recommended: picks the best result from multiple initializations")

Silhouette Score#

$s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}$

$a(i)$ : mean distance to points in the same cluster; $b(i)$ : mean distance to the nearest other cluster. $-1 \leq s(i) \leq 1$ — closer to 1 means better clustering.

Practical applications: Customer segmentation (RFM analysis), image color compression, document clustering, and geographic data clustering. K-Means works well for spherical clusters, but for non-spherical or variable-density clusters, DBSCAN or hierarchical clustering is often a better fit.

7.3 Hierarchical Clustering#

Builds a hierarchical structure based on pairwise distances without requiring $K$ to be specified upfront.

Linkage Method	Formula	Characteristics
Single linkage	$\min\{d(x,y) \mid x \in A, y \in B\}$	Tends to form elongated chain-like clusters
Complete linkage	$\max\{d(x,y) \mid x \in A, y \in B\}$	Produces compact, spherical clusters
Average linkage	$\frac{1}{\lvert A \rvert \lvert B \rvert}\sum_{x \in A}\sum_{y \in B} d(x,y)$	Balanced results
Ward linkage	Minimizes the increase in total within-cluster variance on merge	Most commonly used

The agglomerative (bottom-up) approach is standard: start with each data point as its own cluster, then iteratively merge the closest pair.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

matplotlib.rcParams['font.family'] = 'AppleGothic'
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(42)

X, _ = make_blobs(n_samples=150, centers=3, cluster_std=0.7, random_state=42)
X_sc = StandardScaler().fit_transform(X)

# Compare linkage methods
methods = ['single', 'complete', 'average', 'ward']
method_names = ['Single', 'Complete', 'Average', 'Ward']

fig, axes = plt.subplots(2, 4, figsize=(20, 8))

for i, (method, name) in enumerate(zip(methods, method_names)):
    Z = linkage(X_sc, method=method)
    labels = fcluster(Z, t=3, criterion='maxclust')

    # Dendrogram
    dendrogram(Z, ax=axes[0, i], truncate_mode='lastp', p=12,
               show_contracted=True, no_labels=True, color_threshold=None)
    axes[0, i].set_title(f'{name}\nDendrogram')
    axes[0, i].set_xlabel('Data Points')
    axes[0, i].set_ylabel('Distance')

    # Cluster result
    axes[1, i].scatter(X_sc[:, 0], X_sc[:, 1], c=labels, cmap='tab10', alpha=0.7, s=30)
    sil = silhouette_score(X_sc, labels)
    axes[1, i].set_title(f'Cluster Result\nSilhouette={sil:.3f}')
    axes[1, i].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# Ward linkage — visualizing how to choose K from the dendrogram
Z_ward = linkage(X_sc, method='ward')
plt.figure(figsize=(10, 5))
dend = dendrogram(Z_ward, truncate_mode='lastp', p=15, show_contracted=True)

# Find the optimal cut height (largest jump in merge distances)
last_merges = Z_ward[-10:, 2]  # distances of the last 10 merges
acceleration = np.diff(last_merges, 2)
k_suggested = acceleration.argmax() + 2
print(f"Suggested K by dendrogram acceleration: {k_suggested}")
plt.axhline(last_merges[-k_suggested], color='red', linestyle='--',
            label=f'K={k_suggested} cut line')
plt.xlabel('Clusters (merge history)')
plt.ylabel('Merge Distance (Ward linkage)')
plt.title('Hierarchical Clustering Dendrogram — Cut at red line to determine K')
plt.legend()
plt.show()

Practical applications: Gene expression heatmap visualization, language phylogenetics, and market segmentation. Unlike K-Means, results are visualized as a dendrogram, making the cluster structure easy to interpret. The downside is $O(n^2)$ time complexity, which makes it slow on large datasets.

7.4 Model-Based Clustering (GMM)#

Models clusters as probability distributions. Unlike K-Means, each point can partially belong to multiple clusters.

Term	Description
Mixture model	The overall distribution is a weighted sum of component distributions
EM algorithm	Alternates between E-step (compute responsibilities) and M-step (update parameters) until convergence
Responsibility ( $\gamma_{ik}$ )	The probability that data point $x_i$ belongs to cluster $k$
BIC / AIC	Criteria for selecting the optimal number of clusters (lower is better)

Gaussian Mixture Model (GMM)#

$f(x) = \sum_{k=1}^{K} \pi_k \cdot \mathcal{N}(x \mid \mu_k, \Sigma_k)$

$\pi_k$ : mixing proportion of cluster $k$ ( $\sum \pi_k = 1$ )

Multivariate Normal PDF#

$f(x) = \frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}} \exp\!\left(-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)\right)$

EM Algorithm#

E-step: Compute responsibilities $\gamma_{ik} = \frac{\pi_k f_k(x_i \mid \theta_k)}{\sum_{j=1}^{K} \pi_j f_j(x_i \mid \theta_j)}$

M-step: Update $\pi_k$ , $\mu_k$ , $\Sigma_k$ based on responsibilities → repeat until convergence

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_blobs

matplotlib.rcParams['font.family'] = 'AppleGothic'
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(42)

# Non-spherical cluster data (GMM has an advantage over K-Means here)
X1 = np.random.multivariate_normal([0, 0], [[1, 0.8], [0.8, 1]], 200)
X2 = np.random.multivariate_normal([4, 0], [[1, -0.6], [-0.6, 1]], 200)
X3 = np.random.multivariate_normal([2, 3], [[0.5, 0], [0, 2]], 200)
X = np.vstack([X1, X2, X3])

# Select optimal K using BIC/AIC
k_range = range(1, 9)
bic_scores, aic_scores = [], []

for k in k_range:
    gmm = GaussianMixture(n_components=k, random_state=42, n_init=5)
    gmm.fit(X)
    bic_scores.append(gmm.bic(X))
    aic_scores.append(gmm.aic(X))

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

# BIC / AIC
axes[0].plot(k_range, bic_scores, 'b-o', markersize=6, label='BIC')
axes[0].plot(k_range, aic_scores, 'r-o', markersize=6, label='AIC')
axes[0].axvline(np.argmin(bic_scores) + 1, color='blue', linestyle='--',
                label=f'Optimal K(BIC)={np.argmin(bic_scores)+1}')
axes[0].set_xlabel('K (number of clusters)')
axes[0].set_ylabel('Information Criterion Score (lower is better)')
axes[0].set_title('Selecting Optimal K via BIC / AIC')
axes[0].legend()
axes[0].grid(alpha=0.3)

# GMM vs K-Means comparison
from sklearn.cluster import KMeans

best_k = np.argmin(bic_scores) + 1
gmm_best = GaussianMixture(n_components=best_k, random_state=42, n_init=5)
km_best = KMeans(n_clusters=best_k, random_state=42, n_init=10)

gmm_labels = gmm_best.fit_predict(X)
km_labels  = km_best.fit_predict(X)

for ax, labels, title in [
    (axes[1], km_labels,  f'K-Means (K={best_k})'),
    (axes[2], gmm_labels, f'GMM (K={best_k})'),
]:
    ax.scatter(X[:, 0], X[:, 1], c=labels, cmap='tab10', alpha=0.5, s=15)
    ax.set_title(title)
    ax.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# GMM soft membership probabilities (contrast with K-Means)
probs = gmm_best.predict_proba(X[:5])
print("=== GMM Soft Membership Probabilities (first 5 points) ===")
print(f"{'Point':<6}", end='')
for k in range(best_k):
    print(f"  Cluster{k+1}", end='')
print()
for i, p in enumerate(probs):
    print(f"  {i+1}   ", end='')
    for prob in p:
        print(f"  {prob:.3f}   ", end='')
    print()
print("\n→ K-Means: each point is hard-assigned to exactly one cluster")
print("→ GMM: membership is expressed as a continuous probability over all clusters (soft assignment)")

Practical applications: User behavior segmentation (with soft boundaries), anomaly detection (low-density regions = outliers), and speech recognition (as the observation model in HMMs). GMM provides soft membership in contrast to K-Means' hard assignment, and handles elliptical clusters well.

7.5 Scaling and Categorical Variables#

Many unsupervised learning algorithms are distance-based, making scaling especially important.

Method	Formula	Characteristics
Standardization (Z-score)	$(x - \mu) / \sigma$	Mean 0, std 1 — sensitive to outliers
Min-Max Normalization	$(x - \min) / (\max - \min)$	Range $[0, 1]$ — sensitive to outliers
Robust Scaling	$(x - Q_2) / (Q_3 - Q_1)$	IQR-based — robust to outliers
Gower Distance	Handles mixed numeric + categorical	Suitable for mixed-type data

Min-Max Normalization#

$x'_j = \frac{x_j - \min(x_j)}{\max(x_j) - \min(x_j)}$

Gower Distance#

A distance metric that handles numeric, categorical, and binary variables simultaneously:

$D_{ij} = \frac{1}{p} \sum_{k=1}^{p} d_{ij}^{(k)}, \quad d_{ij}^{(k)} = 1 - s_{ij}^{(k)}$

For categorical variables: $d=0$ if equal, $d=1$ if different. For numeric variables: normalized absolute difference.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

matplotlib.rcParams['font.family'] = 'AppleGothic'
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(42)

# Comparing scaling methods
np.random.seed(42)
n = 200
X_raw = np.column_stack([
    np.random.normal(1000, 500, n),   # Annual income (10k KRW) — large scale
    np.random.normal(35, 10, n),      # Age — medium scale
    np.random.normal(0.5, 0.1, n),    # Purchase rate — small scale
])
# Add outliers
X_raw[[0, 10, 20]] = [[5000, 80, 0.95], [100, 18, 0.1], [4500, 75, 0.9]]

scalers = {
    'No Scaling (raw)':         None,
    'Standardization (Z-score)': StandardScaler(),
    'Min-Max Normalization':     MinMaxScaler(),
    'Robust Scaling':            RobustScaler(),
}

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

for ax, (name, scaler) in zip(axes, scalers.items()):
    X_plot = scaler.fit_transform(X_raw) if scaler else X_raw
    km = KMeans(n_clusters=3, random_state=42, n_init=10)
    labels = km.fit_predict(X_plot)
    sil = silhouette_score(X_plot, labels)

    ax.scatter(X_plot[:, 0], X_plot[:, 1], c=labels, cmap='tab10', alpha=0.6, s=20)
    ax.set_xlabel('Feature1 (Income)' if not scaler else 'Feature1 (scaled)')
    ax.set_ylabel('Feature2 (Age)' if not scaler else 'Feature2 (scaled)')
    ax.set_title(f'{name}\nSilhouette={sil:.3f}')
    ax.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# Gower distance: example with mixed-type data
print("=== Gower Distance (Mixed-Type Data) ===")
df_mixed = pd.DataFrame({
    'Age':    [25, 30, 35, 50],
    'Income': [3000, 4500, 3200, 8000],
    'Gender': ['M', 'F', 'M', 'F'],          # Categorical
    'Tier':   ['Standard', 'Premium', 'Standard', 'VIP'],  # Categorical
})
print(df_mixed)
print("\n→ Numeric features (Age, Income): normalized difference; Categorical (Gender, Tier): 0 if equal, 1 if different")
print("→ Gower distance integrates all variable types on the same scale")

try:
    from gower import gower_matrix
    dist_matrix = gower_matrix(df_mixed)
    print("\nGower Distance Matrix:")
    print(pd.DataFrame(dist_matrix, columns=range(4), index=range(4)).round(3))
except ImportError:
    print("\n※ Install with: pip install gower")

# Scaling method selection guide
print("\n=== Scaling Method Selection Guide ===")
guide = {
    'StandardScaler': 'No outliers, roughly normal distribution — PCA, logistic regression, SVM',
    'MinMaxScaler':   'Need a specific range ([0,1]) — image pixels, neural network inputs',
    'RobustScaler':   'Many outliers — financial data, medical data',
    'Gower Distance': 'Mixed numeric + categorical data — customer profile clustering',
}
for method, desc in guide.items():
    print(f"  {method:<20}: {desc}")

Practical applications: Customer data is often mixed — age (numeric), region (categorical), purchase flag (binary). In these cases, use Gower distance to handle everything uniformly, or encode categoricals and apply Robust scaling. The choice of scaling method has a significant impact on clustering results.

Chapter 7 key takeaways:

PCA: Projects onto eigenvectors of the covariance matrix → dimensionality reduction that preserves variance — standardization required
K-Means: Iteratively updates centroids to form clusters — sensitive to initialization, best for spherical clusters
Choosing optimal K: Use the elbow method (inertia) together with the silhouette score
Hierarchical clustering: Visualizes structure via dendrogram — no need to specify K upfront, but inefficient on large datasets
GMM: Probabilistic soft membership — handles elliptical clusters well, select K via BIC/AIC
Scaling: Essential for distance-based algorithms — choose the method based on the presence of outliers
Gower distance: A unified distance metric for mixed numeric + categorical data

// Related Posts

[Data Science]·03/04/2026·20min

Practical Statistics for Data Scientists — Quiz

Test your understanding across all seven chapters — from descriptive statistics to unsupervised learning — with concept, calculation, code, and scenario problems.

#Statistics#Quiz#Data Science#Interview Prep

[Data Science]·03/04/2026·13min

Statistics Fundamentals Chapter 6: Statistical Machine Learning

A hands-on walkthrough of KNN, decision trees, random forests, AdaBoost, and gradient boosting — the core concepts behind tree-based ensemble models, complete with code examples.

#Statistics#Machine Learning#Random Forest#Boosting#Python