🔖 キーワード索引

このページで扱う主要キーワード（クリックで該当セクションへ）：

🎨 直感で掴む

EM (Expectation-Maximization) アルゴリズムは、 「潜在変数」や「欠損データ」を含む確率モデルの最尤推定を、 2 ステップ反復で行う汎用アルゴリズムです。 1977 年、 Dempster, Laird, Rubin により定式化されました。

直感的には：

• 「もし潜在変数が分かれば、 パラメータを推定するのは簡単」
• 「もしパラメータが分かれば、 潜在変数の分布を計算するのは簡単」
• この 2 つを交互に繰り返す

代表的な応用：

1. ガウス混合モデル（GMM）— データを $K$ 個の正規分布の重ね合わせとして説明
2. HMM（隠れマルコフモデル）— 系列データの潜在状態推定
3. Latent Dirichlet Allocation — 文書のトピック
4. 因子分析・確率的 PCA — 低次元潜在因子
5. 欠損補完 — 欠損値を潜在変数として扱う

🎨 概念図で押さえる（EM アルゴリズムの可視化）

EM アルゴリズムは「潜在変数を含む最尤推定」を E ステップと M ステップの交互更新で行う。 混合分布・クラスタリング・収束過程の 3 枚で、 EM の挙動と GMM への応用を視覚化する。

→ クラスタリング比較（GMM = EM 応用）、 k-means（特殊ケース）、 混合分布（柔軟な確率モデル）。 EM の汎用性と GMM への帰着を理解できる。

✅ 理解度チェック

1. EM アルゴリズムの E ステップと M ステップでそれぞれ何を計算するか？
2. EM が単調に対数尤度を増加させることが保証される数学的理由は？
3. EM が局所最適に捕まる典型的な状況は？ 対策は？
4. k-means と GMM (EM-based) の最大の違いは？
5. scikit-learn で EM を使う代表的なクラス名は？

→ すべて即答できれば、 EM を欠損データ・混合モデルで実践的に使える基礎力は十分。

🧾 発表前の最終確認

EMアルゴリズムでは、 初期値、収束判定、対数尤度の推移、局所解への依存を説明します。 複数初期値で同じ解に近づくかを確認し、 収束した値が大域最適とは限らない点も明示します。

🎮 触って理解する

下のツールは 1 次元・2 成分の混合ガウス分布に対する EM アルゴリズムを、 スライダーとボタンで手動実行できるインタラクティブ教材です。 初期パラメータ（平均 $\mu$・標準偏差 $\sigma$・混合比 $\pi$）をスライダーで決め、 E ステップ（各点の責任度 $\gamma$ を計算＝点の色が変わる）と M ステップ（責任度で重み付けした平均・分散・混合比に更新＝曲線が動く）を交互に押すと、 対数尤度が単調に増えながら収束していく様子が体感できます。 ガウス密度・責任度・重み付き平均分散はすべて数式どおり正確に計算しています。

※ 使用データは 2 つのガウス分布から決定論的に生成した合成デモデータ（60 点）で、 実統計値ではありません。 実データでの計算例は下の「🧮 実値で計算してみる」節（SSDSE-B-2026 の実測値）を参照してください。

🧭 直感 — 何が起きているのか

E ステップでは現在の 2 本の釣鐘曲線を「定規」として、 各点が成分1と成分2のどちらから来たかの確率（責任度 $\gamma$）を計算します。 画面では点の色が青〜赤のグラデーションで表示され、 中間色の点が「どちらとも言えない境界の点」です。 M ステップでは、 その責任度で重み付けした平均・分散・混合比に曲線を更新します。 責任度が高い点ほどその成分の平均を強く引っ張る、 という「加重平均」がすべての正体です。 E と M を繰り返すたびに対数尤度（データへの当てはまりの良さ）が必ず増えるか横這いになり、 やがて曲線が動かなくなれば収束です。

⚠️ よくある落とし穴（触って確かめる）

🚀 発展

ここでは 1 次元・2 成分でしたが、 同じ E→M の反復は多次元・多成分（$K$ 個）、 共分散行列つきの一般的な GMM にそのまま一般化できます（下の Python 実装・SSDSE-B-2026 の 47 都道府県クラスタを参照）。 さらに責任度の計算を「前向き・後ろ向きアルゴリズム」に置き換えれば HMM の Baum-Welch、 E ステップを変分近似に置き換えれば 変分 EM / VAE になります。 「E で潜在変数の分布、 M でパラメータ」という骨格は共通です。

関連ページ： 混合ガウスモデル (GMM) k-means（ハード版） 正規分布 平均 分散

📐 数式・定義

観測 $X$、 潜在 $Z$、 パラメータ $\theta$ とする。 観測尤度 $p(X|\theta) = \sum_Z p(X,Z|\theta)$ を直接最大化するのは難しいので、 補助関数 $Q$ を作って交互に最適化する。

これを単調に繰り返すと、 観測対数尤度 $\log p(X|\theta)$ は単調非減少。 Jensen の不等式を使うと容易に証明できる：

📐 数式: ELBO 分解と EM の単調増加性

EM がなぜ単調増加するか、 ELBO 分解で美しく示せる。 これは EM を理解する上で**最重要の数式**。

🔬 数式を言葉で読み解く — ELBO 分解

任意の分布 $q(Z)$ に対して、 観測対数尤度は次のように分解できる: $$\log p(X|\theta) = \mathcal{L}(q, \theta) + \mathrm{KL}(q(Z) \| p(Z|X,\theta))$$ $$\mathcal{L}(q, \theta) = \sum_Z q(Z) \log \frac{p(X,Z|\theta)}{q(Z)}$$

• $\mathcal{L}(q, \theta)$ — Evidence Lower Bound (ELBO)、 観測尤度の下界
• $\mathrm{KL}(q \| p)$ — 提案分布 $q$ と事後分布 $p(Z|X,\theta)$ の KL ダイバージェンス（$\geq 0$）
• $q(Z) = p(Z|X,\theta^{old})$ と置くと KL=0 となり、 ELBO=対数尤度

E ステップは「$q$ について ELBO を最大化」、 M ステップは「$\theta$ について ELBO を最大化」と統一的に書ける。 ELBO は常に対数尤度の下界なので、 ELBO を増やせば対数尤度も必ず増加（または不変）。 これが **EM の単調収束保証**の正体。

🔬 数式を言葉で読み解く — GMM の E ステップ・M ステップ

GMM では具体的に次の式になる: $$\text{E:} \quad \gamma_{ik} = \frac{\pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)}{\sum_{j=1}^K \pi_j \mathcal{N}(x_i|\mu_j, \Sigma_j)}$$ $$\text{M:} \quad \mu_k = \frac{\sum_i \gamma_{ik} x_i}{N_k}, \quad \Sigma_k = \frac{\sum_i \gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^\top}{N_k}, \quad \pi_k = \frac{N_k}{N}$$ ここで $N_k = \sum_i \gamma_{ik}$ は実効的サンプル数。

E ステップは「観測 $x_i$ がクラスタ $k$ から生成された確率（責任度 $\gamma_{ik}$）」をベイズで計算。 M ステップは「責任度で重み付けした平均・共分散・混合比」を計算。 重み付き標本統計量と完全に同じ形になる。

🔬 数式を言葉で読み解く — Jensen の不等式と E ステップの根拠

$\log$ は凹関数なので Jensen の不等式から: $$\log p(X|\theta) = \log \sum_Z p(X, Z|\theta) = \log \sum_Z q(Z) \frac{p(X,Z|\theta)}{q(Z)} \geq \sum_Z q(Z) \log \frac{p(X,Z|\theta)}{q(Z)}$$ 右辺が ELBO そのもの。 等号成立は $\frac{p(X,Z|\theta)}{q(Z)}$ が $Z$ について定数のとき、 つまり $q(Z) \propto p(X,Z|\theta) = p(X|\theta) p(Z|X,\theta)$、 すなわち $q(Z) = p(Z|X,\theta)$。 これが「E ステップで事後分布を選べ」の数学的根拠。

🔬 数式を言葉で読み解く

GMM (Gaussian Mixture Model) における具体形：

