prml9.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第9章
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{chapter}{8}
\chapter{「混合モデルとEM」の数式の補足}
面倒なので特に紛らわしいと思わない限り$\bm{x}$を$x$と書いたりします. また対数尤度関数を$F$と書くことが多いです.

\newcommand{\wji}{w_{ji}^{(1)}}
\newcommand{\wkj}{w_{kj}^{(2)}}

\section{復習}
よく使ういくつかの式を書いておく. どれも今までに既に示したものである.\\
2章や3章を参照.

\subsection{行列の公式}
\begin{eqnarray*}
&& \quads{A}{x}=\tr\left(Ax\trans{x}\right),\\
&& \dif{A}\log|A|=\trans{(A^{-1})},\\
&& \dif{x}\log|A|=\tr\left(A^{-1}\dif{x}A\right),\\
&& \dif{A}\tr(A^{-1}B) = -\trans{(A^{-1}BA^{-1})}.
\end{eqnarray*}
\vspace{0pt}

\subsection{微分}
関数$f$に対して対数関数の微分は
$$
(\log f)'=\frac{f'}{f}.
$$
よって逆に
$$
f'=f \cdot (\log f)'.
$$
ガウス分布など対数の微分が分かりやすいときによく使う.

\subsection{ガウス分布}
$$
\calN = \calN(x\,|\,\mu,\Sigma)=\frac{1}{(2\pi)^{D/2}}|\Sigma|^{-1/2}\exp\left(-\half\quads{{\Sigma^{-1}}}{(x-\mu)}\right).
$$
期待値と分散について
$$
\EE[x]=\mu, \quad \var[x]=\Sigma, \quad
 \EE[\outp{x}]=\outp{\mu}+\Sigma, \quad \EE[\inp{x}]=\inp{\mu}+\tr(\Sigma).
$$
最後の式は3番目から出る.
$$
\EE[x_i^2]=(\outp{\mu})_{ii}+\Sigma_{ii}=\mu_i^2+\Sigma_{ii}.
$$
よって
$$
\EE[\inp{x}]=\sum_i \EE[x_i^2]=\inp{\mu}+\tr(\Sigma).
$$
\vspace{0pt}

\section{混合ガウス分布}
離散的な潜在変数を用いて混合ガウス分布の定式化を行う.
$K$次元$2$値確率変数$z$を考える（どれか一つの成分のみが$1$であとは$0$）．
つまり
$$
\sum_k z_k = 1.
$$
$z$の種類は$K$個である.
どれか一つの成分のみが1であることを示すのに$K$個の変数$\{z_k\}$で
表現するのはまわりくどい印象を受けるが,
$z$を複数扱う場合にはこの方が記号が簡単になる.

$0 \le \pi_k \le 1$という係数を用いて
$$
p(z_k=1)=\pi_k
$$
という確率分布を与える.
$$
p(z)=\prod_k \pi_k^{z_k}, \quad p(x\,|\,z_k=1)=\calN(x\,|\,\mu_k,\Sigma_k)
$$
なので
$$
p(x\,|\,z)=\prod_k \calN(x\,|\,\mu_k,\Sigma_k)^{z_k}.
$$
これらを合わせて
\begin{eqnarray*}
p(x) &=& \sum_z p(z)\,p(x\,|\,z)
 = \sum_z \prod_k \left(\pi_k\,\calN(x\,|\,\mu_k,\Sigma_k)\right)^{z_k}\\
 && （\text{$z_k$はどれか一つのみが$1$であとは$0$なので}）\\
 &=& \sum_{k} \pi_k\,\calN(x\,|\,\mu_k,\Sigma_k).
\end{eqnarray*}

$x$が与えられたときの$z$の条件付き確率$p(z_k=1\,|\,x)$を$\gamma(z_k)$とする.
$$
\gamma(z_k)=\frac{p(z_k=1)\,p(x\,|\,z_k=1)}{\sum_j p(z_j=1)\,p(x\,|\,z_j=1)}
 = \frac{\pi_k\,\calN(x\,|\,\mu_k,\Sigma_k)}{\sum_j \pi_j\,\calN(x\,|\,\mu_j,\Sigma_j)}.
$$
これを混合要素$k$の観測値$x$に対する負担率という.

