EMアルゴリズムとともだちになろう

イメージをつかむ 数式で理解する 1/34 Rでやってみる 2022年8月15日 (資料公開日：2022年8月16日) ベイズ統計学勉強会 夏’22@福島 EMアルゴリズムと ともだちになろう 武藤 拓之 (Hiroyuki Muto) 京都大学人と社会の未来研究院

イメージをつかむ 数式で理解する 2/34 Rでやってみる この発表のモチベーション ⚫ そろそろ真面目にEMアルゴリズムについて勉強したほうがよい気がした。 ⚫ そんなとき神サイト(↓)を見つけた。全人類読むべし。 ⚫ この内容をさらに咀嚼して直感的に説明できるようになりたい。 ⚫ 自分にとっての分かりやすさを最優先して資料を作成したので， 一般的な記法とは異なるかもしれません。間違いもあるかもです。 https://academ-aid.com/ml/em

• https://academ-aid.com/ml/em

イメージをつかむ 数式で理解する 3/34 Rでやってみる EMアルゴリズムとは ◼ 最尤推定やMAP推定に使えるアルゴリズム。 ◼ ExpectationとMaximizationのステップを交互に実行するのでEM。 (なかなかエモい。emだけに。) ◼ 混合分布モデルとか階層モデルみたいなめんどくさい系のモデルで使える。 ⚫ 対数尤度関数にlog σ 𝑓(・)が出てきちゃうやつ。 e.g., log 𝑝 𝑦1 , 𝑧1,1 |𝜽 + 𝑝 𝑦1 , 𝑧1,2 |𝜽 + log 𝑝 𝑦2 , 𝑧2,1 |𝜽 + 𝑝 𝑦2 , 𝑧2,2 |𝜽 ⚫ logの中に足し算が入ると計算がしんどい。 (足し算じゃなくて掛け算なら分割してそれぞれ微分するだけなのでラクチン。) ◼ 変分ベイズを理解するためには知っておいた方がいい。 (申し訳程度のベイズ要素)

イメージをつかむ 数式で理解する Rでやってみる 4/34 イメージをつかむ 余談： 心的イメージの研究者としては「イメージ」をこういう意味ではあんまり使いたくない。

イメージをつかむ 数式で理解する 5/34 Rでやってみる 今回扱う具体例 母集団分布 𝛼 = 0.3 (未知) 標本 サンプリング N = 1000 推定 𝛼ො = ? 𝑍𝑖 ~ Bernoulli(𝛼) Normal 𝜇1 = 0, 𝜎 = 1 if 𝑍𝑖 = 0 𝑌𝑖 ~ ቊ Normal 𝜇2 = 5, 𝜎 = 1 if 𝑍𝑖 = 1 ◼ 混合分布から生成されたデータを使って， 混合分布のパラメータを推定する。 ◼ 簡単のため，𝑍𝑖 = 1となる確率𝛼だけが未知であるとする。 ◼ このモデルは非常に簡単なのでEMアルゴリズムを使うまでもなく解ける。 (それゆえに答え合わせがしやすくて勉強にはよい。)

6/34 Rでやってみる イメージをつかむ 数式で理解する モデルと対数尤度関数 【確率モデル】 𝑍𝑖 ~ Bernoulli(𝛼) Normal 𝜇1 = 0, 𝜎 = 1 if 𝑍𝑖 = 0 𝑌𝑖 ~ ቊ Normal 𝜇2 = 5, 𝜎 = 1 if 𝑍𝑖 = 1 𝛼 = 0.3 (未知) 噛み砕くと：コインを投げて裏が出たら平均0，表が出たら平均5の 正規分布から観測値が出てくるとき，表が出る確率を推定したい。 ◼ 対数尤度関数 ⚫ log 𝐿 𝛼ො 𝒚 𝑛 = ෍ log 𝑝(𝑦𝑖 |𝛼) ො 𝑖=1 𝑛 2 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 条件付き独立の仮定より 同時分布を周辺化したものとみる 𝑗=1 ⚫ 間に𝑍𝑖 が入っているせいで，計算しにくいlog σ 𝑓(・)が出てきてしまう。 ⚫ なんとかしてこの関数を最大にするような𝛼を求めたい。 ො