E ステップで $\gamma_{ik}$ を計算し、 M ステップで $(\pi_k, \mu_k, \Sigma_k)$ を更新、 この 2 つを収束まで繰り返します。

🧮 実値で計算してみる（SSDSE-B-2026）

🧮 数式に値を入れて手で計算する: SSDSE-B-2026 A4103 の 2 ガウス混合での E-step

SSDSE-B-2026 の A4103 (合計特殊出生率) は 2 つの中心 (都市圏 ≈ 1.1, 沖縄等 ≈ 1.5) を持つ可能性がある。 2 ガウス混合 N1(μ=1.10, σ=0.07) と N2(μ=1.50, σ=0.07) を仮定し、 各都道府県の値が「都市タイプ」か「地方タイプ」かの責任度 γ を計算する。 実データでは東京都 A4103=0.99（全国最小）、 沖縄県 A4103=1.60（全国最大）、 全国中央値は 1.30 である。

Step 1: 初期パラメータと観測

Step 2: 東京都 x=0.99 の責任度 γ

Step 3: 沖縄 x=1.60 の責任度

🐍 Python で再現

このコードでやること: SSDSE-B-2026 から東京 (0.99) と沖縄 (1.60) の A4103 を取り、 2 ガウス混合の責任度 γ を scipy.stats.norm で計算する。

📥 入力データ:

1
2
3
4
5
6
7
8

from scipy.stats import norm
def gamma(x, mu1=1.10, mu2=1.50, sig=0.07):
    f1 = norm.pdf(x, mu1, sig)
    f2 = norm.pdf(x, mu2, sig)
    return round(f1/(f1+f2), 4), round(f2/(f1+f2), 4)
print(f"東京 A4103=0.99: γ={gamma(0.99)}")
print(f"沖縄 A4103=1.60: γ={gamma(1.60)}")
print(f"境界例 A4103=1.30: γ={gamma(1.30)}")

📤 実行結果

💬 手計算 (Step 2/3) と Python 出力が完全一致。 東京 (0.99) はほぼ完全に「都市タイプ」、 沖縄 (1.60) はほぼ完全に「地方タイプ」、 中間値 1.30 では責任度が 50:50 で分かれる。 EM アルゴリズムはこの γ を全 47 県で計算し、 M-step で μ・σ・π を更新する。

🐍 Python 実装

例 1：scikit-learn の GaussianMixture（SSDSE-B-2026 で 47 県クラスタリング）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

import pandas as pd, numpy as np
from sklearn.mixture import GaussianMixture

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
df['高齢化率'] = df['A1303'] / df['A1101']   # 65歳以上人口 / 総人口

X = df[['A1101','高齢化率']].to_numpy().astype(float)
X = (X - X.mean(axis=0)) / X.std(axis=0)
gmm = GaussianMixture(n_components=3, covariance_type='full',
                       n_init=10, max_iter=200, random_state=0,
                       init_params='kmeans')
gmm.fit(X)
df['cluster'] = gmm.predict(X)
proba = gmm.predict_proba(X)
print(df[['Prefecture','cluster']].head(10))
print('mixing weights:', gmm.weights_.round(3))
print('対数尤度:', round(gmm.score(X)*len(X), 1), ' BIC:', round(gmm.bic(X), 1))

例 2：BIC でクラスタ数を選ぶ

1
2
3
4
5
6
7
8

bics = []
for k in range(1, 8):
    g = GaussianMixture(n_components=k, n_init=10, random_state=0).fit(X)
    bics.append((k, g.bic(X)))
for k, b in bics:
    print(f'K={k}  BIC={b:.1f}')
best_k = min(bics, key=lambda t: t[1])[0]
print('最良 K:', best_k)

例 3：手書き EM（1 次元混合ガウス、 2 成分）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import numpy as np
from scipy.stats import norm

# 例 1 の df（cp932, 年度=2023）をそのまま利用。総人口 A1101 の対数値を 2 成分混合で当てはめ
x = np.log(df['A1101'].astype(float).to_numpy())
pi = np.array([0.5, 0.5])
mu = np.array([x.min(), x.max()])
sigma = np.array([x.std(), x.std()])

for it in range(100):
    # E
    r = np.vstack([pi[k]*norm.pdf(x, mu[k], sigma[k]) for k in (0,1)])
    r = r / r.sum(axis=0)
    # M
    nk = r.sum(axis=1)
    pi = nk / nk.sum()
    mu = (r * x).sum(axis=1) / nk
    sigma = np.sqrt(((r * (x - mu[:,None])**2).sum(axis=1)) / nk)
print('pi:', pi.round(3), ' mu:', mu.round(3), ' sigma:', sigma.round(3))
print('exp(mu) [人]:', np.exp(mu).round().astype(int))

例 4：欠損補完への応用（EM-imputation）

1
2
3
4
5
6
7

from sklearn.experimental import enable_iterative_imputer  # noqa
from sklearn.impute import IterativeImputer
# IterativeImputer は内部で「他列を使った回帰を反復」する EM 風アルゴリズム
imp = IterativeImputer(max_iter=20, random_state=0)
num = df.select_dtypes(include='number')          # A1101 など数値列
filled = imp.fit_transform(num)
print('NaN 数（補完前 vs 補完後）:', num.isna().sum().sum(), '→', int(np.isnan(filled).sum()))

📂 ケーススタディ・追加実装例

ケース 1：HMM の Baum-Welch（EM の特殊形）

隠れマルコフモデルの推定で、 E ステップは「前向き・後ろ向きアルゴリズム」、 M ステップは「遷移行列・出力行列の重み付き正規化」。

ケース 2：LDA（潜在ディリクレ配分）

文書 $\to$ トピック $\to$ 単語、 という二段の生成モデル。 変分 EM や Gibbs sampler で推定。 ニュース記事を 10 トピックに分けるなどが典型例。

ケース 3：欠損補完への EM 応用

多変量正規分布を仮定すると、 「欠損値の条件付き期待値」と「パラメータの再推定」を交互に行う EM が組める（Little & Rubin の古典）。

ケース 4：MNAR / MAR / MCAR の前提

ケース 5：GMM の対数尤度の振る舞い

1
2
3
4
5
6
7
8
9

from sklearn.mixture import GaussianMixture
import numpy as np
losses = []
for seed in range(10):
    g = GaussianMixture(n_components=3, n_init=1, random_state=seed,
                         max_iter=500, tol=1e-6).fit(X)
    losses.append(g.score(X) * len(X))
print('対数尤度の最小・最大・中央:', round(np.min(losses),1), round(np.max(losses),1), round(np.median(losses),1))
# 必ず複数 seed で確認

ケース 6：BIC によるクラスタ数選定

BIC = $-2 \ln L + k \ln n$。 $L$＝最大尤度、 $k$＝パラメータ数、 $n$＝データ数。 BIC が最小の $K$ を選ぶ。 AIC は罰則が弱く、 BIC は強い。

🪜 ステップバイステップ チュートリアル

チュートリアル：SSDSE で GMM クラスタリング

ステップ 1：問題設定

47 都道府県を「総人口」「高齢化率」「消費支出（二人以上の世帯）」の 3 特徴量で混合ガウス K=3 にクラスタリング。

ステップ 2：データ準備

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import pandas as pd, numpy as np
from sklearn.mixture import GaussianMixture

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
df['高齢化率'] = df['A1303'] / df['A1101']        # 65歳以上人口 / 総人口

# A1101=総人口, L3221=消費支出（二人以上の世帯）
X = df[['A1101','高齢化率','L3221']].to_numpy().astype(float)
X = (X - X.mean(axis=0)) / X.std(axis=0)

ステップ 3：BIC で K 選定

1
2
3
4
5
6
7

bics = []
for k in range(1, 8):
    g = GaussianMixture(n_components=k, n_init=10, random_state=0).fit(X)
    bics.append(g.bic(X))
best_k = 1 + int(np.argmin(bics))
print('bics:', [round(b,1) for b in bics])
print('BIC 最小 K =', best_k)

ステップ 4：EM 実行と事後確率

1
2
3
4
5
6
7

