混合ガウスモデルと EM アルゴリズムのしくみ：見えないグループを確率で推定する

2 つの品種のリンゴが、品種の記録なしで 1 つの箱に混ざっているとします。重さを量ると、山が 2 つある分布になります。ここから「それぞれの品種の平均の重さとばらつき、混ざっている割合」を推定するにはどうすればよいでしょうか。各リンゴがどちらの品種か分かれば平均を取るだけですが、それが分からないのが問題です。

このような「見えないグループ（潜在変数）があるデータ」を扱う代表的な方法が、混合ガウスモデル と、その推定に使う EM アルゴリズム です。クラスタリングと次元削減のしくみ で説明した k-means を、確率の考え方で一般化したものと見ることもできます。

結論を先に書くと、次のようになります。

• 混合ガウスモデルは、データを「いくつかの正規分布の重ね合わせ」で表す。各分布の割合・平均・ばらつき（共分散）を推定する
• EM アルゴリズムは、2 つの手順をくり返す。E ステップで各データが各グループに属する確率（負担率）を計算し、M ステップでその確率を重みにして、割合・平均・ばらつきを推定し直す。くり返すたびに、尤度は必ず上がる（下がらない）
• k-means との違いは、所属を確率で表すことと、まとまりの形（向き・細長さ）まで学習すること。細長いまとまりのデータでは、k-means の正解率 0.593 に対して、混合ガウスモデルは 0.992 だった
• グループの数は BIC などで選ぶ。また、1 つのグループが 1 点に縮むと尤度がいくらでも大きくなるので、ばらつきに下限を設ける

説明の途中に出てくる計算例は、Python と NumPy で実際に動かした結果です。

環境

計算例を動かした環境です。

項目	バージョン
Python	3.13
NumPy	2.5

混合ガウスモデル

正規分布の重ね合わせ

混合ガウスモデルは、データの分布を次の形で表します。

p(x) = π₁ × N(x | μ₁, σ₁²) + π₂ × N(x | μ₂, σ₂²) + … + π_K × N(x | μ_K, σ_K²)

N(x | μ, σ²) は平均 μ・分散 σ² の正規分布、π_k は k 番目のグループの割合（合計は 1）です。データは「まず割合 π に従ってグループを選び、そのグループの正規分布から値が出てくる」というしくみで作られた、と考えます。選ばれたグループは記録されていないので、潜在変数（見えない変数）と呼びます。

推定したいもの

推定したいのは、割合 π・平均 μ・標準偏差 σ です。ふつうの正規分布なら、ベイズ機械学習のしくみ などで出てきた 最尤推定（データが出てくる確率＝尤度がいちばん高くなるパラメーターを選ぶ）で、平均と標準偏差を計算するだけで求まります。混合ガウスモデルでは、尤度の式の中に「どのグループから出たか」についての足し算が入るので、一度に解く式が作れません。

EM アルゴリズム

2 つの手順のくり返し

EM アルゴリズム（Dempster ら, 1977）は、「グループが分かれば推定は簡単」「推定が分かればグループの見当は付く」という 2 つを交互にくり返します。

手順	名前	やること
E ステップ	期待値（Expectation）	今のパラメーターで、各データが各グループから出てきた確率（負担率）を計算する
M ステップ	最大化（Maximization）	負担率を重みにして、割合・平均・標準偏差を計算し直す

E ステップの負担率は、ベイズの定理で計算します。

r_nk = π_k × N(x_n | μ_k, σ_k²) / Σ_j π_j × N(x_n | μ_j, σ_j²)

M ステップは、負担率を「そのデータがそのグループに属する割合」と考えて、重み付きの平均を取るだけです。

N_k = Σ_n r_nk
π_k = N_k / N
μ_k = Σ_n r_nk × x_n / N_k
σ_k² = Σ_n r_nk × (x_n − μ_k)² / N_k

リンゴの重さで計算する

平均 250 g・標準偏差 20 g の品種を 300 個、平均 320 g・標準偏差 25 g の品種を 200 個混ぜたデータを作り、品種を知らないものとして推定します。初期値は、わざと大きく外して、平均 200 g と 400 g・標準偏差 50 g から始めます。