\section{混合ガウス分布のEMアルゴリズム}
混合ガウス分布において観測したデータ集合を$\trans{X}=(x_1,\ldots, x_N)$,
対応する潜在変数を$\trans{Z}=(z_1,\ldots, z_N)$とする.
$X$は$N \times D$行列で$Z$は$N \times K$行列.

対数尤度関数の最大点の条件をもとめる.
$$
F=\log p(X\,|\,\vpi,\mu,\Sigma)=\sum_{n=1}^N \log\left(\sum_{j=1}^K \pi_j\,\calN(x_n\,|\,\mu_j,\Sigma_j)\right)
$$
とする.
$$
\dif{\mu} \log \calN(x\,|\,\mu,\Sigma)
= \dif{\mu} \left(-\half\quads{{\Sigma^{-1}}}{{(x-\mu)}}\right)
= \Sigma^{-1}(x-\mu)
$$
より
$$
\dif{\mu} \calN = \calN \cdot \left(\dif{\mu} \log \calN\right) = \calN \cdot \Sigma^{-1}(x-\mu).
$$
もちろんガウス分布の微分は普通にそのまましてもよい.
だが今回は対数をとってから微分をとった方が,
微分してでてくる$\calN$が$\gamma(z_{nk})$の一部となることを見通しやすいのでそうしてみた.

さて$\calN_{nk}=\calN(x_n\,|\,\mu_k,\Sigma_k)$とおいて
\begin{eqnarray*}
\dif{\mu_k}F
 &=& \sum_n \frac{\pi_k \dif{\mu_k}\calN_{nk}}{\sum_j \pi_j\,\calN_{nj}}
 = \sum_n \left(\frac{\pi_k\,\calN_{nk}}{\sum_j \pi_j\,\calN_{nj}}\right)\dif{\mu_k}\log \calN_{nk} \\
 &=& \sum_n \gamma(z_{nk}) \dif{\mu_k}\log \calN_{nk}
 = \Sigma_k^{-1}\left(\sum_n \gamma(z_{nk})(x_n-\mu_k)\right)=0.
\end{eqnarray*}
よって
$$
\sum_n \gamma(z_{nk})x_n - \left(\sum_n \gamma(z_{nk})\right)\mu_k=0.
$$
$$
N_k=\sum_n \gamma(z_{nk})
$$
とおくと
$$
\mu_k=\frac{1}{N_k}\sum_n \gamma(z_{nk})x_n.
$$
これは$\mu_k$が$X$の重みつき平均であることを示している.
次に$\Sigma_k$に関する微分を考える.
$$
\calN = \calN(x\,|\,\mu,\Sigma)
$$
のとき
$$
\log \calN = -\frac{D}{2}\log(2\pi)-\half\log |\Sigma| - \half\tr\left(\Sigma^{-1}\outp{(x-\mu)}\right)
$$
なので$\trans{\Sigma}=\Sigma$だから
$$
\dif{\Sigma} (\log \calN)
= -\half(\Sigma^{-1})
  +\half\left(\Sigma^{-1}(x-\mu)\trans{(x-\mu)}\Sigma^{-1}\right).
$$
よって$\mu_k$の微分と同様にして
\begin{eqnarray*}
\dif{\Sigma_k}F
 &=& \sum_n \gamma(z_{nk}) \dif{\Sigma_k} \log \calN_{nk}\\
 &=&\sum_n \gamma(z_{nk})\left(-\half(\Sigma_k^{-1})
        + \half\left(\Sigma_k^{-1}\outp{(x_n-\mu_k)}\Sigma_k^{-1}\right)\right) = 0.
\end{eqnarray*}
よって
$$
\sum_n \gamma(z_{nk})\left(I-\outp{(x_n-\mu_k)}\Sigma_k^{-1}\right)=0.
$$
$$
\Sigma_k=\frac{1}{N_k}\sum_n \gamma(z_{nk})\outp{(x_n-\mu_k)}.
$$