gmm = GaussianMixture(n_components=3, n_init=10, random_state=0).fit(X)  # 解釈性優先で K=3
df['cluster'] = gmm.predict(X)
proba = gmm.predict_proba(X)
print(df[['Prefecture','cluster']].sort_values('cluster').to_string())
print('混合比 π:', gmm.weights_.round(3))
print('平均 μ (標準化後):')
print(gmm.means_.round(2))

ステップ 5：解釈

クラスタごとに代表都道府県と特徴量平均を出し、 「大都市圏 / 地方中核 / 高齢過疎」 のように命名。 境界の県（沖縄・宮城など）は事後確率が分散することを確認。

ステップ 6：手書き EM との比較

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

from scipy.stats import multivariate_normal
def manual_em(X, K, n_iter=200, tol=1e-6):
    n, d = X.shape
    # 初期化 (k-means)
    from sklearn.cluster import KMeans
    km = KMeans(n_clusters=K, n_init=10, random_state=0).fit(X)
    mu = km.cluster_centers_
    sigma = np.array([np.cov(X.T) for _ in range(K)])
    pi = np.ones(K) / K
    ll_prev = -np.inf
    for it in range(n_iter):
        # E
        r = np.zeros((n, K))
        for k in range(K):
            r[:, k] = pi[k] * multivariate_normal.pdf(X, mu[k], sigma[k])
        r /= r.sum(axis=1, keepdims=True)
        # M
        nk = r.sum(axis=0)
        pi = nk / n
        mu = (r.T @ X) / nk[:, None]
        for k in range(K):
            d_ = X - mu[k]
            sigma[k] = (r[:, k][:, None, None] * (d_[:, :, None] @ d_[:, None, :])).sum(axis=0) / nk[k]
        ll = sum(np.log(sum(pi[k]*multivariate_normal.pdf(X, mu[k], sigma[k]) for k in range(K))))
        if abs(ll - ll_prev) < tol: break
        ll_prev = ll
    return pi, mu, sigma

pi, mu, sigma = manual_em(X, 3)
print('手書き EM の π:', pi.round(3))

🚀 現場での応用シナリオ（8 例）

応用 1：顧客セグメンテーション

購買行動データを GMM でクラスタリング。 「ヘビーユーザ」「ライト」「離脱寸前」の確率的所属で施策設計。

応用 2：音声認識の隠れマルコフモデル

Baum-Welch（EM の特殊形）で音素遷移を学習。 古典的 ASR の中核。

応用 3：トピックモデル

LDA：文書から潜在トピックを抽出。 ニュース記事の分類、 推薦システム。

応用 4：欠損データ補完

多変量正規 + EM、 もしくは IterativeImputer。 アンケートやセンサーの欠損対応。

応用 5：因子分析

心理学・教育測定のテスト得点を、 少数の潜在因子で説明。 EM で因子負荷を推定。

応用 6：画像セグメンテーション

ピクセル値を GMM で K クラスに分け、 各クラスを領域とする。 古典的セグメンテーション。

応用 7：バイオインフォマティクス

SNP・遺伝子発現データのクラスタリング、 系統樹推定。

応用 8：信用スコアリング

債務不履行確率を潜在変数モデルで推定。 EM で連続的に更新。

🏋️ 演習問題（8 題）

1. SSDSE 47 県を GMM K=3 でクラスタリングし、 BIC で K を選定せよ。
2. EM の E ステップ・M ステップを手書きで実装せよ。
3. k-means と GMM を同じデータで比較し、 ソフト割当の利点を確認せよ。
4. GMM の共分散タイプを 'full', 'tied', 'diag', 'spherical' で振り、 BIC で比較せよ。
5. Hidden Markov Model を hmmlearn で訓練し、 状態系列を推定せよ。
6. LDA でニュース記事を 10 トピックに分け、 トピックごとの単語頻度を出力せよ。
7. 欠損データを IterativeImputer (MICE 風 EM) で補完せよ。
8. EM の対数尤度が反復ごとに単調非減少することを実験で確認せよ。

⚠️ よくある落とし穴

❓ よくある質問（FAQ）

📜 歴史と背景

歴史と位置づけ：EM アルゴリズムは 1977 年、 Dempster, Laird, Rubin による論文「Maximum Likelihood from Incomplete Data via the EM Algorithm」（Journal of the Royal Statistical Society, Series B）で統一的に定式化されました。 ただし、 個別の問題（混合分布・欠損データ・隠れマルコフ）に対する反復解法はそれ以前から複数の研究者が独立に発見していました（例：Baum-Welch、 1970）。

主な発展：

• 1977 — Dempster-Laird-Rubin が EM を一般化
• 1980 年代 — 単調収束性の証明、 GEM (Generalized EM) への拡張
• 1993 — Neal & Hinton による「EM を変分推論として見直す」視点
• 2000 年代 — 変分 EM、 Stochastic EM、 オンライン EM
• 現代 — 潜在変数モデル全般（VAE、 LDA など）の標準推定手法

EM は「最尤推定を 反復的・分解的 に実行する」普遍的な道具立てとして、 統計・機械学習・信号処理・バイオインフォマティクスに広く浸透しました。

⚠️ EM を実装する前に知っておくべき 7 つの落とし穴

1. 初期値依存で局所最適に捕まる — EM は対数尤度を**局所**最大化するだけで大域最適は保証しない。 sklearn の n_init=10（10 通りの初期値から始めて最良を選ぶ）、 k-means++ 初期化、 焼きなまし併用が標準対策。
2. 共分散の縮退（singular）で尤度が無限大に発散 — 1 クラスタが 1 点に縮退すると共分散がゼロ→対数尤度が +∞。 sklearn の reg_covar=1e-6 で対角に微小値を加える「数値正則化」が必須。
3. Label switching でクラスタ番号が走査の度に入れ替わる — EM のクラスタ番号は順序を持たないので、 並列実行で「クラスタ 1」が指す対象が変わる。 比較したいなら混合比でソート、 もしくはアラインメント手続き（ハンガリアン法）で番号を揃える。
4. K 選定で AIC は過大、 BIC は過小気味 — AIC は罰則 2k、 BIC は罰則 k log n。 n が大きいほど BIC は K を小さく選ぶ。 ICL（Integrated Completed Likelihood）が GMM では推奨されることも。
5. 共分散構造の誤指定で歪んだクラスタが出る — covariance_type='spherical' で実際は楕円形のデータを当てると、 球形に押し込められて変な分割になる。 まず full で試し、 BIC で他構造と比較。
6. 収束判定が緩いと尤度が低い解で止まる — sklearn の tol=1e-3 はやや緩め。 厳しい応用では tol=1e-6 + max_iter=1000 に。 ただし無限ループ気味になることがあるので双方併用。
7. MAR（Missing At Random）でないと EM は無効 — 欠損値補完で EM を使う場合、 欠損メカニズムが MAR でないと最尤推定が偏る。 完全な MNAR（Missing Not At Random）には選択モデルや pattern-mixture モデルが必要。

🔁 EM 派生・拡張アルゴリズム — ECM / GEM / MCEM / Stochastic EM

EM はモジュラーで、 E と M をそれぞれ「近似」したり「分割」することで多様な派生が生まれる。

📊 表: 主要派生アルゴリズム

🐍 ECM 風の実装 — 平均と分散を別ステップで更新

このコードでやること: 1D GMM（K=2）の EM を素朴に手書きし、 E ステップ → 平均更新 → 分散更新と分割。 SSDSE 47 県の人口を 2 つの正規分布の混合と仮定。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import numpy as np
import pandas as pd

# latest はこのあとのブロックで作っているので、ここでも用意しておく
_df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
latest = _df[_df['SSDSE-B-2026'] == 2023].copy()

from scipy.stats import norm

pop = latest['A1101'].values.astype(float)
n = len(pop)

# 初期化
mu = np.array([pop.min(), pop.max()])
sigma = np.array([pop.std(), pop.std()])
pi = np.array([0.5, 0.5])

for it in range(50):
    # E step: 責任度
    r0 = pi[0] * norm.pdf(pop, mu[0], sigma[0])
    r1 = pi[1] * norm.pdf(pop, mu[1], sigma[1])
    gamma = np.column_stack([r0, r1]) / (r0 + r1)[:, None]
    Nk = gamma.sum(axis=0)
    # M step 1: 平均
    mu = (gamma * pop[:, None]).sum(axis=0) / Nk
    # M step 2: 分散
    var = (gamma * (pop[:, None] - mu) ** 2).sum(axis=0) / Nk
    sigma = np.sqrt(var)
    # M step 3: 混合比
    pi = Nk / n

