機械学習の勉強をしていて、初学者がつまづくポイントは多々あると思います。EMアルゴリズムはその典型で、数式だけを追っていても自分の中でしっくりこない人も多いのではないでしょうか。
今回は、EMアルゴリズムの概念を読むだけでわかるように、できるだけ丁寧に書き下します。
本記事の内容
- EMアルゴリズムが解決する問題
- 不完全データと潜在変数の概念
- EステップとMステップの役割
- 混合ガウスモデルでの具体例とPython実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
EMアルゴリズムが解決する問題
通常の最尤推定では、パラメータ $\theta$ について対数尤度を直接最大化します。しかし、データの一部が観測できない（潜在変数がある）場合、対数尤度を直接最大化することが困難になります。
EMアルゴリズムは、このような潜在変数を含むモデルのパラメータを最尤推定するためのアルゴリズムです。
不完全データと潜在変数
具体例で考えましょう。2つの工場AとBで同じ製品を作っているとします。各工場の製品の重さはそれぞれ異なる正規分布に従います。
$$ \text{工場A}: \mathcal{N}(\mu_A, \sigma_A^2), \quad \text{工場B}: \mathcal{N}(\mu_B, \sigma_B^2) $$
もし各製品にどちらの工場で作られたかのラベルが付いていれば（完全データ）、工場ごとにデータを分けて最尤推定すれば済みます。
しかし、ラベルがなく製品の重さだけが分かる場合（不完全データ）、どのデータがどちらの工場由来かが分かりません。この「どちらの工場か」という情報が潜在変数です。
EMアルゴリズムの概要
EMアルゴリズムは以下の2ステップを交互に繰り返します。
Eステップ（Expectation Step）: 現在のパラメータ推定値を使って、潜在変数の期待値（負担率）を計算する
Mステップ（Maximization Step）: Eステップで求めた負担率を使って、パラメータを更新する
これを収束するまで繰り返すことで、パラメータが最尤推定値に収束していきます。
混合ガウスモデルでの定式化
混合ガウスモデル（Gaussian Mixture Model; GMM）は、EMアルゴリズムの最も典型的な応用例です。
$K$ 個のガウス分布の重み付き和として、
$$ \begin{equation} p(x|\theta) = \sum_{k=1}^{K} \pi_k \mathcal{N}(x|\mu_k, \sigma_k^2) \end{equation} $$
ここで $\pi_k$ は混合係数（$\sum_k \pi_k = 1$, $\pi_k \geq 0$）、$\theta = \{\pi_k, \mu_k, \sigma_k^2\}_{k=1}^{K}$ がパラメータです。
Eステップ
データ点 $x_n$ がクラスタ $k$ に属する負担率（responsibility）$\gamma_{nk}$ を計算します。
$$ \begin{equation} \gamma_{nk} = \frac{\pi_k \mathcal{N}(x_n|\mu_k, \sigma_k^2)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(x_n|\mu_j, \sigma_j^2)} \end{equation} $$
$\gamma_{nk}$ は、データ点 $x_n$ がクラスタ $k$ から生成された確率を表しています。
Mステップ
負担率を使ってパラメータを更新します。
$$ \begin{align} N_k &= \sum_{n=1}^{N} \gamma_{nk} \\ \mu_k^{\text{new}} &= \frac{1}{N_k} \sum_{n=1}^{N} \gamma_{nk} x_n \\ (\sigma_k^2)^{\text{new}} &= \frac{1}{N_k} \sum_{n=1}^{N} \gamma_{nk} (x_n – \mu_k^{\text{new}})^2 \\ \pi_k^{\text{new}} &= \frac{N_k}{N} \end{align} $$
ここで $N_k$ はクラスタ $k$ に割り当てられたデータの有効数です。
Python での実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
def em_gmm_1d(data, K, max_iter=100, tol=1e-6):
"""1次元混合ガウスモデルのEMアルゴリズム"""
N = len(data)
# 初期化
np.random.seed(0)
mu = np.random.choice(data, K)
sigma2 = np.ones(K) * np.var(data)
pi = np.ones(K) / K
gamma = np.zeros((N, K))
log_likelihoods = []
for iteration in range(max_iter):
# Eステップ: 負担率の計算
for k in range(K):
gamma[:, k] = pi[k] * norm.pdf(data, mu[k], np.sqrt(sigma2[k]))
gamma /= gamma.sum(axis=1, keepdims=True)
# Mステップ: パラメータ更新
N_k = gamma.sum(axis=0)
for k in range(K):
mu[k] = np.sum(gamma[:, k] * data) / N_k[k]
sigma2[k] = np.sum(gamma[:, k] * (data - mu[k])**2) / N_k[k]
pi[k] = N_k[k] / N
# 対数尤度の計算
ll = 0
for n in range(N):
p = sum(pi[k] * norm.pdf(data[n], mu[k], np.sqrt(sigma2[k])) for k in range(K))
ll += np.log(p + 1e-300)
log_likelihoods.append(ll)
# 収束判定
if iteration > 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) < tol:
break
return mu, sigma2, pi, gamma, log_likelihoods
# データ生成（3つのガウス分布の混合）
np.random.seed(42)
data = np.concatenate([
np.random.normal(1, 0.5, 100),
np.random.normal(4, 0.8, 150),
np.random.normal(7, 0.6, 80)
])
np.random.shuffle(data)
K = 3
mu, sigma2, pi, gamma, log_liks = em_gmm_1d(data, K)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# 左: フィッティング結果
ax1 = axes[0]
x_range = np.linspace(-2, 10, 500)
ax1.hist(data, bins=40, density=True, alpha=0.4, color='gray', label='Data')
total_pdf = np.zeros_like(x_range)
colors = ['blue', 'green', 'red']
for k in range(K):
pdf = pi[k] * norm.pdf(x_range, mu[k], np.sqrt(sigma2[k]))
total_pdf += pdf
ax1.plot(x_range, pdf, color=colors[k], linewidth=1.5,
label=f'k={k+1}: mu={mu[k]:.2f}, sigma={np.sqrt(sigma2[k]):.2f}')
ax1.plot(x_range, total_pdf, 'k--', linewidth=2, label='GMM total')
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('Density', fontsize=12)
ax1.set_title('EM Algorithm: GMM Fitting', fontsize=13)
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)
# 中央: 対数尤度の推移
ax2 = axes[1]
ax2.plot(log_liks, 'b-', linewidth=2)
ax2.set_xlabel('Iteration', fontsize=12)
ax2.set_ylabel('Log-likelihood', fontsize=12)
ax2.set_title('Log-Likelihood Convergence', fontsize=13)
ax2.grid(True, alpha=0.3)
# 右: 負担率の可視化
ax3 = axes[2]
cluster_assignment = gamma.argmax(axis=1)
for k in range(K):
mask = cluster_assignment == k
ax3.scatter(data[mask], gamma[mask, k], s=10, alpha=0.5, color=colors[k], label=f'Cluster {k+1}')
ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('Responsibility', fontsize=12)
ax3.set_title('Cluster Responsibilities', fontsize=13)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフでは3つのガウス分布の混合がうまくフィッティングされていること、中央のグラフでは対数尤度が単調に増加して収束していること、右のグラフでは各データ点のクラスタへの帰属確率が確認できます。
EMアルゴリズムの特徴
- 対数尤度が単調増加する: 各イテレーションで対数尤度は必ず増加する（減少しない）
- 局所最適解に収束する: 大域的最適解への収束は保証されない。初期値に依存する
- 収束が遅い場合がある: 特にクラスタが重なっている場合
まとめ
本記事では、EMアルゴリズムについて解説しました。
- EMアルゴリズムは潜在変数を含むモデルのパラメータを最尤推定するための反復手法
- Eステップで潜在変数の期待値（負担率）を計算し、Mステップでパラメータを更新する
- 混合ガウスモデルはEMアルゴリズムの最も典型的な応用例
- 対数尤度は単調増加するが、局所最適解に収束する可能性がある
次のステップとして、以下の記事も参考にしてください。