最後に$\pi_k$に関する微分を考える.
$\sum_k \pi_k=1$の制約を入れる.
$$
G=F+\lambda\left(\sum_k \pi_k-1\right)
$$
とすると
$$
\dif{\pi_k}{G}
=\sum_n \frac{\calN_{nk}}{\sum_j \pi_j\calN_{nj}}+\lambda=\sum_n \gamma(z_{nk})/\pi_k+\lambda=N_k/\pi_k+\lambda=0.
$$
つまり$N_k = -\lambda \pi_k$.
よって
$$
N=\sum_k N_k=\sum_k (-\lambda \pi_k) = -\lambda.
$$
よって
$$
\pi_k=\frac{N_k}{-\lambda}=\frac{N_k}{N}.
$$
\vspace{0pt}

\section{混合ガウス分布再訪}
$$
p(z)=\prod_k \pi_k^{z_{k}}, \quad p(x\,|\,z)=\prod_k \calN(x\,|\,\mu_k,\Sigma_k)^{z_k}
$$
より
\begin{eqnarray*}
F
 &=& \log p(X,Z\,|\,\mu,\Sigma,\vpi)
 = \log \left(\prod_{n,k} \pi_k^{z_{nk}}\calN(x_n\,|\,\mu_k,\Sigma_k)^{z_{nk}}\right)\\
 &=& \sum_{n,k} z_{nk}(\log \pi_k + \log \calN_{nk}).
\end{eqnarray*}
$z_n$は$(0, 0,\ldots, 1, 0, \ldots, 0)$の形でその成分が$z_{nk}$.
また$\sum_k \pi_k=1$である.
上式の微分を考えると
$$
G=F+\lambda\left(\sum_k \pi_k-1\right)
$$
として
$$
\dif{\pi_k}G=\sum_n z_{nk}\frac{1}{\pi_k}+\lambda=\left(\sum_n z_{nk}\right)/\pi_k + \lambda=0.
$$
よって
$$
\pi_k = -\frac{1}{\lambda}\sum_n z_{nk}.
$$
$$
\sum_k \pi_k=-\frac{1}{\lambda}\sum_{n,k}z_{nk}=-\frac{N}{\lambda}=1.
$$
よって$\lambda=-N$.
つまり
$$
\pi_k=\frac{1}{N}\sum_n z_{nk}.
$$
完全データ集合についての対数尤度関数の最大化は解けるが,
潜在変数が分からない場合の不完全データに関する
対数尤度関数の最大化は困難である.
この場合は潜在変数の事後分布に関する完全データ尤度関数の期待値を考える.
$$
p(Z\,|\,X,\mu,\Sigma,\vpi)
 = \frac{p(X,Z\,|\,\mu,\Sigma,\vpi)}{p(X\,|\,\mu,\Sigma,\vpi)}
 \propto \prod_{n,k} (\pi_k\,\calN_{nk})^{z_{nk}}.
$$
$$
\EE[z_{nk}]=\frac{\sum_{z_n} z_{nk} \prod_j (\pi_j\,\calN_{nj})^{z_{nj}}}
               {\sum_{z_n} \prod_j(\pi_j\,\calN_{nj})^{z_{nj}}}
         =\frac{\pi_k\,\calN_{nk}}{\sum_j \pi_j\,\calN_{nj}} = \gamma(z_{nk}).
$$
よって
$$
F=\EE_Z[\log p(X,Z\,|\,\mu,\Sigma,\vpi)]=\sum_{n,k}\gamma(z_{nk})(\log \pi_k + \log \calN_{nk}).
$$
まずパラメータ$\mu$, $\Sigma$, $\vpi$を適当に決めて
負担率$\gamma(z_{nk})$を求め,
それを固定した範囲で$\mu_k$, $\Sigma_k$, $\pi_k$を動かし$F$を最大化する.
今までと同様にできる.
$F'=F+\lambda\sum_k\left(\pi_k-1\right)$として
$$
\dif{\pi_k}F'=\sum_n \gamma(z_{nk})(1/\pi_k)+\lambda=0
$$
より
$$
\sum_n \gamma(z_{nk})=-\lambda \pi_k.
$$
$$
\sum_{n,k}\gamma(z_{nk})=-\lambda\left(\sum_k \pi_k\right)=-\lambda=N
$$
より
$$
\pi_k=\frac{1}{N}\sum_n \gamma(z_{nk})=\frac{N_k}{N}.
$$
\begin{eqnarray*}
\dif{\mu_k}F
&=&\sum_n\gamma(z_{nk})\left(\Sigma_k^{-1}(x_n-\mu_k)\right)
\\
&=&\Sigma_k^{-1}\left(\sum_n \gamma(z_{nk})x_n-\left(\sum_n \gamma(z_{nk})\right)\mu_k\right)=0.
\end{eqnarray*}
よって
$$
\mu_k=\frac{1}{N_k}\sum_n \gamma(z_{nk})x_n.
$$
$$
\dif{\Sigma_k}F=\sum_n \gamma(z_{nk}) \dif{\Sigma_k}\log \calN_{nk}=0
$$
として同様（流石に略）．

