k-means法とGMM：クラスタリング手法の理論とPython実装

k-means法の収束保証からGMM・EMアルゴリズム、scikit-learnによる実装、クラスタ数選択（エルボー法・シルエット分析・BIC）までを解説します。

2 min read

はじめに

クラスタリングは、ラベルなしデータをグループに分割する教師なし学習の代表的なタスクです。本記事では最も広く使われるk-means法と、確率モデルに基づく**混合ガウスモデル（GMM）**を解説します。

k-means法

アルゴリズム

  1. \(k\) 個のセントロイド \(\boldsymbol{\mu}_1, \ldots, \boldsymbol{\mu}_k\) をランダムに初期化
  2. 割り当てステップ: 各データ \(\mathbf{x}_i\) を最も近いセントロイドのクラスタに割り当て
\[c_i = \arg\min_{j} \|\mathbf{x}_i - \boldsymbol{\mu}_j\|^2 \tag{1}\]
  1. 更新ステップ: 各クラスタのセントロイドを再計算
\[\boldsymbol{\mu}_j = \frac{1}{|C_j|} \sum_{i \in C_j} \mathbf{x}_i \tag{2}\]
  1. 収束するまで 2-3 を繰り返す

目的関数

k-means は以下の目的関数（イナーシャ）を最小化します。

\[J = \sum_{j=1}^{k} \sum_{i \in C_j} \|\mathbf{x}_i - \boldsymbol{\mu}_j\|^2 \tag{3}\]

各ステップで \(J\) は単調に減少するため、有限回で収束することが保証されます。ただしグローバル最適性は保証されず、初期値依存があります。

k-means++

初期セントロイドをデータ点間の距離に比例した確率で逐次選択することで、初期値問題を緩和します。scikit-learn のデフォルトです。

混合ガウスモデル（GMM）

モデル定義

\(k\) 個のガウス分布の混合で確率密度を表現します。

\[p(\mathbf{x}) = \sum_{j=1}^{k} \pi_j \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j) \tag{4}\]

ここで \(\pi_j\) は混合比率（\(\sum_j \pi_j = 1\)）、\(\boldsymbol{\Sigma}_j\) は共分散行列です。

EM アルゴリズム

GMMのパラメータは期待値最大化（EM）アルゴリズムで推定します。

E ステップ: 各データの各クラスタへの帰属確率（責任度）を計算

\[\gamma_{ij} = \frac{\pi_j \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}{\sum_{l=1}^{k} \pi_l \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_l, \boldsymbol{\Sigma}_l)} \tag{5}\]

M ステップ: パラメータを更新

\[\boldsymbol{\mu}_j = \frac{\sum_i \gamma_{ij} \mathbf{x}_i}{\sum_i \gamma_{ij}}, \quad \boldsymbol{\Sigma}_j = \frac{\sum_i \gamma_{ij} (\mathbf{x}_i - \boldsymbol{\mu}_j)(\mathbf{x}_i - \boldsymbol{\mu}_j)^T}{\sum_i \gamma_{ij}} \tag{6}\]\[\pi_j = \frac{\sum_i \gamma_{ij}}{n} \tag{7}\]

k-means との比較

特徴k-meansGMM
割り当てハード（0 or 1）ソフト（確率）
クラスタ形状球状（等方的）楕円体（任意の共分散）
目的関数イナーシャ対数尤度
計算コスト低い高い

Python実装

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

# --- データ生成 ---
np.random.seed(42)
X, y_true = make_blobs(n_samples=500, centers=[[-2, -2], [2, 2], [0, 4]],
                        cluster_std=[0.8, 1.2, 0.6], random_state=42)
scaler = StandardScaler()
X = scaler.fit_transform(X)

# --- k-means ---
kmeans = KMeans(n_clusters=3, init='k-means++', n_init=10, random_state=42)
km_labels = kmeans.fit_predict(X)

# --- GMM ---
gmm = GaussianMixture(n_components=3, covariance_type='full', random_state=42)
gmm.fit(X)
gmm_labels = gmm.predict(X)
gmm_probs = gmm.predict_proba(X)

# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# k-means
axes[0].scatter(X[:, 0], X[:, 1], c=km_labels, cmap='viridis', s=15, alpha=0.6)
axes[0].scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],
                c='red', marker='X', s=200, edgecolors='k')
axes[0].set_title('k-means')

# GMM（ハード割り当て）
axes[1].scatter(X[:, 0], X[:, 1], c=gmm_labels, cmap='viridis', s=15, alpha=0.6)
axes[1].set_title('GMM (hard)')

# GMM（ソフト割り当て：不確実性を透明度で表現）
uncertainty = 1 - gmm_probs.max(axis=1)
axes[2].scatter(X[:, 0], X[:, 1], c=gmm_labels, cmap='viridis',
                s=15, alpha=1 - uncertainty * 0.8)
axes[2].set_title('GMM (soft: uncertainty)')

for ax in axes:
    ax.set_xlabel('Feature 1')
    ax.set_ylabel('Feature 2')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

クラスタ数の選択

エルボー法

イナーシャの減少が緩やかになる「肘」の位置を目視で判断します。

inertias = []
K_range = range(1, 10)
for k in K_range:
    km = KMeans(n_clusters=k, n_init=10, random_state=42)
    km.fit(X)
    inertias.append(km.inertia_)

plt.plot(K_range, inertias, 'bo-')
plt.xlabel('k')
plt.ylabel('Inertia')
plt.title('Elbow Method')
plt.grid(True, alpha=0.3)
plt.show()

シルエット分析

\[s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))} \tag{8}\]

\(a(i)\) は同クラスタ内の平均距離、\(b(i)\) は最も近い他クラスタとの平均距離です。\(s(i) \in [-1, 1]\) で、1に近いほど良いクラスタリングです。

BIC（GMM の場合）

bics = []
K_range = range(1, 10)
for k in K_range:
    gm = GaussianMixture(n_components=k, random_state=42)
    gm.fit(X)
    bics.append(gm.bic(X))

plt.plot(K_range, bics, 'ro-')
plt.xlabel('k')
plt.ylabel('BIC')
plt.title('BIC for GMM')
plt.grid(True, alpha=0.3)
plt.show()

BICが最小となる \(k\) を選びます。

関連記事

参考文献

  • Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. Chapters 9.
  • MacQueen, J. (1967). “Some methods for classification and analysis of multivariate observations”. Proceedings of the 5th Berkeley Symposium.
  • scikit-learn Clustering documentation

関連記事