← 統計検定テキスト 一覧

📊 対象級:1級 | 重要度:C(低頻度)

要点(BLUF)

データに穴(欠測)があるとき、何が起きるかはなぜ穴が空いたかで決まります。穴の空き方は3種類——MCAR(完全にランダム)・MAR(観測値で説明できるランダム)・MNAR(欠測値そのものに依存)——で、解析の難しさが段違いです。

そのうえで本題は EMアルゴリズム。これは「見えていないデータ(潜在変数・欠測値)」を含む状況で最尤推定をやり切るための反復法です。核心は2ステップを繰り返すこと:

  E(期待値)ステップ:Q(θθ(t))=E ⁣[logLc(θ)|xobs,θ(t)]M(最大化)ステップ:θ(t+1)=argmaxθQ(θθ(t))  \boxed{\; \begin{aligned} \text{E(期待値)ステップ:}&\quad Q(\theta\mid\theta^{(t)})=E\!\left[\log L_c(\theta)\,\middle|\,x_{\mathrm{obs}},\theta^{(t)}\right]\\[2pt] \text{M(最大化)ステップ:}&\quad \theta^{(t+1)}=\arg\max_\theta\, Q(\theta\mid\theta^{(t)}) \end{aligned} \;}

要するに「見えない部分を今のパラメータで埋めた『つもり』の対数尤度を作り(E)、それを最大化する(M)。これを繰り返す」。そして本ノートの主役は、この反復が観測データの尤度を決して下げないこと(単調増加性)の証明です。1級(統計数理)ではEステップ・Mステップの定式化、単調増加性の論述、正規混合分布への適用が問われます(出題頻度は高くないが理論問題として狙われる。範囲・配点は要最新確認)。


1. 欠測データのメカニズム(MCAR / MAR / MNAR)

EMに入る前に、「そもそも欠測をどう扱ってよいか」を決める前提を押さえます。これは Rubin が整理した欠測メカニズムの分類で、解析の妥当性を左右します。

記号を決めます。完全データを X=(Xobs,Xmis)X=(X_{\mathrm{obs}},X_{\mathrm{mis}}) に分け(観測された部分 XobsX_{\mathrm{obs}} と欠測した部分 XmisX_{\mathrm{mis}})、各値が欠測したかどうかを表す欠測指標 RR(観測なら1、欠測なら0)を導入します。問題は「RR がどの変数に依存して決まるか」です。

