ガウス混合モデル(GMM) と EMアルゴリズムの可視化(MATLAB実装) #初心者 – Qiita

Gaussian Mixture Model (ガウス混合モデル、GMM) は、データが複数のガウス分布（正規分布）から生成されていると仮定する確率モデルです。自然現象や社会現象の多くはこのガウス分布に従います。以下の図は一次元のガウス分布を示しています。

GMMをまた別の言い方で表現すると、一つのデータ集合を複数のガウス分布の寄せ集めで表現しようというモデルです。以下の図はそのイメージを表しています。緑のような分布のデータが存在するときに、赤と青のガウス分布を重ね合わせであると解釈することができます。

各ガウス分布は一つのクラスタ（データのグループ）を表し、GMMを用いることでデータをクラスタごとに分けたり（クラスタリング）、データの分布そのものを推定したりすることができます。以下の図はGMMによってクラスタリングを行った結果を示しています。

クラスタリングは教師なし学習の一種で、ラベルのないデータを性質の似たグループに分類する手法です。GMMは特にソフトクラスタリングによく使われ、各データ点が複数のクラスタに属する確率を計算します。これは、必ずどれか一つのクラスタにデータを割り当てるハードクラスタリング（例: k-means）とは対照的です。例えばGMMでは、境界付近のデータは複数クラスタにまたがる曖昧な所属確率を持つため、不確実性を表現できます。またGMMは各クラスタに共分散を持たせることで、データの非円形の分布（例えば細長い楕円形状のクラスタ）も表現でき、k-meansより柔軟にクラスタを表現できます。

応用例: GMMはクラスタリング以外にも、確率モデルとしてデータの分布を推定する密度推定（generative model）に使われます。たとえば音声データの特徴抽出や異常検知、画像のピクセル値のクラスタリング（セグメンテーション）など、様々な領域で利用されています。

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

Expectation-Maximization (期待値最大化)アルゴリズムは、GMMのように潜在変数（クラスタの割り当て$z_n$など観測されない隠れた変数）が存在するモデルのパラメータ推定によく用いられる反復型アルゴリズムです。基本的なアイデアは、「もし各データ点がどのガウス成分から来たか（クラスタ所属）が分かっていればパラメータは簡単に推定でき、一方でパラメータが分かっていれば各点のクラスタ所属確率は簡単に計算できる」という点にあります。そこでEステップ（期待値計算）では現在のパラメータで各点のクラスタ所属確率を計算し、Mステップ（最大化）ではその確率をもとにパラメータを更新します。このEとMのステップを交互に繰り返すことで尤度を徐々に高め、パラメータ推定を行います。EMアルゴリズムは各反復で尤度を増大させることが保証されており、収束すると（一般に局所最適な）推定値が得られます。

アルゴリズムの概要

1. 初期化

   各クラスタの平均$\boldsymbol{\mu}_k$, 共分散$\Sigma_k$, 混合係数$\pi_k$に初期値を設定します（例えばデータからランダムに初期中心を選ぶなど）。

2. Eステップ (期待値ステップ)

   現在のパラメータで、各データ点が各クラスタに属する事後確率（責任度$\gamma_{nk}$）を計算します。

3. Mステップ (最大化ステップ)

   Eステップで求めた事後確率を重みに、パラメータ（${\mu}_k, \Sigma_k, \pi_k$）を再計算します。

4. 繰り返し

   パラメータの変化が僅かになるまで、Eステップと Mステップを交互に繰り返します。収束したパラメータが最尤推定解となります。

3.1 Eステップ: 事後確率の計算

\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}

P(\mathbf{x}_n) = \sum_{k=1}^{K} \pi_k \; \mathcal{N}(\mathbf{x}_n \mid \boldsymbol{\mu}_k, \Sigma_k) \tag{1}

\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}

%% ---------- データ合成 ----------
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