\section{$K$-meansとの関連}
PRML式(9.43)は不正確. $E$ではなく$\epsilon E$を考えないとPRML式(9.43)の右辺にはならない.
PRML式(9.40)を$E$とおく.
$$
E=\sum_{n,k} \gamma(z_{nk})\left(\log \pi_k + \log \calN(x_n\,|\,\mu_k,\Sigma_k)\right).
$$
$\epsilon E$に
$$
\calN(x\,|\,\mu_k,\Sigma_k)
=\frac{1}{(2\pi \epsilon)^{D/2}}
 \exp\left(-\frac{1}{2\epsilon}\Norm{x-\mu_k}^2\right)
$$
を代入する.
$$
\epsilon E = \sum_{n,k} \gamma(z_{nk}) \left(
               \epsilon \log \pi_k
               - \frac{D}{2}\epsilon \log (2\pi \epsilon)
               - \half\Norm{x_n-\mu_k}^2
             \right).
$$
$\epsilon \rightarrow 0$で
$$
\gamma(z_{nk}) \rightarrow r_{nk}, \quad
\epsilon \log \pi_k \rightarrow 0, \quad
\epsilon \log (2\pi \epsilon) \rightarrow 0
$$
より
$$
\epsilon E \rightarrow -\half\sum_{n,k} r_{nk} \Norm{x_n-\mu_k}^2 = -J.
$$
よって期待完全データ対数尤度の最大化は$J$の最小化と同等.

\section{混合ベルヌーイ分布}
$x$を$D$次元ベクトル$x:=\trans{(x_1,\ldots,x_D)}$, $0\le x_i \le 1$,
$\mu$も$D$次元ベクトル$\mu:=\trans{(\mu_1, \ldots, \mu_D)}$, $0 < \mu_i < 1$とする.
まず
$$
p(x\,|\,\mu):=\prod_{i=1}^D \mu_i^{x_i}(1-\mu_i)^{(1-x_i)}
$$
について考える. $\EE[x]=\mu$は容易に分かる.
$$
\EE[x_i x_j]=
\begin{cases}
\mu_i \mu_j & (i \ne j),\\
\mu_i & (i = j).
\end{cases}
$$
よって
$$
\var[x]_{ij}=\EE\left[\outp{(x-\mu)}\right]_{ij}=\EE[x_i x_j]-(\outp{\mu})_{ij}=(\mu_i - \mu_i^2)\delta_{ij}
$$
より
$$
\var[x]=\diag(\mu_i(1-\mu_i)).
$$

さて次に$K$個の$D$次元ベクトル$\mu_k$の組$\{\mu_1, \ldots, \mu_K\}$を$\vmu$として
$$
p(x\,|\,\mu_k):=\prod_i \mu_{ki}^{x_i}(1-\mu_{ki})^{(1-x_i)}
$$
の$\vpi:=\{\pi_1, \ldots, \pi_K\}$による混合分布
$$
p(x\,|\,\vmu):=\sum_k \pi_k\, p(x\,|\,\mu_k)
$$
を考えよう.
\begin{eqnarray*}\\&&
\EE[x]=\int x\, p(x\,|\,\vmu)\,dx
= \sum_k \pi_k \int x\, p(x\,|\,\mu_k)\,dx
= \sum_k \pi_k \EE_k[x]=\sum_k \pi_k \mu_k,
\\&&
\EE_k[\outp{x}]=\var_k[x]+\outp{\mu_k}=\Sigma_k+\outp{\mu_k}
\end{eqnarray*}
より
\begin{eqnarray*}
\var[x]
 &=& \EE\left[\outp{(x-\EE[x])}\right]
  = \EE\left[\outp{x}\right]-\outp{\EE[x]}\\
  &=&\sum_k \pi_k\left(\Sigma_k+\outp{\mu_k}\right)-\outp{\EE[x]}.