メカニズム欠測確率 P(RX)P(R\mid X) の依存先直観尤度ベース解析で無視できるか
MCAR(完全無作為欠測)何にも依存しない(P(RX)=P(R)P(R\mid X)=P(R)サイコロで欠測が決まる。観測者の都合で紙が破れた等できる(最も安全)
MAR(無作為欠測)観測値 XobsX_{\mathrm{obs}} のみに依存高齢者ほど所得を未回答。年齢(観測済)で欠測が説明できるできる(無視可能=ignorable)
MNAR(無作為でない欠測)欠測値 XmisX_{\mathrm{mis}} 自身に依存(観測値で条件付けても残る)所得が高い人ほど所得を隠す。隠した所得の値そのものが欠測理由できない(無視不可能)

数式で書くと差は明確です。ψ\psi を欠測メカニズムのパラメータとして、

MCAR:P(RXobs,Xmis,ψ)=P(Rψ)\text{MCAR:}\quad P(R\mid X_{\mathrm{obs}},X_{\mathrm{mis}},\psi)=P(R\mid\psi)

要するに「欠測の有無はデータの中身と完全に無関係」。

MAR:P(RXobs,Xmis,ψ)=P(RXobs,ψ)\text{MAR:}\quad P(R\mid X_{\mathrm{obs}},X_{\mathrm{mis}},\psi)=P(R\mid X_{\mathrm{obs}},\psi)

要するに「観測できている部分で条件付ければ、欠測の有無は欠測値の中身に依存しない」。

MNAR:P(RXobs,Xmis,ψ) が Xmis に真に依存\text{MNAR:}\quad P(R\mid X_{\mathrm{obs}},X_{\mathrm{mis}},\psi)\ \text{が}\ X_{\mathrm{mis}}\ \text{に真に依存}

要するに「観測値でどんなに条件付けても、欠測の理由が欠測値自身に残る」。

1.1 なぜ MAR で欠測を「無視」できるのか(ignorability)

ここは1級で問われうる理論的な勘所です。観測されたデータ (Xobs,R)(X_{\mathrm{obs}},R) の同時尤度は、欠測値を積分消去して

L(θ,ψXobs,R)=P(Xobs,Xmisθ)P(RXobs,Xmis,ψ)dXmisL(\theta,\psi\mid X_{\mathrm{obs}},R)=\int P(X_{\mathrm{obs}},X_{\mathrm{mis}}\mid\theta)\,P(R\mid X_{\mathrm{obs}},X_{\mathrm{mis}},\psi)\,dX_{\mathrm{mis}}

と書けます。要するに「見えなかった XmisX_{\mathrm{mis}} のあらゆる可能性で平均した尤度」です。ここで MAR だと P(R,ψ)=P(RXobs,ψ)P(R\mid\cdot,\psi)=P(R\mid X_{\mathrm{obs}},\psi) が積分の外に出せて、

L=P(RXobs,ψ)θを含まない  ×  P(Xobs,Xmisθ)dXmis=P(Xobsθ)L=\underbrace{P(R\mid X_{\mathrm{obs}},\psi)}_{\theta\,\text{を含まない}}\;\times\;\underbrace{\int P(X_{\mathrm{obs}},X_{\mathrm{mis}}\mid\theta)\,dX_{\mathrm{mis}}}_{=\,P(X_{\mathrm{obs}}\mid\theta)}

積に分離します。要するに「興味のあるパラメータ θ\theta の最尤推定をするには、欠測メカニズムの項 P(RXobs,ψ)P(R\mid X_{\mathrm{obs}},\psi) を無視して P(Xobsθ)P(X_{\mathrm{obs}}\mid\theta) だけ最大化すればよい」。これが**無視可能性(ignorability)**です。MNAR では P(R)P(R\mid\cdot)XmisX_{\mathrm{mis}} に依存して積分の外に出せないため、この分離が成立せず、欠測メカニズムを同時にモデル化しないとバイアスが残ります。

重要な注意:「無視できる」のは尤度ベース(最尤・ベイズ)の解析の話です。後述の完全ケース分析が無バイアスになるのは、より強い MCAR が必要です。MAR は「尤度を正しく書けば(EM等を使えば)大丈夫」というレベルであり、MAR で完全ケース分析をするとバイアスが出ます。MAR と MCAR の混同は最頻出の引っかけです。


2. 素朴な対処はなぜ危険か

EMやその他の正攻法に行く前に、つい手が伸びる「素朴な対処」がどんなバイアスを生むかを押さえます。1級の論述で「なぜ問題か」を説明させられることがあります。

2.1 完全ケース分析(リストワイズ除去)

欠測のある行(ケース)を丸ごと捨て、完全に揃ったケースだけで解析する方法です。最も手軽ですが:

2.2 平均値代入(mean imputation)

欠測を、その変数の観測値の平均で埋める方法です。一見もっともらしいですが、分散を構造的に過小評価するという重大な欠陥があります。理由を式で見ます。

ある変数の観測値が mm 個、欠測が nmn-m 個あり、観測値の平均を xˉobs\bar x_{\mathrm{obs}}、不偏分散を s2s^2 とします。欠測 nmn-m 個をすべて xˉobs\bar x_{\mathrm{obs}} で埋めると、埋めた値はすべて平均そのものなので平均からのズレがゼロ。埋めたあとの標本分散の分子(偏差平方和)は

観測(xixˉobs)2+代入(xˉobsxˉobs)2=0=観測(xixˉobs)2=(m1)s2\sum_{\text{観測}}(x_i-\bar x_{\mathrm{obs}})^2+\sum_{\text{代入}}\underbrace{(\bar x_{\mathrm{obs}}-\bar x_{\mathrm{obs}})^2}_{=\,0}=\sum_{\text{観測}}(x_i-\bar x_{\mathrm{obs}})^2=(m-1)s^2

要するに「代入したぶんは平方和に1ミリも寄与しない」。ところが分母(標本数)は nn に増えるので、埋めたあとの分散はおよそ

Var^代入後(m1)s2n1=m1n1s2  <  s2\widehat{\mathrm{Var}}_{\text{代入後}}\approx\frac{(m-1)s^2}{n-1}=\frac{m-1}{n-1}\,s^2\;<\;s^2

となり、真の分散の m1n1\tfrac{m-1}{n-1} 倍に縮みます。要するに「欠測を平均で埋めると、ばらつきが系統的に小さく見積もられる」。さらに他変数との相関も平均値の一点に潰れるぶん弱まる(相関の希釈)ため、回帰係数や相関係数まで歪みます。点推定(平均)が運よく合っても、ばらつき・関係性の推定が壊れるのが平均値代入の本質的な問題です。

2.3 多重代入(multiple imputation, 概念)

平均値代入の欠陥(不確かさの無視)への正攻法のひとつが多重代入です。考え方は3段階:

  1. 欠測値を「1つの数」ではなく、欠測の予測分布からランダムに何通りも(例:M=550M=5\sim50 通り)埋めた完備データセットを作る。
  2. 各完備データセットで普通に解析し、推定値と標準誤差を得る。
  3. MM 個の結果をRubin の規則で統合する。点推定は平均、分散は「各回内の分散の平均」+「各回の点推定間のばらつき」で合算する。

要するに「埋め方を1通りに決めつけず、欠測の不確かさを複数の埋め方のばらつきとして分散に正しく上乗せする」。これにより平均値代入のような分散の過小評価を避けられます。多重代入は MAR を前提にする手法で、EMと並ぶ欠測対処の主役ですが、本ノートの主題は EM なので概念に留めます。

graph TD
  A["欠測への対処"] --> B["素朴な対処(バイアスあり)"]
  A --> C["正攻法"]
  B --> B1["完全ケース分析<br/>MCAR以外でバイアス"]
  B --> B2["平均値代入<br/>分散を過小評価"]
  C --> C1["EMアルゴリズム<br/>不完全データの最尤推定"]
  C --> C2["多重代入<br/>不確かさを分散に反映"]
  C --> C3["完全情報最尤(FIML)"]

3. EMアルゴリズム(本題)

EMは「見えていないデータを含む最尤推定」の一般手法です。欠測値だけでなく、混合分布の所属ラベル・打ち切り・潜在変数など、「完全データなら最尤推定が簡単なのに、一部が見えないせいで尤度が複雑になる」あらゆる場面に使えます。

3.1 完全データ尤度と観測データ尤度

考え方の出発点は、「見えていたら簡単だったのに」という構造を逆手に取ることです。

要するに「最大化したい観測尤度は積分のせいで厄介。でも積分の中身(完全データ尤度)は素直」。EMはこのギャップを反復で埋めます。最尤法そのものの基礎は 最尤法・モーメント法(推定量の作り方と最尤推定量の漸近論)xmisx_{\mathrm{mis}} を積分消去して周辺化する操作は 同時分布・周辺分布・条件付き分布 の周辺分布です。

3.2 Eステップ:完全データ対数尤度の条件付き期待値

xmisx_{\mathrm{mis}} は見えないので、完全データ対数尤度 logLc(θ)\log L_c(\theta) をそのままは計算できません。そこで、今のパラメータ推定値 θ(t)\theta^{(t)} のもとで xmisx_{\mathrm{mis}} がどうだったかの条件付き分布で平均します。これが QQ 関数です。

Q(θθ(t))=E ⁣[logLc(θ)|xobs,θ(t)]=logp(xobs,xmisθ)  p(xmisxobs,θ(t))dxmisQ(\theta\mid\theta^{(t)}) =E\!\left[\log L_c(\theta)\,\middle|\,x_{\mathrm{obs}},\theta^{(t)}\right] =\int \log p(x_{\mathrm{obs}},x_{\mathrm{mis}}\mid\theta)\;p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta^{(t)})\,dx_{\mathrm{mis}}