import numpy as np

def normal_pdf(x, mu, sd):
    return np.exp(-0.5 * ((x - mu) / sd) ** 2) / (sd * np.sqrt(2 * np.pi))

# 2 つの品種が混ざったリンゴの重さ（g）：品種は記録されていない
rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(250, 20, 300), rng.normal(320, 25, 200)])

pi, mu, sd = np.array([0.5, 0.5]), np.array([200.0, 400.0]), np.array([50.0, 50.0])  # わざと外れた初期値
for it in range(51):
    dens = pi * normal_pdf(x[:, None], mu, sd)              # 各品種から出てくる確率の高さ × 割合
    loglik = np.log(dens.sum(1)).sum()
    if it in (0, 1, 2, 5, 10, 50):
        print(f"{it:2d} 回目: 割合 {np.round(pi, 3)}, 平均 {np.round(mu, 1)}, 標準偏差 {np.round(sd, 1)}, 対数尤度 {loglik:.1f}")
    r = dens / dens.sum(1, keepdims=True)                    # E ステップ：どちらの品種らしいか（負担率）
    nk = r.sum(0)                                            # M ステップ：負担率で重み付けして推定し直す
    pi = nk / len(x)
    mu = (r * x[:, None]).sum(0) / nk
    sd = np.sqrt((r * (x[:, None] - mu) ** 2).sum(0) / nk)

for w in (240, 280, 290, 300, 340):
    d = pi * normal_pdf(w, mu, sd)
    print(f"{w} g のリンゴが品種 1 である確率 {d[0] / d.sum():.3f}")

実行結果:

 0 回目: 割合 [0.5 0.5], 平均 [200. 400.], 標準偏差 [50. 50.], 対数尤度 -3106.2
 1 回目: 割合 [0.678 0.322], 平均 [255.6 323.5], 標準偏差 [26.8 24.9], 対数尤度 -2529.1
 2 回目: 割合 [0.668 0.332], 平均 [254.1 324.5], 標準偏差 [25.  22.8], 対数尤度 -2523.3
 5 回目: 割合 [0.646 0.354], 平均 [251.8 324.2], 標準偏差 [22.6 21.5], 対数尤度 -2519.4
10 回目: 割合 [0.623 0.377], 平均 [250.3 322.3], 標準偏差 [21.4 22.4], 対数尤度 -2518.1
50 回目: 割合 [0.587 0.413], 平均 [248.3 318.9], 標準偏差 [19.9 24.5], 対数尤度 -2517.4
240 g のリンゴが品種 1 である確率 0.997
280 g のリンゴが品種 1 である確率 0.634
290 g のリンゴが品種 1 である確率 0.281
300 g のリンゴが品種 1 である確率 0.075
340 g のリンゴが品種 1 である確率 0.000

• 1 回のくり返しで、平均は 255.6 g と 323.5 g まで近づいた。初期値の 200 g と 400 g からは大きく動いています
• 50 回で、割合 0.587 / 0.413・平均 248.3 / 318.9 g・標準偏差 19.9 / 24.5 g になった。データを作ったときの値（0.6 / 0.4、250 / 320 g、20 / 25 g）にかなり近く、品種の記録がなくても推定できています
• 対数尤度は、くり返すたびに上がっている（−3106.2 → −2517.4）。これは EM アルゴリズムの一般的な性質で、どの回でも尤度が下がることはありません
• 品種の判定は確率で出る。240 g ならほぼ確実に品種 1、280 g なら 0.634 で、どちらとも言い切れない、という答え方ができます

なぜ尤度が下がらないのか

EM アルゴリズムは、本当の対数尤度の「下から接する近似の式」を作り（E ステップ）、その近似の式を最大にする（M ステップ）、という操作になっています。近似の式は今のパラメーターでは本当の対数尤度とぴったり一致し、ほかの場所では下回るので、近似の式を上げれば、本当の対数尤度も同じだけ以上に上がります。

ただし、保証されるのは「下がらない」ことだけで、いちばん良い答えにたどり着く保証はありません。初期値によっては、途中の山（局所的な最大）で止まります。実際には、初期値を変えて何回か実行し、尤度のいちばん高い結果を選びます。