イメージをつかむ 数式で理解する 7/34 Rでやってみる データのイメージ 𝑖 𝑦𝑖 𝑧𝑖 1 ＋5.26 1 2 ＋6.18 1 3 ＋1.36 0 4 ＋0.19 0 5 －1.06 0 6 ＋2.98 0 ⋮ ⋮ ⋮ ↑ 観測される ↑ 実際には 観測されない

イメージをつかむ 数式で理解する 8/34 Rでやってみる 最尤推定の基本 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 𝑗=1 𝛼ො = 0.3 1. 対数尤度関数を数式で表す。 2. 適当な方法で，対数尤度を最大にするパラメータを求める。 (解析的 or 数値的に解く) 3. おしまい。

イメージをつかむ 数式で理解する 9/34 Rでやってみる EMアルゴリズム log σ 𝑓(・)があって 超だるい…… 1. 初期値の設定 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 1. 適当な初期値を決める。(e.g., 𝛼ො0 = 0.8) 𝑗=1

イメージをつかむ 数式で理解する 10/34 Rでやってみる EMアルゴリズム log σ 𝑓(・)があって 超だるい…… 2. Eステップ 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො0 = 0.8 𝑛 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.8 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 𝑖=1 𝑗=1 赤線は黒線を 超えない！ log σ 𝑓(・)がない！ 簡単！ 1. 適当な初期値を決める。(e.g., 𝛼ො0 = 0.8) 2. 真の対数尤度関数を超えないがこの点だけで交わってくれる， 「仮の」対数尤度関数 (i.e., 下界の1つ) を簡単な式で表せる。 (実際には定数部分の計算を省略する。) (log σ 𝑓(・)が出てこない形)

イメージをつかむ 数式で理解する 11/34 Rでやってみる EMアルゴリズム 3. Mステップ 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො0 = 0.8 𝑛 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.8 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 𝑖=1 𝑗=1 赤線は黒線を 超えない！ 𝛼ො1 = 0.330 1. 適当な初期値を決める。(e.g., 𝛼ො0 = 0.8) 2. 真の対数尤度関数を超えないがこの点だけで交わってくれる， 「仮の」対数尤度関数 (i.e., 下界の1つ) を簡単な式で表せる。 3. この「仮の」対数尤度関数を最大にするパラメータ𝛼ො1 を求める。

イメージをつかむ 数式で理解する 12/34 Rでやってみる EMアルゴリズム 4. これを繰り返すと真の対数尤度関数に近づいていく！ 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො0 = 0.8 𝑛 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.8 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො 0 = 0.330 𝑛 赤線と青線は黒線を 超えない！ 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.330 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 𝑖=1 𝑗=1 𝛼ො2 = 0.317 1. 適当な初期値を決める。(e.g., 𝛼ො0 = 0.8) 2. 真の対数尤度関数を超えないがこの点だけで交わってくれる， 「仮の」対数尤度関数 (i.e., 下界の1つ) を簡単な式で表せる。 3. この「仮の」対数尤度関数を最大にするパラメータ𝛼ො1 を求める。 収束まで 繰り返す

イメージをつかむ 数式で理解する 13/34 Rでやってみる EMアルゴリズム 4. これを繰り返すと真の対数尤度関数に近づいていく！ 𝑛 2 log 𝐿 𝛼ො 𝒚 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼) ො 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො0 = 0.8 𝑛 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.8 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 赤線と青線は黒線を 超えない！ 対数尤度関数そのものを最大化する代わりに， 𝛼ො = 0.317 対数尤度関数の「下界」となる関数を 何度も最大化していくのがEMアルゴリズム！ 𝑖=1 𝑗=1 ℒ 𝛼|𝒚, ො 𝛼ො 0 = 0.330 𝑛 2 = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 0.330 log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝛼ො + const 𝑖=1 𝑗=1 2 1. 適当な初期値を決める。(e.g., 𝛼ො0 = 0.8) 2. 真の対数尤度関数を超えないがこの点だけで交わってくれる， 「仮の」対数尤度関数 (i.e., 下界の1つ) を簡単な式で表せる。 3. この「仮の」対数尤度関数を最大にするパラメータ𝛼ො1 を求める。 収束まで 繰り返す

イメージをつかむ 数式で理解する Rでやってみる 14/34 数式で理解する