要するに「見えない xmisx_{\mathrm{mis}} を、今の推定値で予想した分布で平均して、完全データ対数尤度の『期待値版』を作る」。ここで θ\theta(最大化したい引数)と θ(t)\theta^{(t)}(期待値を取るための現在値)は別物であることが決定的に重要です。期待値(積分)に使うのは θ(t)\theta^{(t)} で固定し、QQ は残った θ\theta の関数になります。

3.3 Mステップ:QQ の最大化

作った Q(θθ(t))Q(\theta\mid\theta^{(t)})θ\theta について最大化し、次の推定値とします。

θ(t+1)=argmaxθQ(θθ(t))\theta^{(t+1)}=\arg\max_\theta\, Q(\theta\mid\theta^{(t)})

要するに「期待値版の対数尤度を、普通の最尤推定と同じように最大化する」。QQ は完全データ対数尤度の期待値なので、完全データ尤度が指数型分布族なら QQ の最大化も閉形式で解け、複雑な観測尤度の直接最大化を、易しい問題の繰り返しに置き換えられるのがEMの旨味です。

flowchart TD
  S["初期値 θ⁽⁰⁾ を決める"] --> E
  E["E ステップ:<br/>Q(θ &#124; θ⁽ᵗ⁾) = E[ log Lc(θ) &#124; x_obs, θ⁽ᵗ⁾ ] を計算<br/>(見えない部分を現在値で期待値化)"] --> M
  M["M ステップ:<br/>θ⁽ᵗ⁺¹⁾ = argmax_θ Q(θ &#124; θ⁽ᵗ⁾)<br/>(期待値版の対数尤度を最大化)"] --> C{"収束したか?<br/>対数尤度の増分が小"}
  C -- いいえ --> E
  C -- はい --> D["局所最適 θ̂ を出力"]

