EMアルゴリズム(Expectation-Maximization algorithm)は、潜在変数を含む確率モデルのパラメータを最尤推定する手法です。特に混合ガウスモデル(Gaussian Mixture Model, GMM)のパラメータ推定で広く利用されています。
通常の最尤推定では、対数尤度関数を微分して解析的にパラメータを求めますが、潜在変数が存在する場合、対数の中に和が入る形になり、解析的な解が得られません。EMアルゴリズムは、この問題をEステップ(期待値計算)とMステップ(最大化)の反復で解決します。
本記事の内容
- 混合ガウスモデルの定義と確率密度関数
- 対数尤度関数の導出と直接最適化が困難な理由
- EMアルゴリズムのEステップとMステップの導出
- Pythonでのスクラッチ実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 正規分布とは？定義やグラフの見方を徹底解説
- 最尤推定法の基本的な考え方
混合ガウスモデルとは
混合ガウスモデル(GMM)は、データが複数のガウス分布(正規分布)の重ね合わせから生成されると仮定するモデルです。
イメージとしては、データの背後に $K$ 個のクラスタがあり、それぞれのクラスタがガウス分布で表現されると考えます。各データ点は、いずれかのクラスタから生成されたものとみなします。
混合ガウスモデルの確率密度関数
$K$ 個のガウス分布からなる混合ガウスモデルの確率密度関数は以下のように定義されます。
$$ p(\bm{x} | \bm{\pi}, \bm{\mu}, \bm{\Sigma}) = \sum_{k=1}^{K} \pi_k \mathcal{N}(\bm{x} | \bm{\mu}_k, \bm{\Sigma}_k) $$
ここで、
- $\pi_k$: 第 $k$ 成分の混合比率（$\sum_{k=1}^{K} \pi_k = 1$, $\pi_k \geq 0$）
- $\bm{\mu}_k$: 第 $k$ 成分の平均ベクトル
- $\bm{\Sigma}_k$: 第 $k$ 成分の共分散行列
- $\mathcal{N}(\bm{x} | \bm{\mu}_k, \bm{\Sigma}_k)$: 多変量正規分布の確率密度関数
多変量正規分布の確率密度関数は、
$$ \mathcal{N}(\bm{x} | \bm{\mu}_k, \bm{\Sigma}_k) = \frac{1}{(2\pi)^{d/2} |\bm{\Sigma}_k|^{1/2}} \exp\left( -\frac{1}{2} (\bm{x} – \bm{\mu}_k)^T \bm{\Sigma}_k^{-1} (\bm{x} – \bm{\mu}_k) \right) $$
です。
対数尤度関数の導出
データ $\bm{X} = \{ \bm{x}_1, \bm{x}_2, \dots, \bm{x}_N \}$ が独立に混合ガウスモデルから生成されたとすると、尤度関数は、
$$ p(\bm{X} | \bm{\pi}, \bm{\mu}, \bm{\Sigma}) = \prod_{n=1}^{N} \left( \sum_{k=1}^{K} \pi_k \mathcal{N}(\bm{x}_n | \bm{\mu}_k, \bm{\Sigma}_k) \right) $$
対数をとると、
$$ \ln p(\bm{X} | \bm{\pi}, \bm{\mu}, \bm{\Sigma}) = \sum_{n=1}^{N} \ln \left( \sum_{k=1}^{K} \pi_k \mathcal{N}(\bm{x}_n | \bm{\mu}_k, \bm{\Sigma}_k) \right) $$
ここで問題となるのは、$\ln$ の中に $\sum$ が入っていることです。単一のガウス分布であれば対数と積が消し合って解析解が得られますが、混合モデルではこれが困難になります。
EMアルゴリズムの導出
潜在変数の導入
各データ点 $\bm{x}_n$ がどのクラスタから生成されたかを表す潜在変数 $z_n \in \{1, 2, \dots, K\}$ を導入します。$z_n = k$ は $\bm{x}_n$ が第 $k$ クラスタから生成されたことを意味します。
Eステップ: 負担率の計算
Eステップでは、現在のパラメータ推定値の下で、各データ点が各クラスタに属する事後確率（負担率, responsibility）を計算します。
ベイズの定理を用いると、
$$ \begin{align} \gamma(z_{nk}) &= p(z_n = k | \bm{x}_n) \\ &= \frac{p(z_n = k) \cdot p(\bm{x}_n | z_n = k)}{p(\bm{x}_n)} \\ &= \frac{\pi_k \mathcal{N}(\bm{x}_n | \bm{\mu}_k, \bm{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(\bm{x}_n | \bm{\mu}_j, \bm{\Sigma}_j)} \end{align} $$
$\gamma(z_{nk})$ は「データ点 $\bm{x}_n$ が第 $k$ クラスタから生成された確率」を表し、負担率と呼ばれます。
Mステップ: パラメータの更新
Mステップでは、Eステップで計算した負担率を用いて、各パラメータを更新します。
まず、第 $k$ クラスタの実効的なデータ数を $N_k$ とします。
$$ N_k = \sum_{n=1}^{N} \gamma(z_{nk}) $$
各パラメータの更新式は以下の通りです。
平均ベクトルの更新:
$$ \bm{\mu}_k^{\text{new}} = \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) \bm{x}_n $$
共分散行列の更新:
$$ \bm{\Sigma}_k^{\text{new}} = \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) (\bm{x}_n – \bm{\mu}_k^{\text{new}})(\bm{x}_n – \bm{\mu}_k^{\text{new}})^T $$
混合比率の更新:
$$ \pi_k^{\text{new}} = \frac{N_k}{N} $$
EステップとMステップを交互に繰り返すことで、対数尤度が単調増加し、局所最適解に収束します。
Pythonでの実装
EMアルゴリズムを用いたGMMをPythonでスクラッチ実装します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
np.random.seed(42)
# --- データ生成（3つのクラスタ） ---
n_samples = 300
true_means = [np.array([0, 0]), np.array([5, 5]), np.array([5, 0])]
true_covs = [np.array([[1, 0.5], [0.5, 1]]),
np.array([[1, -0.3], [-0.3, 0.8]]),
np.array([[0.8, 0], [0, 1.2]])]
true_weights = [0.3, 0.4, 0.3]
X = []
true_labels = []
for k in range(3):
n_k = int(n_samples * true_weights[k])
X.append(np.random.multivariate_normal(true_means[k], true_covs[k], n_k))
true_labels.extend([k] * n_k)
X = np.vstack(X)
true_labels = np.array(true_labels)
# --- EMアルゴリズムの実装 ---
class GaussianMixtureEM:
def __init__(self, n_components, max_iter=100, tol=1e-6):
self.K = n_components
self.max_iter = max_iter
self.tol = tol
def fit(self, X):
N, d = X.shape
# 初期化: K-means的にランダムにデータ点を選ぶ
indices = np.random.choice(N, self.K, replace=False)
self.means = X[indices].copy()
self.covs = [np.eye(d) for _ in range(self.K)]
self.weights = np.ones(self.K) / self.K
self.log_likelihoods = []
for iteration in range(self.max_iter):
# Eステップ: 負担率の計算
gamma = self._e_step(X)
# Mステップ: パラメータの更新
self._m_step(X, gamma)
# 対数尤度の計算
ll = self._log_likelihood(X)
self.log_likelihoods.append(ll)
# 収束判定
if iteration > 0 and abs(ll - self.log_likelihoods[-2]) < self.tol:
print(f"収束しました (iteration={iteration+1})")
break
return self
def _e_step(self, X):
N = X.shape[0]
gamma = np.zeros((N, self.K))
for k in range(self.K):
gamma[:, k] = self.weights[k] * multivariate_normal.pdf(
X, mean=self.means[k], cov=self.covs[k])
# 正規化
gamma_sum = gamma.sum(axis=1, keepdims=True)
gamma = gamma / gamma_sum
return gamma
def _m_step(self, X, gamma):
N, d = X.shape
for k in range(self.K):
N_k = gamma[:, k].sum()
# 平均の更新
self.means[k] = (gamma[:, k] @ X) / N_k
# 共分散行列の更新
diff = X - self.means[k]
self.covs[k] = (diff.T @ np.diag(gamma[:, k]) @ diff) / N_k
# 混合比率の更新
self.weights[k] = N_k / N
def _log_likelihood(self, X):
N = X.shape[0]
ll = 0
for n in range(N):
s = 0
for k in range(self.K):
s += self.weights[k] * multivariate_normal.pdf(
X[n], mean=self.means[k], cov=self.covs[k])
ll += np.log(s)
return ll
def predict(self, X):
gamma = self._e_step(X)
return np.argmax(gamma, axis=1)
# --- 実行 ---
gmm = GaussianMixtureEM(n_components=3, max_iter=100)
gmm.fit(X)
pred_labels = gmm.predict(X)
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 真のラベル
ax0 = axes[0]
for k in range(3):
mask = true_labels == k
ax0.scatter(X[mask, 0], X[mask, 1], alpha=0.5, s=15, label=f'Cluster {k}')
ax0.set_title('True Labels')
ax0.legend()
ax0.grid(True, alpha=0.3)
# 推定ラベル
ax1 = axes[1]
for k in range(3):
mask = pred_labels == k
ax1.scatter(X[mask, 0], X[mask, 1], alpha=0.5, s=15, label=f'Cluster {k}')
# 推定された平均を表示
for k in range(3):
ax1.scatter(*gmm.means[k], c='black', marker='x', s=200, linewidths=3)
ax1.set_title('Estimated Labels (EM)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 対数尤度の推移
ax2 = axes[2]
ax2.plot(gmm.log_likelihoods, marker='o', markersize=3)
ax2.set_xlabel('Iteration')
ax2.set_ylabel('Log-Likelihood')
ax2.set_title('Log-Likelihood Convergence')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"推定された混合比率: {gmm.weights}")
for k in range(3):
print(f"クラスタ {k}: 平均 = {gmm.means[k].round(2)}")
まとめ
本記事では、混合ガウスモデル(GMM)とEMアルゴリズムについて解説しました。
- 混合ガウスモデルは複数のガウス分布の重み付き和でデータを表現するモデル
- 対数尤度関数の直接最適化が困難なため、EMアルゴリズムを用いて反復的にパラメータを推定する
- Eステップで負担率（各データ点が各クラスタに属する確率）を計算し、Mステップでパラメータを更新する
- EMアルゴリズムは対数尤度の単調増加を保証するが、初期値に依存して局所最適解に収束する場合がある
次のステップとして、以下の記事も参考にしてください。