k-means法とGMM（混合ガウスモデル）の理論・比較・Python実装：sklearn.cluster.KMeans/sklearn.mixture.GaussianMixture/EMアルゴリズム

はじめに

クラスタリングは、ラベルなしデータをグループに分割する教師なし学習の代表的なタスクです。本記事では最も広く使われるk-means法と、確率モデルに基づく**混合ガウスモデル（GMM）**を解説します。単に更新式を提示するだけでなく、目的関数の収束証明・EMアルゴリズムの導出・両手法の理論的な関係・特異共分散や局所解といった実務上の注意点までを、実際に実行したPythonコードの数値とともに一貫して扱います。

k-means法

アルゴリズム

1. \(k\) 個のセントロイド \(\boldsymbol{\mu}_1, \ldots, \boldsymbol{\mu}_k\) をランダムに初期化
2. 割り当てステップ: 各データ \(\mathbf{x}_i\) を最も近いセントロイドのクラスタに割り当て

3. 更新ステップ: 各クラスタのセントロイドを再計算

4. 収束するまで 2-3 を繰り返す

目的関数と収束証明

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

データ \(i\) の割り当て先クラスタを \(c_i \in \{1,\ldots,k\}\) と書けば、\(J = \sum_{i=1}^n \|\mathbf{x}_i - \boldsymbol{\mu}_{c_i}\|^2\) とも表せます。k-means のアルゴリズムは、この \(J\) を割り当て \(\{c_i\}\) とセントロイド \(\{\boldsymbol{\mu}_j\}\) という2種類の変数について交互に最小化する**座標降下法（block coordinate descent）**とみなせます。

割り当てステップが \(J\) を最小化する証明：\(\{\boldsymbol{\mu}_j\}\) を固定すると、\(J = \sum_i \|\mathbf{x}_i - \boldsymbol{\mu}_{c_i}\|^2\) は各データ点について独立な項の和になります。したがって全体を最小化するには各項 \(\|\mathbf{x}_i - \boldsymbol{\mu}_{c_i}\|^2\) を個別に最小化すればよく、これは最も近いセントロイドを選ぶこと（式(1)）そのものです。つまり割り当てステップは \(\{\boldsymbol{\mu}_j\}\) を固定した下での \(J\) の厳密な最小化です。

更新ステップが \(J\) を最小化する証明：割り当て \(\{c_i\}\) を固定すると、\(J\) はクラスタごとに分離できます。

各 \(J_j\) は \(\boldsymbol{\mu}_j\) について凸な二次関数なので、勾配をゼロと置けば大域最小が求まります。

これは式(2)の更新式そのものであり（ヘッセ行列 \(2|C_j|I\) は正定値なので最小性も保証されます）、更新ステップも厳密な最小化です。

単調非増加性と有限回収束：割り当て・更新の両ステップとも \(J\) を減少させることはあっても増加させることはないため、反復列 \(J^{(0)} \ge J^{(1)} \ge \cdots \ge 0\) は単調非増加かつ下に有界です。さらに重要なのは、これが有限回で厳密に停止する点です。\(n\) 個の点を \(k\) 個のクラスタに割り当てるパターンは高々 \(k^n\) 通りで有限個です。割り当てが変化しない限りセントロイドも変化せずアルゴリズムは停止します。逆に割り当てが変化した場合（同点でのタイブレークを除く）、少なくとも1点がより近いセントロイドへ移動したことになるため \(J\) は厳密に減少します。よって同一の割り当てパターンが2度現れることはあり得ず、高々 \(k^n\) 回のパターンを一度ずつしか通過できないため、アルゴリズムは有限回の反復で必ず収束します（Bishop, 2006, 9.1節; Lloyd, 1982）。ただし収束先は \(J\) の局所最小に過ぎず、大域最小とは限りません（次節の「k-means++」も参照）。

実行検証：この証明を確認するため、後述のPython実装で使う3クラスタのデータに対し、割り当てステップ後・更新ステップ後それぞれの \(J\) を記録する手動実装を用意しました。ランダムに選んだ3点を初期セントロイドとして開始すると、\(J\) は \(672.854 \to 308.005 \to 301.931 \to 301.290 \to 300.958 \to 300.548 \to \cdots \to 300.037\) と単調に減少し続け、8回の反復（16回のJ評価）で割り当てが変化しなくなり収束しました。反復間の差分の最大値は \(0\) （増加は一度も観測されず）であり、単調非増加性が数値的に確認できました。