イメージをつかむ 数式で理解する 15/34 Rでやってみる 下界を求めるための道具 ◼ Jensenの不等式 ⚫ 𝑓(𝑥)は凹関数 (i.e., 上に凸) で，𝜆𝑗 は正の実数かつσ∞ 𝑗=1 𝜆𝑗 = 1であるとする。 ⚫ このとき，任意の𝑥1 , 𝑥2 , … , 𝑥𝑚 および上記条件を満たす任意の𝜆1 , 𝜆2 , … , 𝜆𝑚 について， 𝑚 𝑓 ෍ 𝜆𝑗 𝑥𝑗 𝑗=1 𝑚 ≥ ෍ 𝜆𝑗 𝑓 𝑥𝑗 𝑗=1 が成り立つ。 ⚫ 対数関数の場合は， 𝑚 𝑚 log ෍ 𝜆𝑗 𝑥𝑗 ≥ ෍ 𝜆𝑗 log 𝑥𝑗 𝑗=1 𝑗=1 となる。 ⚫ この𝜆𝑗 を確率測度と考えると後の説明が分かりやすい。 これを使えばlogの中のΣを外に出せそう！

イメージをつかむ 数式で理解する 16/34 Rでやってみる 最大対数尤度とその下界を計算する パラメータの推定値のベクトル ෡ = (ෞ e.g., 𝜽 𝜇1 , 𝜇 ෞ2 , 𝜎, ො 𝛼) ො ◼ 最大対数尤度とその下界 𝑛 ෡𝒚 ⚫ log 𝐿 𝜽 ෡ = ෍ log 𝑝(𝑦𝑖 |𝜽) 𝑖=1 𝑛 2 ←条件付き独立の仮定より 𝑧𝑖𝑗 が取りうる値の数 (分かりやすさのために具体的に2とした。) ※連続変数の場合は総和Σを定積分∫にすればOK ෡ = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝑖=1 𝑛 テクい部分 ←同時分布を周辺化したものとみる 𝑗=1 2 ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) = ෍ log ෍ 𝜆𝑖𝑗 𝜆𝑖𝑗 ←Jensenの不等式を使う準備 𝑖=1 𝑗=1 𝑛 2 ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) ≥ ෍ ෍ 𝜆𝑖𝑗 log 𝜆𝑖𝑗 𝑖=1 𝑗=1 ෡ 𝜦 = ℒ 𝜽|𝒚, ※ 𝜦は各𝜆𝑖𝑗 を要素に持つ𝑛 × 2の行列 ←Jensenの不等式より Jensenの不等式： 𝑚 𝑚 log ෍ 𝜆𝑗 𝑥𝑗 ≥ ෍ 𝜆𝑗 log 𝑥𝑗 𝑗=1 𝑗=1 ෡ 𝜦 を最大化するような𝜦を見つければ この「仮の」対数尤度関数ℒ 𝜽|𝒚, ෡ 𝒚 と一致する。 「本物の」対数尤度関数log 𝐿 𝜽

イメージをつかむ 数式で理解する 17/34 Rでやってみる 対数尤度－下界≥0 ◼ 要するに差が0になればよいので…… 𝑛 ෡ 𝒚 − ℒ 𝜽|𝒚, ෡ 𝑸 ⚫ log 𝐿 𝜽 = 1を掛ける→ 𝑗=1 𝑛 𝑛 𝑖=1 𝑗=1 2 𝑖=1 𝑗=1 ෡ − 𝜆𝑖𝑗 log = ෍ ෍ 𝜆𝑖𝑗 log 𝑝 𝑦𝑖 𝜽 = ෍ ෍ 𝜆𝑖𝑗 = ෍ ෍ 𝜆𝑖𝑗 𝑖=1 𝑗=1 𝑛 2 対数の性質を利用→ ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝜆𝑖𝑗 ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝜆𝑖𝑗 ෡ − log log 𝑝 𝑦𝑖 𝜽 ෡ 𝑦𝑖 𝜽 ෡ 𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽)𝑝 𝜆𝑖𝑗 ෡ 𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽) = − ෍ ෍ 𝜆𝑖𝑗 log 𝜆𝑖𝑗 𝑖=1 𝑗=1 ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝜆𝑖𝑗 ෡ − log log 𝑝 𝑦𝑖 𝜽 𝑖=1 𝑗=1 𝑛 2 条件付き確率の定義式より→ 2 ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) − ෍ ෍ 𝜆𝑖𝑗 log 𝜆𝑖𝑗 𝑖=1 𝑗=1 𝑛 2 𝜆𝑖𝑗 でくくる→ ෡ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝜆𝑖𝑗 ෡ − ෍ ෍ 𝜆𝑖𝑗 log = ෍ ෍ 𝜆𝑖𝑗 log 𝑝(𝑦𝑖 |𝜽) 𝑖=1 𝑗=1 𝑛 2 シグマを統合→ 𝑖=1 𝑗=1 2 ෡ ෍ 𝜆𝑖𝑗 = ෍ log 𝑝 𝑦𝑖 𝜽 𝑖=1 𝑛 2 足し算の順番を変える→ 2 ෡ − ෍ ෍ 𝜆𝑖𝑗 log = ෍ log 𝑝(𝑦𝑖 |𝜽) 𝑖=1 𝑛 σ2𝑗=1 𝜆𝑖𝑗 𝑛 KL情報量！