k-means との違い

k-means は「硬い」EM

クラスタリングと次元削減のしくみ で説明した k-means は、「各点をいちばん近い中心に割り当てる」「割り当てた点の平均を中心にする」をくり返しました。これは、混合ガウスモデルの EM で、次のように決めたものと同じです。

	k-means	混合ガウスモデル
所属	どれか 1 つ（0 か 1）	確率（負担率）
まとまりの形	どれも同じ大きさの丸	向き・細長さ・大きさを、グループごとに学習する（共分散行列）
まとまりの大きさの違い	扱えない	割合 π で扱える
計算	速い	遅い。次元が多いと共分散の推定が難しい

細長いまとまりで比べる

向きの違う細長いまとまりが 3 つあるデータで、k-means と混合ガウスモデルを比べます。2 次元では、標準偏差のかわりに 共分散行列 を推定します。共分散行列は、まとまりが「どちらの向きに、どれだけ伸びているか」を表します。

from itertools import permutations

def mvn_logpdf(X, mu, S):
    d = X.shape[1]
    L = np.linalg.cholesky(S)
    z = np.linalg.solve(L, (X - mu).T)
    return -0.5 * (z**2).sum(0) - np.log(np.diag(L)).sum() - 0.5 * d * np.log(2 * np.pi)

def gmm_em(X, K, rng, iters=200, reg=1e-6):
    n, d = X.shape
    mu = X[rng.choice(n, K, replace=False)]
    S = np.array([np.cov(X.T) for _ in range(K)]); pi = np.full(K, 1 / K)
    for _ in range(iters):
        logp = np.column_stack([np.log(pi[k]) + mvn_logpdf(X, mu[k], S[k]) for k in range(K)])
        m = logp.max(1, keepdims=True)
        loglik = (m[:, 0] + np.log(np.exp(logp - m).sum(1))).sum()
        r = np.exp(logp - m); r /= r.sum(1, keepdims=True)   # E ステップ
        nk = r.sum(0)                                         # M ステップ
        pi = nk / n
        mu = (r.T @ X) / nk[:, None]
        S = np.array([((r[:, k, None] * (X - mu[k])).T @ (X - mu[k])) / nk[k] + reg * np.eye(d) for k in range(K)])
    return pi, mu, S, r, loglik

def kmeans(X, K, rng, iters=100):
    c = X[rng.choice(len(X), K, replace=False)]
    for _ in range(iters):
        lab = ((X[:, None] - c[None]) ** 2).sum(-1).argmin(1)
        c = np.array([X[lab == k].mean(0) for k in range(K)])
    return lab

def match_acc(y, lab, K):
    return max(np.mean(np.array(p)[lab] == y) for p in permutations(range(K)))

# 細長い 3 つのまとまり（向きがそれぞれ違う）
rng = np.random.default_rng(0)
def blob(n, center, angle, long=3.0, short=0.3):
    R = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
    return rng.normal(size=(n, 2)) * [long, short] @ R.T + center
X = np.vstack([blob(300, [0, 0], 0.3), blob(300, [0, 2.5], 0.3), blob(300, [5, -3], 1.4)])
y = np.repeat([0, 1, 2], 300)

best = max((gmm_em(X, 3, np.random.default_rng(s)) for s in range(5)), key=lambda t: t[-1])   # 初期値を変えて 5 回試し、尤度の高いものを選ぶ
print("k-means の正解率", round(float(match_acc(y, kmeans(X, 3, np.random.default_rng(0)), 3)), 3))
print("混合ガウスの正解率", round(float(match_acc(y, best[3].argmax(1), 3)), 3))
print("混合ガウスの割合", np.round(best[0], 3))

対数の確率の足し算は、そのまま指数を取ると桁があふれるので、いちばん大きい値 m を引いてから計算しています（log-sum-exp の工夫）。reg は、共分散行列が 0 に近づいて計算が壊れないように足す、小さな値です。クラスタリングの結果のグループの番号は、正解の番号と対応していないので、match_acc で番号の対応をすべて試し、いちばん合うもので正解率を計算しています。