k-means++と局所解の問題

前節で見た通り、k-means が収束する先は \(J\) の局所最小に過ぎず、初期セントロイドの選び方によって結果は大きく変わり得ます。特にクラスタの大きさや密度が大きく異なるデータでは、ランダム初期化だと複数の初期セントロイドが同一の大きなクラスタ内に集中し、小さなクラスタが見逃されるという典型的な失敗パターンが起こります。

k-means++（Arthur & Vassilvitskii, 2007）はこれを緩和する初期化法です。1個目のセントロイドをデータ点からランダムに選び、以降は各データ点 \(\mathbf{x}\) について既に選ばれたセントロイド集合への最短距離の2乗 \(D(\mathbf{x})^2\) を計算し、

に比例する確率で次のセントロイドを選びます。既存セントロイドから遠い点ほど選ばれやすくなるため、クラスタ全体をカバーしやすい初期配置が得られます。scikit-learn の KMeans はデフォルトで init='k-means++' を使用します。

実行検証：局所解の問題を実際に確認するため、大きく広がった1個の主クラスタ（300点、標準偏差3.0）と、離れた場所に密集した4個の小クラスタ（各20点、標準偏差0.3）からなる、意図的に偏った5クラスタのデータを作成し、init='random' と init='k-means++' をそれぞれ200回ずつ（n_init=1）実行して最終イナーシャを比較しました。参照値として n_init=200 のk-means++で得た大域最適の近似値は \(J_{\text{best}} \approx 1898.748\) でした。

ランダム初期化では200回中7回（3.5%）が大域最適の1.5倍を超える明らかな局所解（小クラスタを1個以上見逃す配置）に陥ったのに対し、k-means++では1度も発生しませんでした。標準偏差も \(361.994\) から \(91.664\) へと約1/4に縮小しており、k-means++が最悪ケースを効果的に排除しつつ結果の安定性を大きく高めていることが定量的に確認できます。一方、大域最適に近い解（大域最適の0.1%以内）に到達した割合は random が14.0%（200回中28回）、k-means++ が12.0%（200回中24回）とほぼ同水準であり、k-means++は「最良の解に到達する確率」よりも「最悪の解を避ける効果」に大きく寄与するという点は誤解しやすいので注意が必要です。

混合ガウスモデル（GMM）

モデル定義

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

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

対数尤度最大化の難しさとEMアルゴリズムの導出

GMMのパラメータ \(\theta = \{\pi_j, \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j\}\) は、観測データの対数尤度

を最大化することで推定します。しかし線形回帰の尤度と異なり、この式は対数の中に和がある（\(\log\sum\) ）ため、\(\theta\) について微分してゼロと置いても閉じた形の解が得られません（\(\boldsymbol{\mu}_j\) の微分でも他のクラスタの項が分母に絡み合ったまま残ります）。そこで、各データ点がどのクラスタに属するかを表す潜在変数 \(z_i \in \{1,\ldots,k\}\) を導入し、任意の分布 \(q_i(j) \ge 0,\ \sum_j q_i(j)=1\) を使って各項を書き換えます。

対数関数 \(\log\) は凹関数なので、Jensenの不等式 \(\log \mathbb{E}[Y] \ge \mathbb{E}[\log Y]\) を \(q_i(j)\) を重みとする期待値に適用すると、

が成り立ちます。両辺を \(i\) について総和すると、対数尤度の下界（ELBO: Evidence Lower Bound）が得られます。

Jensenの不等式で等号が成立するのは、期待値の対象 \(Y=\pi_j \mathcal{N}(\mathbf{x}_i|\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)/q_i(j)\) が \(j\) によらず一定になる場合です。これを解くと等号成立条件は \(q_i(j) \propto \pi_j \mathcal{N}(\mathbf{x}_i|\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)\) 、すなわち正規化して

となります。これはまさに次式(5)の責任度であり、ベイズの定理による潜在変数の事後分布 \(p(z_i=j|\mathbf{x}_i,\theta)\) に他なりません。つまりEステップとは、現在のパラメータ \(\theta^{\text{old}}\) の下でELBOを最大化する \(q\) （＝潜在変数の事後分布）を求める操作です。この \(q=\gamma\) を選ぶと、ELBOは現在のパラメータでの対数尤度に一致します（\(\mathcal{L}(\gamma,\theta^{\text{old}}) = \ell(\theta^{\text{old}})\) 、下界が「タイト」になる）。

