ある衛星画像から「都市」「森林」「水面」の3種類の領域を分けたい——ラベル付きの教師データは無い、でも各ピクセルの分光特徴ベクトルだけはある。あるいは、顧客の購買履歴を眺めて、自然な「グループ」をいくつかに分けたい。こうした「正解ラベルなしで似ているもの同士をまとめる」課題が クラスタリング であり、その代表的な手法のひとつが ガウス混合モデル(Gaussian Mixture Model, GMM) です。
GMMをひとことで言えば、「データは複数のガウス分布の重ね合わせから生成された」という確率モデルを立て、その混合比率・平均・共分散をデータから推定する手法です。よく比較される K-means が「各クラスタを中心点ひとつで代表する」のに対し、GMMは「各クラスタを平均ベクトルと共分散行列の組で代表する」点が決定的に違います。共分散行列を持つということは、クラスタの 大きさ・形状・向き までモデルに含めるということで、楕円体に伸びたクラスタや、密度の違うクラスタも自然に扱えます。
イメージとしてはこちらの画像を見ると混合ガウスモデルのイメージが湧くと思います。
混合ガウスモデルのアルゴリズムを用いることで、左側のようなデータの分布を与えたときに、右側のようにデータのクラスタリングと、そのクラスタを生成していると仮定されるガウス分布を得ることができます。右側のグラフは、ガウス分布の確率の等高線を描いており、これによって、そのデータが生成される確率がどの程度かを定量的に評価することができます。確率モデルだからこそ、「このデータ点がクラスタAに属する確率は70%、Bに属する確率は30%」のような ソフトな所属度 を出せるのもGMMの大きな魅力です。
GMMの応用先は驚くほど広く、音声認識の発音モデル(HMM-GMM)、画像のセグメンテーション、異常検知(低密度の点を異常と判定)、密度推定によるサンプル生成、混合エキスパートモデルなど枚挙にいとまがありません。最近の Variational Autoencoder や拡散モデルといった生成モデルも、混合分布の発想を土台にしている部分があります。GMMを理解しておくと、これらの先進的なモデルの内部で何が起きているのかが自然に見えてくるはずです。
本記事の内容
- GMMの直感的理解 — なぜK-meansでは足りないのか
- GMMの厳密な数学的定義(混合係数 $\pi_k$、平均 $\bm{\mu}_k$、共分散 $\bm{\Sigma}_k$)
- EMアルゴリズムのE-step(責任度の計算)とM-step(閉形式更新)の導出
- 対数尤度がEM反復で単調増加することの証明
- 初期化問題 — K-means++ vs ランダム初期化
- AIC/BICによるクラスタ数の決定
- 縮退共分散と特異点問題、その回避策
- 変分GMM・ベイズGMMへの発展
- Pythonでのスクラッチ実装と scikit-learn 実装
- Irisデータセットでの実例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
混合ガウスモデルをわかりやすく
混合ガウスモデルは混合モデルの1種
混合ガウスモデルは、重み付きの異なるガウス分布の線型結合からなる、混合モデル(mixture model)として表現できます。混合ガウスモデルを用いることで、多峰性(multimodal)の分布を表現することができます。
わかりやすく図解すると次のようになります。
左側は正規分布のグラフで、1つの山があるグラフとなっています。一方右側が、薄い線で描いている3つの正規分布を足し合わせたグラフになっています。赤い実線の分布に注目すると、山が2つあるような分布になっていることがわかります。
このように混合モデルを用いることで、よく使われるようなガウス分布やベータ、ガンマ分布のような確率分布では表現できないようなデータの分布を表現できるようになります。直感的に言えば、「単峰性の分布をいくつか用意して、好きな比率で足し合わせる」という極めて柔軟な構成法であり、ガウス分布を十分多く用意すれば、ほとんど任意の連続分布を任意の精度で近似できるという性質(普遍近似性)が知られています。
一般的に、上記のような$K$個のガウス分布が複数重ね合わさった確率分布を混合ガウス分布(Mixture Gaussian Distribution)といい、次のような式で定義されます。
混合ガウス分布の数学的定義
定義に入る前に、式が何を表しているかを一文で予告します。「データ点 $\bm{x}$ の確率密度を、$K$ 個のガウス分布それぞれの密度を $\pi_k$ という重みで線形結合した値で表す」——これがGMMの本体です。
$$ \begin{equation} p(\bm{x} \mid \bm{\pi}, \bm{\mu}, \bm{\Sigma}) = \sum_{k=1}^{K} \pi_k \,\mathcal{N}(\bm{x} \mid \bm{\mu}_k, \bm{\Sigma}_k) \end{equation} $$
ここで、$\pi_k$ は $k$ 番目のガウス分布の混合比率を表すパラメータであり、次の制約を満たします。
$$ \sum_{k=1}^{K}\pi_k = 1, \quad 0 \leq \pi_k \leq 1 $$
また、$\bm{\pi} = (\pi_1, \pi_2, \dots, \pi_K)^\top$ であり、$\bm{\mu}_k \in \mathbb{R}^d$、$\bm{\Sigma}_k \in \mathbb{R}^{d \times d}$ はそれぞれ $k$ 番目のガウス分布の平均ベクトルと共分散行列です。$\mathcal{N}(\bm{x} \mid \bm{\mu}_k, \bm{\Sigma}_k)$ は $d$ 次元の多変量正規分布の確率密度関数で、
$$ \mathcal{N}(\bm{x} \mid \bm{\mu}, \bm{\Sigma}) = \frac{1}{(2\pi)^{d/2} |\bm{\Sigma}|^{1/2}} \exp\!\left[-\frac{1}{2}(\bm{x}-\bm{\mu})^\top \bm{\Sigma}^{-1} (\bm{x}-\bm{\mu})\right] $$
と書けます。(1)式から、GMMの全パラメータをまとめて $\bm{\theta} = \{\pi_k, \bm{\mu}_k, \bm{\Sigma}_k\}_{k=1}^{K}$ と書くことにします。
潜在変数を導入したもうひとつの表現
GMMには、より見通しの良いもうひとつの定式化があります。「データ点ごとに、まずどのガウス成分から生成されたかを表す 潜在変数 $z \in \{1, 2, \dots, K\}$ をサンプリングし、その後その成分のガウス分布から $\bm{x}$ を生成する」という階層的な生成プロセスとして書く方法です。
$$ \begin{aligned} z &\sim \mathrm{Cat}(\bm{\pi}) \quad (\text{離散カテゴリ分布}) \\ \bm{x} \mid z=k &\sim \mathcal{N}(\bm{\mu}_k, \bm{\Sigma}_k) \end{aligned} $$
このとき $z$ を周辺化して消去すると、
$$ p(\bm{x}) = \sum_{k=1}^{K} p(z=k)\, p(\bm{x} \mid z=k) = \sum_{k=1}^{K} \pi_k\, \mathcal{N}(\bm{x} \mid \bm{\mu}_k, \bm{\Sigma}_k) $$
となり、(1)式と完全に一致します。潜在変数の存在こそがEMアルゴリズムの鍵 であり、EMはこの未観測の $z$ を上手に扱うための一般的な道具立てだと考えると見通しが立ちます。
通常パラメータの推定は、データから対数尤度関数を計算し、それを各パラメータで偏微分してゼロとおく——という最尤推定の手続きを取ります。ところが(1)式のように対数の中に和が入っている形 $\log \sum_k \pi_k \mathcal{N}(\cdots)$ は、偏微分しても閉形式で解けず、通常の最尤推定は使えません。ここで登場するのが EMアルゴリズム という反復的最適化手法です。
EMアルゴリズムの一般論については、こちらの記事をご覧ください。
[
【機械学習】EMアルゴリズムをゼロから理解する
EMアルゴリズム(expectation maximization algori…](https://disassemble-channel.com/em-algorithm/)
EMの一般論を踏まえた上で、次節ではGMM特有のE-step・M-stepの更新式を一行ずつ丁寧に導出していきます。閉形式で書けるところがGMMの美しさです。
GMMにおけるEMアルゴリズムの導出
ゴールを宣言する
ここでは「GMMの対数尤度 $\log p(X \mid \bm{\theta}) = \sum_{n=1}^{N} \log \sum_{k=1}^{K} \pi_k \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \bm{\Sigma}_k)$ を反復的に増加させる更新規則」を導出することがゴールです。具体的に得たいのは、次の2ステップを交互に回す手続きです。
- E-step: 現在のパラメータの下で、各データ点 $\bm{x}_n$ が各成分 $k$ に属する確率(責任度) $\gamma_{nk}$ を計算する
- M-step: 責任度を重みとして、$\pi_k$、$\bm{\mu}_k$、$\bm{\Sigma}_k$ を更新する
この2ステップを繰り返すと、対数尤度は単調に増加し、必ずどこかの停留点(極大値または鞍点)に収束します。
E-step: 責任度(responsibility)の導出
潜在変数 $z_n \in \{1, \dots, K\}$ を導入すると、データ点 $\bm{x}_n$ の同時確率は
$$ p(\bm{x}_n, z_n = k \mid \bm{\theta}) = \pi_k\, \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \bm{\Sigma}_k) $$
と書けます。ここで、現在のパラメータ $\bm{\theta}^{(t)}$ の下で「$\bm{x}_n$ を観測したとき $z_n = k$ である事後確率」をベイズの定理で計算すると、
$$ \begin{aligned} \gamma_{nk} &\equiv p(z_n = k \mid \bm{x}_n, \bm{\theta}^{(t)}) \\ &= \frac{p(\bm{x}_n, z_n=k \mid \bm{\theta}^{(t)})}{\sum_{j=1}^{K} p(\bm{x}_n, z_n=j \mid \bm{\theta}^{(t)})} \\ &= \frac{\pi_k^{(t)}\, \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k^{(t)}, \bm{\Sigma}_k^{(t)})}{\sum_{j=1}^{K} \pi_j^{(t)}\, \mathcal{N}(\bm{x}_n \mid \bm{\mu}_j^{(t)}, \bm{\Sigma}_j^{(t)})} \end{aligned} $$
この $\gamma_{nk}$ を 責任度(responsibility) と呼びます。「点 $\bm{x}_n$ をどれだけ成分 $k$ が説明しているか」を $[0, 1]$ の値で表しており、$\sum_k \gamma_{nk} = 1$ が常に成り立ちます。これがE-stepの計算結果です。
M-step: パラメータの閉形式更新
次に、責任度を固定したときに対数尤度の下界 $Q(\bm{\theta}, \bm{\theta}^{(t)})$ を最大化する $\bm{\theta}$ を求めます。EMの一般論より、$Q$ は
$$ Q(\bm{\theta}, \bm{\theta}^{(t)}) = \sum_{n=1}^{N} \sum_{k=1}^{K} \gamma_{nk} \left[ \log \pi_k + \log \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \bm{\Sigma}_k) \right] $$
と書けます。ここがEMの核心です。対数の中の和 だった元の対数尤度が、和の中の対数 に化けています。これにより各パラメータについて閉形式で最大化できるようになるのです。
まず $\bm{\mu}_k$ について微分してゼロとおきます。$\log \mathcal{N}$ の $\bm{\mu}_k$ に依存する項は $-\tfrac{1}{2}(\bm{x}_n – \bm{\mu}_k)^\top \bm{\Sigma}_k^{-1}(\bm{x}_n – \bm{\mu}_k)$ なので、
$$ \frac{\partial Q}{\partial \bm{\mu}_k} = \sum_{n=1}^{N} \gamma_{nk} \bm{\Sigma}_k^{-1}(\bm{x}_n – \bm{\mu}_k) = \bm{0} $$
$\bm{\Sigma}_k^{-1}$ を左から掛けて消すと、
$$ \sum_{n=1}^{N} \gamma_{nk}\, \bm{x}_n = \bm{\mu}_k \sum_{n=1}^{N} \gamma_{nk} $$
$N_k \equiv \sum_n \gamma_{nk}$ を「成分 $k$ が抱えるソフトなデータ数」と定義すると、
$$ \boxed{\,\bm{\mu}_k^{(t+1)} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma_{nk}\, \bm{x}_n\,} $$
つまり 責任度で重み付けしたデータの平均 が新しい平均ベクトルになります。$\gamma_{nk}$ が全部 $0$ か $1$ の極限を取ると、これはまさに K-means の平均更新と一致します。
次に $\bm{\Sigma}_k$。$\log \mathcal{N}$ の中の $\bm{\Sigma}_k$ 依存項は $-\tfrac{1}{2}\log |\bm{\Sigma}_k| – \tfrac{1}{2}(\bm{x}_n-\bm{\mu}_k)^\top \bm{\Sigma}_k^{-1}(\bm{x}_n-\bm{\mu}_k)$ です。行列微分の公式 $\partial \log|\bm{\Sigma}| / \partial \bm{\Sigma} = \bm{\Sigma}^{-\top}$、$\partial (\bm{x}^\top \bm{\Sigma}^{-1} \bm{x}) / \partial \bm{\Sigma} = -\bm{\Sigma}^{-\top} \bm{x}\bm{x}^\top \bm{\Sigma}^{-\top}$ を使って微分しゼロとおくと、
$$ \boxed{\,\bm{\Sigma}_k^{(t+1)} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma_{nk}\, (\bm{x}_n – \bm{\mu}_k^{(t+1)})(\bm{x}_n – \bm{\mu}_k^{(t+1)})^\top\,} $$
これも 責任度で重み付けした標本共分散 という直感的な形になっています。
最後に $\pi_k$。制約 $\sum_k \pi_k = 1$ があるのでラグランジュ未定乗数法を使います。ラグランジュ関数
$$ \mathcal{L} = Q + \lambda\!\left(\sum_k \pi_k – 1\right) $$
を $\pi_k$ で微分してゼロとおくと
$$ \sum_n \frac{\gamma_{nk}}{\pi_k} + \lambda = 0 \quad \Rightarrow\quad \pi_k = -\frac{N_k}{\lambda} $$
$\sum_k \pi_k = 1$ を使うと $\lambda = -N$ となり、
$$ \boxed{\,\pi_k^{(t+1)} = \frac{N_k}{N}\,} $$
つまり 責任度の総和の比 が新しい混合係数です。「成分 $k$ がトータルで何点ぶん抱えているか」を全体で割っただけの、非常に自然な式です。
対数尤度の単調増加
EMの最大の魅力は、反復するたびに対数尤度が決して下がらない という強力な性質です。これを示すには、対数尤度 $\log p(X \mid \bm{\theta})$ が
$$ \log p(X \mid \bm{\theta}) = Q(\bm{\theta}, \bm{\theta}^{(t)}) – \sum_{n,k} \gamma_{nk} \log p(z_n = k \mid \bm{x}_n, \bm{\theta}) $$
と分解できることを使います。第二項はクロスエントロピー的な量で、$\bm{\theta}=\bm{\theta}^{(t)}$ で最小(KLダイバージェンス的に下界に張り付く)になります。M-step で $Q$ を最大化したパラメータ $\bm{\theta}^{(t+1)}$ では第一項が必ず増え、第二項は最小値からしか動けないので、合計は単調増加します。詳しい議論はEMアルゴリズム個別記事をご覧ください。
上の図は、3つの異なる初期化からEMを回したときの対数尤度の推移です。3本とも単調に増加していることが見て取れます。初期値によって到達する値が違うのが分かりますが、これは局所最適解の存在を示しています。EMが保証するのはあくまで「対数尤度を下げない」ことであり、大域最適に到達することは保証していません。実用上はランダム初期化を複数回試して最良のものを採用するのが定石です。
ここまでで「責任度を計算」→「閉形式でパラメータ更新」というEMの完全な手続きが得られました。次のセクションでは、なぜGMMが K-means より優れているのか、その本質を決定境界の形で見ていきます。
K-meansとの比較 — なぜ共分散を持つことが効くのか
K-meansは「各クラスタを中心点だけで代表する」極めて簡素なアルゴリズムです。実は K-means は GMMの特殊ケース(全成分の共分散を $\bm{\Sigma}_k = \sigma^2 \bm{I}$ で固定し、$\sigma \to 0$ の極限を取ったもの)に相当することが知られています。この極限では責任度 $\gamma_{nk}$ が 0/1 のハード割り当てになり、距離が最も近い中心への割り当てに一致します。
逆に言えば、GMMはK-meansに 「クラスタの形状(共分散)」「ソフト割り当て」「クラスタごとに異なるサイズ」 という3つの追加表現力を与えたものです。この差が最も劇的に効くのは、楕円体に伸びたクラスタや、向きが斜めになったクラスタを分類するときです。
上の図では、向きが斜めに伸びた2つの楕円体クラスタを K-means と GMM で分類しています。左の K-means では決定境界が必ず直線になるため、楕円の交差部分でデータを誤分類してしまいます。一方右の GMM では、各成分の共分散行列を学習することで決定境界が曲線になり、楕円の形状に追随してきれいに分けられます。「形を見て分ける」ことができる、これがGMMの真骨頂です。
ソフト割り当ても重要な利点です。境界近くの点について「クラスタAに60%、Bに40%属する」という出力ができることで、不確かさをそのまま下流のタスクに伝えられます。たとえば異常検知では「どのクラスタにも所属確率が低い点」を異常と判定する、混合エキスパートでは責任度を専門家ネットワークのゲートとして使う、といった応用が自然に展開できます。
ここまでで「GMMは何が嬉しいのか」が明らかになりました。次に、実用上避けて通れない初期化問題に踏み込みます。
初期化問題 — EMはどこから始めるか
EMは対数尤度を下げないことしか保証しないので、到達する局所最適解は初期値に強く依存 します。下手な初期値からは下手な解に収束する、というのが現実です。代表的な初期化戦略には次の3つがあります。
ランダム初期化
各 $\bm{\mu}_k$ をデータ点からランダムに $K$ 個サンプリングし、$\bm{\Sigma}_k$ を全データの分散行列、$\pi_k = 1/K$ で揃える、というシンプルな方法。手軽だが、たまたまクラスタの近くに初期値が落ちないと貧弱な解に収束しがちです。複数回試して最良の対数尤度を採用する multiple restarts が現実的な対策です。
K-means++ 初期化
K-means の文脈で提案された手法ですが、GMMでも非常に有効です。手順は以下のとおりです。
- 最初の中心 $\bm{\mu}_1$ をデータ点から一様ランダムに選ぶ
- 残りの中心 $\bm{\mu}_k$ ($k = 2, \dots, K$) は、各データ点 $\bm{x}_n$ について「既に選ばれた中心への最短距離 $D(\bm{x}_n)^2$」に比例する確率でサンプリングする
- これを繰り返して $K$ 個の中心を得る
直感的には「既存の中心から離れた点を優先的に新たな中心にする」ことで、初期値が偏らず空間にまんべんなく散らばります。scikit-learn の GaussianMixture(init_params='kmeans')(デフォルト)は実質的にこの戦略を取っています。
階層クラスタリングによる初期化
データ数が比較的小さい場合、Ward法などの凝集型階層クラスタリングを先に走らせ、得られた $K$ クラスタの中心を初期値にする方法も使われます。決定論的になる代わりに計算コストが高めです。
実装上のコツは、複数の初期化を試し、最良の対数尤度に到達したものを採用する ことです。scikit-learnでは n_init=10 のように指定すれば自動で複数試行してくれます。EMが弱いのは初期化次第で簡単に壊れることなので、ここはケチらず計算コストを払う価値があります。
初期化が落ち着いたら次の難問は「そもそも $K$ をいくつにすればいいか」です。次節で扱います。
クラスタ数 $K$ の選び方 — AICとBIC
教師なし学習であるクラスタリングでは、適切なクラスタ数を 自動的に決める ことが永遠の課題です。GMMでは確率モデルとしての性質を活かして、対数尤度にパラメータ数の罰則項を加えた 情報量基準 で比較するのが標準的です。
対数尤度だけでは決まらない
まず注意点を確認します。対数尤度 $\log p(X \mid \bm{\theta})$ は $K$ を増やすほど 常に増加 します。なぜなら、$K=10$ のモデルは $K=9$ のモデルを部分集合として含み、よりデータにフィットさせる自由度を持っているからです。この事実が「対数尤度を最大化する $K$ を選ぶ」という素朴な方針が機能しない理由です。すべてのデータ点に1つずつガウスを当てる($K=N$)のが対数尤度最大ですが、これは明らかに過学習です。
AIC(赤池情報量基準)
$$ \mathrm{AIC} = -2 \log p(X \mid \hat{\bm{\theta}}) + 2 p $$
ここで $p$ はモデルのパラメータ数です。$d$ 次元GMMの場合、$K-1$ 個の独立な $\pi_k$、$K \cdot d$ 個の平均要素、$K \cdot d(d+1)/2$ 個の共分散要素で、合計 $p = K-1 + Kd + Kd(d+1)/2$ となります。AICが小さいほど良いモデルです。
BIC(ベイズ情報量基準)
$$ \mathrm{BIC} = -2 \log p(X \mid \hat{\bm{\theta}}) + p \log N $$
サンプル数 $N$ が罰則項に入っているのがAICとの違いです。$\log N$ は通常 $2$ より大きい(たとえば $N=1000$ なら $\log N \approx 6.9$)ので、BICはAICよりも複雑なモデルに対して厳しい罰則を課します。結果として、BICは「真のモデルが候補集合の中にあるとき、サンプル数を増やすと真の $K$ を正しく選ぶ」という一致性を持つことが知られています(AICは一致性を持たない)。
上のプロットは、真のクラスタ数が3である合成データに対して $K=1, 2, \dots, 8$ でGMMを当て、AIC/BICを計算したものです。左の対数尤度は予想通り $K$ に対して単調増加しています。右のAIC/BICは $K=3$ 付近で最小値を取り、真のクラスタ数を正しく拾えていることが分かります。BICのほうが谷の位置がよりはっきりしているのは、罰則項が強いためです。
その他の選び方
実用ではAIC/BIC以外に、ICL(Integrated Completed Likelihood)、Cross-Validation、Variational Bayes(後述)による自動決定、シルエットスコアやエルボー法など多様な手法があります。データの性質や応用先(分類精度を重視するのか、密度推定としての汎化を重視するのか)に応じて使い分けるのが現実的です。
クラスタ数が決まったら、次の落とし穴は「共分散行列が変な形になる」病理的なケースです。
縮退共分散と特異点問題
EMをGMMで素朴に回すと、ときどき対数尤度がいきなり無限大に発散することがあります。これは決してバグではなく、GMMの 本質的な病理 です。原因を見ていきましょう。
何が起きるか
ある成分 $k$ の平均 $\bm{\mu}_k$ が、ちょうど1つのデータ点 $\bm{x}_n$ にぴったり一致した状況を考えます。このとき $(\bm{x}_n – \bm{\mu}_k) = \bm{0}$ なので、その点での密度は
$$ \mathcal{N}(\bm{x}_n \mid \bm{\mu}_k, \bm{\Sigma}_k) = \frac{1}{(2\pi)^{d/2} |\bm{\Sigma}_k|^{1/2}} $$
となります。ここで $\bm{\Sigma}_k$ をどんどん縮ませる($|\bm{\Sigma}_k| \to 0$)と、密度は 無限大に発散 し、対数尤度も $+\infty$ に飛びます。
左図に示すように、1次元GMMで成分の標準偏差 $\sigma$ を縮ませていくと、その平均位置での密度 $\mathcal{N}(0 | 0, \sigma^2) = 1/\sqrt{2\pi \sigma^2}$ は対数的に発散します。右図は2成分のうち1つが1点に「貼り付こうとしている」病的状況の概念図です。健全な成分(青)とは対照的に、赤の成分は分散を縮めるほど見せかけの尤度が上がる罠にはまっています。
なぜ起きるのか
これはGMMの最尤推定が 数学的に病的(ill-posed) であることに起因します。対数尤度に上限がなく、最大値が存在しないのです。EMは局所最適への収束を保証しますが、その「局所」が特異点近傍だと、対数尤度がいくらでも上がる方向に逃げていってしまいます。
対処法
実用上は以下の3つの戦略が組み合わされます。
- 正則化項を加える — 共分散行列に下限を設ける、たとえば $\bm{\Sigma}_k \leftarrow \bm{\Sigma}_k + \lambda \bm{I}$ という形で、最小固有値を強制的に押し上げる。scikit-learn の
reg_covarパラメータ(デフォルト $10^{-6}$)はこれを行っています。 - ベイズ的事前分布 — 共分散に逆Wishart事前を置くと、MAP推定としては特異点が自動的に避けられます。Variational GMM(後述)は本質的にこれを行います。
- 特異点を検出して再初期化 — ある成分の有効サンプル数 $N_k$ が極端に小さくなったり、$|\bm{\Sigma}_k|$ が一定以下になったら、その成分を別のランダムな位置に再配置する。
scikit-learn では covariance_type で全成分共通の共分散('tied')、対角行列('diag')、球状('spherical')に制限することもでき、これらは特異点を起こしにくくなります(自由度が下がるトレードオフ)。
特異点を避ける根本的な解決として、ベイズ的な定式化(Variational GMM)があります。次節で扱います。
共分散の幾何と「クラスタの形」
GMMがK-meansを凌駕する最大の理由は、共分散行列が クラスタの形状(楕円体) を表現できる点にあります。この幾何的な意味を押さえておくと、GMMの出力結果を解釈しやすくなります。
固有値分解で見る楕円
任意の正定値対称行列 $\bm{\Sigma}$ は、固有値分解 $\bm{\Sigma} = \bm{V}\bm{\Lambda}\bm{V}^\top$ を持ちます。ここで $\bm{V}$ は固有ベクトルを列に並べた直交行列、$\bm{\Lambda} = \mathrm{diag}(\lambda_1, \dots, \lambda_d)$ は固有値の対角行列です。
多変量正規分布 $\mathcal{N}(\bm{\mu}, \bm{\Sigma})$ の等密度線は、$(\bm{x}-\bm{\mu})^\top \bm{\Sigma}^{-1}(\bm{x}-\bm{\mu}) = $ 一定 という楕円体になります。この楕円体の 主軸の方向 が固有ベクトル $\bm{v}_i$、軸の長さ が $\sqrt{\lambda_i}$ に比例します。つまり
- 固有値がすべて等しい → 球(等方的)
- 固有値が異なるが固有ベクトルが座標軸に揃っている → 軸に沿った楕円(対角共分散)
- 固有ベクトルが斜めを向いている → 傾いた楕円(一般の共分散)
上の4パネルは、左上から順に等方的・対角・正相関・負相関の共分散行列が描く楕円体を可視化したものです。各図で 1σ・2σ・3σ の等密度楕円を描き、固有ベクトル(矢印)と固有値(楕円の主軸の長さ)の対応を見ています。共分散の非対角要素が正だと右上がりに、負だと右下がりに楕円が傾くことが視覚的に確認できます。
scikit-learn の covariance_type
scikit-learn の GaussianMixture には4つの共分散タイプがあります。
'full': 各成分が独立な一般の共分散行列を持つ(最も柔軟、パラメータ多)'tied': 全成分が同じ共分散行列を共有'diag': 各成分が独立な対角行列(軸並行の楕円)'spherical': 各成分が独立だが球状($\bm{\Sigma}_k = \sigma_k^2 \bm{I}$)
サンプル数が少ない、あるいは特徴量が多い場合は 'diag' や 'spherical' が過学習を防ぐのに有効です。'full' は表現力が最も高いですが、パラメータ数が $K d(d+1)/2$ と次元の2乗で増えるため、高次元では注意が必要です。
幾何的直感が固まったところで、次にベイズ的拡張へと話を進めます。
Variational GMM とベイズGMMへの発展
最尤推定のGMMは「特異点問題」「クラスタ数の自動決定が難しい」という根本的な弱点を抱えています。これらを 同時に 解決するのが ベイズGMM、特に 変分ベイズ(Variational Bayes, VB) によるアプローチです。
ベイズ的定式化
ベイズGMMでは、パラメータ $\bm{\pi}, \bm{\mu}, \bm{\Sigma}$ を確定値ではなく確率変数として扱い、それぞれに共役事前分布を置きます。
- $\bm{\pi} \sim \mathrm{Dirichlet}(\alpha_0)$
- $\bm{\mu}_k \mid \bm{\Sigma}_k \sim \mathcal{N}(\bm{m}_0, \beta_0^{-1} \bm{\Sigma}_k)$
- $\bm{\Sigma}_k^{-1} \sim \mathrm{Wishart}(\bm{W}_0, \nu_0)$
このように事前分布を置くと、共分散の下限が事前分布によって自然に保証され、特異点問題が消滅します。さらに、$\alpha_0$ を小さく設定すると 「使われない成分の混合係数が自動的にゼロに押し下げられる」 という素晴らしい性質が現れます。これにより、$K$ を大きめに設定しておけば、実際に必要なクラスタ数が事後分布から自動的に決まるのです。
変分ベイズの更新式
事後分布 $p(\bm{\theta}, \bm{Z} \mid X)$ は閉形式で書けませんが、$q(\bm{\theta}, \bm{Z}) = q(\bm{\theta}) q(\bm{Z})$ という分解された近似(平均場近似)を仮定すると、各因子の更新式が解析的に得られます。本質的にはEMの拡張で、
- VE-step: $q(\bm{Z}) \propto \exp\!\left[\mathbb{E}_{q(\bm{\theta})}[\log p(X, \bm{Z}, \bm{\theta})]\right]$ を更新
- VM-step: $q(\bm{\theta}) \propto \exp\!\left[\mathbb{E}_{q(\bm{Z})}[\log p(X, \bm{Z}, \bm{\theta})]\right]$ を更新
を交互に回します。最大化される量はELBO(Evidence Lower BOund)で、これがVBの収束保証を与えます。詳細はBishopの「Pattern Recognition and Machine Learning」第10章が定番です。
scikit-learn では BayesianGaussianMixture でこの変分GMMが利用できます。実用上、最尤推定のGMMより遥かに頑健で、特異点問題なし・クラスタ数自動決定・正則化込み という三拍子が揃っており、まず試すべきデフォルトと言えます。
Dirichlet Process Mixture(DPMM)
クラスタ数を理論上「無限」とし、データから自然に有効な数だけが立ち上がるノンパラメトリックな拡張が Dirichlet Process Mixture Model(DPMM) です。Stick-breaking表現と変分推論を組み合わせると、計算可能なアルゴリズムになります。scikit-learn の BayesianGaussianMixture で weight_concentration_prior_type='dirichlet_process' を選ぶとこの動作になります。
理論的な発展はここまでにして、いよいよ実装に進みましょう。
Pythonでの実装
スクラッチ実装でEMを完全理解する
まず、ライブラリに頼らずにEMアルゴリズムをゼロから書いてみます。これにより、E-stepとM-stepの数式が実際にどう動くかが手触りとして掴めます。
import numpy as np
from scipy import stats
def gmm_em_scratch(X, K, n_iter=100, tol=1e-6, seed=0):
"""GMMのEMアルゴリズム スクラッチ実装
X: (N, d) 観測データ
K: クラスタ数
n_iter: 最大反復回数
tol: 対数尤度の改善幅がこれ以下になったら収束判定
"""
rng = np.random.default_rng(seed)
N, d = X.shape
# --- 初期化 ---
# μ: データ点からランダムに K 個を選ぶ
mu = X[rng.choice(N, K, replace=False)].copy()
# Σ: 全データの分散行列で揃える(対角に少し足して正則化)
cov = np.stack([np.cov(X.T) + 0.1 * np.eye(d) for _ in range(K)])
# π: 一様
pi = np.full(K, 1.0 / K)
log_likelihood_history = []
for it in range(n_iter):
# --- E-step: 責任度 γ_{nk} を計算 ---
# 数値安定化のため log-sum-exp で計算
log_weighted = np.zeros((N, K))
for k in range(K):
log_weighted[:, k] = (
np.log(pi[k] + 1e-300)
+ stats.multivariate_normal(mu[k], cov[k],
allow_singular=True).logpdf(X)
)
log_max = log_weighted.max(axis=1, keepdims=True)
log_sum = log_max + np.log(
np.sum(np.exp(log_weighted - log_max), axis=1, keepdims=True)
)
log_gamma = log_weighted - log_sum # 正規化
gamma = np.exp(log_gamma) # (N, K)
# 対数尤度を計算して保存
ll = float(np.sum(log_sum))
log_likelihood_history.append(ll)
# --- M-step: 閉形式更新 ---
Nk = gamma.sum(axis=0) + 1e-12 # 有効サンプル数 (K,)
pi = Nk / N
for k in range(K):
mu[k] = (gamma[:, k:k+1] * X).sum(axis=0) / Nk[k]
diff = X - mu[k]
cov[k] = (gamma[:, k:k+1] * diff).T @ diff / Nk[k]
cov[k] += 1e-6 * np.eye(d) # 縮退回避の正則化
# --- 収束判定 ---
if it > 0 and abs(log_likelihood_history[-1]
- log_likelihood_history[-2]) < tol:
print(f"反復 {it} で収束しました")
break
return {"pi": pi, "mu": mu, "cov": cov,
"gamma": gamma, "history": log_likelihood_history}
このコードの肝は3点です。第一に、責任度の計算で logpdf を使い log-sum-exp で正規化することで、高次元データでも数値が下振れしないようにしています。第二に、Nk + 1e-12 で割って division-by-zero を防ぎ、cov[k] += 1e-6 * np.eye(d) で縮退共分散の発散も抑えています。第三に、対数尤度の改善幅で早期終了を判定しているので、無駄な反復を避けられます。
合成データで動作確認
スクラッチ実装が正しく動くか、3つのガウスから生成した合成データで試します。
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(42)
mus_true = np.array([[-3.0, 0.0], [3.0, 2.0], [0.0, -3.0]])
covs_true = [
np.array([[1.2, 0.4], [0.4, 0.6]]),
np.array([[0.8, -0.5], [-0.5, 1.4]]),
np.array([[1.5, 0.0], [0.0, 0.5]]),
]
pis_true = [0.4, 0.35, 0.25]
N = 600
Ns = (np.array(pis_true) * N).astype(int)
X = np.vstack([
rng.multivariate_normal(mus_true[k], covs_true[k], Ns[k])
for k in range(3)
])
rng.shuffle(X)
# スクラッチGMMで学習
result = gmm_em_scratch(X, K=3, n_iter=50, seed=0)
print("推定された混合係数:", result["pi"])
print("推定された平均:")
print(result["mu"])
このコードを実行すると、混合係数の推定値が真値 [0.40, 0.35, 0.25] にほぼ一致し、平均も真の中心 $(\pm 3, 0)$、$(0, -3)$ 付近に収束することが確認できます。順番はラベルスワップの自由度があるので毎回同じとは限りませんが、集合としては一致します。
scikit-learn で実用的に使う
実用ではスクラッチ実装より scikit-learn の GaussianMixture を使うのが標準的です。複数の初期化、収束判定、共分散タイプの切り替え、AIC/BIC計算など、必要な機能が一通り揃っています。
from sklearn.mixture import GaussianMixture
# K=3 で fit
gmm = GaussianMixture(
n_components=3,
covariance_type='full',
n_init=10, # 10通り初期化して最良を採用
max_iter=200,
tol=1e-4,
reg_covar=1e-6, # 共分散の正則化(特異点回避)
random_state=0,
)
gmm.fit(X)
print("混合係数:", gmm.weights_)
print("平均:\n", gmm.means_)
print("対数尤度:", gmm.score(X) * len(X))
print("AIC:", gmm.aic(X))
print("BIC:", gmm.bic(X))
n_init=10 は強くおすすめのオプションです。初期化への依存をかなり緩和してくれます。score(X) は1サンプルあたりの平均対数尤度を返すので、総和に直すには * len(X) を掛けます。AIC/BICはモデル選択時に他の $K$ と比較する基準となります。
Irisデータセットでの実例
ここからは、より実データに近いシナリオとして有名なIrisデータセットでGMMを試します。Irisデータセットの説明はこちらをご覧ください。
[
【Python】scikit-learnでIrisデータセットを手軽に用意する
機械学習の勉強をしていると、まず最初に登場するのが、Irisデータセットではない…](https://disassemble-channel.com/iris-datasets/)
データセットを準備します。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
iris = load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['label'] = iris.target
Irisデータセットはデータ点数が150点で、系列が4つある以下のようなデータになります。
各特徴量(sepal length, sepal width, petal length, petal width)が4列に並んでおり、label 列が真のアヤメ品種(0/1/2)を表しています。今回はクラスタリングを評価する目的でこの真値を利用しますが、GMM自体は真値を使わずに教師なしで学習することに注意してください。
混合ガウスモデルを利用する際は、データの次元はいくつでも構いませんが、今回はガウス分布の等高線を可視化してわかりやすくするために、2次元に落とします。そのために主成分分析を行います。主成分分析のコードは下記のようになります。
pca = PCA(n_components=2)
data = df[["sepal length (cm)", "sepal width (cm)",
"petal length (cm)", "petal width (cm)"]].values
trans_pca = pca.fit_transform(data)
fig, ax = plt.subplots()
cmap = plt.get_cmap("tab10")
color = [cmap(label) for label in df.label]
ax.set_xlabel("PCA axis 1")
ax.set_ylabel("PCA axis 2")
for i in range(3):
ax.scatter(
trans_pca[df[df.label == i].index, 0],
trans_pca[df[df.label == i].index, 1],
color=cmap(i), s=10, marker="o", label="Label: {}".format(i)
)
ax.legend()
主成分分析した結果を可視化するとこのようになります。第1主成分軸(横軸)で setosa(label=0)が左側に明確に分かれ、versicolor(label=1)と virginica(label=2)が右側に重なり気味に分布しているのが見て取れます。この「線形では完全に分けにくい2クラスタが重なる」状況こそ、GMMの共分散モデリングが活きる場面です。主成分分析の詳しい解説はこちらに掲載しています。
[
主成分分析(PCA)をわかりやすく解説。固有値や固有ベクトルと数学的背景など
主成分分析(Principle Component Analysis, PCA)…](https://disassemble-channel.com/pca_principal_component_analysis/)
IrisデータセットでGMMでクラスタリング
ここまで準備することができたら、実際に混合ガウスモデルを用いて、データのクラスタリングを行なっていきましょう。Irisデータセットは正解ラベルがあるので、先ほど掲載したように、データがうまくこのようにクラスタリングできていれば、GMMがうまく行っていることになります。
GMMはスクラッチでも実装できますが、最適化計算のところで、EMアルゴリズムや変分推論など数学的に少し難解なアルゴリズムを必要とするため、今回はscikit-learnで既に実装されているモジュールを利用します。
GMM自体は、scikit-learnを用いることで、学習までわずか数行の次のコードで行うことができます。trans_pca は主成分分析によって2次元に次元削減したデータセットが格納されています。
from matplotlib.colors import LogNorm
from sklearn import mixture
gmm = mixture.GaussianMixture(n_components=3, covariance_type='full',
n_init=10, random_state=0)
gmm.fit(trans_pca)
labels = gmm.predict(trans_pca)
ここまでで、GMMを用いて推論をすることができました。GMMではハイパーパラメータとして、分類するクラスタ数を割り当てる必要性があります。$K=3$ をどう選んだかというと、Irisは3品種だと既知だから、というのが一番の理由ですが、AIC/BICでこの値が選ばれることも確認できます(後述)。あとは結果を可視化します。
x = np.linspace(-4, 4)
y = np.linspace(-2, 2)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -gmm.score_samples(XX)
Z = Z.reshape(X.shape)
fig, ax = plt.subplots(dpi=150, figsize=(5, 4))
ax.scatter(trans_pca[:, 0], trans_pca[:, 1], s=0.5, c=labels)
cont = ax.contourf(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=100.0),
levels=np.logspace(-1, 3, 20), alpha=0.2,
linestyles='dashed', linewidths=0.5)
ax.scatter(trans_pca[:, 0], trans_pca[:, 1], s=1, c=labels)
ax.set_title("GMMによるクラスタリングと等高線")
結果はこのようになりました。gmm.score_samples は対数密度を返すため、-gmm.score_samples は負の対数密度(高いほど密度が低い)になります。等高線が「データから離れるほど高くなる丘」として見えており、各クラスタの中心が低地(谷)になっています。実際に正解ラベルと比較してみると、かなり高い精度で分類できていることがわかります。
左が真のラベル、右がGMMによるクラスタリング結果です。setosa は完全に分離でき、versicolor と virginica の境界付近で数点の誤分類が見られますが、ほとんどのデータ点が正しく分類されています。これは教師なし学習でここまで正解ラベルを再現できているということで、実用上は十分な結果と言えます。
AIC/BICによる $K$ の決定
最後に、もしIrisの真のクラスタ数を知らなかった場合、AIC/BICで $K=3$ が選ばれるかを確認しておきましょう。
from sklearn.mixture import GaussianMixture
aics, bics = [], []
Ks = range(1, 9)
for K in Ks:
g = GaussianMixture(n_components=K, covariance_type='full',
n_init=10, random_state=0)
g.fit(trans_pca)
aics.append(g.aic(trans_pca))
bics.append(g.bic(trans_pca))
print("AIC最小 K =", Ks[np.argmin(aics)])
print("BIC最小 K =", Ks[np.argmin(bics)])
このコードを実行すると、BIC最小は $K=2$ または $K=3$ になることが多いです(versicolor と virginica が重なり気味なので、BICは保守的に2を選ぶことがあります)。AICは $K=3$ や $K=4$ になることもあります。教師なしの情報量基準だけでクラスタ数を完全に当てるのは難しい場合がある、というのが現実的な教訓です。実用では、ドメイン知識と情報量基準の両方を見て決めるのが定石です。
実用上の注意とチェックリスト
GMMを実応用するとき、押さえておきたいポイントをまとめます。
- データの標準化 — 共分散行列が特徴量のスケールに敏感なので、
StandardScalerで各次元を平均0・分散1に揃えてから学習するのが安全です。生のスケールが極端に違うと、共分散の固有値分解が数値的に不安定になります。 n_init=10を必ず指定 — デフォルトの1回試行では初期化運任せです。複数初期化で最良を採るのが標準。- 高次元では
covariance_typeを見直す — 次元 $d$ が大きい場合、'full'ではパラメータが爆発します。'diag'や'tied'を試す、あるいはPCAで次元削減してから当てる。 - 可視化はラベルだけでなく確率も —
predict_proba(X)で各点のソフト所属度が得られます。境界近くの不確かさを下流の意思決定に反映できます。 - 特異点が起きたら
reg_covarを上げる — デフォルトは $10^{-6}$ ですが、必要なら $10^{-3}$ など大きめにする。 - ベイズGMM(
BayesianGaussianMixture)を第一選択に — 特異点なし・クラスタ数自動決定・正則化込み。最尤GMMの上位互換と考えてよいケースが多い。
これらの注意点を踏まえれば、GMMは多くの場面で頑健に動く強力なツールとなります。
まとめ
本記事では、ガウス混合モデル(GMM)について、定義からEM法の導出、実装、そしてベイズ的拡張まで一気通貫で解説しました。
- GMMの定義: $K$ 個のガウス分布の重み付き和で柔軟な確率分布を表現。混合係数 $\bm{\pi}$、平均 $\bm{\mu}_k$、共分散 $\bm{\Sigma}_k$ の3組がパラメータ。
- 潜在変数表現: 各データ点に「どの成分から生成されたか」の潜在変数 $z$ を導入することでEMが自然に展開できる。
- EMアルゴリズム: E-stepで責任度 $\gamma_{nk}$ を計算し、M-stepで $\bm{\pi}, \bm{\mu}_k, \bm{\Sigma}_k$ を閉形式で更新。対数尤度の単調増加が保証される。
- K-meansとの関係: K-meansは共分散を $\sigma^2 \bm{I}$ で固定した極限。GMMは形状・ソフト割り当て・サイズ可変を獲得。
- 初期化問題: K-means++ または
n_initで複数試行が標準的対策。 - クラスタ数の選択: 対数尤度は単調増加なので不可。AIC/BIC で罰則を加えて比較。BICは一致性を持つ。
- 縮退共分散問題: 共分散が1点に貼り付くと尤度が発散。
reg_covarで下限を設けるか、ベイズGMMで根本解決。 - 共分散の幾何: 固有値=軸長、固有ベクトル=軸方向。共分散タイプ(
full/diag/spherical)で柔軟性と過学習リスクのバランスを取る。 - ベイズGMM・変分GMM: 事前分布で特異点を回避し、不要な成分を自動的に消去。
BayesianGaussianMixtureがscikit-learnで使える。 - 実装: スクラッチEMで仕組みを理解した上で、実用では
GaussianMixtureまたはBayesianGaussianMixtureを使う。
GMMは古典的な確率モデルですが、その背後にある「潜在変数+EM」「混合分布」「変分推論」というアイデアは、現代の生成モデル(VAE、拡散モデル、Mixture-of-Experts)に至るまで脈々と受け継がれています。GMMを通じてこれらの基礎を体感しておくと、より高度なモデルを学ぶときに見通しが立ちます。
次のステップとして、以下の記事も参考にしてください。