4. 単調増加性の証明(1級の核心)

EMがなぜ正当なのか——反復ごとに観測データの対数尤度が決して減らないこと

  logp(xobsθ(t+1))  logp(xobsθ(t))  \boxed{\;\log p(x_{\mathrm{obs}}\mid\theta^{(t+1)})\ \ge\ \log p(x_{\mathrm{obs}}\mid\theta^{(t)})\;}

を省略なく導きます。1級ではこの論述が問われます。

4.1 QHQ-H 分解

任意の θ\theta について、完全データの同時分布を条件付き分布で割って書き直します。確率の連鎖律から

p(xobs,xmisθ)=p(xobsθ)p(xmisxobs,θ)p(x_{\mathrm{obs}},x_{\mathrm{mis}}\mid\theta)=p(x_{\mathrm{obs}}\mid\theta)\,p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta)

なので、両辺の対数を取って観測尤度について解くと

logp(xobsθ)=logp(xobs,xmisθ)logp(xmisxobs,θ)\log p(x_{\mathrm{obs}}\mid\theta)=\log p(x_{\mathrm{obs}},x_{\mathrm{mis}}\mid\theta)-\log p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta)

要するに「観測尤度 = 完全データ尤度 −(欠測部分の条件付き尤度)」。この式の左辺は xmisx_{\mathrm{mis}} を含みません。そこで両辺を p(xmisxobs,θ(t))p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta^{(t)}) で期待値化します(左辺は xmisx_{\mathrm{mis}} を含まないので期待値を取っても変わらない点がポイント)。

logp(xobsθ)=E ⁣[logp(xobs,xmisθ)|xobs,θ(t)]=Q(θθ(t))E ⁣[logp(xmisxobs,θ)|xobs,θ(t)]=H(θθ(t))\log p(x_{\mathrm{obs}}\mid\theta) =\underbrace{E\!\left[\log p(x_{\mathrm{obs}},x_{\mathrm{mis}}\mid\theta)\,\middle|\,x_{\mathrm{obs}},\theta^{(t)}\right]}_{=\,Q(\theta\mid\theta^{(t)})} -\underbrace{E\!\left[\log p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta)\,\middle|\,x_{\mathrm{obs}},\theta^{(t)}\right]}_{=\,H(\theta\mid\theta^{(t)})}