Eステップで \(q=\gamma\) を固定した後、Mステップはこの \(q\) の下でELBO \(\mathcal{L}(\gamma,\theta)\) を \(\theta\) について最大化する操作です。例えば \(\boldsymbol{\mu}_j\) については、\(\mathcal{L}(\gamma,\theta)\) の中で \(\boldsymbol{\mu}_j\) に依存する項は \(\sum_i \gamma_{ij}\log\mathcal{N}(\mathbf{x}_i|\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j) = -\frac12\sum_i \gamma_{ij}(\mathbf{x}_i-\boldsymbol{\mu}_j)^T\boldsymbol{\Sigma}_j^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_j) + \text{const}\) だけです。これを \(\boldsymbol{\mu}_j\) で微分してゼロと置くと

となり、後述の式(6)の \(\boldsymbol{\mu}_j\) の更新式と一致します。\(\boldsymbol{\Sigma}_j\) や \(\pi_j\) （制約 \(\sum_j\pi_j=1\) 付きのラグランジュ関数を微分）についても同様の手順で式(6)(7)が導かれます。\(\log\mathcal{N}\) がガウス分布の二次形式であることから、責任度 \(\gamma_{ij}\) による重み付き最尤推定（重み付き平均・重み付き共分散）という直感的な形になっています。

EM アルゴリズム

GMMのパラメータは期待値最大化（EM）アルゴリズムで推定します。上記の通り、EステップとMステップはいずれもELBOの最大化として統一的に理解できます。

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

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

EMアルゴリズムの単調性証明

EMアルゴリズムの各反復で対数尤度 \(\ell(\theta)\) が単調に非減少であることを証明します。反復 \(t\) でのパラメータを \(\theta^{(t)}\) 、対応するEステップの責任度を \(\gamma^{(t)}\) （ELBOの \(q\) ）とすると、

が成り立ちます。各不等式・等式の根拠は次の通りです。

1. 左の不等号：ELBOは任意の \(q\) に対して対数尤度の下界なので（Jensenの不等式）、\(q=\gamma^{(t)}\) とパラメータ \(\theta^{(t+1)}\) の組でも \(\ell(\theta^{(t+1)}) \ge \mathcal{L}(\gamma^{(t)},\theta^{(t+1)})\) が成り立ちます。
2. 中央の不等号：Mステップは \(\gamma^{(t)}\) を固定してELBOを \(\theta\) について最大化する操作であり、\(\theta^{(t+1)}\) はその最大化の結果なので、\(\theta^{(t)}\) を含む任意の \(\theta\) における値以上になります。
3. 右の等号：Eステップで選んだ \(\gamma^{(t)}=p(z|\mathbf{x},\theta^{(t)})\) はJensenの不等式の等号成立条件を満たすため、ELBOは \(\theta^{(t)}\) においてタイトです。

これらを繋げると \(\ell(\theta^{(t+1)}) \ge \ell(\theta^{(t)})\) が任意の反復で成立し、EMアルゴリズムは対数尤度を単調に非減少にします（EMは一般にMM法＝Majorize-Maximizationアルゴリズムの一種として説明されます）。ただしk-meansと同様、収束先は対数尤度の局所最大に過ぎず、大域最適性は保証されません。

実行検証：この単調性を確認するため、後述のPython実装と同じ3クラスタデータに対しEMアルゴリズムを手動実装しました（共分散に \(10^{-6}\) の正則化項を追加、初期値はランダムに選んだ3点を平均とし単位行列を共分散として開始）。各反復での対数尤度を記録した結果、初期値の対数尤度は \(-1442.438\) でしたが、反復ごとに単調に増加し（\(-1139.789 \to -1030.377 \to -1012.214 \to \cdots\) ）、30回の反復後には \(-911.193\) に収束しました。反復間の差分の最小値は \(3.60 \times 10^{-6}\) （常に非負）であり、単調非減少性が数値的に確認できました。総増加量は \(531.245\) でした。

k-meansとの関係：等方分散→0の退化極限