\end{eqnarray*}
データ集合$X=\{x_1, \ldots, x_N\}$が与えられたとき, 対数尤度関数は
$$
\log p(X\,|\,\vmu,\vpi)=\sum_n \log\left(\sum_k \pi_k\, p(x_n\,|\,\mu_k)\right).
$$
対数の中に和があるので解析的に最尤解をもとめられないためEMアルゴリズムを使う.
$x$に対応する潜在変数を$z=\trans{(z_1,\ldots,z_K)}$を導入する.
どれか一つのみ$1$でその他は$0$のベクトルである.
$z$の事前分布を
$$
p(z\,|\,\pi)=\prod_k \pi_k^{z_k}
$$
とする. $z$が与えられたときの条件付き確率は
$$
p(x\,|\,z,\vmu)=\prod_k p(x\,|\,\mu_k)^{z_k}.
$$
$$
p(x,z\,|\,\vmu,\vpi)=p(x\,|\,z,\vmu)\,p(z\,|\,\vpi)=\prod_k(\pi_k\, p(x\,|\,\mu_k))^{z_k}.
$$
よって
$$
p(x\,|\,\vmu,\vpi)=\sum_{z} p(x,z\,|\,\vmu,\vpi)=\sum_k \pi_k\, p(x\,|\,\mu_k).
$$
完全データ対数尤度関数は, $X=\trans{(x_1,\ldots,x_N)}\!$,
$Z=\trans{(z_1,\ldots,z_N)}\!$として
\begin{eqnarray*}
\log p(X,Z\,|\,\vmu,\vpi)
 &=& \sum_{n,k} z_{nk}\underbrace{\left(\log \pi_k + \sum_i x_{ni} \log \mu_{ki}+(1-x_{ni})\log (1-\mu_{ki})\right)}_{=:A_{nk}}\\
 &=& \sum_{n,k} z_{nk} A_{nk}.
\end{eqnarray*}
\begin{eqnarray*}
\EE[z_{nk}]
 &=& \frac{\sum_{z_n} z_{nk} \prod_j (\pi_j\, p(x_n\,|\,\mu_j))^{z_{nj}}}
          {\sum_{z_n} \prod_j (\pi_j\, p(x_n\,|\,\mu_j))^{z_{nj}}}\\
 && （z_{nk}=1 \text{となるものだけとるので}）\\
 &=& \frac{\pi_k\, p(x_n\,|\,\mu_k)}{\sum_j \pi_j\, p(x_n\,|\,\mu_j)}
\end{eqnarray*}
を$\gamma(z_{nk})$とおく. すると
$$
\EE_Z[\log p(X,Z\,|\,\vmu,\vpi)]=\sum_{n,k} \gamma(z_{nk}) A_{nk}.
$$
$$
N_k=\sum_n \gamma(z_{nk}), \quad \bar{x}_k = \frac{1}{N_k}\sum_n \gamma(z_{nk})x_n
$$
とおく.
\begin{eqnarray*}
F
 &:=& \EE_Z[\log p(X,Z\,|\,\vmu,\vpi)]\\
 &=& \sum_k (\log \pi_k)\left(\sum_n \gamma(z_{nk})\right)+\sum_{k,i} \log \mu_{ki}\left(\sum_n \gamma(z_{nk})x_{ni}\right)\\
 &&+ \sum_{k,i} \log (1-\mu_{ki})\left(\sum_n \gamma(z_{nk})(1-x_{ni})\right)\\
 &=& \sum_k N_k \log \pi_k + \sum_{k,i} N_k \bar{x}_{ki} \log \mu_{ki}
   + \sum_{k,i} \log (1-\mu_{ki}) N_k(1-\bar{x}_{ki}).
\end{eqnarray*}
$\mu_{ki}$に関する最大化.
\begin{eqnarray*}
\dif{\mu_{ki}}F
 &=& N_k \bar{x}_{ki} \frac{1}{\mu_{ki}}+\frac{-1}{1-\mu_{ki}} N_k (1-\bar{x}_{ki})
 \\
 &=& \frac{N_{k}}{\mu_{ki}(1-\mu_{ki})}(\bar{x}_{ki}(1-\mu_{ki})-(1-\bar{x}_{ki})\mu_{ki}) = 0.