これが QHQ-H 分解です。

  logp(xobsθ)=Q(θθ(t))H(θθ(t))  \boxed{\;\log p(x_{\mathrm{obs}}\mid\theta)=Q(\theta\mid\theta^{(t)})-H(\theta\mid\theta^{(t)})\;}

要するに「観測対数尤度を、Eステップで作る QQ と、もう一つの項 HH の差にきれいに分けられる」。QQ はすでに見た期待完全データ対数尤度、HH は条件付き尤度の期待値です。

4.2 増分を2つに分ける

θ=θ(t+1)\theta=\theta^{(t+1)}θ=θ(t)\theta=\theta^{(t)} で上式を書き、引き算します。

logp(xobsθ(t+1))logp(xobsθ(t))=[Q(θ(t+1)θ(t))Q(θ(t)θ(t))](I) Q の増分[H(θ(t+1)θ(t))H(θ(t)θ(t))](II) H の増分\log p(x_{\mathrm{obs}}\mid\theta^{(t+1)})-\log p(x_{\mathrm{obs}}\mid\theta^{(t)}) =\underbrace{\big[Q(\theta^{(t+1)}\mid\theta^{(t)})-Q(\theta^{(t)}\mid\theta^{(t)})\big]}_{(\text{I})\ Q\ \text{の増分}} -\underbrace{\big[H(\theta^{(t+1)}\mid\theta^{(t)})-H(\theta^{(t)}\mid\theta^{(t)})\big]}_{(\text{II})\ H\ \text{の増分}}

観測尤度が増える(0\ge0)ことを言うには、(I) 0\ge0 かつ (II) 0\le0 を示せば十分です(増分 = (I) − (II) なので、(I) が非負で (II) が非正なら全体は非負)。

4.3 第1項 (I):Mステップから 0\ge0

これはMステップの定義そのものから出ます。Mステップは θ(t+1)=argmaxθQ(θθ(t))\theta^{(t+1)}=\arg\max_\theta Q(\theta\mid\theta^{(t)}) なので、θ(t+1)\theta^{(t+1)}θ(t)\theta^{(t)} も含むあらゆる候補の中で QQ を最大にします。したがって

Q(θ(t+1)θ(t))  Q(θ(t)θ(t))(I)  0Q(\theta^{(t+1)}\mid\theta^{(t)})\ \ge\ Q(\theta^{(t)}\mid\theta^{(t)}) \quad\Longrightarrow\quad (\text{I})\ \ge\ 0

要するに「Mステップは QQ を最大化するのだから、少なくとも今の値より QQ は下がらない」。実は厳密な最大化でなくても、QQ を「増やしさえすれば」(Q(θ(t+1)θ(t))Q(θ(t)θ(t))Q(\theta^{(t+1)}\mid\theta^{(t)})\ge Q(\theta^{(t)}\mid\theta^{(t)}) であれば)(I)0\ge0 は保たれます。これが一般化EM(GEM)が成立する根拠です。

4.4 第2項 (II):Gibbsの不等式から 0\le0

ここが証明の山です。HHθ=θ(t)\theta=\theta^{(t)} で最大になること、すなわち任意の θ\theta に対して H(θθ(t))H(θ(t)θ(t))H(\theta\mid\theta^{(t)})\le H(\theta^{(t)}\mid\theta^{(t)}) を示します。これが言えれば θ=θ(t+1)\theta=\theta^{(t+1)} を入れて (II)0\le0 が出ます。

差を取って書き下します。表記簡略のため pθ:=p(xmisxobs,θ)p_\theta:=p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta)pt:=p(xmisxobs,θ(t))p_t:=p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta^{(t)}) とおくと、期待値は ptp_t で取るので