print(f"μ = {(mu/1e4).round(1)} (万人)")
print(f"σ = {(sigma/1e4).round(1)} (万人)")
print(f"π = {pi.round(3)}")

📤 実行例:

💬 47 県の人口は「平均 136 万人の小規模県 75%」と「平均 639 万人の大規模県 25%」の混合と解釈できる。 ECM 風に M を分解しても sklearn の同時最大化と同じ結果に収束する。

🛠 実装チェックリスト — EM を初めて書く前に

• ✅ 入力データに NaN/inf が無いことを確認（np.isfinite(X).all()）
• ✅ 特徴量を標準化（StandardScaler）— 共分散構造の収束安定化
• ✅ n_init=10 以上で初期値依存を緩和
• ✅ reg_covar=1e-6 で共分散縮退を防止
• ✅ K 候補を 1..10 で BIC 比較、 elbow と最小値を確認
• ✅ random_state 固定で再現性確保
• ✅ 収束後、 predict_proba で責任度を確認（0/1 にしか割り当てられていない → クラスタが分離しすぎ）

🌐 関連手法・派生・対概念

🌐 EM の応用領域 — GMM 以外で何に効くか

• HMM（Baum-Welch）: 観測系列から状態遷移行列・出力分布を推定。 音声認識・遺伝子配列。
• LDA トピックモデル: 文書を潜在トピックの混合とみなし、 単語の生成確率を EM/変分 EM で推定。
• 因子分析 (Factor Analysis): 観測変数を低次元の連続潜在因子 + 雑音で説明。 EM で因子負荷量を推定。
• 確率的 PCA: PCA を生成モデルとして書き直し、 欠損値があっても EM で訓練可能に。
• 欠損値補完 (MICE): 連鎖式で各変数を他変数の関数として補完。 E 風に補完、 M 風に回帰係数更新。
• Tobit 回帰: 打ち切りデータ（左/右）に対する最尤推定。 観測切断を潜在変数として EM 化。
• 共参照解析 (NLP): 同一エンティティを指す表現を潜在変数とみなして EM。
• 画像セグメンテーション: 画素を GMM クラスタにソフト割り当て。 医療画像で標準。
• RFM 顧客セグメンテーション: マーケティングで顧客 3D 特徴の GMM クラスタリング。

🔗 関連用語（前提・並列・発展）

📘 前提となる用語

🤝 並列・関連する用語

🚀 発展・応用

📚 関連グループ教材

このページの内容を、より広い文脈で学びたいときは：

• 📚 機械学習の基礎 — 教師なし学習と潜在変数モデル
• 📚 統計推定 — 最尤推定の応用としての EM
• 📚 クラスタリング — ハード vs ソフトクラスタの比較

📚 参考文献・出典

• Dempster, A. P., Laird, N. M., Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS B, 39(1), 1-38.
• Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Chapter 9. Springer.
• Neal, R., Hinton, G. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. Learning in Graphical Models.
• Meng, X.-L., Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm. Biometrika, 80(2), 267-278.
• Liu, C., Rubin, D. B., Wu, Y. N. (1998). Parameter expansion to accelerate EM. Biometrika, 85(4), 755-770.
• Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. The Annals of Statistics, 11(1), 95-103.
• Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective, Chapter 11. MIT Press.
• Cappé, O., Moulines, E. (2009). On-line expectation-maximization algorithm for latent data models. JRSS B, 71(3), 593-613.
• Anandkumar, A. et al. (2014). Tensor decompositions for learning latent variable models. JMLR, 15, 2773-2832.
• Baum, L. E., Welch, L. R. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. Annals of Math. Statistics, 41(1), 164-171.
• Blei, D. M., Ng, A. Y., Jordan, M. I. (2003). Latent Dirichlet Allocation. JMLR, 3, 993-1022.
• Kingma, D. P., Welling, M. (2014). Auto-Encoding Variational Bayes. ICLR.
• 政府統計の総合窓口 e-Stat. SSDSE-B-2026 都道府県別社会経済データ集. 統計センター.
• scikit-learn documentation. Gaussian mixture models.
• hmmlearn documentation. Hidden Markov Models in Python.

🧪 SSDSE-B-2026 で 47 都道府県を Gaussian Mixture でクラスタリング

EM アルゴリズムの代表的応用が **GMM（混合ガウス分布）** によるソフトクラスタリング。 ここでは SSDSE-B-2026 の社会経済指標 5 次元で 47 都道府県を K=3 クラスタに分割し、 EM の挙動と BIC による K 選定までを一気通貫で扱う。

📥 入力データ: 47 県 × 5 次元の特徴行列

SSDSE-B-2026 から最新年の 5 列を抽出: A1101 総人口, A4101 出生数, L3221 消費支出, G7101 延べ宿泊者数, I5102 一般診療所数。 5 次元なので可視化は PCA で 2D に落とす。

🐍 GMM を EM で訓練（K=3 固定）

このコードでやること: sklearn.mixture.GaussianMixture で K=3 の GMM を最尤推定（内部は EM）。 各県の責任度 γ を出して「どのクラスタに何 % 所属するか」を可視化。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

import pandas as pd
import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
latest = df[df['SSDSE-B-2026'] == 2023].copy()
cols = ['A1101', 'A4101', 'L3221', 'G7101', 'I5102']
X = latest[cols].values
Xs = StandardScaler().fit_transform(X)

gmm = GaussianMixture(n_components=3, covariance_type='full',
                      random_state=42, max_iter=200).fit(Xs)
labels = gmm.predict(Xs)
gamma = gmm.predict_proba(Xs)

print(f"対数尤度 = {gmm.score(Xs) * len(Xs):.2f}")
print(f"混合比 π = {gmm.weights_.round(3)}")
print(f"反復回数 = {gmm.n_iter_}")
print(f"\n責任度の例 (最初の 3 県):")
for i in range(3):
    print(f"  {latest['Prefecture'].iloc[i]}: γ = {gamma[i].round(3)}")

📤 実行例:

💬 結果の読み方: K=3 は「小規模 23 県」（混合比 46.8%、 平均人口 約 102 万人。 青森・岩手・秋田など）、 「中規模 15 県」（34.0%、 約 215 万人。 宮城・福島・茨城など）、 「大規模 9 県」（19.2%、 約 763 万人。 東京・神奈川・大阪・愛知に加えて面積の大きい北海道）に分かれた。 責任度 γ を見ると、 北海道と青森県はそれぞれ 1 つのクラスタに γ=1.0 でほぼ確実に所属する一方、 岩手県は γ=[0.936, 0, 0.064] と 2 つのクラスタにまたがっている。 この「どっちつかず」を 0/1 ではなく確率で表せるのが、 k-means にはない EM（混合モデル）の利点である。 収束までの反復は 12 回。 混合比が 0.47/0.34/0.19 とどれも十分な大きさで、 縮退した（ほぼ空の）クラスタが生じていないことも確認できる。

🐍 BIC で K=1..8 を比較してクラスタ数選定

このコードでやること: K を 1 から 8 まで動かして BIC を計算。 BIC 最小の K がデータが支持する最適クラスタ数。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

bic_scores, aic_scores = [], []
for k in range(1, 9):
    g = GaussianMixture(n_components=k, covariance_type='full',
                        random_state=42, max_iter=200).fit(Xs)
    bic_scores.append(g.bic(Xs))
    aic_scores.append(g.aic(Xs))

import pandas as pd
res = pd.DataFrame({'K': range(1, 9), 'BIC': bic_scores, 'AIC': aic_scores})
print(res.to_string(index=False))
print(f"\nBIC 最小 K = {res.loc[res['BIC'].idxmin(), 'K']}")
print(f"AIC 最小 K = {res.loc[res['AIC'].idxmin(), 'K']}")

📤 実行例:

