1. GMMとは
Gaussian Mixture Model (ガウス混合モデル、GMM) は、データが複数のガウス分布（正規分布）から生成されていると仮定する確率モデルです。自然現象や社会現象の多くはこのガウス分布に従います。以下の図は一次元のガウス分布を示しています。
GMMをまた別の言い方で表現すると、一つのデータ集合を複数のガウス分布の寄せ集めで表現しようというモデルです。以下の図はそのイメージを表しています。緑のような分布のデータが存在するときに、赤と青のガウス分布を重ね合わせであると解釈することができます。
各ガウス分布は一つのクラスタ（データのグループ）を表し、GMMを用いることでデータをクラスタごとに分けたり（クラスタリング）、データの分布そのものを推定したりすることができます。以下の図はGMMによってクラスタリングを行った結果を示しています。
クラスタリングは教師なし学習の一種で、ラベルのないデータを性質の似たグループに分類する手法です。GMMは特にソフトクラスタリングによく使われ、各データ点が複数のクラスタに属する確率を計算します。これは、必ずどれか一つのクラスタにデータを割り当てるハードクラスタリング（例: k-means）とは対照的です。例えばGMMでは、境界付近のデータは複数クラスタにまたがる曖昧な所属確率を持つため、不確実性を表現できます。またGMMは各クラスタに共分散を持たせることで、データの非円形の分布（例えば細長い楕円形状のクラスタ）も表現でき、k-meansより柔軟にクラスタを表現できます。
応用例: GMMはクラスタリング以外にも、確率モデルとしてデータの分布を推定する密度推定（generative model）に使われます。たとえば音声データの特徴抽出や異常検知、画像のピクセル値のクラスタリング（セグメンテーション）など、様々な領域で利用されています。
2. GMMの数式表現とパラメータ
GMMでは、データの確率分布を複数のガウス分布の重ね合わせ（混合）として定義します。数式で表すと、データ点$\mathbf{x}_n$の確率密度は次のようになります:
P(\mathbf{x}_n) = \sum_{k=1}^{K} \pi_k \; \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k) \tag{1}
ここで$K$はガウス分布の混合成分の個数（クラスタ数に相当）です。各項$\mathcal{N}(\mathbf{x}|\boldsymbol{\mu_k}, \Sigma_k)$は平均${\mu}_k$、共分散行列$\Sigma_k$を持つガウス分布（正規分布）の確率密度です。$\pi_k$は各成分$k$の混合係数（重み）で、$\pi_k \ge 0$かつ$\sum_{k=1}^{K}\pi_k = 1$を満たします。この混合係数$\pi_k$は「データ点が成分$k$から生成される事前確率」を表します。式(1)は「各ガウス分布に$\pi_k$の確率で従いデータが生成される」というモデルを意味しています。
パラメータの意味
GMMは以下のパラメータから構成されます:
- 平均 $\boldsymbol{\mu}_k$ – 第$k$成分ガウス分布の平均ベクトル。クラスタの中心を表します。
- 共分散 $\Sigma_k$ – 第$k$成分ガウス分布の共分散行列。クラスタ内のデータの広がりや形状（楕円の方向やサイズ）を表します。1次元の場合は分散$\sigma_k^2$になります。
- 混合係数 $\pi_k$ – 第$k$成分の混合比率。全データ中でクラスタ$k$が占める割合（事前確率）を表します。全$\pi_k$の総和は1になります。
以上のパラメータ集合 ${\pi_k, \boldsymbol{\mu}_k, \Sigma_k}、{k=1..K}$ が定まれば、式(1)によってデータ分布（混合ガウス分布）が定義されます。なお、各ガウス分布$\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \Sigma)$とは、平均${\mu}$と共分散$\Sigma$によって形が決まる「鐘形」の連続確率分布です（1次元ならおなじみの正規分布の形）。
共分散(分散)の理解を深めるために、一次元のガウス分布において共分散を変化さえた時の様子を以下に示します。共分散の増加し従い、データが広がっている様子が確認できます。
以下は混合係数の理解を深めるために、一次元のガウス分布において混合係数を変化させた時の、データの分布を以下に示します。実際のデータの分布と緑の曲線が近くなるように混合係数を推定します。
尤度と最尤推定
モデルのパラメータが決まれば、データ集合$X={\mathbf{x}_1,...,\mathbf{x}_N}$が得られる確率（尤度）を計算できます。その尤度は各データ点の確率を掛け合わせたものになります。例えばパラメータセット$\theta=({\pi_k},{\boldsymbol{\mu}_k},{\Sigma_k})$のもとでの尤度関数は以下のように表されます。
L(\theta) = \prod_{n=1}^N P(\mathbf{x}_n \mid \theta) = \prod_{n=1}^{N}\sum_{k=1}^{K} \pi_k\, \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k)
モデルをデータにフィットさせるには、この尤度$L(\theta)$が最大となるようパラメータ$\theta$を調整します（最尤推定）。尤度 (likelihood) とはパラメータのもとでデータが観測される確率のことです。対数を取って微分による解析的な最大化が行われますが、混合分布の尤度はガウス分布の和を含むためそのような解析が困難です。この最適化問題を解くために用いられるのがExpectation-Maximization (EM)アルゴリズムです。
$$
\log L(\theta) = \sum_{n=1}^{N} \log \left( \sum_{k=1}^{K} \pi_k, \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k) \right)
$$
3. EMアルゴリズムによるパラメータ推定
Expectation-Maximization (期待値最大化)アルゴリズムは、GMMのように潜在変数（クラスタの割り当て$z_n$など観測されない隠れた変数）が存在するモデルのパラメータ推定によく用いられる反復型アルゴリズムです。基本的なアイデアは、「もし各データ点がどのガウス成分から来たか（クラスタ所属）が分かっていればパラメータは簡単に推定でき、一方でパラメータが分かっていれば各点のクラスタ所属確率は簡単に計算できる」という点にあります。そこでEステップ（期待値計算）では現在のパラメータで各点のクラスタ所属確率を計算し、Mステップ（最大化）ではその確率をもとにパラメータを更新します。このEとMのステップを交互に繰り返すことで尤度を徐々に高め、パラメータ推定を行います。EMアルゴリズムは各反復で尤度を増大させることが保証されており、収束すると（一般に局所最適な）推定値が得られます。
アルゴリズムの概要
-
初期化
各クラスタの平均$\boldsymbol{\mu}_k$, 共分散$\Sigma_k$, 混合係数$\pi_k$に初期値を設定します（例えばデータからランダムに初期中心を選ぶなど）。
-
Eステップ (期待値ステップ)
現在のパラメータで、各データ点が各クラスタに属する事後確率（責任度$\gamma_{nk}$）を計算します。
-
Mステップ (最大化ステップ)
Eステップで求めた事後確率を重みに、パラメータ（${\mu}_k, \Sigma_k, \pi_k$）を再計算します。
-
繰り返し
パラメータの変化が僅かになるまで、Eステップと Mステップを交互に繰り返します。収束したパラメータが最尤推定解となります。
3.1 Eステップ: 事後確率の計算
Eステップでは、現在のパラメータ値を固定して各データ点 $\mathbf{x}_n$がクラスタ$k$に属する事後確率（責任度とも呼ぶ）$\gamma_{nk}$を計算します。ベイズの定理に基づき、この確率は以下のようになります。
\gamma_{nk} \;=\; P(z_n = k \mid \mathbf{x}_n) \;=\; \frac{\pi_k \; \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \; \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_j, \Sigma_j)} \tag{2}
ここで $z_n$ はデータ点$\mathbf{x}_n$のクラスタラベル（潜在変数）です。$\gamma_{nk}$ は「データ点 $\mathbf{x}_n$ がクラスタ $k$ に属する確率」を表し、$\sum_{k=1}^K \gamma_{nk} = 1$ になります。この式の分子$\pi_k \mathcal{N}(\mathbf{x}_n|\mu_k,\Sigma_k)$は「クラスタ$k$の事前確率$\pi_k$」と「$\mathbf{x}_n$がクラスタ$k$のガウス分布から生成される尤度」の積であり、分母は全クラスタについてその積を足し合わせたものです。したがって$\gamma_{nk}$は事後確率＝「$\mathbf{x}_n$が与えられた下で$z_n=k$である確率」を表現しています（クラスタ$k$の事後確率）。各データ点についてこの計算を行うと、$N\times K$の責任度行列$\Gamma=(\gamma_{nk})$が得られます。
3.2. Eステップ具体例
ここではEステップの計算を実際の数値での例をもとに確認していきます。一次元の以下のデータを考えます。
P(\mathbf{x}_n) = \sum_{k=1}^{K} \pi_k \; \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k) \tag{1}
またガウス分布は2つとして、混合係数、平均、分散の初期値を以下のように与えます。
$$
\pi^{(0)} = (0.4,;0.6), \quad\mu^{(0)} = \bigl(-2,;2\bigr), \quad\sigma^2{}^{(0)} = \bigl(1^2,;1^2\bigr)
$$
それぞれの点に対応する確率密度の値を2つのガウス分布ごとに計算します。
| $x_n$ | ${pdf}_1$ | $pdf_2$ | $\gamma_{n1}$ | $\gamma_{n2}$ |
|---|---|---|---|---|
| -3 | 0.241971 | 1e-06 | 0.999991 | 9e-06 |
| -2 | 0.398942 | 0.000134 | 0.999497 | 0.000503 |
| -1 | 0.241971 | 0.004432 | 0.973261 | 0.026739 |
| 1 | 0.004432 | 0.241971 | 0.012063 | 0.987937 |
| 2 | 0.000134 | 0.398942 | 0.000224 | 0.999776 |
| 3 | 1e-06 | 0.241971 | 4e-06 | 0.999996 |
3.3. Mステップ: パラメータの更新
Mステップでは、Eステップで求めた$\gamma_{nk}$を重みとしてパラメータを再推定します。直感的には、各データ点を実際にクラスタ$k$に属する「重み$\gamma_{nk}$」だけ属するとみなし、その重み付き平均や分散を計算するイメージです。具体的なクラスタkに対する更新式は以下のとおりです:
-
混合係数の更新
$$
\pi_k^{\text{new}} = \frac{1}{N}\sum_{n=1}^{N} \gamma_{nk} \tag{3a}
$$各クラスタに属するデータの割合（平均的な責任度）で$\pi_k$を更新します。
-
平均ベクトルの更新
\boldsymbol{\mu}_k^{\text{new}} = \frac{\sum_{n=1}^{N} \gamma_{nk}\, \mathbf{x}_n}{\sum_{n=1}^{N} \gamma_{nk}} \tag{3b}クラスタ$k$の責任度を重みとしたデータ点の加重平均で$\mu_k$を更新します。各クラスタの「重心」を計算し直すイメージです。
-
共分散行列の更新:
\Sigma_k^{\text{new}} = \frac{\sum_{n=1}^{N} \gamma_{nk}\, (\mathbf{x}_n - \boldsymbol{\mu}_k^{\text{new}})(\mathbf{x}_n - \boldsymbol{\mu}_k^{\text{new}})^\top}{\sum_{n=1}^{N} \gamma_{nk}} \tag{3c}責任度を重みとした分散共分散で$\Sigma_k$を更新します。各クラスタ内のデータの広がりを現在の割り当てに合わせて再計算します。
以上により、新しいパラメータ${\pi_k^{new}, \mu_k^{new}, \Sigma_k^{new}}$が得られます。更新後は再びEステップに戻り、$\gamma_{nk}$を再計算→パラメータ再更新…と繰り返していきます。この反復によってモデルの対数尤度は徐々に改善し、十分繰り返すとパラメータ推定が収束します。ただし、EMアルゴリズムは一般に局所最適に収束することがあるため、異なる初期値で複数回試行してベストの結果を選ぶこともよく行われます。
3.4. Mステップ具体例
先ほどの計算結果を以下にまた示します。
| $x_n$ | $\gamma_{n1}$ | $\gamma_{n2}$ |
|---|---|---|
| -3 | 0.999991 | 9e-06 |
| -2 | 0.999497 | 0.000503 |
| -1 | 0.973261 | 0.026739 |
| 1 | 0.012063 | 0.987937 |
| 2 | 0.000224 | 0.999776 |
| 3 | 4e-06 | 0.999996 |
先ほどのEステップでの計算結果を用いて、M ステップで重み付き平均・分散を計算し，次のように更新します。
\begin{aligned}
N_k &= \sum_{n=1}^{6} \gamma_{nk}
&= (2.99,\;3.02)\\[2pt]
\pi_k^{\text{new}} &= \frac{N_k}{6}
&= (0.50,\;0.50) \\[2pt]
\mu_k^{\text{new}} &= \frac{\sum_n^6 \gamma_{nk}x_n}{N_k}
&= (-2.00,\; 1.98) \\[2pt]
\sigma_k^{2\,\text{new}} &= \frac{\sum_n^6 \gamma_{nk}(x_n-\mu_k^{\text{new}})^2}{N_k}
&= (0.70,\;0.74)
\end{aligned}
4. MATLABによる2次元クラスタリング
シンプルな二次元データにGMMを適用し、EMアルゴリズムの手順を具体的に見てみましょう。ここでは、3つのグループに分かれた点群データを想定します（例えば中心が異なる3つのガウス分布からサンプルしたデータ）。最初にクラスタの平均などパラメータを適当に初期化し、EMアルゴリズムによってそれらがどのように収束していくかを観察します。
以下の動画はEMアルゴリズムが繰り返されるにつれてクラスタリングされていく様子を示しています。
×印がランダムに選んだクラスタ中心の初期値を示しています。データ点は各初期クラスタに属する確率が高い方の色でプロットされています。反復回数の少ない初期の段階では2つの中心が近いため正しく分離できていません。反復が増加するにつれて赤いクラスターの中心が右上に移動していくことが確認できます。
アルゴリズム収束後で、それぞれの中心が3つの真のデータ分布の中心に一致し、データは正しく3クラスタに分離されています。それぞれのクラスタは楕円形のガウス分布でモデル化されており、境界付近の点は色が混合され所属確率が分散している様子が確認できます。
このように、Eステップでは各点の事後確率$\gamma_{nk}$が計算され（図中の色の割当てに相当）、Mステップではそれに基づきクラスタ中心（平均$\mu_k$）が移動していく様子が確認できます。最終的に混合ガウス2成分のパラメータがデータから適切に推定されています。
以下の図は GMM の 3 次元サーフェスプロットを示しています。これは、平面上の各点 $(x,y)$ における確率密度 $p(x,y)$を高さとして可視化したものであり、各山はモデルが学習した混合ガウス成分の平均 $\mu_k$ を中心に、混合比 $\pi_k$ や共分散 $\Sigma_k$ に応じた広がりとともに現れています。成分間の重なりや分離具合、モデルが捉えた密度構造を直感的かつ立体的に理解できます
5.コード
最後に今回MATLABでEMアルゴリズムを実行した際のコードを以下に示します
%% ---------- データ合成 ----------
N = 500;
mu = [0 0; 1.5 1.2; -1.2 1.8];
sigma = cat(3, ...
[2.0 0.7; 0.7 1.5], ...
[1.8 -0.4; -0.4 2.0], ...
[2.2 0.1; 0.1 1.6]);
X = []; trueIdx = [];
for k = 1:3
Xk = mvnrnd(mu(k,:), sigma(:,:,k), N);
X = [X; Xk];
trueIdx = [trueIdx; k*ones(N,1)];
end
%% ---------- EM 初期化 ----------
K = 3; [N, D] = size(X);
idx0 = randi(K, N, 1); % 初期クラスタ割り当て
pi_k = rand(K,1); pi_k = pi_k / sum(pi_k);
mu_k = X(randsample(N,K), :); % 初期平均
Sigma = zeros(D,D,K);
for k = 1:K, v = 0.5+rand(1,D); Sigma(:,:,k) = diag(v); end
maxIter = 500;
%% ---------- EM 反復 ----------
for it = 1:maxIter
% === E-step ===
resp = zeros(N, K);
for k = 1:K
resp(:,k) = pi_k(k) * mvnpdf(X, mu_k(k,:), Sigma(:,:,k));
end
resp = resp ./ sum(resp, 2); % 責任度正規化
% === M-step ===
Nk = sum(resp, 1); % 各成分の重み
pi_k = Nk / N;
for k = 1:K
mu_k(k,:) = (resp(:,k)' * X) / Nk(k);
Xc = X - mu_k(k,:);
Sigma(:,:,k) = (Xc' * (Xc .* resp(:,k))) / Nk(k) + 1e-6 * eye(D);
end
end
6. 参考
Gaussian Mixture Model | Brilliant Math & Science Wiki
In Depth: Gaussian Mixture Models | Python Data Science Handbook