イメージをつかむ 数式で理解する 18/34 Rでやってみる ここまでの要点のまとめと次のポイント ◼ 対数尤度関数とその下界 𝑛 ෡𝒚 ⚫ 対数尤度関数：log 𝐿 𝜽 𝑛 ෡ 𝜦 ⚫ その下界：ℒ 𝜽|𝒚, 2 2 ෡ = ෍ log ෍ 𝑝(𝑦𝑖 , |𝑧𝑖𝑗 , 𝜽)𝑝(𝑧 ෡ ෡ ෡ = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝜽) 𝑖𝑗 |𝜽) ≥ ℒ 𝜽|𝒚, 𝜦 𝑖=1 2 𝑗=1 ෡ 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 = ෍ ෍ 𝜆𝑖𝑗 log 𝜆𝑖𝑗 𝑖=1 𝑗=1 ෡ 𝒚 − ℒ 𝜽|𝒚, ෡ 𝜦 ⚫ 差：log 𝐿 𝜽 𝑛 𝑛 𝑖=1 𝑗=1 ←これを最大化すれば最大対数尤度関数と一致する。 2 ෡ 𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽) ෡ ≥0 = − ෍ ෍ 𝜆𝑖𝑗 log = KL 𝜆𝑖𝑗 ||𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽) 𝜆𝑖𝑗 𝑖=1 𝑗=1 ◼ ポイント ෡ 𝜦 は， 𝑛 × 2の行列𝜦を適当に与えてやると「仮の」尤度関数になる。 ⚫ 下界ℒ 𝜽|𝒚, (各要素は正の実数で行和は1) → 関数の形を決めるパラメータとして𝜦を考えると分かりやすい。 ෡ (0) を適当に決めてやればその条件下で理想的な𝜦が決まる。 ⚫ パラメータの初期値𝜽 ෡ (0) から𝜦を決める方法を考える。 次に，初期値𝜽

イメージをつかむ 数式で理解する 19/34 Rでやってみる ෡ 0 からΛを決める方法 初期値𝜽 ◼ KL情報量が0になればよいから…… ⚫ 全く同じ2つの確率分布のKL情報量は0だから， ෡ ෡ = 0となる。 𝜆𝑖𝑗 = 𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽)のときKL 𝜆𝑖𝑗 ||𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽) ෡ ෡ 0 を使って， ⚫ 望ましい𝜆𝑖𝑗 の関数形は分かったが𝜽は未知なので，適当な初期値 𝜽 ෡ (0) 𝑝 𝑧𝑖𝑗 𝑝 𝑦 𝑧 , 𝜽 𝑖 𝑖𝑗 (0) ෡ (0) = 𝜆𝑖𝑗 = 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 ෡ (0) 𝑝 𝑧𝑖1 + 𝑝 𝑦𝑖 𝑧𝑖2 , 𝜽 ෡ (0) 𝑝 𝑧𝑖2 𝑝 𝑦𝑖 𝑧𝑖1 , 𝜽 (0) と𝜆𝑖𝑗 の初期値𝜆𝑖𝑗 を定めてやる。(計算のためにベイズ定理を利用) ⚫ 先の混合分布モデルの場合，尤度を使って次の式で𝜦(0) の値を決められる。 (0) 𝜆𝑖1 = (0) 𝜆𝑖2 = ෡ (0) 𝛼ො (0) 𝑝 𝑦𝑖 𝑧𝑖 = 0, 𝜽 ෡ (0) (1 − 𝛼ො (0) ) + 𝑝 𝑦𝑖 𝑧𝑖 = 1, 𝜽 ෡ 0 𝛼ො (0) 𝑝 𝑦𝑖 𝑧𝑖 = 0, 𝜽 ෡ (0) 𝛼ො (0) 𝑝 𝑦𝑖 𝑧𝑖 = 1, 𝜽 ෡ (0) (1 − 𝛼ො (0) ) + 𝑝 𝑦𝑖 𝑧𝑖 = 1, 𝜽 ෡ (0) 𝛼ො (0) 𝑝 𝑦𝑖 𝑧𝑖 = 0, 𝜽