💬 結果の読み方: 実測では BIC は K=5（12.33）が最小、 AIC は K=8（−230.85）が最小と、 2 つの基準が食い違った。 これは偶然ではなく両者の設計の違いがそのまま出たもので、 罰則項が BIC は $k\ln n$、 AIC は $2k$。 n=47 では $\ln 47 \approx 3.85 > 2$ なので BIC のほうがパラメータ数に厳しく、 小さい K を選ぶ。 5 次元・full 共分散だと K を 1 増やすごとにパラメータが 21 個も増えるため、 この差が効いてくる。 さらに BIC の並び（287.0 → 90.9 → 172.5 → 137.2 → 12.3 → 62.4 → 123.9 → 78.1）は単調でも凸でもなく上下に振れている。 n=47 に対してパラメータが多すぎ、 EM が K ごとに違う局所解へ落ちているためで、 「最小値の K」を機械的に採用するのは危うい。 実務では共分散を 'diag' に制限して自由度を落とすか、 前の実行で見た「小規模・中規模・大規模」の 3 クラスタのように解釈可能性を優先して K を決めるのが妥当である。

🐍 共分散構造の比較 (spherical / diag / tied / full)

このコードでやること: GMM の共分散構造を 4 種類変えて BIC を比較。 パラメータ数と汎化性能のトレードオフを観察。

1
2
3
4
5
6
7
8
9

for cov_type in ['spherical', 'diag', 'tied', 'full']:
    g = GaussianMixture(n_components=3, covariance_type=cov_type,
                        random_state=42, max_iter=200).fit(Xs)
    n_params = (g.weights_.size - 1) + g.means_.size + (
        3 if cov_type == 'spherical' else
        3 * 5 if cov_type == 'diag' else
        15 if cov_type == 'tied' else
        3 * 15)
    print(f"{cov_type:10s} BIC={g.bic(Xs):7.2f} params={n_params}")

📤 実行例:

💬 full が最良 BIC だが、 パラメータ数も最多。 47 県のような小サンプルでは diag や spherical で十分なことも多い。 共分散構造は事前知識（特徴間の相関の有無）で決めるのが王道。

🗺 学習ロードマップ

🗺 学習ロードマップ

1. レベル 1 — 潜在変数の概念。 観測尤度と完全データ尤度の区別。
2. レベル 2 — E ステップ・M ステップの 1 反復を手書きで（1 次元混合ガウス）。
3. レベル 3 — sklearn の GaussianMixture で SSDSE クラスタリング。 BIC で K 選定。
4. レベル 4 — HMM、 LDA、 因子分析、 確率的 PCA — 各種潜在変数モデル。
5. レベル 5 — 変分 EM、 オンライン EM、 Stochastic EM。
6. レベル 6 — VAE、 GAN、 拡散モデルでの潜在表現と EM の関係。

📊 比較表（兄弟手法・選択肢）

クラスタリング手法の比較

📖 用語ミニ辞典

🍳 コードレシピ（コピペ用 15 連発）

1
2

from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(3).fit(X)

labels = gmm.predict(X)

proba = gmm.predict_proba(X)

gmm.bic(X)

gmm.aic(X)

GaussianMixture(3, covariance_type='diag')

GaussianMixture(3, init_params='kmeans', n_init=10)

GaussianMixture(3, tol=1e-6, max_iter=500)

gmm.score(X)  # 1 サンプル当たりの平均対数尤度

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import numpy as np
import pandas as pd

# X（時系列として扱う観測値）を用意する。hmmlearn は 2 次元を求める
_d = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
_tokyo = _d[_d['Prefecture'] == '東京都'].sort_values('SSDSE-B-2026')
X = _tokyo[['A1101', 'A4101']].astype(float).to_numpy()

from hmmlearn.hmm import GaussianHMM
hmm = GaussianHMM(n_components=3).fit(X)

1
2

from sklearn.decomposition import LatentDirichletAllocation
lda = LatentDirichletAllocation(n_components=5).fit(X)

1
2
3

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
filled = IterativeImputer().fit_transform(X)

gmm.lower_bound_  # ELBO 風

gmm.means_

gmm.weights_

📊 SSDSE-B-2026 で EM を回す詳細レシピ

SSDSE-B-2026（都道府県別社会経済データ集 2026）を題材に、 EM アルゴリズムを段階的に適用するレシピを用意しました。 単一のコードブロックではなく、 探索→ 前処理 → モデル選択 → 解釈の 4 段階に分けて理解することで、 実務的な使い方が身につきます。

1. 探索的データ分析（EDA）: df.describe() で各列の平均・分散・歪度を確認。 SSDSE-B-2026 の A1101（人口）は東京が極端な外れ値で歪度が高いため、 対数変換が有効。
2. 標準化: StandardScaler で各列を平均 0・分散 1 に。 EM/GMM はスケールに敏感なので必須。
3. K の選定: K=2〜10 で BIC を計算し最小点を選ぶ。 都道府県データでは K=3〜5 が解釈しやすいことが多い。
4. 初期値多試行: n_init=10 以上で局所最適を回避。 init_params='kmeans' で k-means の解を初期値に使うと収束が速い。
5. 解釈: 各クラスタの平均ベクトルから「これは都市部」「これは地方」とラベル付け。 シルエット係数や ARI で他手法と比較。

SSDSE-B-2026 を使うと、 数値結果と地理的直感が一致するかどうかで EM の妥当性を検証できます。 例えば 東京・大阪・愛知が同一クラスタに入らない結果は、 特徴量の選択を見直すサインです。

この 4 段階レシピは、 EM/GMM だけでなく k-means、 階層クラスタリング、 DBSCAN など他のクラスタリング手法にも汎用的に応用できます。 SSDSE-B-2026 を題材にすれば、 数値結果の妥当性を「東京は都市部」「秋田は地方」など地理感覚と照合して直感的に検証できます。 これが他のデータセットでは難しい SSDSE-B-2026 ならではの利点です。

EM アルゴリズムが解く問題群の俯瞰：

                        【最尤推定】
                             │
                  ┌──────────┴──────────┐
            観測のみ                潜在変数あり
                 │                       │
            微分で解析的解        EM アルゴリズム
                                          │
                      ┌──────────┬───────┴──────┬──────────┐
                     GMM        HMM            LDA       欠損補完
                  (混合ガウス) (系列の潜在)  (トピック)   (IterativeImputer)

🔍 理論深掘り：EM が単調収束する理由（Jensen → Q 関数）

EM が「観測対数尤度を単調非減少にする」ことは、 Jensen の不等式と KL ダイバージェンスの非負性で示せます。 ここでは Bishop (2006) の議論に従って、 SSDSE-B-2026 を念頭にステップごとに丁寧に追います。

ステップ 1：観測尤度の下界（ELBO）を作る

観測 $X$、 潜在 $Z$、 パラメータ $\theta$。 任意の分布 $q(Z)$ について、 観測対数尤度を分解できます：

ここで $\mathcal{L}(q,\theta) = \mathbb{E}_q[\log p(X,Z|\theta)] - \mathbb{E}_q[\log q(Z)]$ は ELBO（Evidence Lower Bound）。 KL は常に非負なので、 $\log p(X|\theta) \ge \mathcal{L}(q,\theta)$ となります。

ステップ 2：E ステップで KL を 0 にする

$q(Z) = p(Z|X,\theta^{(t)})$ と置けば KL は 0 になり、 ELBO は観測対数尤度と一致。 これが「現パラメータでの潜在の事後分布」を求める E ステップの意味。

ステップ 3：M ステップで ELBO を最大化

$\theta^{(t+1)} = \arg\max_\theta \mathcal{L}(q,\theta) = \arg\max_\theta \mathbb{E}_{p(Z|X,\theta^{(t)})}[\log p(X,Z|\theta)]$。 これが Q 関数最大化。

ステップ 4：単調性の確認

$\log p(X|\theta^{(t+1)}) \ge \mathcal{L}(q^{(t)}, \theta^{(t+1)}) \ge \mathcal{L}(q^{(t)}, \theta^{(t)}) = \log p(X|\theta^{(t)})$。 これで「反復ごとに観測尤度が単調非減少」が示されました。

SSDSE-B-2026 での実証