H(θθ(t))H(θ(t)θ(t))=E ⁣[logpθlogpt|xobs,θ(t)]=ptlogpθptdxmisH(\theta\mid\theta^{(t)})-H(\theta^{(t)}\mid\theta^{(t)}) =E\!\left[\log p_\theta-\log p_t\,\middle|\,x_{\mathrm{obs}},\theta^{(t)}\right] =\int p_t\,\log\frac{p_\theta}{p_t}\,dx_{\mathrm{mis}}

要するに「HH の増分は、ptp_t から見た比 pθ/ptp_\theta/p_t の対数の平均」。ここでJensenの不等式を使います。log\log は上に凸(concave)なので、確率測度 ptp_t に関する期待値 Ept[]E_{p_t}[\cdot] について E[log()]logE[]E[\log(\cdot)]\le\log E[\cdot]。これを適用すると

ptlogpθptdxmis  log ⁣(ptpθptdxmis)=log ⁣(pθdxmis)=log1=0\int p_t\,\log\frac{p_\theta}{p_t}\,dx_{\mathrm{mis}} \ \le\ \log\!\left(\int p_t\cdot\frac{p_\theta}{p_t}\,dx_{\mathrm{mis}}\right) =\log\!\left(\int p_\theta\,dx_{\mathrm{mis}}\right) =\log 1=0

最後は pθ=p(xmisxobs,θ)p_\theta=p(x_{\mathrm{mis}}\mid x_{\mathrm{obs}},\theta) が(xmisx_{\mathrm{mis}} についての)確率分布なので全積分が1であることを使いました。よって