\end{eqnarray*}
よって
$$
\bar{x}_{ki}-\bar{x}_{ki} \mu_{ki} - \mu_{ki}+\bar{x}_{ki} \mu_{ki}=\bar{x}_{ki}- \mu_{ki} = 0.
$$
よって
$$
\mu_k = \bar{x}_k.
$$
$\pi_k$に関する最適化.
$G=F+\lambda(\sum_k \pi_k-1)$とすると
$$
\dif{\pi_k}G = \frac{N_k}{\pi_k}+\lambda=0.
$$
よって
$$
N_k=-\lambda \pi_k, \quad N=\sum_k N_k = -\lambda \sum_k \pi_k = -\lambda.
$$
つまり$\lambda=-N$となり
$$
\pi_k=\frac{N_k}{N}.
$$
$0 \le p(x_n\,|\,\mu_k) \le 1$より
$$
\log p(X\,|\,\vmu,\vpi)=\sum_n \log \left(\sum_k \pi_k\, p(x_n\,|\,\mu_k)\right) \le \sum \log \left(\sum_k \pi_k\right)=0.
$$
よって尤度関数が発散することはない.

\section{ベイズ線形回帰に関するEMアルゴリズム}
EMアルゴリズムに基づいてベイズ線形回帰を考えてみる.
$w$を潜在関数と見なしてそれを最大化する方針を採る.
$$
p(w\,|\,t)=\calN(w\,|\,m_N,S_N)
$$
で$w$の事後分布が求まっているとする.
$$
p(t\,|\,w,\beta)=\prod_n \calN(t_n\,|\,\trans{w}\phi(x_n), \beta^{-1}), \quad
p(w\,|\,\alpha)=\calN(w\,|\,0, \alpha^{-1}I)
$$
であった.
このとき完全データ対数尤度関数は
$$
\log p(t,w\,|\,\alpha,\beta)=\log p(t\,|\,w,\beta) + \log p(w\,|\,\alpha).
$$
なので
\begin{eqnarray*}
F &=& \EE[\log p(t,w\,|\,\alpha,\beta)]\\
 &=& \EE\left[\sum_n \left(\half \log \left(\frac{\beta}{2\pi}\right)
     -\frac{\beta}{2}(t_n-\trans{w}\phi_n)^2\right)
     +\frac{M}{2}\log \left(\frac{\alpha}{2\pi}\right) - \frac{\alpha}{2}\inp{w}\right]\\
 &=& \frac{M}{2}\log \left(\frac{\alpha}{2\pi}\right)-\frac{\alpha}{2}\EE\left[\inp{w}\right]
     +\frac{N}{2}\log \left(\frac{\beta}{2\pi}\right)-\frac{\beta}{2}\sum_n \EE\left[(t_n-\trans{w}\phi_n)^2\right].
\end{eqnarray*}
$\alpha$に関する最大化
$$
\dif{\alpha}F=\frac{M}{2}\frac{1}{\alpha}-\half \EE[\inp{w}]=0.
$$
よって
$$
\alpha=\frac{M}{\EE[\inp{w}]}=\frac{M}{\inp{m_N}+\tr(S_N)}.
$$
$\beta$に関する最大化
$$
\dif{\beta}F=\frac{N}{2}\frac{1}{\beta}-\half\sum_n \EE\left[(t_n-\trans{w}\phi_n)^2\right]=0.
$$
よって
$$
\frac{1}{\beta}=\frac{1}{N}\sum_n \EE\left[(t_n-\trans{w}\phi_n)^2\right].
$$
\vspace{0pt}

\section{一般のEMアルゴリズム}
潜在変数をもつ確率モデルの最尤解を求めるための一般的手法.
$X$を確率変数, $Z$を潜在変数, $\theta$をパラメータとする.
目的は, 完全データ対数尤度関数$\log p(X,Z\,|\,\theta)$の最適化が
容易であるという仮定の下で$p(X\,|\,\theta)=\sum_Z p(X,Z\,|\,\theta)$を
最大化することである.