実行結果:

k-means の正解率 0.593
混合ガウスの正解率 0.992
混合ガウスの割合 [0.335 0.333 0.332]

k-means は、どのまとまりも「丸い」と考えるので、細長く平行に並んだ 2 つのまとまりを、横に切り分けてしまいます。混合ガウスモデルは、それぞれのまとまりの向きと細長さを共分散行列として学習するので、ほぼ正しく分けられました。割合も、正解の 1/3 ずつにほぼ一致しています。

グループの数の選び方

尤度だけでは決められない

グループの数 K を増やせば、データにより細かく合わせられるので、尤度はほぼ必ず上がります。そこで、パラメーターの数が多いことへの罰則 を足した指標で比べます。よく使われるのが BIC（ベイズ情報量規準, Schwarz, 1978）です。

BIC = −2 × 対数尤度 + パラメーターの数 × log(データの数)

BIC が小さいほど良いモデルです。2 次元の混合ガウスモデルのパラメーターの数は、グループごとに平均 2 個・共分散行列 3 個、それに割合 K − 1 個です。

# まとまりの数を BIC で選ぶ
n, d = X.shape
for K in range(1, 7):
    loglik = max(gmm_em(X, K, np.random.default_rng(s))[-1] for s in range(5))
    n_params = K * (d + d * (d + 1) / 2) + (K - 1)            # 平均・共分散・割合の数
    print(K, "対数尤度", round(float(loglik), 1), "BIC", round(float(-2 * loglik + n_params * np.log(n)), 1))

実行結果:

1 対数尤度 -4579.7 BIC 9193.3
2 対数尤度 -3856.5 BIC 7787.8
3 対数尤度 -3419.7 BIC 6955.0
4 対数尤度 -3413.0 BIC 6982.4
5 対数尤度 -3402.1 BIC 7001.5
6 対数尤度 -3396.1 BIC 7030.4

対数尤度は K を増やすほど上がり続けますが、3 を超えると上がり方がわずかになります。BIC は K = 3 でいちばん小さくなり、正しいまとまりの数を選べました。

落とし穴：尤度がいくらでも大きくなる

混合ガウスモデルの尤度には、困った性質があります。1 つのグループの平均を 1 個のデータの上 に置き、標準偏差を 0 に近づけると、そのデータの確率の高さが無限に大きくなり、尤度もいくらでも大きくなってしまうのです。

リンゴのデータに、1 個のリンゴの上に置いた 3 つ目の山を足して、その幅を細くしていきます。

# 3 つ目の山を、1 個のデータの上に置いて、幅を細くしていく
base_pi, base_mu, base_sd = np.array([0.6, 0.4]), np.array([250.0, 320.0]), np.array([20.0, 25.0])
for s in (10.0, 1.0, 1e-3, 1e-6):
    pi = np.append(base_pi * 0.99, 0.01); mu = np.append(base_mu, x[0]); sd = np.append(base_sd, s)
    print(f"3 つ目の山の標準偏差 {s:g}: 対数尤度 {np.log((pi * normal_pdf(x[:, None], mu, sd)).sum(1)).sum():.1f}")

実行結果:

3 つ目の山の標準偏差 10: 対数尤度 -2518.7
3 つ目の山の標準偏差 1: 対数尤度 -2519.3
3 つ目の山の標準偏差 0.001: 対数尤度 -2517.6
3 つ目の山の標準偏差 1e-06: 対数尤度 -2510.7

標準偏差を 1000 分の 1 にするたびに、対数尤度は約 log(1000) ≒ 6.9 ずつ上がり、0 に近づけるほど限りなく大きくなります。もちろん、これは「1 個のリンゴだけの品種がある」という意味のない答えです。EM アルゴリズムの途中でこの状態に落ち込むことがあるので、実際には次のような対策をします。

• 分散（共分散行列の対角）に小さな値を足して、下限を設ける（上のコードの reg、scikit-learn の reg_covar）
• 負担率の合計 N_k が小さくなりすぎたグループは、初期化し直す
• パラメーターに事前分布を置くベイズ的な推定（変分ベイズの混合ガウスモデル）にする