H(θθ(t))H(θ(t)θ(t))  0(任意の θH(\theta\mid\theta^{(t)})-H(\theta^{(t)}\mid\theta^{(t)})\ \le\ 0\quad\text{(任意の }\theta\text{)}

要するに「HH は現在値 θ(t)\theta^{(t)} で最大。だからどこへ動いても HH は増えない」。θ=θ(t+1)\theta=\theta^{(t+1)} を入れれば (II)0\le0

なお ptlogptpθdxmis=KL(ptpθ)0\displaystyle\int p_t\log\frac{p_t}{p_\theta}\,dx_{\mathrm{mis}}=\mathrm{KL}(p_t\,\|\,p_\theta)\ge0 は2分布間の Kullback–Leibler ダイバージェンスで、これが非負であること(Gibbsの不等式)と同じことを言っています。符号を反転すれば上の (II)0\le0 そのものです。

4.5 結論:観測尤度は単調非減少

(I)0\ge0 かつ (II)0\le0 なので、

logp(xobsθ(t+1))logp(xobsθ(t))=(I)(II)  0\log p(x_{\mathrm{obs}}\mid\theta^{(t+1)})-\log p(x_{\mathrm{obs}}\mid\theta^{(t)})=(\text{I})-(\text{II})\ \ge\ 0

すなわち EMの各反復で観測データの対数尤度は決して減らない(単調非減少)。証明完了。

4.6 何が言えて、何が言えないか(収束の限界)

単調非減少なので、対数尤度が上に有界(普通そうである)なら数列 {logp(xobsθ(t))}\{\log p(x_{\mathrm{obs}}\mid\theta^{(t)})\} は収束します。ただし——

要するに「EMは尤度を着実に登るが、必ずしも一番高い山には登らない」。これがEMの本質的限界で、実務では初期値を変えて複数回走らせ、最も尤度が高い解を採るのが定石です。


5. 具体例:2成分正規混合分布のEM

EMの典型例です。「各データがどちらの成分から生まれたか」が見えない(=潜在変数)状況で、混合分布のパラメータを推定します。混合分布の構成要素である正規分布は 正規分布(標準正規・標準化) を参照。

5.1 モデルと潜在変数

データ x1,,xnx_1,\dots,x_n が、2つの正規分布の混合から生成されたとします。混合比 π\pi(成分1の確率)と各成分のパラメータ (μ1,σ12),(μ2,σ22)(\mu_1,\sigma_1^2),(\mu_2,\sigma_2^2) で、1点の密度は

p(xiθ)=πN(xiμ1,σ12)+(1π)N(xiμ2,σ22)p(x_i\mid\theta)=\pi\,\mathcal N(x_i\mid\mu_1,\sigma_1^2)+(1-\pi)\,\mathcal N(x_i\mid\mu_2,\sigma_2^2)

要するに「各点は確率 π\pi で成分1から、1π1-\pi で成分2から出た、という2つの山の重ね合わせ」。観測対数尤度 ilog[πN()+(1π)N()]\sum_i\log\big[\pi\mathcal N(\cdot)+(1-\pi)\mathcal N(\cdot)\big]log\log の中に和があるため直接微分しても閉形式で解けません。これがEMの出番です。

ここで潜在変数 zi{1,2}z_i\in\{1,2\} を「点 ii がどちらの成分から出たか」とします。ziz_i が見えていれば(完全データなら)各成分ごとに点を仕分けて普通に正規分布の最尤推定をするだけ——という単純さが鍵です。潜在変数の条件付き期待値という発想は 同時分布・周辺分布・条件付き分布 の条件付き分布そのものです。

5.2 Eステップ:負担率(responsibility)

潜在変数 ziz_i は見えないので、現在値 θ(t)\theta^{(t)} のもとで「点 ii が成分1由来である事後確率」を計算します。これを負担率 γi\gamma_i と呼びます。ベイズの定理から

γi=P(zi=1xi,θ(t))=π(t)N(xiμ1(t),σ12(t))π(t)N(xiμ1(t),σ12(t))+(1π(t))N(xiμ2(t),σ22(t))\gamma_i=P(z_i=1\mid x_i,\theta^{(t)}) =\frac{\pi^{(t)}\,\mathcal N(x_i\mid\mu_1^{(t)},\sigma_1^{2(t)})} {\pi^{(t)}\,\mathcal N(x_i\mid\mu_1^{(t)},\sigma_1^{2(t)})+(1-\pi^{(t)})\,\mathcal N(x_i\mid\mu_2^{(t)},\sigma_2^{2(t)})}

要するに「今のパラメータで見て、この点はどれくらい成分1っぽいか(0〜1のソフトな所属度)」。完全データ対数尤度の期待値 QQ を取ると、ziz_i の指示関数の期待値がちょうどこの γi\gamma_i になります(これがEステップの実体)。

5.3 Mステップ:負担率で重みづけた最尤推定

QQθ\theta で最大化すると、各パラメータは負担率を重みとした更新式になります(有効標本数 N1=iγiN_1=\sum_i\gamma_iN2=i(1γi)N_2=\sum_i(1-\gamma_i))。

μ1(t+1)=iγixiiγi,σ12(t+1)=iγi(xiμ1(t+1))2iγi,π(t+1)=iγin=N1n\mu_1^{(t+1)}=\frac{\sum_i \gamma_i\,x_i}{\sum_i \gamma_i},\qquad \sigma_1^{2(t+1)}=\frac{\sum_i \gamma_i\,(x_i-\mu_1^{(t+1)})^2}{\sum_i \gamma_i},\qquad \pi^{(t+1)}=\frac{\sum_i \gamma_i}{n}=\frac{N_1}{n}

(成分2は γi\gamma_i1γi1-\gamma_i に置き換えた同形の式。)要するに「普通の標本平均・標本分散・割合の式そのもの。ただし各点を『所属度 γi\gamma_i で重みづけ』しているだけ」。所属が確実(γi{0,1}\gamma_i\in\{0,1\})なら、これは点を2群に振り分けてそれぞれ最尤推定するのと同じ式に戻ります。Eで γi\gamma_i を更新→Mでパラメータを更新→またEへ、と反復し、対数尤度の増分が小さくなったら停止します。

この「ソフトな所属度で重みづけ」という構造は、潜在変数を介して観測を説明する潜在変数モデルに共通します。因子分析(因子分析)も観測変数の背後に潜在因子を置くモデルで、そのパラメータ推定にEMが使えます。


6. 引っかけ・頻出論点


試験での問われ方(1級)

1級(統計数理)でのEMは、出題頻度こそ高くないものの、出れば記述式で深く問われます(範囲・配点は改訂されうるため要最新確認)。典型的な問われ方:

前提として最尤法(最尤法・モーメント法(推定量の作り方と最尤推定量の漸近論))の対数尤度の最大化、条件付き期待値(同時分布・周辺分布・条件付き分布)が完全に使えることが要求されます。


よくある疑問(Q&A)

Q1. EMアルゴリズムは結局、観測データの尤度を最大化しているのですか? それとも別の量?

最終的に最大化したいのは観測データの対数尤度 logp(xobsθ)\log p(x_{\mathrm{obs}}\mid\theta) です。ただしそれを直接いじるのではなく、各反復では完全データ対数尤度の期待値 QQ を最大化します。QHQ-H 分解により「QQ を増やせば観測尤度も減らない」ことが保証されるので、易しい QQ の最大化を繰り返すことで、難しい観測尤度を間接的に登れる——という二段構えです。だから「EMはEステップで作った QQ を最大化する手続きだが、結果として観測尤度が単調に増える」が正確な答えです。

Q2. MAR なのに完全ケース分析(欠測行を捨てる)をしてはいけないのですか?

避けるべきです。MAR で「欠測を無視できる」のは尤度を正しく書いて最尤・ベイズ・EM・多重代入で解く場合の話です。完全ケース分析が無バイアスでいられるのは、より強い MCAR のときだけ。MAR では残った完全ケースが(観測変数を通じて)偏っているため、単純に捨てるとバイアスが出ます。MAR のときは欠測行を捨てずに、EMや多重代入で観測情報を活かして推定するのが正しい対処です。

Q3. 平均値代入はなぜダメなのですか? 平均で埋めれば平均は変わらないのでよさそうに見えます。

平均そのものは(MCARなら)大きく狂いませんが、ばらつきと関係性が壊れます。欠測をすべて平均で埋めると、埋めた点は偏差ゼロなので分散への寄与がゼロ。なのに標本数は増えるので、分散の推定値が真の値の m1n1\tfrac{m-1}{n-1} 倍に系統的に縮みます(mm=観測数、nn=全数)。分散が小さく見えれば標準誤差も小さく出て、信頼区間は狭すぎ・検定は甘すぎになります。さらに他変数との相関も平均一点に潰れて弱まります。「点推定が合えばよい」のでなく、ばらつきの推定が壊れるのが本質的な害です。

Q4. EMは必ず最尤推定値(大域最適)にたどり着きますか?

いいえ。EMが保証するのは「反復ごとに観測対数尤度が減らない」ことだけで、収束先は**停留点(局所最適または鞍点)**です。尤度が多峰の場合(正規混合が典型)、初期値が悪いと低い山に登って止まります。対策は、異なる初期値から複数回EMを走らせ、最も対数尤度が高くなった解を採用すること。k-means の結果を初期値にする、混合比や中心を散らして初期化する、などが実務でよく使われます。

Q5. Eステップの「負担率」は、各点をどちらかの成分に割り当てているのですか?

割り当て(0か1のハードなラベル付け)ではありません。負担率 γi\gamma_i は「点 ii が成分1由来である事後確率(0〜1の連続値)」です。例えば γi=0.7\gamma_i=0.7 なら「7割がた成分1っぽい」というソフトな所属度で、Mステップではこの 0.70.7 を重みとして各成分のパラメータ更新に寄与します。0か1に丸めてから更新するのは ハードEM(k-means に近い近似)で、EM本来の更新ではありません。ソフトに重みづけるからこそ、QHQ-H 分解による単調増加の理論保証が成り立ちます。

Q6. EステップとMステップ、どちらが計算的に重いことが多いですか?

モデルによりますが、典型的にはEステップは閉形式(負担率や条件付き期待値を式に代入するだけ)で軽く、Mステップは最大化問題です。完全データ尤度が指数型分布族なら、Mステップも閉形式(重みつき平均・分散など)で軽く解けます。Mステップが閉形式で解けない場合は、数値最適化を1回だけ回す・QQ を増やしさえすればよいという**一般化EM(GEM)**で済ませることもできます(それでも単調増加性は保たれる)。


まとめ


関連ノート