Phân cụm — Clustering

k-Means chia dữ liệu thành \(K\) cụm sao cho tổng bình phương khoảng cách trong cụm (within-cluster sum of squares, WCSS / SSE) là nhỏ nhất. Mỗi cụm \(C_k\) được đại diện bởi một trọng tâm \(\mu_k\):

$$J=\sum_{k=1}^{K}\sum_{x\in C_k}\lVert x-\mu_k\rVert^2$$

Đây là bài toán NP-khó nếu giải tối ưu toàn cục, nên trên thực tế ta dùng thuật toán Lloyd — lặp hai bước cho tới khi hội tụ:

• Bước GÁN (assignment): gán mỗi điểm \(x_i\) vào cụm có trọng tâm gần nhất, tức \(C_k=\{x_i:\lVert x_i-\mu_k\rVert^2\le\lVert x_i-\mu_j\rVert^2\ \forall j\}\).
• Bước CẬP NHẬT (update): tính lại mỗi trọng tâm bằng trung bình các điểm trong cụm: \(\mu_k=\dfrac{1}{|C_k|}\sum_{x\in C_k}x\).

Mỗi vòng lặp đều làm \(J\) giảm (hoặc giữ nguyên), nên thuật toán bảo đảm hội tụ — nhưng chỉ tới cực tiểu cục bộ, phụ thuộc vào điểm khởi tạo.

Khởi tạo k-means++: thay vì chọn ngẫu nhiên, k-means++ chọn các tâm ban đầu sao cho chúng phân tán nhau (điểm càng xa các tâm đã chọn thì xác suất được chọn càng cao). Cách này giảm mạnh nguy cơ kẹt ở cực tiểu cục bộ tồi.

Chọn số cụm \(K\):

• Elbow method: vẽ \(J\) theo \(K\), chọn điểm "khuỷu tay" nơi \(J\) bắt đầu giảm chậm lại.
• Silhouette: đo mức độ một điểm "hợp" với cụm của nó so với cụm gần nhất; điểm số càng gần \(1\) càng tốt.

k-Medians là biến thể của k-Means nhưng dùng trung vị (median) làm tâm cụm và đo khoảng cách theo chuẩn Manhattan (L1) thay vì Euclid bình phương:

$$J=\sum_{k=1}^{K}\sum_{x\in C_k}\lVert x-\text{med}_k\rVert_1$$

Trong đó \(\text{med}_k\) là trung vị (tính theo từng chiều) của các điểm trong cụm \(C_k\). Vì trung vị ít nhạy với giá trị cực đoan hơn trung bình, k-Medians bền với outlier hơn k-Means: một vài điểm nằm rất xa sẽ không kéo lệch tâm cụm như khi dùng trung bình.

Khác với k-Means/k-Medians (gán cứng — mỗi điểm thuộc đúng một cụm), GMM thực hiện phân cụm mềm (soft): mỗi điểm thuộc về nhiều cụm cùng lúc, với một xác suất tương ứng.

Mô hình: giả định dữ liệu được sinh ra từ một hỗn hợp của \(K\) phân phối Gaussian:

$$p(x)=\sum_{k=1}^{K}\pi_k\,\mathcal{N}(x\mid\mu_k,\Sigma_k),\qquad \sum_k\pi_k=1,\ \pi_k\ge 0$$

trong đó \(\pi_k\) là trọng số (tỉ lệ) của thành phần thứ \(k\), còn \(\mu_k,\Sigma_k\) là kỳ vọng và ma trận hiệp phương sai. Ta ước lượng các tham số này bằng thuật toán EM, lặp hai bước:

• E-step (Expectation): với tham số hiện tại, tính "trách nhiệm" (responsibility) — xác suất hậu nghiệm mỗi điểm thuộc về cụm \(k\): $$\gamma_{ik}=P(k\mid x_i)=\frac{\pi_k\,\mathcal{N}(x_i\mid\mu_k,\Sigma_k)}{\sum_{j=1}^{K}\pi_j\,\mathcal{N}(x_i\mid\mu_j,\Sigma_j)}$$
• M-step (Maximization): dùng \(\gamma\) làm trọng số để cập nhật lại tham số. Đặt \(N_k=\sum_i\gamma_{ik}\): $$\mu_k=\frac{1}{N_k}\sum_i\gamma_{ik}x_i,\quad \Sigma_k=\frac{1}{N_k}\sum_i\gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^\top,\quad \pi_k=\frac{N_k}{N}$$

EM lặp tới khi hội tụ và bảo đảm cực đại hoá log-likelihood của dữ liệu (mỗi vòng lặp likelihood không giảm). Ưu điểm lớn là cụm có thể có hình elip bất kỳ (nhờ \(\Sigma_k\) đầy đủ), không bị bó buộc dạng cầu.

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.metrics import silhouette_score

# Chuẩn hoá trước khi phân cụm (rất quan trọng)
X = StandardScaler().fit_transform(X_raw)

# --- k-Means (gán cứng) ---
km = KMeans(n_clusters=4, init="k-means++", n_init=10, random_state=42)
labels = km.fit_predict(X)
print("inertia (SSE):", km.inertia_)
print("silhouette:", silhouette_score(X, labels))

# --- Elbow: quét K để chọn số cụm ---
for k in range(2, 9):
    m = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
    print(k, m.inertia_)

# --- Gaussian Mixture (gán mềm, EM) ---
gmm = GaussianMixture(n_components=4, covariance_type="full", random_state=42)
gmm.fit(X)
proba = gmm.predict_proba(X)   # xác suất thuộc từng cụm (responsibility)
hard  = gmm.predict(X)         # nhãn cứng nếu cần
print("log-likelihood TB:", gmm.score(X))