EM アルゴリズムの広がり

EM アルゴリズムは、混合ガウスモデル専用の方法ではなく、見えない変数があるモデル全般 の推定に使われます。

モデル	見えない変数	使われる場面
混合ガウスモデル	どのグループから出たか	クラスタリング、密度の推定、異常検知（異常検知のしくみ）
隠れマルコフモデル（バウム・ウェルチ法）	各時刻の見えない状態	昔の音声認識（音声認識と音声合成のしくみ）、品詞の付与、遺伝子の配列の解析
トピックモデル	各単語がどの話題から出たか	文書の分類、話題の抽出
欠損値の推定	抜けている値	欠損のある表のデータから、平均や相関を推定する
因子分析	観測できない要因	心理学のアンケートの分析

深層学習の VAE（画像生成 AI のしくみ）も、潜在変数の分布を近似しながら尤度の下限を最大にする、という点で、EM アルゴリズムと同じ考え方の上にあります。

混合ガウスモデルと EM の流れ

年	手法・出来事	ポイント
1894	正規分布の混合の推定（ピアソン）	カニの甲羅の大きさのデータを、2 つの正規分布の混合として分析した
1970	バウム・ウェルチ法	隠れマルコフモデルのパラメーターを推定する。EM アルゴリズムの特別な場合
1977	EM アルゴリズム（Dempster・Laird・Rubin）	見えない変数があるモデルの最尤推定を、一般的な形にまとめた
1978	BIC（Schwarz）	パラメーターの数の罰則を付けて、モデルの複雑さを選ぶ
1983	EM の収束の性質（Wu）	尤度が下がらないこと、局所的な最大に収束しうることを示した
1998	EM の見直し（Neal と Hinton）	E ステップも M ステップも、1 つの関数を交互に最大にする操作として整理した。変分推論につながる
2003	潜在ディリクレ配分法（LDA）	文書の話題を見えない変数とするトピックモデル
2013	VAE	潜在変数のモデルを、ニューラルネットワークで近似して学習する

まとめ

• 混合ガウスモデルは、データを割合 π・平均 μ・ばらつきの違う正規分布の重ね合わせで表す。どのグループから出たかは見えない変数として扱う
• EM アルゴリズムは、負担率を計算する E ステップと、負担率で重み付けして推定し直す M ステップをくり返す。尤度は下がらないが、局所的な最大で止まることがあるので、初期値を変えて何回か試す
• k-means は、所属を 0 か 1 にし、まとまりを丸と決めた、混合ガウスモデルの特別な場合。細長いまとまりでは、混合ガウスモデルの方がずっと正しく分けられた
• グループの数は、対数尤度にパラメーターの数の罰則を足した BIC などで選ぶ
• 1 つのグループが 1 点に縮むと尤度がいくらでも大きくなるので、ばらつきに下限を設ける

参考

• K. Pearson, "Contributions to the Mathematical Theory of Evolution", Philosophical Transactions of the Royal Society of London A 1894
• L. E. Baum et al., "A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains", The Annals of Mathematical Statistics 1970
• A. P. Dempster, N. M. Laird and D. B. Rubin, "Maximum Likelihood from Incomplete Data via the EM Algorithm", Journal of the Royal Statistical Society Series B 1977
• G. Schwarz, "Estimating the Dimension of a Model", The Annals of Statistics 1978
• C. F. J. Wu, "On the Convergence Properties of the EM Algorithm", The Annals of Statistics 1983
• R. M. Neal and G. E. Hinton, "A View of the EM Algorithm that Justifies Incremental, Sparse, and other Variants", Learning in Graphical Models 1998
• D. M. Blei, A. Y. Ng and M. I. Jordan, "Latent Dirichlet Allocation", Journal of Machine Learning Research 2003
• C. M. Bishop, "Pattern Recognition and Machine Learning", Springer 2006（9 章）
• D. P. Kingma and M. Welling, "Auto-Encoding Variational Bayes", ICLR 2014（arXiv:1312.6114）
• scikit-learn User Guide: Gaussian mixture models