47 都道府県の総人口・65 歳以上人口（A1101・A1303）を K=3 の GMM に当てはめると、 数反復で対数尤度（ELBO）が単調に増加して収束することが観察できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

import pandas as pd
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
X = df[['A1101','A1303']].astype(float).to_numpy()
X = StandardScaler().fit_transform(X)
gmm = GaussianMixture(n_components=3, n_init=10, max_iter=200,
                       random_state=0, warm_start=False).fit(X)
print('対数尤度の下界:', gmm.lower_bound_)
print('反復回数:', gmm.n_iter_)
print('収束:', gmm.converged_)

🧪 数値計算上の注意点（実装で必ず踏む罠）

1. log-sum-exp トリック

E ステップで $\gamma_{ik} = \pi_k \mathcal{N}(x_i|\mu_k,\Sigma_k) / \sum_l \pi_l \mathcal{N}(x_i|\mu_l,\Sigma_l)$ を計算するとき、 高次元では分子・分母ともに極端に小さく、 オーバーフロー／アンダーフローを起こします。 ログ空間で計算して log-sum-exp を使うのが定石。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

import numpy as np
import pandas as pd

# log_normal_pdf / X / pi / mu / sigma を用意する（1 次元 2 成分の混合）
def log_normal_pdf(x, mu, sigma):
    x = np.asarray(x).reshape(-1, 1)
    return -0.5 * np.log(2 * np.pi * sigma ** 2) - (x - mu) ** 2 / (2 * sigma ** 2)

_d = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
_d = _d[_d['SSDSE-B-2026'] == 2023]
X = np.log(_d['A1101'].astype(float).to_numpy()).reshape(-1, 1)   # sklearn は 2 次元を求める
pi = np.array([0.5, 0.5])
mu = np.array([X.min(), X.max()])
sigma = np.array([X.std(), X.std()])

import numpy as np
from scipy.special import logsumexp
# log γ_ik = log π_k + log N(x|μ_k,Σ_k) - logsumexp over k
log_w = np.log(pi)[None, :] + log_normal_pdf(X, mu, sigma)  # (n,K)
log_gamma = log_w - logsumexp(log_w, axis=1, keepdims=True)
gamma = np.exp(log_gamma)

2. 共分散行列の正則化

1 つのクラスタにデータがほぼ 1 点しか入らないと、 $\Sigma_k$ が特異になり対数尤度が発散します（縮退）。 sklearn の reg_covar=1e-6 がデフォルトで対角に微小値を加える正則化。

1
2
3
4

from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=5, covariance_type='full',
                       reg_covar=1e-4,  # 既定 1e-6 では小さすぎる場合
                       n_init=20, max_iter=500, random_state=0).fit(X)

3. 初期化の重要性

EM は局所最適に陥る。 sklearn は既定で init_params='kmeans'（k-means の解で初期化）。 n_init=10 で 10 回試行し最良を採用する。

4. 収束判定の閾値

ELBO の改善量 $\Delta \mathcal{L} < \epsilon$ で停止。 既定は $\epsilon = 10^{-3}$。 精度が要るなら $10^{-6}$ まで小さく。

🌀 EM の変種・拡張

1. Generalized EM (GEM)

M ステップで Q を完全に最大化せず、 増やすだけでも単調性は保たれる。 計算が軽い。

2. Stochastic EM (SEM)

E ステップで潜在変数を分布からサンプリング。 ノイズが入る代わりに局所最適から抜けやすい。

3. Online EM

大規模データで全データを保持できないとき、 ミニバッチで増分的に更新。

4. 変分 EM (VBEM)

E ステップで $q(Z)$ を近似分布族に制限。 ベイズ的に事前分布も入れられる。 LDA や VAE の基礎。

5. ECM / ECME

Expectation Conditional Maximization。 M ステップを複数の条件付き最大化に分解。 Meng & Rubin (1993)。

6. PX-EM (Parameter Expansion)

パラメータ空間を拡大して収束を加速。 Liu, Rubin, Wu (1998)。

7. SAEM (Stochastic Approximation EM)

SEM とロビンス・モンロー型確率近似の融合。 強い理論保証。

8. MCEM (Monte Carlo EM)

E ステップで期待値を MC サンプリングで近似。 複雑な潜在モデルで活躍。

📈 経験的振る舞い：SSDSE-B-2026 での K 選定

SSDSE-B-2026 の 47 都道府県を様々な特徴量集合で GMM クラスタリングしたときの BIC を比較します。

特徴量セットによって BIC 最小の K は 2〜5 と変わり、 都市部・地方中核・過疎地域・特殊県のような構造が現れます。 なお「全数値列」（112 列）は n=47 に対して次元が大きすぎ、 full 共分散の GMM が過剰適合して BIC が K とともに発散します（次元の呪い）。 実務では少数の解釈可能な特徴量に絞るのが鉄則です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

import pandas as pd
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
cols = ['A1101','A1303','A4101','A4200']  # 人口・65 歳以上・出生・死亡
X = StandardScaler().fit_transform(df[cols].astype(float))

bics = []
for k in range(1, 9):
    g = GaussianMixture(n_components=k, n_init=20, max_iter=500,
                         random_state=0).fit(X)
    bics.append((k, g.bic(X), g.aic(X)))
for k, b, a in bics:
    print(f'K={k}  BIC={b:.1f}  AIC={a:.1f}')
best_k = min(bics, key=lambda t: t[1])[0]
print('best K (BIC):', best_k)

📉 収束モニタリング：対数尤度の軌跡を見る

EM は単調収束しますが、 収束速度はパラメータの初期化と問題の凸性で大きく変わります。 SSDSE-B-2026 で 47 都道府県を GMM K=3 にしたときの対数尤度履歴を可視化します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
X = StandardScaler().fit_transform(df[['A1101','A1303']].astype(float))

histories = []
for seed in range(5):
    history = []
    gmm = GaussianMixture(n_components=3, n_init=1, max_iter=1,
                          warm_start=True, random_state=seed)
    for it in range(50):
        gmm.fit(X)
        history.append(gmm.score(X) * len(X))
    histories.append(history)

fig, ax = plt.subplots(figsize=(8, 5))
for seed, h in enumerate(histories):
    ax.plot(h, label=f'seed={seed}')
ax.set_xlabel('反復回数')
ax.set_ylabel('対数尤度')
ax.set_title('SSDSE 47 都道府県 GMM K=3 — 初期値ごとの収束')
ax.legend()
plt.savefig('em_convergence.png', dpi=100)

⚔️ GMM vs k-means：SSDSE-B-2026 で比較

同じデータで GMM と k-means を実行し、 ARI（Adjusted Rand Index）で結果の一致度を測ります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import pandas as pd
from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
X = StandardScaler().fit_transform(df[['A1101','A1303']].astype(float))

gmm = GaussianMixture(n_components=3, n_init=20, random_state=0).fit(X)
km = KMeans(n_clusters=3, n_init=20, random_state=0).fit(X)
labels_g = gmm.predict(X)
labels_k = km.labels_

print('ARI:', adjusted_rand_score(labels_g, labels_k))
print('GMM シルエット:', silhouette_score(X, labels_g))
print('k-means シルエット:', silhouette_score(X, labels_k))
print('GMM 事後確率の最大（境界の県）:')
import numpy as np
margin = gmm.predict_proba(X).max(axis=1)
border_idx = np.argsort(margin)[:5]
print(df.iloc[border_idx][['Prefecture']].values.ravel())
print('確率:', margin[border_idx])

🎼 GMM 共分散タイプの比較

GaussianMixture には 4 種類の共分散構造があります。 SSDSE-B-2026 47 都道府県で比較。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

import pandas as pd
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
X = StandardScaler().fit_transform(
    df[['A1101','A1303','A4101','A4200']].astype(float))

for ct in ['spherical','diag','tied','full']:
    g = GaussianMixture(n_components=3, covariance_type=ct,
                         n_init=20, random_state=0).fit(X)
    print(f'{ct:10s}  BIC={g.bic(X):.1f}  AIC={g.aic(X):.1f}')

🧭 ベイズ的視点：EM ＝ MAP の特殊形