GMMとk-meansは無関係な2つの手法ではなく、GMMの特殊な極限としてk-meansが得られるという関係にあります。全クラスタの共分散を等方的かつ共通の値 \(\boldsymbol{\Sigma}_j = \epsilon I\) に固定し（\(\epsilon\) は推定せず定数とする）、混合比も均等 \(\pi_j = 1/k\) とします。このとき式(5)の責任度は

という、温度パラメータ \(\epsilon\) を持つソフトマックス（Gibbs分布）の形になります。\(\epsilon \to 0\) の極限では、指数の中の \(-\|\mathbf{x}_i-\boldsymbol{\mu}_j\|^2/(2\epsilon)\) は、最も近いセントロイド（\(\|\mathbf{x}_i-\boldsymbol{\mu}_j\|^2\) が最小のクラスタ）以外では \(-\infty\) に発散するため、そのクラスタの責任度だけが1に近づき他は0に収束します。

これはまさに式(1)のk-meansの割り当てステップ（ハード割り当て）です。さらにMステップの \(\boldsymbol{\mu}_j\) 更新式もこの極限で \(\gamma_{ij} \to \mathbb{1}[i \in C_j]\) となるため、式(2)の単純平均に一致します。つまりk-meansは、共分散が等方的かつ分散をゼロに近づけた極限でのGMM-EMアルゴリズムに他ならないという関係が成り立ちます（Bishop, 2006, 9.3.2節）。この意味で、k-meansの「ハード割り当て」はGMMの「ソフト割り当て」が退化した特殊ケースであると理解できます。

実行検証：この極限を確認するため、上記のk-means実装で得られたセントロイド（3クラスタ）を固定し、共通の等方分散 \(\epsilon\) を \(1.0\) から \(0.0001\) まで段階的に小さくしながら責任度 \(\gamma_{ij}\) を計算しました。

\(\epsilon\) を小さくするほど最大責任度が1に近づき（\(\epsilon=0.0001\) では実質的に \(1.0000\) ）、\(\arg\max_j \gamma_{ij}\) によるハード割り当てはすべての \(\epsilon\) でk-meansの割り当てと100%一致しました。セントロイドの位置自体はk-meansの解に固定しているため割り当ては変化しませんが、責任度の「柔らかさ」が \(\epsilon \to 0\) でハード割り当てへ収束していく様子が定量的に確認できます。

k-means との比較

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()

上記コードを実際に実行した結果が上図です。左からk-means、GMM（ハード割り当て）、GMM（ソフト割り当て：クラスタ境界付近の不確実な点ほど透明度が高い）を並べています。境界付近の点で色の濃淡が変わっている様子から、GMMが単純な多数決ではなく確率的な帰属度を保持していることが視覚的に確認できます。

真のクラスタラベルとの一致度をsklearn.metrics.adjusted_rand_score（ARI、1に近いほど真のラベルと一致）で比較すると、k-meansが0.861、GMMが0.902という結果になりました。本例ではクラスタごとに標準偏差（0.8/1.2/0.6）が異なる楕円形に近い分布を使っているため、球状クラスタを仮定するk-meansより、共分散行列を個別に推定できるGMMの方が真の構造を正確に捉えられることが定量的に確認できます。

エッジケースと注意点

特異共分散行列と尤度の発散

GMMの対数尤度には、教科書ではあまり強調されない重大な病理があります。あるクラスタ \(j\) に割り当てられる（責任度で重み付けされた）データ点が実質的に1点しかない場合、あるいは複数点が線形従属（同一直線上、あるいは完全に重複）である場合、共分散行列

が特異行列（行列式ゼロ、逆行列が存在しない）に近づきます。ガウス分布の密度 \(\mathcal{N}(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}) \propto |\boldsymbol{\Sigma}|^{-1/2}\exp(\cdots)\) は \(|\boldsymbol{\Sigma}| \to 0\) で発散するため、GMMの対数尤度の大域的な上限は実は \(+\infty\) であり、1つの成分がデータ点1つに収縮する退化解がその発散を引き起こします。これはGMMのMLEが原理的に抱える既知の病理であり（Day, 1969; McLachlan & Peel, 2000）、EM反復で偶然そのような割り当てに近づくと、共分散が特異に近づき数値的に不安定になったり、以降の反復で破綻します。