$Z$に対する分布を$q(Z)$とする
$$
p(X,Z\,|\,\theta)=p(Z\,|\,X,\theta)\,p(X\,|\,\theta).
$$
$$
\calL(q,\theta)=\sum_Z q(Z) \log \frac{p(X,Z\,|\,\theta)}{q(Z)}, \quad
\KL(q\,\|\,p)=-\sum_Z q(Z) \log \frac{p(Z\,|\,X,\theta)}{q(Z)}
$$
とおく.
$\KL(q\,\|\,p)$は$q(Z)$と事後分布$p(Z\,|\,X,\theta)$との距離なので常に$0$以上
（3章のカルバック距離を参照）．
$$
\calL(q,\theta)+\KL(q\,\|\,p)
 = \sum_Z q(Z) \log \frac{p(X,Z\,|\,\theta)}{p(Z\,|\,X,\theta)}
 = \sum_Z q(Z) \log p(X\,|\,\theta)
 = \log p(X\,|\,\theta).
$$
よって
$$
\log p(X\,|\,\theta)
 = \calL(q,\theta)+\KL(q\,\|\,p)
 \ge \calL(q,\theta).
$$
したがって$\calL(q,\theta)$は$\log p(X\,|\,\theta)$の下界.
パラメータの現在の値が$\theta^o$だったときに
\begin{itemize}
\item[] Eステップでは$\theta^o$を固定して$\calL(q,\theta)$を$q(Z)$について最大化する.
$\log p(X\,|\,\theta)$は$q$によらないのでそれは$\KL=0$のとき, つまり
$$
q(Z)=p(Z\,|\,X,\theta^o)
$$
のときである.
\item[] Mステップでは$q(Z)$を固定して$\calL(q,\theta)$を$\theta$について最大化する.
その$\theta$を$\theta^n$とする.
最大値になっていなければ, 必ず$\calL$が増加し, $\log p(X\,|\,\theta)$も増える.
このときの$\KL(q\,\|\,p)$は$\theta^o$を使って計算されていた（そして値は$0$）ので新しい$\theta^n$を使って計算し直すと通常正となる.
\end{itemize}
$q(Z)=p(Z\,|\,X,\theta^o)$より
\begin{eqnarray*}
\calL(q,\theta) &=& \sum_ Z q(Z) \log \frac{p(X,Z\,|\,\theta)}{q(Z)}\\
     &=& \sum_Z p(Z\,|\,X,\theta^o) \log p(X,Z\,|\,\theta) - \sum_Z p(Z\,|\,X,\theta^o)\log p(Z\,|\,X,\theta^o)\\
     && （\calQ(\theta,\theta^o)=\sum_Z p(Z\,|\,X,\theta^o) \log p(X,Z\,|\,\theta){\text とおいて}）\\
     &=& \calQ(\theta,\theta^o) + \theta{\text {に非依存}}.
\end{eqnarray*}
つまり$\calL(q,\theta)$の最大化は$\calQ(\theta,\theta^o)$の最大化に等しい.

\section{混合ガウス分布のオンライン版EMアルゴリズム}
各EMのステップで一つのデータ点のみの更新を行うことを考える. これは$m$番目のデータ以外を潜在変数とするEMアルゴリズムとみなすことができる.