EM は最尤推定（MLE）を反復で実行する手法ですが、 事前分布を加えると MAP（Maximum A Posteriori）推定に拡張できます。 これを「Bayesian EM」または「MAP EM」と呼びます。 SSDSE-B-2026 のように小標本（n=47）では、 事前分布で正則化すると安定します。

BayesianGaussianMixture との関係

scikit-learn の BayesianGaussianMixture は、 混合比に Dirichlet 事前、 平均・共分散に Normal-Wishart 事前を置いた変分 EM。 「データに支持されない余分なクラスタの混合比を自動的に 0 に近づける」性質があり、 K 選定が緩くなります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

import pandas as pd
from sklearn.mixture import BayesianGaussianMixture
from sklearn.preprocessing import StandardScaler

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
X = StandardScaler().fit_transform(
    df[['A1101','A1303','A4101','A4200']].astype(float))

bgm = BayesianGaussianMixture(n_components=10,  # 大きめに設定
                                weight_concentration_prior_type='dirichlet_process',
                                weight_concentration_prior=0.1,  # 小さいほど少数クラスタを好む
                                n_init=10, max_iter=500,
                                random_state=0).fit(X)
print('実効クラスタ数（重み > 0.01）:', (bgm.weights_ > 0.01).sum())
print('weights:', bgm.weights_.round(3))

🤖 EM と VAE（変分オートエンコーダ）の橋渡し

近年深層学習で広く使われる VAE（Variational Autoencoder）は、 EM の精神を引き継ぎ、 ニューラルネットワークでパラメータ化したものと見なせます。

SSDSE-B-2026 のように小標本では GMM/EM が圧勝。 画像・テキストのような大規模・高次元では VAE が優位。 「EM の発想 → 変分推論 → VAE」というロードマップは、 現代統計＆ML の中核。

🧩 欠損データ補完への EM 応用（詳細）

SSDSE-B-2026 はほぼ完全データですが、 一部の年度・項目に欠損があります。 EM で多変量正規分布を仮定して補完する古典的手法を紹介。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import pandas as pd
import numpy as np
from sklearn.experimental import enable_iterative_imputer  # noqa
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
cols = ['A1101','A1303','A4101','A4200','L3221']  # 数値 5 列（人口・65歳以上・出生・死亡・消費支出）
sub = df[cols].astype(float)

# 5% を人工的に欠損させる（実データ忠実なので「無作為に 5% 削除」のみ）
sub_missing = sub.copy()
mask = np.tri(len(sub), len(cols), dtype=bool)  # 三角マスクで代用
# 注: 実 SSDSE は元来欠損がほぼなく、 ここでは方法論のデモ
sub_missing.iloc[10:15, 2] = np.nan

# IterativeImputer は EM 風の反復補完（各列を他列で予測 → 繰り返し）
imp = IterativeImputer(estimator=BayesianRidge(), max_iter=30, tol=1e-4)
filled = imp.fit_transform(sub_missing)
print('補完前 NaN 数:', sub_missing.isna().sum().sum())
print('補完後 NaN 数:', np.isnan(filled).sum())
# 真値との比較（行 10-14、 列 2）
print('真値:', sub.iloc[10:15, 2].to_numpy())
print('補完値:', filled[10:15, 2])

⛓ HMM の Baum-Welch（EM の系列版）

HMM（Hidden Markov Model）における EM 推定アルゴリズムは「Baum-Welch」と呼ばれ、 EM の有名な特殊形です。 SSDSE-B-2026 のような年度別パネル時系列（2023, 2022, 2021…）にも応用できます。

HMM の構造

• 状態 $S_t \in \{1,...,K\}$：時刻 $t$ の潜在状態（例：景気「拡大」「停滞」「収縮」）
• 観測 $X_t$：状態に依存する確率分布から生成（例：消費支出の変化率）
• 遷移確率 $A_{ij} = P(S_{t+1}=j|S_t=i)$
• 出力確率 $B_i(x) = P(X_t=x|S_t=i)$

Baum-Welch の流れ

1. 前向き $\alpha_t(i)$：$P(X_{1:t}, S_t=i|\lambda)$
2. 後ろ向き $\beta_t(i)$：$P(X_{t+1:T}|S_t=i,\lambda)$
3. E ステップ：$\gamma_t(i) = \alpha_t(i)\beta_t(i) / \sum_j \alpha_t(j)\beta_t(j)$、 $\xi_t(i,j) = \alpha_t(i) A_{ij} B_j(X_{t+1}) \beta_{t+1}(j) / Z$
4. M ステップ：$\hat A_{ij} = \sum_t \xi_t(i,j) / \sum_t \gamma_t(i)$、 $\hat\pi_i = \gamma_1(i)$、 $\hat B_i$ は重み付き再推定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import pandas as pd
import numpy as np
from hmmlearn.hmm import GaussianHMM

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
# 東京都（Code='R13000'）の年度別データ（2012〜2023 の 12 年分）を抽出
tokyo = df[df['Code']=='R13000'].sort_values('SSDSE-B-2026')
X = tokyo[['A1101','A1303']].astype(float).to_numpy()

# 3 状態 HMM
hmm = GaussianHMM(n_components=3, covariance_type='full',
                   n_iter=100, random_state=0)
hmm.fit(X)
print('遷移行列 A:')
print(hmm.transmat_.round(3))
print('状態系列:', hmm.predict(X))
print('対数尤度:', hmm.score(X))

🎓 上級者向け FAQ（10 問）

🛠 ハンズオン演習：SSDSE-B-2026 で EM 完全攻略

以下は、 SSDSE-B-2026 を題材にした 6 段階のハンズオン演習です。 各段階で実行可能なコードと期待される出力を示します。

段階 1：データのロードと探索

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import pandas as pd
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
print('行数:', len(df), ' 列数:', df.shape[1])
print('年度の範囲:', df['SSDSE-B-2026'].min(), '〜', df['SSDSE-B-2026'].max())
print('都道府県数:', df['Prefecture'].nunique())

# 2023 年だけに絞る
d23 = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
print('2023 年データ:', len(d23), '件')
print(d23[['Prefecture','A1101','A1303','A4101','A4200']].head())

段階 2：特徴量エンジニアリング

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import pandas as pd
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d23 = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)

# 派生特徴量
d23['aging_rate'] = d23['A1303'].astype(float) / d23['A1101'].astype(float)
d23['birth_rate'] = d23['A4101'].astype(float) / d23['A1101'].astype(float) * 1000
d23['death_rate'] = d23['A4200'].astype(float) / d23['A1101'].astype(float) * 1000
print(d23[['Prefecture','aging_rate','birth_rate','death_rate']].head())
print(d23[['aging_rate','birth_rate','death_rate']].describe().round(3))

段階 3：標準化と GMM 訓練

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d23 = df[df['SSDSE-B-2026']==2023].reset_index(drop=True)
d23['aging_rate'] = d23['A1303'].astype(float) / d23['A1101'].astype(float)
d23['birth_rate'] = d23['A4101'].astype(float) / d23['A1101'].astype(float) * 1000
d23['death_rate'] = d23['A4200'].astype(float) / d23['A1101'].astype(float) * 1000

scaler = StandardScaler()
X = scaler.fit_transform(d23[['aging_rate','birth_rate','death_rate']])
gmm = GaussianMixture(n_components=4, covariance_type='full',
                       n_init=20, max_iter=500, random_state=42).fit(X)
d23['cluster'] = gmm.predict(X)
d23['confidence'] = gmm.predict_proba(X).max(axis=1)
for k in range(4):
    members = d23[d23['cluster']==k]['Prefecture'].tolist()
    print(f'クラスタ {k}: {len(members)} 県 — {members[:5]} ...')

段階 4：クラスタ可視化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

import matplotlib.pyplot as plt

# 2D に PCA で射影してプロット
from sklearn.decomposition import PCA
pca = PCA(n_components=2).fit(X)
Xp = pca.transform(X)
fig, ax = plt.subplots(figsize=(9, 7))
for k in range(4):
    mask = d23['cluster']==k
    ax.scatter(Xp[mask, 0], Xp[mask, 1], s=80, label=f'C{k}', alpha=0.7)
    for i, prefname in d23[mask][['Prefecture']].iterrows():
        ax.annotate(prefname['Prefecture'], (Xp[i, 0], Xp[i, 1]),
                    fontsize=8, alpha=0.7)
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.legend()
ax.set_title('SSDSE-B-2026 47 都道府県 GMM K=4 クラスタ（PCA 投影）')
plt.savefig('em_clusters.png', dpi=120, bbox_inches='tight')

