$\gdef \vec#1{\boldsymbol{#1}} \\ \gdef \rank {\mathrm{rank}} \\ \gdef \det {\mathrm{det}} \\ \gdef \Bern {\mathrm{Bern}} \\ \gdef \Bin {\mathrm{Bin}} \\ \gdef \Mn {\mathrm{Mn}} \\ \gdef \Cov {\mathrm{Cov}} \\ \gdef \Po {\mathrm{Po}} \\ \gdef \HG {\mathrm{HG}} \\ \gdef \Geo {\mathrm{Geo}}\\ \gdef \N {\mathrm{N}} \\ \gdef \LN {\mathrm{LN}} \\ \gdef \U {\mathrm{U}} \\ \gdef \t {\mathrm{t}} \\ \gdef \F {\mathrm{F}} \\ \gdef \Exp {\mathrm{Exp}} \\ \gdef \Ga {\mathrm{Ga}} \\ \gdef \Be {\mathrm{Be}} \\ \gdef \NB {\mathrm{NB}} \\ \gdef \indep {\mathop{\perp\!\!\!\!\perp}} \\ \gdef \tr {\mathrm{tr}}$
今回はデータに欠損/潜在変数がある時に最尤解を求める手法である『EMアルゴリズム』を紹介します。
例えば被説明変数$y_i$、説明変数$x_i~\small{(i=1,\cdots,n)}$が得られた時、単回帰をするとします。この時、$(x_i,y_i) ~ \small{(i=1,\cdots,n)}$が全て観測されていれば、背景の分布に基づく最尤解が求められましたね。
その通りです。
「データに欠損がある時に最尤解を求める」とは、上記の$x_1,y_2$が欠損するなど一部のデータに欠損があっても最尤解が求められるということでしょうか？
もちろん全ての場合に求められるわけではないですが、イメージとしてはその通りです。
はい。では、潜在変数がある場合とはどの様な場合でしょうか？
潜在変数について説明します。
例えば東京都において1000人をランダムサンプルして身長を測定するという状況を考えます。
手元にある身長のデータを確認したところ、異なるグループから構成されたデータであると解析者は考えました。そこで、解析者はグループの違いを考慮したモデリングをしたいとします。グループ1、グループ2の身長の元となる分布をそれぞれ$\N(\mu_1,\sigma^2_1)$、$\N(\mu_2,\sigma^2_2)$とし、$Y_i$がグループを表す確率変数とした時、被験者$i$の身長$X_i$を以下の様にモデリングすることができます。
$$\begin{aligned} X_i | Y_i=1（グループ1） \overset{i.i.d} \sim \N(\mu_1,\sigma_1^2) \\ X_i | Y_i=2（グループ2） \overset{i.i.d} \sim \N(\mu_2,\sigma_2^2) \end{aligned}$$
この様に、収集はしてないが背景にあると考えられる$Y_i$の様な変数を『潜在変数』と言います。 なお、$Y_i=1$となる確率を$\pi$とする時、$X_i$が従う分布を書くとすれば、$X_i \overset{i.i.d} \sim \pi \N(\mu_1,\sigma_1^2) + (1-\pi) \N(\mu_2,\sigma_2^2) ~~{\small (i=1, \ldots, n)}$、の様に混合正規分布となります。 通常の最尤推定ではこの混合正規分布の尤度関数を最大化するパラメータを計算することになります。ただその直接的な計算が難しいため、EMアルゴリズムでは潜在変数を用いた最適化によってその様なパラメータを求めます。
1. EMアルゴリズムを使う場面
EMアルゴリズムは、データに欠損/潜在変数がある時に最尤解を得るために使用されます。
確率変数$X,Z$について、
・$X$：観測データ
・$Z$：欠損データ（ないし潜在変数データ）
・$(X,Z) \sim p(x,z;\theta)$
とした時、本来は対数尤度$\log p(x;\theta) = \log \{ \int p(x,z;\theta) dz \}$を最大化する最尤解を解析的に求められれば良いわけです。
しかしながら、それが困難であるために、上記対数尤度を最大化する最尤解を数値的に求める手法であるEMアルゴリズムが使用されます。
2. EMアルゴリズムの実際
全体像としては、E (Expectation) stepでQ関数を計算し、M (Maximization) stepでQ関数を最大化（する$\theta$を求める）、を繰り返す、ことになります。
つまり、M stepのたびに$\theta$が更新されていくということです。
○E step: Q関数の計算
以下の$Q(\theta, \theta^*)$を計算します。
EMアルゴリズムではstepが複数回繰り返されますが、$\theta^*$は現stepでのfixされた値、$\theta$はfixされてない値、と見ます。
$$\begin{aligned} Q(\theta,\theta^*) &= E_{Z | X, \theta^*} [\log p(X, Z; \theta)] \\[10px] &= \int \log p(x,z;\theta) \cdot p(z | x, \theta^*) dz \end{aligned}$$
Q関数は完全データの尤度を（$Z$の）事後分布で期待値をとったものになります。
なお、この式は$\int$計算において$z$が消えるため、$(x,\theta,\theta^*)$を含む式になりますが、この内fixされてないのは$\theta$のみです。
○M step: Q関数を最大化（する$\theta$を求める）
$$\begin{aligned} \theta^{new} = \underset{\theta}{\argmax} Q(\theta,\theta^*) \end{aligned}$$を求め、
・$| \theta^{new}-\theta^* | \lt$（事前に定めたある小さな値）→収束したとみなしてstepを終了
・$| \theta^{new}-\theta^* | \geqq$（事前に定めたある小さな値）→収束しきってないとみなして$\theta^* \leftarrow \theta^{new}$と置換して次のstep(E step)へ
とします。
補足.
- EMアルゴリズムにおいては$\theta$の初期値$\theta_0$が必要となるため、最初にそれを与えます。stepの進行毎に$| \theta^{new}-\theta^* |$が小さくなり、それが（事前に定めたある小さな値）を下回れば終了となります。ただし、使用する統計ソフト・パッケージによってはパラメータではなく尤度関数（対数尤度関数）の変化幅がある小さな値を下回ったときにstepを終了するとしてるものもあります。
- EMアルゴリズムにおいては、各M stepにおいて尤度は単調増加となります。（下の「EMアルゴリズムでやっていることの意味」で詳しく説明します）つまり、stepを回せば回すほど尤度は大きくなることが保証されており、尤度がほぼ変化しなくなった段階で終了となるわけです。
EMアルゴリズムが初見の方はイメージが湧いてないと思いますので、以下の例（データに欠損がある時）を通じてEMアルゴリズムのイメージを確立してください。なお、計算は結構大変なので全て理解できなくてもおおよその流れを理解すれば問題ないかと思います。
例題1.
$X_i \overset{i.i.d} \sim \N(\mu,\sigma^2) ~~ \small{(i=1,\cdots,n)}$なる$X_i$のデータを収集する予定であったが、実際にはある値$c ~~ \small{(\gt 0)}$以上の$(n-m)$個の$X_i ~~ \small{(i=m+1,\cdots,n)}$についてはデータが収集されなかった。
つまり、$M_i$を、
$$\begin{aligned} M_i = \begin{cases} 0 ~~ (X_iが収集された) \\ 1 ~~ (X_iが収集されなかった（欠損）) \end{cases} \end{aligned}$$とすると、データ$(X_i,M_i)$は、$(X_1,0),\cdots,(X_m,0),(X_{m+1},1),\cdots,(X_n,1)$という構造をもつ。
この時、EMアルゴリズムを用いて最尤推定値$\hat{\mu}$を求めよ。
（注意：$\sigma^2$は既知とする）
解答.
Q関数を求める前にQ関数の計算に必要な、
a. 完全データの尤度
b. 事後分布
を求めておく。
（なお、$(X_1,\cdots,X_n),(M_1,\cdots,M_n)$をそれぞれ$X_{1:n},M_{1:n}$と略記する）
まず、
$$\begin{aligned} Pr\{ M_i = m_i | X_i \} &= m_i \cdot \vec 1\{ X_i \geqq c \} + (1-m_i) \cdot \vec 1\{ X_i \lt c \} \\[10px] Pr\{ X_i | \mu \} &= \N(X_i | \mu, \sigma^2) \end{aligned}$$であるから
（注意：$\sigma^2$は既知である。またここでは”$\N(X_i|\mu,\sigma^2)$”はその確率密度関数を表すものとする。）、
a.完全データの尤度$\prod_{i=1}^n Pr\{ X_i, \mu_i | \mu \}$は、
$$\begin{aligned} \prod_{i=1}^n Pr\{ X_i, M_i | \mu \} &= \prod_{i=1}^n [ M_i \cdot \vec 1\{ X_i \geqq c \} + (1-M_i) \cdot \vec 1\{ X_i \lt c \} ] \times \N(X_i | \mu, \sigma^2) ] \\[10px] &= \prod_{i=1}^{\textcolor{red}{m}} \vec 1\{ X_i \lt c \} \cdot \N(X_i | \mu, \sigma^2) \times \prod_{i=\textcolor{red}{m+1}}^{\textcolor{red}{n}} \vec 1\{ X_i \geqq c \} \cdot \N(X_i | \mu, \sigma^2) \\ &{\scriptsize(場合わけした)} \\[10px] &= \prod_{i=1}^m \N(X_i | \mu, \sigma^2) \times \prod_{i=m+1}^{n} \N(X_i | \mu, \sigma^2) ~~~~~ \mathrm{(A)} \end{aligned}$$
また、b. 事後分布$Pr\{ X_{m+1:n} | X_{1:m},M_{1:n},\mu^{*} \}$は、
$$\begin{aligned} Pr\{ X_{m+1:n} | X_{1:m},M_{1:n},\mu^{*} \} &= \prod_{i=m+1}^{n} Pr\{ X_i | M_i=1,\ \mu^* \} \\[10px] &= \prod_{i=m+1}^{n} \frac{ Pr\{ X_i,M_i=1 | \mu^* \} }{ Pr\{ M_i=1 | \mu^* \} } \\ &{\scriptsize(ベイズの定理より)} \\[10px] &= \prod_{i=m+1}^{n} \frac{ \N(X_i | \mu^*,\sigma^2 ) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } ~~~~~ \mathrm{(B)} \\ &{\scriptsize(今後\Phi,\phiはそれぞれ\N(0,1)の分布関数、密度関数を表すものとする)} \end{aligned}$$となる。
以上よりQ関数$Q(\mu | \mu^*)$は、
$$\begin{aligned} Q(\mu | \mu^*) &= E_{X_{m+1:n | X_{1:m}, M_{1:n}, \mu^*}} [\log \mathrm{(A)}] \\[10px] &= E_{X_{m+1:n | X_{1:m}, M_{1:n}, \mu^*}} [ \sum_{i=1}^m \{ -\frac{1}{2} \log (2 \pi \sigma^2)\underbrace{-\frac{(X_i-\mu)^2}{2 \sigma^2}}_{①} \} + \sum_{i=m+1}^{n} \{ -\frac{1}{2} \log (2 \pi \sigma^2)\underbrace{-\frac{(X_i-\mu)^2}{2 \sigma^2}}_{②} \} ] \end{aligned}$$となる。
Q関数というのは、完全データの尤度$\mathrm{(A)}$を事後分布$\mathrm{(B)}$で期待値をとったものに相当しました。
M stepではQ関数を最大化する$\mu$を求めるが、具体的には$\frac{\partial Q}{\partial \mu}=0$なる$\mu$を求めることになる。
$\frac{\partial Q}{\partial \mu}$を計算するにあたり上記$Q(\mu | \mu^*)$の中で$\mu$に関する部分は①②のみであり、計算が大変な②を以下先に計算しておく。
$$\begin{aligned} ② &= \int_{X_{m+1:n}} \sum_{i=m+1}^n -\frac{(X_i-\mu)^2}{2 \sigma^2} \cdot \prod_{\textcolor{red}{j}=m+1}^n \frac{ \N(X_{\textcolor{red}{j}} | \mu^*,\sigma^2) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } dX_{m+1:n} \\[10px] &= \sum_{i=m+1}^n \int_{X_{m+1:n}} -\frac{(X_i-\mu)^2}{2 \sigma^2} \cdot \prod_{\textcolor{red}{j}=m+1}^n \frac{ \N(X_{\textcolor{red}{j}} | \mu^*,\sigma^2) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } dX_{m+1:n} \\ &{\scriptsize(\intと\sumを入れ替えた)} \\[10px] &= \sum_{i=m+1}^n \int_{X_{i}} -\frac{(X_i-\mu)^2}{2 \sigma^2} \cdot \frac{ \N(X_{\textcolor{red}{i}} | \mu^*,\sigma^2) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } dX_{i} \\ &{\scriptsize(X_{m+1},\cdots,X_nの独立性より)} \\[10px] &= \sum_{i=m+1}^n \int_c^{\infty} -\frac{(X_i-\mu)^2}{2 \sigma^2} \cdot \frac{ \frac{1}{\sqrt{2 \pi \sigma^2}} \exp[-\frac{(x_i-\mu^*)^2}{2 \sigma^2}] }{ \Phi(-\frac{c-\mu^*}{\sigma}) } dX_{i} \\ &{\scriptsize(以下、Y_i=\frac{X_i-\mu^*}{\sigma}, Z=\frac{\mu-\mu^*}{\sigma}、と置いて計算する)} \\[10px] &= \sum_{i=m+1}^n \int_{\frac{c-\mu^*}{\sigma}}^{\infty} -\frac{1}{2} (Y_i-Z)^2 \cdot \frac{ \frac{1}{\sqrt{2 \pi}} \exp[ -\frac{Y_i^2}{2} ] }{ \Phi(-\frac{c-\mu^*}{\sigma}) } dY_{i} \\[10px] &= \sum_{i=m+1}^n -\frac{1}{2 \sqrt{2 \pi}} \cdot \frac{1}{\Phi(-\frac{c-\mu^*}{\sigma})} \cdot \{ \underbrace{ \int_{\frac{c-\mu^*}{\sigma}}^{\infty} Y_i^2 \cdot \exp[-\frac{Y_i^2}{2}] dY_i }_{\muに関係ない項} + \int_{\frac{c-\mu^*}{\sigma}}^{\infty} (-2Y_i Z) \cdot \exp[-\frac{Y_i^2}{2}] dY_i + \int_{\frac{c-\mu^*}{\sigma}}^{\infty} (Z^2) \cdot \exp[-\frac{Y_i^2}{2}] dY_i \} \\[10px] &= \sum_{i=m+1}^n -\frac{1}{2 \sqrt{2 \pi}} \cdot \frac{1}{\Phi(-\frac{c-\mu^*}{\sigma})} \cdot \{ (\muに関係ない項)+[ 2 Z \exp[-\frac{Y_i^2}{2}] ]_{\frac{c-\mu^*}{\sigma}}^{\infty} + Z^2 \int_{\frac{c-\mu^*}{\sigma}}^{\infty} \exp [-\frac{Y_i^2}{2}] dY_i \} \\[10px] &= \sum_{i=m+1}^n -\frac{1}{2 \sqrt{2 \pi}} \cdot \frac{1}{\Phi(-\frac{c-\mu^*}{\sigma})} \cdot \{ (\muに関係ない項) -2Z \cdot \exp[-\frac{(c-\mu^*)^2}{2 \sigma^2}] + Z^2 \cdot \sqrt{2 \pi} \int_{\frac{c-\mu^*}{\sigma}}^{\infty} \underbrace{ \frac{1}{\sqrt{2 \pi}} \exp [-\frac{Y_i^2}{2}] }_{\phi(Y_i)である} dY_i \\[10px] &= \sum_{i=m+1}^n -\frac{1}{2 \sqrt{2 \pi}} \cdot \frac{1}{\Phi(-\frac{c-\mu^*}{\sigma})} \cdot \{ (\muに関係ない項) -2Z \cdot \sqrt{2 \pi} \cdot \underbrace{ \frac{1}{\sqrt{2 \pi}} \cdot \exp[-\frac{(c-\mu^*)^2}{2 \sigma^2}]}_{=\phi(\frac{c-\mu^*}{\sigma})} + Z^2 \cdot \sqrt{2 \pi} \cdot \Phi(\textcolor{red}{-}\frac{c-\mu^*}{\sigma}) \} \\[10px] &= \sum_{i=m+1}^{\infty} -\{ (\muに関係ない項) -\frac{\mu-\mu^*}{\sigma} \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } + \frac{(\mu-\mu^*)^2}{2 \sigma^2} \} \\[10px] &= -(n-m) \{ (\muに関係ない項) -\frac{\mu-\mu^*}{\sigma} \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } + \frac{(\mu-\mu^*)^2}{2 \sigma^2} \} \end{aligned}$$
以上より$Q(\mu | \mu^*)$は、
$$\begin{aligned} Q(\mu | \mu^*) &= (\muに関係ない項) – \sum_{i=1}^m \frac{(X_i-\mu)^2}{2 \sigma^2} – (n-m) \cdot \{ -\frac{\mu-\mu^*}{\sigma} \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } + \frac{(\mu-\mu^*)^2}{2 \sigma^2} \} \\ &= (\muに関係ない項) -\frac{1}{2 \sigma^2} [ \underbrace{ \sum_{i=1}^m (X_i-\mu)^2 + (n-m) \{ -2 \sigma (\mu-\mu^*) \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } + (\mu-\mu^*)^2 \} }_{Q'(\mu | \mu^*)とする} ] \end{aligned}$$となる。
これより、
$$\begin{aligned} \frac{\partial Q’}{\partial \mu} &= -2 \sum_{i=1}^m (X_i-\mu) + (n-m) \{ -2 \sigma \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) } + 2(\mu-\mu^*) \} \end{aligned}$$となるので、
$$\begin{aligned} \frac{\partial Q’}{\partial \mu} &= 0 \\ \Rightarrow \mu &= \frac{1}{n} \{ \sum_{i=1}^m X_i + (n-m) [\mu^* + \sigma \cdot \frac{ \phi(\frac{c-\mu^*}{\sigma}) }{ \Phi(-\frac{c-\mu^*}{\sigma}) }] \} \end{aligned}$$
3. EMアルゴリズムでやっていることの意味
EMアルゴリズムの目的は、欠損/潜在変数（$z$）がある時に$\log p(x;\theta)$を最大化する最尤解$\hat{\theta}$を数値的に求めることでした。
ここではEMアルゴリズムでやってきたことの意味を確認してみます。
まず、EMアルゴリズムの根幹をなす式$\mathrm{(A)}$を準備します。
$$\begin{aligned} \log p(x;\theta) &= \int \log p(x; \theta) q(z) dz \\ &{\scriptsize(\int q(z) dz = 1なるqに対して)} \\[10px] &= \int \{ (\log p(x;\theta) + \textcolor{red}{\log p(z|x,\theta)}) \textcolor{red}{-\log p(z|x,\theta)} \} q(z) dz \\[10px] &= \int \{ \log p(x,z;\theta)-\log p(z|x,\theta) \} q(z) dz \\ &{\scriptsize(\log p(x,z;\theta) = \log p(x;\theta) + \log p(z|x,\theta))} \\[10px] &= \int [ \{ \log p(x,z;\theta) \textcolor{red}{-\log q(z)} \} + \{ \textcolor{red}{\log q(z)}-\log p(z|x,\theta) \} ] q(z) dz \\[10px] &= \underbrace{\int \{ \log (\frac{p(x,z;\theta)}{q(z)}) \cdot q(z) \} dz}_{\textcolor{red}{L(q,\theta)}} + \underbrace{\int \log (\frac{q(z)}{p(z|x,\theta)}) \cdot q(z) dz}_{\textcolor{red}{KL(q||p_{z|x,\theta})}} ~~~~~ \mathrm{(A)} \end{aligned}$$
$\mathrm{(A)}$から、$\log p(x;\theta)$は$L(q,\theta)$と$KL(q||p_{z|x,\theta})$という2つの項に分解されることがわかります。
ここでのポイントは式$\mathrm{(A)}$の第2項はKL情報量であり必ず0以上であることが保証されているため（参照：＜カルバック・ライブラー情報量＞）、$\log p(x;\theta)$は少なくとも$L(q,\theta)$以上であるということです。
適当な分布$q,p$に対して、常に$KL(q||p) \geqq 0$は成立しました。
また、$KL(q||p_{z|x,\theta})=0$となるのは$q=p$の時でした。
これは下記でも使うので要チェックです。
KL情報量についてはとりあえずこれだけわかっていればokですが、復習したい方は＜カルバック・ライブラー情報量＞を参照ください。
さて、以上を抑えた上でステップを具体的にみてみましょう。
EMのstep iでの対数尤度を$\log p(x | \theta_i)$とします。
上記関係式
$$\begin{aligned} \log p(x | \theta) = L(q,\theta) + KL(q||p_{z|x,\theta}) \end{aligned}$$から、$L(q,\theta_i)$は$\log p(x | \theta_i)$の下限であり、
$$\begin{aligned} q = q_{i+1} = p (z | x, \theta_i) \end{aligned}$$とすることで、下限$L(q,\theta_i)$は現時点での対数尤度は$\log p(x | \theta_i)$と一致します：
$$\begin{aligned} \log p(x | \theta_i) = L(q,\theta_i) \end{aligned}$$
ここで、$L(q_{i+1},\theta^*) \gt L(q_{i+1},\theta_i)$なる$\theta^*$が見つかったとします。
すると上記関係式から、
$$\begin{aligned} \log p(x | \theta^*) &= L(q_{i+1}, \theta^*) + KL(q_{i+1} ||p_{z|x,\theta^*}) \\ &\geqq L(q_{i+1},\theta^*) \\ &\gt L(q_{i+1}, \theta_i) = \log p(x | \theta_i) \end{aligned}$$となり、対数尤度が増加します。
したがって、例えば$L(q_{i+1}, \theta)$を最大化する様な$\theta$を求められれば、その様な$\theta$の対数尤度は$\theta_i$の対数尤度よりも大きくなることがわかります。 ここで問題になるのは、その様な最大化をどの様にするかということです。
ここで、$q_{i+1} = p(z | x, \theta_i)$、とした時、
$$\begin{aligned} L(q_{i+1}, \theta) &= \int \log \frac{ p(x, z ; \theta) }{ p(z | x, \theta_i) } \cdot p(z | x, \theta_i) dz \\ &= Q(\theta, \theta_i)-\underbrace{\int \log p(z | x, \theta_i) \cdot p(z | x, \theta_i) dz}_{\theta_iのみに依存する項} \\ &{\scriptsize(Q(\theta,\theta_i)は上記で定義したものを参照のこと)} \end{aligned}$$となるため、$Q(\theta,\theta_i)$を最大化する$\theta$を求めることができれば、$L(q_{i+1},\theta)$を最大化する$\theta$を求めることができるのです。
$Q(\theta,\theta_i)$なんてものを持ち出さなくても$L(q_{i+1},\theta)$を最大化する$\theta$を求めてもよいのです。
しかし、$Q(\theta,\theta_i)$を準備してこれを最大化する$\theta$を求める方が計算が楽なので、$Q(\theta,\theta_i)$をわざわざ持ち出したという次第です。
まとめ.
- EMアルゴリズムはデータに欠損/潜在変数がある時にパラメータの最尤解を得る手法である。
- E step、M stepを繰り返し、パラメータの更新幅が事前に定めた小さな値を下回れば更新を終了する。
- 各M stepにおいて尤度（対数尤度）は単調増加であることが保証されている。