Eステップでは分布$p(Z\,|\,X,\theta^o)$を求める必要があるが, Mステップで必要な$\mu_k$, $\Sigma_k$, $\pi_k$の更新式の右辺を見ると必要なデータは$\gamma(z_{nk})$のみであることが分かる. つまりそれらの差分さえ分かればアルゴリズムを書き下すことができる.
$$N_k=\sum_n \gamma(z_{nk})$$
を$m$番目の値だけ更新する. 新しい値を$N_k'$とすると
$$
N_k'=\sum_{n \ne m} \gamma(z_{nk}) + \gamma'(z_{mk})=N_k + \gamma'(z_{mk})-\gamma(z_{mk}).
$$
$d:=\gamma'(z_{mk})-\gamma(z_{mk})$とおくと$N_k'=N_k + d$.
$\pi_k=N_k/N$なので
$$
\pi_k'=\frac{N_k'}{N}=\frac{N_k+d}{N_k}=\pi_k + \frac{d}{N}.
$$
$\mu_k=(1/N_k)\sum_n \gamma(z_{nk})x_n$より
$$
\mu_k'=\frac{1}{N_k'}\left(\sum_{n \ne m}\gamma(z_{nk})x_n+\gamma'(z_{mk})x_m\right).
$$
よって
$$
N_k' \mu_k'=N_k \mu_k - \gamma(z_{mk})x_m + \gamma'(z_{mk})x_m=(N_k'-d)\mu_k+dx_m=N_k'\mu_k+d(x_m-\mu_k).
$$
よって
$$
\mu_k'=\mu_k+\frac{d}{N_k'}(x_m-\mu_k).
$$
$S:=\Sigma_k=(1/N_k)\sum_n \gamma(z_{nk})\outp{(x_n-\mu_k)}$より$S':=\Sigma_k'$とすると
$$
N_k'S'=\sum_{n \ne m}\gamma(z_{nk})\outp{(x_n-\mu_k')}+\gamma'(z_{mk})\outp{(x_m-\mu_k')}.
$$
以下式変形をひたすら行う.
$$
x_m-\mu_k'=x_m-\mu_k-\frac{d}{N_k'}(x_m-\mu_k)=\left(1-\frac{d}{N_k'}\right)(x_m-\mu_k)=\frac{N_k}{N_k'}(x_m-\mu_k).
$$
$$
x_n-\mu_k'=(x_n-\mu_k)-\frac{d}{N_k'}(x_m-\mu_k).
$$
$A:=\outp{(x_m-\mu_k)}$とおくと
\begin{eqnarray*}
&&\outp{(x_n-\mu_k')}=\outp{(x_n-\mu_k)}-2\frac{d}{N_k'}(x_n-\mu_k)\trans{(x_m-\mu_k)}+\frac{d^2}{N_k'^2}A,
\\
&&\sum_{n \ne m}\gamma(z_{nk})\outp{(x_n-\mu_k)}
=N_k S - \gamma(z_{mk})\outp{(x_m-\mu_k)}
\\
&&\hphantom{\sum_{n \ne m}\gamma(z_{nk})\outp{(x_n-\mu_k)}}
= N_k S -\gamma(z_{mk})A,
\\&&
\sum_{n \ne m}\gamma(z_{nk})(x_n-\mu_k)
 = \sum_{n \ne m}\gamma(z_{nk})x_n - \sum_{n \ne m}\gamma(z_{nk})\mu_k\\
&&\hphantom{\sum_{n \ne m}\gamma(z_{nk})(x_n-\mu_k)}
 = N_k \mu_k - \gamma(z_{mk})x_m-(N_k-\gamma(z_{mk}))\mu_k\\
&&\hphantom{\sum_{n \ne m}\gamma(z_{nk})(x_n-\mu_k)}
 = -\gamma(z_{mk})(x_m-\mu_k).
\end{eqnarray*}
よって$\gamma:=\gamma(z_{mk})$とおくと
\begin{eqnarray*}
N_k' S'
 &=& N_k S - \gamma A + 2\frac{d}{N_k'}\gamma A + (N_k-\gamma)\frac{d^2}{N_k'^2}A + (\gamma+d)A \frac{N_k^2}{N_k'^2}\\
 &=& N_k S + \frac{A}{N_k'^2}\left(-\gamma(N_k+d)^2+2d\gamma(N_k+d)+(N_k-\gamma)d^2+(\gamma+d)N_k^2\right)\\
 &=& N_k S + \frac{A}{N_k'^2}\left(-\gamma N_k^2-2d\gamma N_k - \gamma d^2 + 2d \gamma N_k  \right.\\
 & & \hphantom{N_k S + \frac{A}{N_k'^2}}\quad \left. + 2d^2 \gamma + N_k d^2-\gamma d^2 + \gamma N_k^2 + d N_k^2\right)\\
 &=& N_k S + \frac{A}{N_k'^2}N_k d(N_k+d) = N_k S + \frac{A N_k d}{N_k'}.
\end{eqnarray*}
よって
$$
S'=\frac{N_k}{N_k'}\left(S + \frac{d}{N_k'}\outp{(x_m-\mu_k)}\right).
$$