段階 5：クラスタ意味づけと指標

1
2
3
4
5
6
7
8
9

import pandas as pd
# クラスタごとの平均で意味づけ
cluster_stats = d23.groupby('cluster')[['aging_rate','birth_rate','death_rate']].agg(['mean','std']).round(3)
print(cluster_stats)
print()
print('クラスタごとの代表的都道府県（事後確率最大の 3 件）:')
for k in range(4):
    sub = d23[d23['cluster']==k].sort_values('confidence', ascending=False).head(3)
    print(f'  C{k}: {sub["Prefecture"].tolist()}')

段階 6：時系列での追跡（2014→2023）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])

years = [2014, 2018, 2023]
for yr in years:
    d = df[df['SSDSE-B-2026']==yr].reset_index(drop=True)
    d['aging'] = d['A1303'].astype(float) / d['A1101'].astype(float)
    d['birth'] = d['A4101'].astype(float) / d['A1101'].astype(float) * 1000
    X = StandardScaler().fit_transform(d[['aging','birth']])
    gmm = GaussianMixture(n_components=3, n_init=10, random_state=42).fit(X)
    labels = gmm.predict(X)
    print(f'年度 {yr}: 平均高齢化率 {d["aging"].mean():.3f}')
    for k in range(3):
        n_in_cluster = (labels==k).sum()
        mean_aging = d.loc[labels==k, 'aging'].mean()
        print(f'  C{k}: {n_in_cluster} 県, 高齢化率 {mean_aging:.3f}')

📖 拡張用語集（EM 周辺の専門用語 30 件）

🔗 隣接手法への橋渡し

EM アルゴリズムは単独の最適化手法ではなく、 潜在変数モデル (混合正規 ・隠れマルコフ ・LDA) を統一的に扱う一般枠組みである。 MAP 推定 ・変分ベイズ ・MCMC との位置関係を意識して使う。

• 上流（前処理・データ整備）: 最尤推定 (ロジスティック回帰) / 確率分布 / 欠損値処理
• 並列（代替・補完手法）: k-means / 変分推論 (ベイズ) / MCMC (ベイズ推論)
• 下流（評価・解釈・応用）: ガウス混合モデル / 時系列モデル (ARIMA) / クラスタリング

EM アルゴリズムは「潜在変数モデルの最尤推定を反復で解く」一般手法で、 上流のモデル仮定 (混合正規・隠れマルコフ等) で E/M 式が決まり、 並列の MAP・変分ベイズと比較され、 下流のクラスタリング・欠測値補完で広く使われる。

🌳 手法選択フロー

EM を採るかは「潜在変数の存在・分布の解析的扱いやすさ・他手法 (変分・MCMC) との比較」の 3 軸で決まる。 共役な指数族なら EM が高速、 複雑分布なら変分推論や MCMC が必要になる。

1. 欠測や隠れ変数があるか
   あるなら EM の出番。 混合分布の所属クラスタ、 欠測値、 潜在因子などが典型。 全部観測できているなら、 素直に最尤推定すればよい。
2. 初期値をどう決めるか
   EM は局所最適に落ちるので、 初期値で結果が変わる。 複数の初期値から走らせ、 対数尤度が最大になったものを採る。
3. 収束をどう判定するか
   対数尤度の増分が閾値を下回ったら止める。 尤度は単調増加するので、 減ったら実装の誤り。
4. 結果をどう検証するか
   クラスタ数のような設定は EM 自体では決まらない。 AIC・BIC や交差検証で外側から選ぶ。

SSDSE-B-2026 の消費支出（L3221）を 2-3 個の正規分布の混合とみなしてクラスタリングするとき、 まず初期化を 2 通り試して同じ解に収束するか確認し、 BIC で混合数を決めるのが安定した運用手順である。

📝 補足：EM をもう一歩深く理解する

※ 本節は既存の解説を壊さずに追記した「深掘り補足」です。 上の各節と重複する話題もありますが、 直感・落とし穴・発展を一望できるまとめとして活用してください。

🎨 直感：「鶏と卵」を反復で断ち切る

EM の本質は、 潜在変数がある最尤推定でぶつかる「鶏と卵」問題を、 交互更新で解くところにあります。 「各データがどのクラスタから来たか（＝潜在変数 $Z$）が分かればパラメータ $\theta$ は簡単に推定できる」。 逆に「$\theta$ が分かれば $Z$ の分布は簡単に計算できる」。 どちらも相手が分からないと決まらない — これが鶏と卵です。

EM はこれを、 E（期待）ステップで「今の $\theta$ のもとで潜在変数の期待値（責任度 $\gamma$）を埋める」→ M（最大化）ステップで「その期待値を真値とみなして $\theta$ を更新する」という反復で断ち切ります。 重要なのは、 この 2 ステップを繰り返すたびに観測データの対数尤度が単調に増える（または横這い）ことが数学的に保証されている点です（上の「理論深掘り」節の ELBO 分解を参照）。 だから「毎回良くなるのに、 なぜ止まると分かるのか」に悩まなくてよい — 増える一方で上に有界なので、 必ずどこかに収束します。

• E ステップ＝相手（潜在変数）の「今のところの見込み」を確率で埋める（ソフトな穴埋め）
• M ステップ＝埋めた見込みを使ってパラメータを最尤更新する（重み付き平均・分散・混合比）
• この往復が、 欠損補完・混合分布・隠れマルコフに共通する「EM の骨格」

⚠️ 落とし穴（重要）：EM が「うまくいったように見えて外している」典型

🚀 発展：EM の骨格が見えると視界が広がる

• 混合ガウス（GMM）への適用 — EM の代表的応用。 E で責任度 $\gamma_{ik}$、 M で重み付き平均・共分散・混合比を更新。 各点を確率で複数クラスタに割り振るソフトクラスタリングが得られる（上の SSDSE 実測例を参照）。
• k-means は EM の「ハード版」 — 各点をどれか 1 つのクラスタに割り切る（責任度が 0/1）、 かつ全成分の共分散を等方・等分散に固定した GMM が k-means に一致します。 EM はより柔軟な確率モデルを許す一般化、 k-means は軽くて速い特殊ケース、 という関係です。
• 下界最大化（ELBO）の視点 — EM は「観測対数尤度を直接上げる」のではなく、 その下界 ELBO を E と M で交互に押し上げるアルゴリズムと読み替えられます。 E ステップ＝$q(Z)$ について ELBO 最大化（KL を 0 に）、 M ステップ＝$\theta$ について ELBO 最大化。 この見方が変分推論への橋になります。
• 変分ベイズ／EM の一般化 — E ステップの事後分布 $p(Z|X,\theta)$ が解析的に求まらないとき、 近似分布族 $q(Z)$ に制限して ELBO を最大化するのが変分 EM（VBEM）。 さらにニューラルネットで $q$ を表現すれば VAE、 サンプリングで期待値を近似すれば MCEM／Stochastic EM になります（上の「変種・拡張」節を参照）。
• 欠測データへの適用 — EM はそもそも Dempster–Laird–Rubin (1977) が欠測データ解析のために定式化した枠組みです。 欠測値を潜在変数とみなし、 E で条件付き期待値を補完、 M でパラメータを再推定。 MAR 仮定のもとで妥当な最尤推定が得られ、 多重代入法（MICE／IterativeImputer）の理論的基礎にもなっています（MNAR では欠測機序をモデルに含める必要あり）。

🔗 関連ページ（前提・並列・発展）

EM の理解を広げる関連用語ページ：

※ 「最尤推定（の尤度そのもの）」「潜在変数」の独立した用語ページは本用語集には未整備のため、 前提概念としてはテキストで補足します：尤度はパラメータをデータで評価する関数、 潜在変数は観測されない隠れ変数（クラスタ所属・欠測値など）で、 EM はこの潜在変数を含む尤度の最大化を反復で解く手法です。