実務的な対処法は、共分散行列の対角成分に小さな正の値を加える正則化（sklearn.mixture.GaussianMixtureのreg_covarパラメータ、デフォルト \(10^{-6}\) ）です。これにより \(\boldsymbol{\Sigma}_j + \text{reg\_covar}\cdot I\) を使うことで、行列式が厳密に0になることを防ぎます。

実行検証：この病理を実際に再現するため、200点の主クラスタ（正規分布）に加えて、完全に同一の座標 \((10,10)\) を持つ点を5個複製したデータセットを作り、n_components=3のGMMをreg_covar=0（正則化なし）でフィットしました。結果、scikit-learnは実際に例外を送出しました。

ValueError: Fitting the mixture model failed because some components have
ill-defined empirical covariance (for instance caused by singleton or
collapsed samples). Try to decrease the number of components, increase
reg_covar, or scale the input data.

これはまさに理論通り、複製点に1成分が収縮しようとして共分散が特異になったために発生したエラーです。一方、デフォルトのreg_covar=1e-6で同じデータをフィットすると正常に収束し（全データでの対数尤度は \(-524.729\) ）、複製点に対応する成分の共分散行列式は \(1.0\times 10^{-12}\) （他の2成分は \(0.496\) 、\(0.597\) ）と、他の成分より8桁以上小さい値になりました。この成分の混合比は \(0.0244\) で、これは複製点数5個を全データ数205個で割った値（\(5/205=0.0244\) ）と正確に一致しており、GMMがこの5点だけに1成分をほぼ割り当てて退化した解に近づいていることが定量的に確認できます。正則化がなければこの成分の共分散はさらに収縮を続け、尤度は理論上無限大に発散していたはずです。

局所解と初期値依存性

k-meansの局所解問題とk-means++による緩和については前述の通りです（「k-means++と局所解の問題」節を参照）。ここで見落とされがちな実務上の注意点があります。scikit-learn の KMeans の n_init パラメータのデフォルトは 'auto' ですが、これは init='k-means++'（デフォルト）の場合は1回だけ実行に解決されます（init='random' の場合のみ10回に解決されます。sklearn.cluster._kmeans._BaseKMeans._check_params_vs_input の実装で確認できます）。つまり「k-means++を使っているから自動的に複数回試行される」という思い込みは誤りで、明示的に n_init を指定しない限り scikit-learn は k-means++ 初期化を1回しか行いません。本記事のPython実装例（後述）では n_init=10 を明示的に指定していますが、これは前節の実験で見た局所解リスク（ランダム初期化では200回中7回、3.5%が大域最適の1.5倍を超える悪い局所解に陥ったのに対し、k-means++では0回だった）を踏まえた意図的な選択です。

GMMの GaussianMixture も同様に対数尤度の局所最大にしか収束を保証せず、n_init パラメータのデフォルトは明示的に 1 です。KMeansの 'auto' のような賢い解決ロジックはないため、複数のクラスタ構造が複雑な場合は明示的に n_init を増やす必要があります。

クラスタ数の選択

k-means・GMMはいずれもクラスタ数 \(k\) を事前に指定する必要があります。ここでは代表的な3つの選択基準を、理論的根拠と限界とともに解説します。

エルボー法

k-meansのイナーシャ \(J\) は \(k\) を増やすほど単調に減少します（\(k=n\) で \(J=0\) ）。イナーシャの減少が緩やかになる「肘」の位置を目視で判断するのがエルボー法です。理論的な最適性の保証はなく、あくまで経験則的なヒューリスティックです。減少カーブが滑らかで肘が不明瞭な場合（後述の実験のように、クラスタが近接している場合など）は判断が難しく、シルエット分析やBICと併用するのが実務上望ましいとされます。

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()

シルエット分析

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

シルエット分析の落とし穴：シルエットスコアは「クラスタ内の凝集度」と「クラスタ間の分離度」のトレードオフを最大化する \(k\) を選びますが、これは必ずしも真のクラスタ数と一致しません。特にクラスタが互いに近接している場合、隣接する2クラスタを1つに統合した方が見かけ上のクラスタ間分離が改善し、シルエットスコアの上ではより小さい \(k\) が有利に評価されることがあります。

BIC（GMM の場合）