20/34 Rでやってみる イメージをつかむ 数式で理解する Eステップ (Expectation) ◼ 用意した𝜦(0) で下界のひとつを計算すると…… 𝑛 2 ෡ 𝜦(0) = ෍ ෍ 𝜆(0) log ⚫ ℒ 𝜽|𝒚, 𝑖𝑗 ෡ 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 (0) 𝜆𝑖𝑗 𝑖=1 𝑗=1 𝑛 2 ෡ (0) = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 𝑖=1 𝑗=1 𝑛 2 ෡ 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 log ෡ (0) 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 ←Jensenの不等式から得た式 (0) ←𝜆𝑖𝑗 に先の値を代入 ෡ (0) log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 ෡ − log 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 ෡ (0) = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 ←対数の性質より 𝑖=1 𝑗=1 𝑛 2 ෡ (0) log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 ෡ + const ←𝜽を含まない定数部分を分離 ෡ = ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 𝑖=1 𝑗=1 𝑧𝑖𝑗 についての対数尤度の期待値の形になっている。 ෡ 𝜽 ෡ (0) ) + const = 𝒬(𝜽|𝒚, ෡ ←𝜽の関数部分 (期待対数尤度) を𝒬と置く。 ෡ 𝜽 ෡ (0) )を求める。 Eステップ：対数尤度の期待値すなわち𝒬(𝜽|𝒚,

イメージをつかむ 数式で理解する 21/34 Rでやってみる Mステップ (Maximization) ◼ この下界を最大にするパラメータを推定する ෡ 𝜦(0) ෡ (1) = arg max ℒ 𝜽|𝒚, ⚫ 𝜽 ෡ 𝜽 ෡ 𝜽 ෡ (0) ) = arg max𝒬(𝜽|𝒚, ෡ 𝜽 𝑛 ෡ ←「仮の」対数尤度関数を最大化する𝜽を採用 ෡ 𝜽 ෡ (0) )を最大化 ←定数部分は影響しないので𝒬(𝜽|𝒚, 2 ෡ (0) log 𝑝 𝑦𝑖 , 𝑧𝑖𝑗 𝜽 ෡ = arg max ෍ ෍ 𝑝 𝑧𝑖𝑗 𝑦𝑖 , 𝜽 ෡ 𝜽 ←ちゃんと書くとこうなる 𝑖=1 𝑗=1 ⚫ この式は簡単なので，偏微分を使って最大化できる。 (よくある最適化問題) ෡ (1) は，𝜽 ෡ (0) のときよりも対数尤度を高くしてくれることが保証される。 ⚫ この𝜽 ෡ (1) を計算する。 Mステップ：期待対数尤度を最大化する𝜽 ෡ (1) を使ってまた「Eステップ→Mステップ」を実行する。これを気の済むまで繰り返す。 → この𝜽

イメージをつかむ 数式で理解する 22/34 Rでやってみる まとめ 1. 対数尤度関数を求めたいが，log σ 𝑓(・)が出てきて計算がしんどい。 2. そこで，対数尤度関数よりも下に来る下界 (「仮の」対数尤度関数) に注目する。 ෡ 𝜦 を表せる。 3. Jensenの不等式を使うと簡単な式で下界ℒ 𝜽|𝒚, 4. 無数にある下界のうち1つを選ぶには𝜦を決める必要がある。 5. 対数尤度関数と下界の差がKL情報量になることに注目すると， ෡ とりあえず𝜆𝑖𝑗 = 𝑝(𝑧𝑖𝑗 |𝑦𝑖 , 𝜽)という関数形にする必要があることが分かる。 ෡ ෡ (0) を使って𝜆(0) をひとつに決める。 6. まだ正確な𝜽は分からないので，適当な初期値 𝜽 𝑖𝑗 ෡ 𝜦(0) を得ることができたので，こいつの最大化を目指す。 7. こうしてひとつの下界ℒ 𝜽|𝒚, (こいつの最大値は最大対数尤度を絶対に超えないから安心。) ෡ 𝜽 ෡ (0) )＋定数になる。 8. 計算しやすいように式を変形すると，期待対数尤度𝒬(𝜽|𝒚, (＝Eステップ) ෡ 𝜽 ෡ (0) )を最大にする𝜽 ෡ (1) を計算する。これは𝜽 ෡ (0) のときよりも対数尤度を高くしてくれる。 9. この𝒬(𝜽, (＝Mステップ) ෡ (𝑘) が最尤推定値に近づいていく。 10. EステップとMステップを繰り返せば𝜽

イメージをつかむ 数式で理解する Rでやってみる Rでやってみる 23/34

24/34 Rでやってみる イメージをつかむ 数式で理解する 設定 【確率モデル】 𝑍𝑖 ~ Bernoulli(𝛼) Normal 𝜇1 , 𝜎 if 𝑍𝑖 = 0 𝑌𝑖 ~ ቊ Normal 𝜇2 , 𝜎 if 𝑍𝑖 = 1 未知の母数 (推定対象) • 𝛼 = 0.3 • 𝜇1 = 0.0 • 𝜇2 = 5.0 • 𝜎 = 1.0 ◼ データと推定対象 ⚫ 上記モデルから1000個の乱数を生成 (要するにN = 1000)。 ⚫ このデータから4つのパラメータを最尤推定してみる。 ◼ 対数尤度関数 1000 ෡𝒚 ⚫ log 𝐿 𝜽 = ෍ log 𝑝(𝑦𝑖 |𝛼, 𝜇1 , 𝜇2 , 𝜎) 𝑖=1 1000 2 = ෍ log ෍ 𝑝(𝑦𝑖 , 𝑧𝑖𝑗 |𝛼, 𝜇1 , 𝜇2 , 𝜎) 𝑖=1 条件付き独立の仮定より 𝑗=1 同時分布を周辺化したものとみる

イメージをつかむ 数式で理解する Rでやってみる データを生成するコード 25/34

イメージをつかむ 数式で理解する Rでやってみる ෡ 𝜽 ෡ (0) )を定義する 関数𝒬(𝜽|𝒚, 26/34

イメージをつかむ 数式で理解する Rでやってみる 定数部分(const)を定義する ※本当は不要だけど 確認用として作った。 (定数なしで5m50s， 定数ありで8m41sかかった。) 27/34

イメージをつかむ 数式で理解する 28/34 Rでやってみる 目的関数を定義する パラメータの範囲を制約する ために，exp()とpnorm() を使っている。

イメージをつかむ 数式で理解する Rでやってみる EMアルゴリズムの準備 29/34

イメージをつかむ 数式で理解する 30/34 Rでやってみる EMアルゴリズムの実行 本当は偏微分して解くほうがよいが 今回は省略してoptim()で数値的に解く。

イメージをつかむ 数式で理解する 31/34 Rでやってみる 結果の可視化 真の母数 • 𝛼 = 0.3 • 𝜇1 = 0.0 • 𝜇2 = 5.0 • 𝜎 = 1.0 推定値 • 𝛼ො (50) = 0.316 • 𝜇 ෞ1 (50) = 0.006 • 𝜇 ෞ2 (50) = 5.040 • 𝜎ො (50) = 1.040 最終的な評価値 ෡ 50 𝒚, 𝜽 ෡ 49 = −2083.837 • 𝒬 𝜽 ෡ 50 |𝒚, 𝜦(49) = −2062.462 • ℒ 𝜽 ※数値解を直接求めたときの 最大対数尤度も−2062.462

イメージをつかむ 数式で理解する Rでやってみる おまけ：直接数値的に解くコード 32/34

イメージをつかむ 数式で理解する 33/34 Rでやってみる 留意点 ◼ ラベルスイッチングに注意。 → 順序制約は必須 ◼ 混合分布モデルは局所解に陥りやすい。 → 初期値をいろいろ変えないと怖い (今回もかなり試行錯誤したけど対症療法しかできなかった……) ◼ 比較的単純な混合分布モデルなら， MCMCで事後分布を近似したほうが実践上はよさそう。 (少なくとも局所解を大局解と思い込んだりする事態は避けられそう。) ◼ 今回はEMアルゴリズムの仕組みを理解することが主目的 だったので，効率性とか初期値の決め方とか諸々は 他の文献・資料を参照してください。

イメージをつかむ 数式で理解する Rでやってみる EMjoy! 34/34