BIC（ベイズ情報量規準）は、モデルのエビデンス（周辺尤度）\(p(D|k) = \int p(D|\theta,k)p(\theta|k)\,d\theta\) を近似することで導かれます。この積分は高次元では計算困難なため、最尤推定値 \(\hat\theta\) の周りでラプラス近似（対数エビデンスの2次テイラー展開による正規近似）を適用します。

ここで \(p\) はモデルの自由パラメータ数、\(\mathbf{H}\) は負の対数尤度のヘッセ行列（観測フィッシャー情報行列）です。独立同分布データでは \(\mathbf{H} \approx n\bar{\mathbf{H}}\) （\(\bar{\mathbf{H}}\) はデータ数に依存しない行列）とみなせるため、\(\log|\mathbf{H}| \approx p\log n + \log|\bar{\mathbf{H}}|\) となります。\(n\to\infty\) で支配的なのは \(p\log n\) の項だけなので、\(n\) に依存しない定数項を無視すると

が得られます（Schwarz, 1978）。これに \(-2\) を掛けて最小化問題の形にしたものが、慣例的なBICの定義です。

GMM（共分散full）の自由パラメータ数 \(p\) は、混合比が \(k-1\) 個（制約 \(\sum_j\pi_j=1\) のため独立なのは \(k-1\) 個）、平均が \(kd\) 個、共分散行列が対称行列 \(k\) 個分で各 \(d(d+1)/2\) 個より

となります。AIC（赤池情報量規準）は同じ対数尤度項に対し罰則項が \(2p\) である点が異なります（\(\text{AIC} = -2\ell(\hat\theta) + 2p\) ）。\(p\log n\) の罰則を持つBICは、\(n\) が大きいほどAICより強くモデルの複雑さを罰するため、一般にBICの方が小さい \(k\) を選ぶ傾向があります。

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\) を選びます。

実行検証：エルボー・シルエット・BIC・AICの比較

本記事のPython実装で使った3クラスタデータ（\(k=3\) が真の値）に対し、\(k=2\) から \(9\) までイナーシャ・シルエットスコア・BIC・AICを計算し、また \(d=2,\,k=3\) の自由パラメータ数を式から計算した値（\(p=(3-1)+3\cdot2+3\cdot3=17\) ）が GaussianMixture(n_components=3)._n_parameters() の返す値と一致するかも確認しました。

自由パラメータ数は式の計算通り \(p=17\) で、_n_parameters() の返り値と一致しました。シルエットスコアが最大になったのは \(k=2\) （スコア \(0.6331\) ）で、これは真の \(k=3\) と一致しません。本データは3クラスタのうち2つがやや近接して生成されている（標準偏差1.2の大きなクラスタと0.8のクラスタが隣接）ため、2クラスタに統合した方が見かけ上の分離度が高く評価されたと考えられます。一方、BIC・AICはいずれも \(k=3\) で最小値（\(\text{BIC}=1928.32\) 、\(\text{AIC}=1856.67\) ）を取り、真のクラスタ数を正しく当てました。この結果は、シルエット分析を唯一の基準にせず、密度モデルであるBIC/AICや目視によるエルボー法と併用すべきことを実証的に示しています。

関連記事

• アンサンブル学習の手法と比較 - 教師あり学習の代表的手法。クラスタリング結果を特徴量として活用できます。
• Self-Attentionの仕組みとPython実装 - Attentionの重み付け平均とGMMの混合モデルの類似性を理解できます。
• サポートベクターマシン（SVM） - 教師あり分類との比較。SVMは分類境界、クラスタリングはデータ構造の発見です。
• MCMC入門：メトロポリス・ヘイスティングス法とギブスサンプリング - GMMのパラメータ推定をMCMCで行う手法もあります。
• ベイズ線形回帰の基礎 - GMMの背景にあるベイズ的な考え方を解説しています。
• SGDからAdamまで：勾配降下法の理論と実装 - EMアルゴリズムと勾配法の関係を理解する参考になります。
• 機械学習による時系列予測・分類・異常検知ハブ - k-means/GMM・アンサンブル・LSTM・カルマンの使い分けを 3 軸 × 9 シナリオで整理したハブ記事。
• 時系列データの異常検知：統計的手法からカルマンフィルタまで - クラスタリングを使った教師なし異常検知の発展形。クラスタ密度の低い点や混合成分の責任度を異常スコアとして使えます。

参考文献

• 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