メニュー3 - critter_eng @ ウィキ

\chapter{近似推論法}

\section{変分推論}
モデルパラメータ（および9章での潜在変数）を$\bm{Z}$とし、観測変数を$\bm{X}$と書く。
確率モデルによって同時分布$p(\bm{X},\bm{Z})$が定められ、周辺分布の対数は任意の分布$q(\bm{Z})$を用いて
\begin{eqnarray}
\ln p(\bm{X}) &=& \mathcal{L}(q) + \mathrm{KL}(q||p) \notag \\
\mathcal{L}(q) &=& \int q(\bm{Z}) \ln \left \{ \frac{p(\bm{X},\bm{Z})}{q(\bm{Z})}\right \} d\bm{Z} \notag \\
\mathrm{KL}(q||p) &=& - \int q(\bm{Z}) \ln \left \{ \frac{p(\bm{Z}|\bm{X})}{q(\bm{Z})} \right \} d\bm{Z}
\end{eqnarray}
と分解することができる。

\subsection{分布の分解}
ここでは$\bm{Z}$の要素をいくつかの排反なグループに分割し、$\bm{Z}_{i}\ (i=1,\cdots,M)$と書き
\begin{eqnarray}
q(\bm{Z}) = \prod_{i=1}^{M}q_{i}(\bm{Z}_{i})
\end{eqnarray}
と分解できると仮定する。
$\mathcal{L}(q)$に対する一つの因子に注目するため、$q_{j}(\bm{Z}_{j})$に対する依存項を取り出してみると、$q_{j}(\bm{Z}_{j})$を$q_{j}$と書き
\begin{eqnarray}
\mathcal{L}(q) &=& \int \prod_{i}q_{i} \left \{ \ln p(\bm{X},\bm{Z}) - \sum_{i'}\ln q_{i'} \right \} d\bm{Z} \notag \\
&=& \int q_{j} \left \{ \int \ln p(\bm{X},\bm{Z}) \prod_{i\neq j}q_{i} d\bm{Z}_{i} \right \} d\bm{Z}_{j} - \int q_{j} \ln q_{j}d\bm{Z}_{j}
- \sum_{i\neq j}\int q_{i} \ln q_{i} d\bm{Z}_{i} \notag \\
&=& \int q_{j} \ln \tilde{p}(\bm{X},\bm{Z}_{j})d\bm{Z}_{j} - \int q_{j} \ln q_{j} d\bm{Z}_{j} + \mathrm{const}
\end{eqnarray}
を得る。ここで
\begin{eqnarray}
\ln \tilde{p}(\bm{X},\bm{Z}_{j}) &=& \mathbb{E}_{i\neq j}[\ln p(\bm{X},\bm{Z})] + \mathrm{const} \notag \\
\mathbb{E}_{i\neq j}[\ln p(\bm{X},\bm{Z})] &=& \int \ln p(\bm{X},\bm{Z}) \prod_{i\neq j}q_{i} d\bm{Z}_{i}
\end{eqnarray}
である。
\textcolor{blue}{$\mathrm{const}$は$\tilde{p}(\bm{X},\bm{Z}_{j})$が確率分布として規格化されていることによる。}
$\mathcal{L}(q)$を$q_{j}(\bm{Z}_{j})$について最大化することを考えると、これは$\mathcal{L}(q)$が$q_{j}(\bm{Z}_{j})$と$\tilde{p}(\bm{X},\bm{Z}_{j})$の間の負のKLダイバージェンスであることに注目すると、$q_{j}(\bm{Z}_{j}) = \tilde{p}(\bm{X},\bm{Z}_{j})$が得られることがわかる。よって最適解$q_{j}^{\star}(\bm{Z}_{j})$は一般に
\begin{eqnarray}
\ln q_{j}^{\star}(\bm{Z}_{j}) = \mathbb{E}_{i\neq j}[\ln p(\bm{X},\bm{Z})] + \mathrm{const}
\end{eqnarray}
を満たす。定数項は$q_{j}^{\star}(\bm{Z}_{j})$の正規化に起因する。すなわち
\begin{eqnarray}
q_{j}^{\star}(\bm{Z}_{j}) = \frac{\exp(\mathbb{E}_{i\neq j}[\ln p(\bm{X},\bm{Z})])}{\int \exp(\mathbb{E}_{i\neq j}[\ln p(\bm{X},\bm{Z})]) d\bm{Z}_{j}}
\end{eqnarray}
である。

\subsection{分解による近似のもつ性質}
省略

\textcolor{blue}{
逆のKLダイバージェンス$KL(p||q)$を最小化する話があるが、どうしてそれを考えるのかいまひとつよくわからない。
$\ln p(\bm{X})$と関係のある量なのか？}

\subsection{例：一変数ガウス分布}
省略

\subsection{モデル比較}
隠れ変数$\bm{Z}$だけでなく、変数$m$で表され、事前確率$p(m)$を持つ複数のモデルの候補を比較する場合を考える。
この場合、$q(\bm{Z},m) = q(\bm{Z}|m)q(m)$を考える必要があり、
\begin{eqnarray}
\ln p(\bm{X}) &=& \mathcal{L} - \sum_{m}\sum_{\bm{Z}}q(\bm{Z}|m)q(m) \ln \left \{ \frac{p(\bm{Z},m|\bm{X})}{q(\bm{Z}|m)q(m)} \right \} \notag \\
\mathcal{L} &=& \sum_{m}\sum_{\bm{Z}}q(\bm{Z}|m)q(m) \ln \left \{ \frac{p(\bm{Z},\bm{X},m)}{q(\bm{Z}|m)q(m)} \right \}
\end{eqnarray}
と分解する。$\mathcal{L}$を$q(m)$について最大化すると
\begin{eqnarray}
q(m) &\propto& p(m)\exp\{ \mathcal{L}_{m} \} \notag \\
\mathcal{L}_{m} &=& \sum_{\bm{Z}}q(\bm{Z}|m) \ln \left \{ \frac{p(\bm{Z},\bm{X}|m)}{q(\bm{Z}|m)} \right \}
\end{eqnarray}
を得る。
\textcolor{blue}{
これは
\begin{eqnarray}
f = \sum_{\bm{Z}}q(\bm{Z}|m) \ln \left \{ \frac{p(\bm{Z},\bm{X},m)}{q(\bm{Z}|m)q(m)} \right \} + \lambda \left[ \sum_{m}q(m)-1 \right]
\end{eqnarray}
として$q(m)$による微分を考えると
\begin{eqnarray}
\frac{\partial f}{\partial q(m)} &=& \sum_{\bm{Z}}q(\bm{Z}|m) \ln \left \{ \frac{p(\bm{Z},\bm{X},m)}{q(\bm{Z}|m)q(m)} \right \} -\sum_{\bm{Z}}q(\bm{Z}|m) + \lambda \notag \\
&=& \ln p(m) - \ln q(m) + \sum_{\bm{Z}} q(\bm{Z}|m) \ln \left \{ \frac{p(\bm{Z},\bm{X}|m)}{q(\bm{Z}|m)} \right \} - 1 + \lambda = 0 \notag \\
\end{eqnarray}
となることからわかる。
}

\section{例：変分混合ガウス分布}
ここでは変分推論法を癌す混合モデルに適用することを考える。
観測データの集合を$\bm{X}=\{\bm{x}_{1},\cdots,\bm{x}_{N}\}$とし、対応する潜在変数を$\bm{Z}=\{\bm{z}_{1},\cdots,\bm{z}_{N}\}$とする。
ここで$\bm{z}_{n}$は$(k=1,\cdots,K)$の中に1が一つだけある二値ベクトルである。混合比$\bm{\pi}$が与えられたときの$\bm{Z}$の条件付き分布は
\begin{eqnarray}
p(\bm{Z}|\bm{\pi}) = \prod_{n=1}^{N}\prod_{k=1}^{K} \pi_{k}^{z_{nk}}
\end{eqnarray}
で与えられる。
また、潜在変数と混合要素のパラメータが与えられたときの観測データベクトルの条件付き分布は
\begin{eqnarray}
p(\bm{X}|\bm{Z},\bm{\mu},\bm{\Lambda}) = \prod_{n=1}^{N}\prod_{k=1}^{K} \mathcal{N}(\bm{x}_{n}|\bm{\mu}_{k},\bm{\Lambda}_{k}^{-1})^{z_{nk}}
\end{eqnarray}
となる。

次にパラメータ$\bm{\mu},\bm{\Lambda},\bm{\pi}$の事前分布を導入する。
混合比$\bm{\pi}$にはディレクレ分布を、平均と精度にはガウス‐ウィシャート事前分布を用いると
\begin{eqnarray}
p(\bm{\pi}) &=& \mathrm{Dir}(\bm{\pi}|\bm{\alpha_{0}}) = C(\bm{\alpha}_{0})\prod_{k=1}^{K}\pi_{k}^{\alpha_{0}-1} \notag \\
p(\bm{\mu},\bm{\Lambda}) &=& p(\bm{\mu}|\bm{\Lambda})p(\bm{\Lambda}) \notag \\
&=& \prod_{k=1}^{K}\mathcal{N}(\bm{\mu_{k}}|\bm{m}_{0},(\beta_{0}\bm{\Lambda})^{-1}) \mathcal{W}(\bm{\Lambda}_{k}|\bm{W}_{0},\nu_{0})
\end{eqnarray}
となる。通常$\bm{m}_{0}=\bm{0}$とおく。

\subsection{変分事後分布}
変分ベイズ法で扱うための、全ての確率変数の同時分布は
\begin{eqnarray}
p(\bm{X},\bm{Z},\bm{\pi},\bm{\mu},\bm{\Lambda}) = p(\bm{X}|\bm{Z},\bm{\mu},\bm{\Lambda}) p(\bm{Z}|\bm{\pi})p(\bm{\pi})p(\bm{\mu}|\bm{\Lambda})p(\bm{\Lambda})
\end{eqnarray}
となる。

ここで潜在変数とパラメータに分解した変分近似
\begin{eqnarray}
q(\bm{Z},\bm{\pi},\bm{\mu},\bm{\Lambda}) = q(\bm{Z})q(\bm{\pi},\bm{\mu},\bm{\Lambda})
\end{eqnarray}
を考える。本文(10.9)式を用いると
\begin{eqnarray}
\ln q^{*}(\bm{Z}) &=& \mathbb{E}_{\bm{\pi},\bm{\mu},\bm{\Lambda}}[\ln p(\bm{X},\bm{Z},\bm{\pi},\bm{\mu},\bm{\Lambda}) ] + \mathrm{const} \notag \\
&=& \mathbb{E}_{\bm{\pi}}[\ln p(\bm{Z}|\bm{\pi})] + \mathbb{E}_{\bm{\mu},\bm{\Lambda}}[ \ln p(\bm{X}|\bm{Z},\bm{\mu},\bm{\Lambda})] + \mathrm{const}
\end{eqnarray}
となり、具体的な関数形を代入すると
\begin{eqnarray}
\ln q^{*}(\bm{Z}) &=& \sum_{n=1}^{N}\sum_{k=1}^{K}z_{nk}\ln \rho_{nk} + \mathrm{const} \notag \\
\ln \rho_{nk} &=& \mathbb{E}[\ln \pi_{k}] + \frac{1}{2}\mathbb{E}[\ln |\bm{\Lambda}_{k}|] - \frac{D}{2}\ln(2\pi) \notag \\
&-& \frac{1}{2}\mathbb{E}_{\bm{\mu}_{k},\bm{\Lambda}_{k}}[(\bm{x}_{n}-\bm{\mu}_{k})^{T}\bm{\Lambda}_{k}(\bm{x}_{n}-\bm{\mu}_{k})]
\end{eqnarray}
となるため、指数をとって
\begin{eqnarray}
q^{*}(\bm{Z}) \propto \prod_{n=1}^{N}\prod_{k=1}^{K}\rho_{nk}^{z_{nk}}
\end{eqnarray}
となり、規格化を行うと
\begin{eqnarray}
q^{*}(\bm{Z}) &=& \prod_{n=1}^{N}\prod_{k=1}^{K}r_{nk}^{z_{nk}} \notag \\
r_{nk} &=& \frac{\rho_{nk}}{\sum_{j=1}^{K}\rho_{nj}}
\end{eqnarray}
となる。
\textcolor{blue}{
これは
\begin{eqnarray}
\sum_{\bm{Z}} \prod_{n=1}^{N}\prod_{k=1}^{K}\rho_{nk}^{z_{nk}} = \prod_{n=1}^{N} \left( \sum_{k=1}^{K}\rho_{nk} \right)
\end{eqnarray}
による。}
これより$q^{*}(\bm{Z})$については
\begin{eqnarray}
\mathbb{E}[z_{nk}] = r_{nk}
\end{eqnarray}
が成り立つ。
ここで今後のために
\begin{eqnarray}
N_{k} &=& \sum_{n=1}^{N}r_{nk} \notag \\
\bar{\bm{x}}_{k} &=& \frac{1}{N_{k}} \sum_{n=1}^{N}r_{nk}\bm{x}_{n} \notag \\
\bm{S}_{k} &=& \frac{1}{N_{k}} \sum_{n=1}^{N}r_{nk}(\bm{x}_{n}-\bar{\bm{x}}_{k})(\bm{x}_{n}-\bar{\bm{x}}_{k})^{T}
\end{eqnarray}
を定義する。

次に$q(\bm{\pi},\bm{\mu},\bm{\Lambda})$について考える。
本文(10.9)を用いると
\begin{eqnarray}
\ln q^{*}(\bm{\pi},\bm{\mu},\bm{\Lambda}) &=& \ln p(\bm{\pi}) + \sum_{k=1}^{K}\ln p(\bm{\mu}_{k},\bm{\Lambda}_{k}) + \mathbb{E}_{\bm{Z}}[\ln p(\bm{Z}|\bm{\pi})]
\notag \\
&+& \sum_{k=1}^{K}\sum_{n=1}^{N}\mathbb{E}[z_{nk}]\ln \mathcal{N}(\bm{x}_{n}|\bm{\mu}_{k},\bm{\Lambda}_{k}^{-1}) + \mathrm{const}
\end{eqnarray}
を得る。これは
\begin{eqnarray}
q^{*}(\bm{\pi},\bm{\mu},\bm{\Lambda}) = q^{*}(\bm{\pi}) \prod_{k=1}^{K}q^{*}(\bm{\mu}_{k},\bm{\Lambda}_{k})
\end{eqnarray}
と分解できることを意味する。変分事後分布の$\bm{\pi}$に依存する項を取り出すと
\begin{eqnarray}
\ln q^{*}(\bm{\pi}) = (\alpha_{0}-1)\sum_{k=1}^{K}\ln \pi_{k} + \sum_{k=1}^{K}\sum_{n=1}^{N}r_{nk}\ln\pi_{k} + \mathrm{const}
\end{eqnarray}
となり、指数をとることで
\begin{eqnarray}
q^{*}(\bm{\pi}) = \mathrm{Dir}(\bm{\pi}|\bm{\alpha})
\end{eqnarray}
とディリクレ分布になる。ここで$\bm{\alpha}$の要素は
\begin{eqnarray}
\alpha_{k} = \alpha_{0} + N_{k}
\end{eqnarray}
である。最後に、$q^{*}(\bm{\mu}_{k},\bm{\Lambda}_{k})$について考えると、これは
\begin{eqnarray}
q^{*}(\bm{\mu}_{k},\bm{\Lambda}_{k}) = \mathcal{N}(\bm{\mu_{k}}|\bm{m}_{k},(\beta_{k}\bm{\Lambda}_{k})^{-1})\mathrm{W}(\bm{\Lambda}_{k}|\bm{W}_{k},\nu_{k})
\end{eqnarray}
となる。ただし、
\begin{eqnarray}
\beta_{k} &=& \beta_{0} + N_{k} \notag \\
\bm{m}_{k} &=& \frac{1}{\beta_{k}}(\beta_{0}\bm{m}_{0}+N_{k}\bar{\bm{x}}_{k}) \notag \\
\bm{W}_{k}^{-1} &=& \bm{W}_{0}^{-1} + N_{k}\bm{S}_{k} + \frac{\beta_{0}N_{k}}{\beta_{0}+N_{k}}(\bar{\bm{x}}_{k}-\bm{m}_{0})(\bar{\bm{x}}_{k}-\bm{m}_{0})^{T} \notag \\
\nu_{k} &=& \nu_{0} + N_{k}
\end{eqnarray}
である。これを用いると、$\mathbb{E}[z_{nk}]=r_{nk}$の計算に必要な量が計算できて
\begin{eqnarray}
\mathbb{E}_{\bm{\mu}_{k},\bm{\Lambda}_{k}}[(\bm{x}_{n}-\bm{\mu}_{k})^{T}\bm{\Lambda}_{k}(\bm{x}_{n}-\bm{\mu}_{k})]
&=& D\beta_{k}^{-1} + \nu_{k}(\bm{x}_{n}-\bm{m}_{k})^{T}\bm{W}_{k}(\bm{x}_{n}-\bm{m}_{k}) \notag \\
\ln \tilde{\Lambda}_{k} \equiv \mathbb{E}[\ln |\bm{\Gamma}_{k}|] &=& \sum_{i=1}^{D}\psi\left(\frac{\nu_{k}+1-i}{2} \right) + D\ln 2 + \ln|\bm{W}_{k}| \notag \\
\ln \tilde{\pi}_{k} \equiv \mathbb{E}[\ln \pi_{k}] &=& \psi(\alpha_{k}) - \psi\left(\sum_{k}\alpha_{k}\right)
\end{eqnarray}
を得る。実際にはこれらを交互に繰り返すことになる。

\subsection{変分下界}
省略

\subsection{予測分布}
混合ガウスモデルの新しい観測値$\hat{\bm{x}}$の予測分布を考える。
この観測値には、対応する潜在変数$\hat{\bm{z}}$が存在し、
\begin{eqnarray}
p(\hat{\bm{x}}|\bm{X}) = \sum_{\hat{\bm{z}}}\int\int\int p(\hat{\bm{x}}|\hat{\bm{z}},\bm{\mu},\bm{\Lambda}) p(\hat{\bm{z}}|\bm{\pi})
p(\bm{\pi},\bm{\mu},\bm{\Lambda}|\bm{X})d\bm{\pi}d\bm{\mu}d\bm{\Lambda}
\end{eqnarray}
で与えられ、モデルの具体形を用いると
\begin{eqnarray}
p(\hat{\bm{x}}|\bm{X}) = \sum_{k=1}^{K}\int\int\int \pi_{k}\mathcal{N}(\hat{\bm{x}}|\bm{\mu}_{k},\bm{\Lambda}^{-1})p(\bm{\pi},\bm{\mu},\bm{\Lambda}|\bm{X})d\bm{\pi}d\bm{\mu}d\bm{\Lambda}
\end{eqnarray}
となる。事後分布$p(\bm{\pi},\bm{\mu},\bm{\Lambda}|\bm{X})$を変分近似$q(\bm{\pi})q(\bm{\mu},\bm{\Lambda})$で置き換えると
\begin{eqnarray}
p(\hat{\bm{x}}|\bm{X}) &\sim&
\sum_{k=1}^{K}\int\int\int \pi_{k} \mathcal{N}(\hat{\bm{x}}|\bm{\mu}_{k},\bm{\Lambda}^{-1}) q(\bm{\pi})q(\bm{\mu}_{k},\bm{\Lambda}_{k})
d\bm{\pi}d\bm{\mu}_{k}d\bm{\Lambda}_{k} \notag \\
&\sim& \frac{1}{\hat{\alpha}}\sum_{k=1}^{K}\alpha_{k}\mathrm{St}(\hat{\bm{x}}|\bm{m}_{k},\bm{L}_{k},\nu_{k}+1-D)
\end{eqnarray}
と混合スチューデント$t$分布となる。ここで
\begin{eqnarray}
\bm{L}_{k} = \frac{(\nu_{k}+1-D)\beta_{k}}{1+\beta_{k}}\bm{W}_{k}
\end{eqnarray}
である。

\subsection{混合要素数の決定}
省略

\subsection{導出された分解}
省略

\section{変分線形回帰}
省略

\section{指数型分布族}
省略

\section{局所変分推論法}
$f(x)=\exp(-x)$は$x$について凸関数であり、$x=\xi$での接線は
\begin{eqnarray}
y(x) = f(\xi) + f'(\xi)(x-\xi)
\end{eqnarray}
である。今の場合
\begin{eqnarray}
y(x) = \exp(-\xi) - \exp(-\xi)(x-\xi)
\end{eqnarray}
の形となっている。$\eta = -\exp(-\xi)$と定義すると
\begin{eqnarray}
y(x,\eta) = \eta x - \eta + \eta \ln (-\eta)
\end{eqnarray}
となる。これらは異なる$\eta$が異なる接線に対応していて、$f(x)\geq y(x,\eta)$となっているため
\begin{eqnarray}
f(x) = \max_{\eta}\{ \eta x - \eta + \eta \ln (-\eta) \}
\end{eqnarray}
が成り立つ。

このような法則はより一般に凸相対性として成り立つ。
すなわち下に凸な関数の接線の方程式を$y=\eta x -g(\eta)$と書くと、
\begin{eqnarray}
g(\eta) = \max_{x}\{ \eta x -f(x)\}
\end{eqnarray}
を得る。また、ある$x$での$y$座標は、その点が接点となるときに最大化されるから、
\begin{eqnarray}
f(x) = \max_{\eta}\{\eta x -g(\eta)\}
\end{eqnarray}
が成り立つ。

以下省略

\section{変分ロジスティック回帰}
\subsection{変分事後分布}
ここでは、変分法の枠組みでロジスティック回帰モデルを考える。データ$\bm{t}$が与えられた場合の尤度関数は
\begin{eqnarray}
p(\bm{t}) = \int p(\bm{t}|\bm{w})p(\bm{w}) d\bm{w} = \int \left[ \prod_{n=1}^{N}p(t_{n}|\bm{w}) \right] p(\bm{w})d\bm{w}
\end{eqnarray}
となる。ここで
\begin{eqnarray}
p(\bm{t}|\bm{w}) &=& \sigma(a)^{t}\{ 1-\sigma(a) \}^{1-t} \notag \\
&=& \left( \frac{1}{1+e^{-a}} \right) \left( 1 - \frac{1}{1+e^{-a}} \right)^{1-t} \notag \\
&=& e^{at}\frac{e^{-a}}{1+e^{-a}} = e^{at}\sigma(-a)
\end{eqnarray}
であり、$a=\bm{w}^{T}\bm{\phi}$である。
シグモイド関数の一般的な性質
\begin{eqnarray}
\sigma(z) &\geq& \sigma(\xi) \exp \{ (z-\xi)/2 - \lambda(\xi)(z^{2}-xi^{2}) \} \notag \\
\lambda(\xi) &=& \frac{1}{2\xi} \left[ \sigma(\xi) - \frac{1}{2} \right]
\end{eqnarray}
を用いると
\begin{eqnarray}
p(t|\bm{w}) = e^{at}\sigma(-a) \geq e^{at}\sigma(\xi)\exp\{-(a+\xi)/2 - \lambda(\xi)(a^{2}-\xi^{2}) \}
\end{eqnarray}
となる。したがって、同時分布についての下界
\begin{eqnarray}
p(\bm{t},\bm{w}) = p(\bm{t}|\bm{w})p(\bm{w}) \geq h(\bm{w},\bm{xi})p(\bm{w})
\end{eqnarray}
を得る。ただし、$\bm{\xi}$は変分パラメータの集合$\{\xi_{n}\}$を意味し
\begin{eqnarray}
h(\bm{w},\bm{\xi}) = \prod_{n=1}^{N}\sigma(\xi_{n})\exp\{\bm{w}^{T}\bm{\phi}_{n}t_{n}-(\bm{w}^{T}\bm{\phi}_{n}+\xi_{n})/2
- \lambda(\xi_{n})(|\bm{w}^{T}\bm{\phi}_{n}|^2- \xi_{n}^{2}) \} \notag \\
\end{eqnarray}
である。対数を考えると
\begin{eqnarray}
\ln\{ p(\bm{t}|\bm{w})p(\bm{w}) \} \geq \ln p(\bm{w}) + \sum_{n=1}^{N} \{ \ln \sigma(\xi_{n}) + \bm{w}^{T}\bm{\phi}_{n}t_{n} \notag \\
-(\bm{w}^{T}\bm{\phi}_{n}+\xi_{n})/2 - \lambda(\xi_{n})(|\bm{w}^{T}\bm{\phi}_{n}|^2- \xi_{n}^{2}) \}
\end{eqnarray}
となり、事前分布$p(\bm{w})$の値を考えると$\bm{w}$の関数として
\begin{eqnarray}
&-&\frac{1}{2}(\bm{w}-\bm{m}_{0})^{T}\bm{S}_{0}^{-1}(\bm{w}-\bm{m}_{0}) \notag \\
&+&\sum_{n=1}^{N} \left\{ \bm{w}^{T}\bm{\phi}_{n}(t_{n}-1/2) - \lambda(\xi_{n})\bm{w}^{T}(\bm{\phi}_{n}\bm{\phi}_{n}^{T})\bm{w}  \right\} + \mathrm{const}
\end{eqnarray}
となる。これは$\bm{w}$の二次関数であり、変分近似は
\begin{eqnarray}
q(\bm{w}) &=& \mathcal{N}(\bm{w}|\bm{m}_{N},\bm{S}_{N}) \notag \\
\bm{m}_{N} &=& \bm{S}_{N} \left( \bm{S}_{0}^{-1}\bm{m}_{0} + \sum_{n=1}^{N}(t_{n}-1/2)\bm{\phi}_{n} \right) \notag \\
\bm{S}_{N}^{-1} &=& \bm{S}_{0}^{-1} + 2\sum_{n=1}^{N}\lambda(\xi_{n})\bm{\phi}_{n}\bm{\phi}_{n}^{T}
\end{eqnarray}
となる。
\textcolor{blue}{
本文に、「正規化して変分事後分布$q(\bm{w})$にしてしまうと、もう下界ではなくなってしまう」
とあるが、その観点からしてこの節の取り扱いは大丈夫なのだろうか」
}

\section{変分パラメータの最適化}
ここでは変分パラメータ$\{\xi_{n}\}$の決め方について考える。
周辺尤度については
\begin{eqnarray}
\ln p(\bm{t}) = \ln \int p(\bm{t}|\bm{w})p(\bm{w}) \geq \ln \int h(\bm{w},\bm{\xi})p(\bm{w}) d\bm{w} = \mathcal{L}(\bm{\xi})
\end{eqnarray}
が成り立つ。
最初にこれをEMアルゴリズムを用いて最大化することを考える。
その場合、$\bm{\xi}^{\mathrm{old}}$に対して、
\begin{eqnarray}
Q(\bm{\xi},\bm{\xi}^{\mathrm{old}}) = \mathrm{E}[\ln\{h(\bm{w},\bm{\xi})p(\bm{w})\}]
\end{eqnarray}
を最大化することを考える。これは
\begin{eqnarray}
Q(\bm{\xi},\bm{\xi}^{\mathrm{old}}) = \sum_{n=1}^{N} \left\{ \ln\sigma(\xi_{n})-\xi_{n}/2-\lambda(\xi_{n})(\bm{\phi}_{n}^{T}\mathbb{E}[\bm{w}\bm{w}^{T}]\bm{\phi}_{n} - \xi_{n}^{2}) \right\} + \mathrm{const}
\end{eqnarray}
で与えられる。ここでconst項は$\bm{\xi}$に依存しない項を表す。
これを$\xi_{n}$について微分しそれを$0$とおくと
\begin{eqnarray}
0 = \lambda'(\xi_{n})(\bm{\phi}^{T}\mathbb{E}[\bm{w}\bm{w}^{T}]\bm{\phi}_{n}-\xi_{n}^{2})
\end{eqnarray}
を得る。$\lambda(\xi)$が単調関数であることから$\lambda'(\xi)\neq 0$であり、
\begin{eqnarray}
(\xi_{n}^{\mathrm{new}})^2 = \bm{\phi}^{T}\mathbb{E}[\bm{w}\bm{w}^{T}]\bm{\phi}_{n} = \bm{\phi}^{T}(\bm{S}_{N} + \bm{m}_{N}\bm{m}_{N}^{T})\bm{\phi}_{n}
\end{eqnarray}
を得る。

一方で$L(\bm{\xi})$を直接計算することも可能で
\begin{eqnarray}
L(\bm{\xi}) &=& \frac{1}{2} \ln \frac{|\bm{S}_{N}|}{|\bm{S}_{0}|} + \frac{1}{2}\bm{m}_{N}^{T}\bm{S}_{N}^{-1}\bm{m}_{N} - \frac{1}{2}\bm{m}_{0}^{T}\bm{S}_{0}^{-1}\bm{m}_{0} \notag \\
&+& \sum_{n=1}^{N} \left\{ \ln\sigma(\xi_{n}) - \frac{1}{2}\xi_{n} + \lambda(\xi_{n})\xi_{n}^{2}\right\}
\end{eqnarray}
を得る。

\subsection{超パラメータの推論}
パラメータ$\bm{w}$の事前分布を
\begin{eqnarray}
p(\bm{w}|\alpha) = \mathcal{N}(\bm{w}|\bm{0},\alpha^{-1}\bm{I})
\end{eqnarray}
し、$\alpha$に対する共役超事前分布を
\begin{eqnarray}
p(\alpha) = \mathrm{Gam}(\alpha|a_{0},b_{0})
\end{eqnarray}
とする。このモデルの周辺尤度は
\begin{eqnarray}
p(\bm{t}) &=& \int\int p(\bm{w},\alpha,\bm{t})d\bm{w}d\alpha \notag \\
p(\bm{w},\alpha,\bm{t}) &=& p(\bm{t}|\bm{w})p(\bm{w}|\alpha)p(\alpha)
\end{eqnarray}
で与えられる。
変分分布$q(\bm{w},\alpha)$を導入すると
\begin{eqnarray}
\ln p(\bm{t}) &=& \mathcal{L}(q) + \mathrm{KL}(q||p) \notag \\
\mathrm{L}(q) &=& \int\int q(\bm{w},\alpha) \ln \left\{ \frac{p(\bm{w},\alpha,\bm{t})}{q(\bm{w},\alpha)} \right\} d\bm{w}d\alpha \notag \\
\mathrm{KL}(q||p) &=& - \int\int q(\bm{w},\alpha) \ln \left\{ \frac{p(\bm{w},\alpha|\bm{t})}{q(\bm{w},\alpha)} \right\} d\bm{w}d\alpha
\end{eqnarray}
となる。本文(10.152)より
\begin{eqnarray}
\ln p(\bm{t}) &\geq& \mathcal{L}(q) \geq \tilde{\mathcal{L}}(q,\bm{\xi}) \notag \\
&=& q(\bm{w},\alpha) \ln \left\{ \frac{h(\bm{w},\bm{\xi})p(\bm{w}|\alpha)p(\alpha)}{q(\bm{w},\alpha)} \right\} d\bm{w}d\alpha
\end{eqnarray}
が成り立つ。変分分布が
\begin{eqnarray}
q(\bm{w},\alpha) = q(\bm{w})q(\alpha)
\end{eqnarray}
と分解できる場合を考える。本文(10.9)を用いると最適な$q(\bm{w})$は
\begin{eqnarray}
\ln q(\bm{w}) &=& \mathbb{E}_{\alpha} [\ln \{h(\bm{w},\bm{\xi})p(\bm{w}|\alpha)p(\alpha)\}] + \mathrm{const} \notag \\
&=& \ln h(\bm{w},\bm{\xi}) + \mathbb{E}_{\alpha}[\ln p(\bm{w}|\alpha)] + \mathrm{const}
\end{eqnarray}
となる。具体的な値を代入すると
\begin{eqnarray}
\ln q(\bm{w}) = - \frac{\mathbb{E}[\alpha]}{2}\bm{w}^{T}\bm{w} + \sum_{n=1}^{N}\left\{(t_{n}-1/2)\bm{w}^{T}\bm{\phi}_{n} - \lambda(\xi_{n})\bm{w}^{T}\bm{\phi}_{n}\bm{\phi}_{n}^{T}\bm{w} \right\} + \mathrm{const}
\end{eqnarray}
となるため
\begin{eqnarray}
q(\bm{w}) &=& \mathcal{N}(\bm{w}|\bm{\mu}_{N},\bm{\Sigma}_{N}) \notag \\
\bm{\Sigma}_{N}^{-1}\bm{\mu}_{N} &=& \sum_{n=1}^{N}(t_{n}-1/2)\bm{\phi}_{n} \notag \\
\bm{\Sigma}_{N}^{-1} &=& \mathbb{E}[\alpha]\bm{I} + 2\sum_{n=1}^{N}\lambda(\xi_{n})\bm{\phi}_{n}\bm{\phi}_{n}^{T}
\end{eqnarray}
を得る。同様に$q(\alpha)$についても
\begin{eqnarray}
\ln q(\alpha) &=& \mathbb{E}_{\bm{w}}[\ln p(\bm{w}|\alpha)] + \ln p(\alpha) + \mathrm{const} \notag \\
&=& \frac{M}{2}\ln\alpha - \frac{\alpha}{2}\mathbb{E}[\bm{w}^{T}\bm{w}] + (a_{0}-1)\ln\alpha - b_{0}\alpha + \mathrm{const}
\end{eqnarray}
となるため
\begin{eqnarray}
q(\alpha) &=& \mathrm{Gam}(\alpha|a_{N},b_{N}) \notag \\
a_{N} &=& a_{0} + \frac{M}{2} \notag \\
b_{N} &=& b_{0} + \frac{1}{2}\mathbb{E}_{\bm{w}}[\bm{w}^{T}\bm{w}]
\end{eqnarray}
を得る。
最後に$\xi_{n}$について考えると、$\tilde{\mathcal{L}}(q,\bm{\xi}) $を最大化することで行うが、$\bm{\xi}$に依存しない項を除き$\alpha$について積分すると
\begin{eqnarray}
\tilde{\mathcal{L}}(q,\bm{\xi})  = \int q(\bm{w}) \ln h(\bm{w},\bm{\xi})d\bm{w} + \mathrm{const}
\end{eqnarray}
となるため、前節の結果を使うことができて
\begin{eqnarray}
(\xi_{n}^{\mathrm{new}})^{2} = \bm{\phi}_{n}^{T}(\bm{\Sigma}_{N}+\bm{\mu}_{N}\bm{\mu}_{N}^{T})\bm{\phi}_{n}
\end{eqnarray}
となる。これにより、三つの量$q(\bm{w}),q(\alpha),\bm{\xi}$を再推定する方程式が得られたことになる。この時に必要となるモーメントは
\begin{eqnarray}
\mathbb{E}[\alpha] &=& \frac{a_{N}}{b_{N}} \notag \\
\mathbb{E}[\bm{w}\bm{w}^{T}] &=& \bm{\Sigma}_{N} + \bm{\mu}_{N}\bm{\mu}_{N}^{T}
\end{eqnarray}
で与えられる。

\section{EP法}
まず$p(\bm{z})$を固定された確率分布としたとき、$KL(p||q)$を$q(\bm{z})$について最小化する問題を考える。
$q(\bm{z})$を指数型分布族とすると
\begin{eqnarray}
q(\bm{z}) = h(\bm{z})g(\bm{\eta})\exp\{ \bm{\eta}^{T}\bm{u}(\bm{z}) \}
\end{eqnarray}
と書くことができる。このとき$\bm{\eta}$の関数としてのKLダイバージェンスは
\begin{eqnarray}
KL(p||q) = -\ln g(\bm{\eta}) - \bm{\eta}^{T}\mathbb{E}_{p(\bm{z})}[\bm{u}(\bm{z})] + \mathrm{const}
\end{eqnarray}
となる。
$\bm{\eta}$についての最小化は上の式の勾配を$0$とおいて
\begin{eqnarray}
-\nabla \ln g(\bm{\eta}) = \mathbb{E}_{p(\bm{z})}[\bm{u}(\bm{z})]
\end{eqnarray}
によって得られる。本文(2.226)の結果から左辺を変形すると
\begin{eqnarray}
\mathbb{E}_{q(\bm{z})}[\bm{u}(\bm{z})] = \mathbb{E}_{p(\bm{z})}[\bm{u}(\bm{z})]
\end{eqnarray}
となる。例えば$q(\bm{z})$がガウス分布$\mathcal{N}(\bm{z}|\bm{\mu},\bm{\Sigma})$の場合、平均$\bm{\mu}$と分散$\bm{\Sigma}$を分布$p(\bm{z})$と一致させることでKLダイバージェンスを最小化することができる。

この結果を用いて近似推論のアルゴリズムの導出を考える。
多くの確率モデルにおいて、データとパラメータの同時分布は
\begin{eqnarray}
p(\mathcal{D},\bm{\theta}) = \prod_{i}f_{i}(\theta)
\end{eqnarray}
と因子の積の形でかける。
\textcolor{blue}{
($f_{i}(\mathcal{D},\theta)$と書いた方が誤解が少ないように思えるが。)
}
この場合、事後分布は
\begin{eqnarray}
p(\bm{\theta}|\mathcal{D}) = \frac{1}{p(\mathcal{D})} \prod_{i}f_{i}(\bm{\theta})
\end{eqnarray}
であり、モデルエビデンスは
\begin{eqnarray}
p(\mathcal{D}) = \int \prod_{i}f_{i}(\bm{\theta})d\bm{\theta}
\end{eqnarray}
となる。

EP法では、このモデルに対して
\begin{eqnarray}
q(\bm{\theta}) = \frac{1}{Z}\prod_{i}\tilde{f}_{i}(\bm{\theta})
\end{eqnarray}
をにより近似することを考える。そして収束するまで以下を繰り返す。
\begin{enumerate}
\item 改良したい因子$\tilde{f}_{j}(\bm{\theta})$を選ぶ
\item
\begin{eqnarray}
q^{\backslash j}(\bm{\theta}) &=& \frac{q(\bm{\theta})}{\tilde{f}_{j}(\bm{\theta})} \notag \\
Z_{j} &=& \int q^{\backslash j}(\bm{\theta}) f_{j}(\bm{\theta}) d\bm{\theta}
\end{eqnarray}
とし、$q^{\mathrm{new}}(\bm{\theta})$を
\begin{eqnarray}
\mathrm{KL}\left( \frac{f_{j}(\bm{\theta})q^{\backslash j}(\bm{\theta})}{Z_{j}} \middle|\middle| q^{\mathrm{new}}(\theta) \right)
\end{eqnarray}
最小化するように定義。
\textcolor{blue}{（本文にはこれが容易に計算できる操作だと仮定するとあるが$f_{j}(\bm{\theta})$がわからないのでなぜこのように仮定できるかがいまひとつ不明）}
\item 新しい因子を
\begin{eqnarray}
\tilde{f}_{j}(\bm{\theta}) = Z_{j}\frac{q^{\mathrm{new}}(\bm{\theta})}{ q^{\backslash j}(\bm{\theta}) }
\end{eqnarray}
と定義する。
\end{enumerate}

\subsection{例：雑音データ問題}
省略

\subsection{グラフィカルモデルとEP法}
省略

\chapter{サンプリング法}
確率分布$p(\bm{z})$に対して、関数$f(\bm{z})$の期待値は
\begin{eqnarray}
\mathbb{E}[f] = \int f(\bm{z})p(\bm{z})d\bm{z}
\end{eqnarray}
で与えられるが、分布$p(\bm{z})$から独立に抽出されたサンプルの集合$\bm{z}^{(l)}(l=1,\cdots,L)$を用いた有限和で
\begin{eqnarray}
\hat{f} = \frac{1}{L}\sum_{l=1}^{L}f(\bm{z}^{(l)})
\end{eqnarray}
と近似することができる。このときその期待値と分散については
\begin{eqnarray}
\hat{f} &=& \frac{1}{L}\sum_{l=1}^{L}f(\bm{z}^{(l)}) \notag \\
\mathrm{var}[\hat{f}] &=& \frac{1}{L}\mathrm{E} \left[ (f-\mathbb{E}[f])^2 \right]
\end{eqnarray}
が成り立つ。

\section{基本的なサンプリングアルゴリズム}

\subsection{標準的な分布}
ここでは区間$(0,1)$で一様に分布する変数$z$を$y=f(z)$と変換した場合の$y$の分布を考える。
$y$の分布は
\begin{eqnarray}
p(y) = p(z) \left| \frac{dz}{dy} \right|
\end{eqnarray}
に従う。このとき、$y$と$z$の関係は
\begin{eqnarray}
z = h(y) \equiv \int^{y}_{-\infty} p(\hat{y})d\hat{y}
\end{eqnarray}
となる。このことから確率分布$p(y)$に従う変数$y$を得たい場合は$y=h^{-1}(z)$とすればよいことがわかる。

\subsection{棄却サンプリング}
標準的でない分布$p(\bm{z})$に従うサンプルを得たいが、直接$p(\bm{z})$からサンプリングすることが困難な場合を考える。
また与えられた$\bm{z}$の値について$p(\bm{z})$の値を求めることが正規化定数$Z$を除いて容易だとする。
すなわち
\begin{eqnarray}
p(z) = \frac{1}{Z_{p}}\tilde{p}(z)
\end{eqnarray}
において$\tilde{p}(z)$はすぐに求められるが$Z_{p}$がわからないとする。
棄却サンプリング法とは以下の手順でサンプリングを行うことをいう。
\begin{enumerate}
\item 簡単なサンプリング分布$q(z)$を準備する
\item 定数$k$を決めて全ての$z$に対して$kq(z) \geq \tilde{p}(z)$が成り立つようにその値を定める。
\item 乱数$z_{0}$を分布$q(z)$から生成する。
\item 乱数$u_{0}$を区間$[0,kq(z_{0})]$上の一様分布から生成する。
\item $u_{0}\leq\tilde{p}(z_{0})$ならばサンプルを保持し、そうでなければ棄却する。
\end{enumerate}
これによって、$p(z)$に従うサンプルが生成される。

\subsection{適応的棄却サンプリング}
省略

\subsection{重点サンプリング}
ここでは与えられた任意の$\bm{z}$の値について$p(\bm{z})$は計算できると仮定する。
極めて単純な戦略の１つは$\bm{z}$空間を均一なグリッドで離散化し、積分を
\begin{eqnarray}
\mathrm{E}[f] \approx \frac{1}{L}\sum_{l=1}^{L}p(\bm{z}^{(l)})f(\bm{z}^{(l)})
\end{eqnarray}
で計算することであるが、これは和を取る項の数が$\bm{z}$の次元に対して指数的に大きくなる。
重点サンプリングではサンプリングが容易な提案分布$q(\bm{z})$を用いて
\begin{eqnarray}
\mathbb{E}[f] &=& \int f(\bm{z})p(\bm{z})d\bm{z} \notag \\
&=& \int f(\bm{z})\frac{p(\bm{z})}{q(\bm{z})}q(\bm{z})d\bm{z} \notag \\
&\approx& \frac{1}{L} \sum_{l=1}^{L}\frac{p(\bm{z}^{(l)})}{q(\bm{z}^{(l)})}f(\bm{z}^{(l)})
\end{eqnarray}
と計算する。

また、$p(\bm{z})=\tilde{p(\bm{z})}Z_{p}$のように分布が正規化定数を除いてしか評価できない場合を考える。
さらに、用いたい重点サンプリングの分布も$q(\bm{z})=\tilde{q(\bm{z})}Z_{q}$のように正規化定数を除いてしか評価できないものを用いたいとする。
この場合は
\begin{eqnarray}
\mathbb{E}[f] &=& \int f(\bm{z})p(\bm{z})d\bm{z} \notag \\
&=& \frac{Z_{q}}{Z_{p}} \int f(\bm{z})\frac{\tilde{p}(\bm{z})}{q(\bm{z})}\tilde{q}(\bm{z})d\bm{z} \notag \\
&\approx& \frac{Z_{q}}{Z_{p}} \frac{1}{L} \sum_{l=1}^{L} \tilde{r}_{l}f(\bm{z}^{(l)})
\end{eqnarray}
となる。
\textcolor{blue}{
この議論はどの分布に基づいてサンプルを生成することを想定しているのかよくわからない。
}

\subsection{SIR}
省略

\subsection{サンプリングとEMアルゴリズム}
EMアルゴリズムにおいて、Mステップで$\bm{\theta}$に関して最適化される関数は完全データの対数尤度の期待値
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}^{\mathrm{old}} = \int p(\bm{Z}|\bm{X},\bm{\theta}^{\mathrm{old}} \ln p(\bm{Z},\bm{X}|\bm{\theta}) d\bm{Z}
\end{eqnarray}
で与えられるが、これを事後分布の現在の推定$p(\bm{Z}|\bm{X},\bm{\theta}^{old}$から抽出されたサンプル$\{\bm{Z}^{(l)}\}$を用いて
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}^{\mathrm{old}}) \approx \frac{1}{L}\sum_{l=1}^{L}\ln p(\bm{Z}^{(l)},\bm{X}|\bm{\theta})
\end{eqnarray}
のように近似するアルゴリズムをモンテカルロEMアルゴリズムという。

また、$p(\bm{Z}|\bm{X})$からサンプリングを行いたいが直接はできず、$p(\bm{\theta}|\bm{Z},\bm{X})$からサンプリングすることが比較的容易な場合に以下のアルゴリズムをIPアルゴリズムという。
\begin{itemize}
\item まず$p(\bm{\theta})$に対する現在の推定に基づいてサンプル$\bm{\theta}^{(l)}$を抽出し、それを用いて$p(\bm{Z}|\bm{\theta}^{(l)},\bm{X})$からサンプル$\bm{Z}^{(l)}$を抽出する。(Iステップ)
\item $\bm{\theta}$の事後分布に対する推定を
\begin{eqnarray}
p(\bm{\theta}|\bm{X}) \approx \frac{1}{L}\sum_{l=1}^{L}p(\bm{\theta}|\bm{Z}^{(l)},\bm{X})
\end{eqnarray}
によって計算し更新する。(Pステップ)
\end{itemize}

\section{マルコフ連鎖モンテカルロ}
確率$p(\bm{z}) = \tilde{p}(\bm{z})/Z_{p}$に従うサンプリングを行う以下のアルゴリズムをMetropolis法とよぶ
\begin{enumerate}
\item $\bm{z}^{(0)}$を適当に選ぶ
\item $\bm{z}^{(\tau)}$が与えられた場合、提案分布$q(\bm{z}^{*}|\bm{z}^{(\tau)})$に基づいて$\bm{z}^{(\tau+1)}$の候補を提案。
このとき$q(\bm{z}_{A}|\bm{z}_{B}) = q(\bm{z}_{B}|\bm{z}_{A})$であるとする。
\item 確率
\begin{eqnarray}
A(\bm{z}^{*},\bm{z}^{(\tau)}) = \min \left(1, \frac{\tilde{p}(\bm{z}^{*})}{\tilde{p}(\bm{z}^{(\tau)})} \right)
\end{eqnarray}
で$\bm{z}^{*}$を$\bm{z}^{(\tau+1)}$として受理し、棄却された場合は$\bm{z}^{(\tau+1)} = \bm{z}^{(\tau)}$とする。
\end{enumerate}

\subsection{マルコフ連鎖}
ここではマルコフ連鎖の一般的性質を考える。
１次マルコフ連鎖とは確率変数の系列であって
\begin{eqnarray}
p(\bm{z}^{(m+1)}|\bm{z}^{(1)},\cdots,\bm{z}^{(m)}) = p(\bm{z}^{(m+1)}|\bm{z}^{(m)})
\end{eqnarray}
が成立するものとして定義される。
よって、初期変数の確率分布$p(\bm{z}^{(0)})$と、遷移確率$T_{m}(\bm{z}^{(m)},\bm{z}^{(m+1)})\equiv p(\bm{z}^{(m+1)}|p(\bm{z}^{(m)})$を与えることで
マルコフ連鎖を指定することができる。全ての$m$について遷移確率が同じマルコフ連鎖は、均一マルコフ連鎖と呼ばれる。

また、特定の変数の周辺確率は、ひとつ前の周辺確率を用いて
\begin{eqnarray}
p(\bm{z}^{(m+1)}) = \sum_{\bm{z}^{(m)}} p(\bm{z}^{(m+1)}|\bm{z}^{(m)}) p(\bm{z}^{(m)})
\end{eqnarray}
と表すことができる。
よって均一なマルコフ連鎖において
\begin{eqnarray}
p^{*}(\bm{z}) = \sum_{\bm{z}'}T(\bm{z}',\bm{z})p^{*}(\bm{z}')
\end{eqnarray}
が成り立つ場合、分布$p^{*}(\bm{z})$は不変である。なお、与えられたマルコフ連鎖が２つ以上の不変分布を持ち得る可能性がある。

分布$p^{*}(\bm{z})$と遷移確率に関する
\begin{eqnarray}
p^{*}(\bm{z})T(\bm{z},\bm{z}') = p^{*}(\bm{z}')T(\bm{z}',\bm{z})
\end{eqnarray}
で表される等式を、詳細つり合い条件というが、これを満たす場合$p^{*}(\bm{z})$は不変分布になる。
なぜなら
\begin{eqnarray}
\sum_{\bm{z}'}p^{*}(\bm{z}')T(\bm{z}',\bm{z}) = \sum_{\bm{z}'}p^{*}(\bm{z})T(\bm{z},\bm{z}') = p^{*}(\bm{z})\sum_{\bm{z}'}p(\bm{z}'|\bm{z}) = p^{*}(\bm{z})
\end{eqnarray}
となるからである。

\subsection{Metropolis-Hastingsアルゴリズム}
確率$p(\bm{z}) = \tilde{p}(\bm{z})/Z_{p}$に従うサンプリングを行う以下のアルゴリズムを法とよぶ
\begin{enumerate}
\item $\bm{z}^{(0)}$を適当に選ぶ
\item $\bm{z}^{(\tau)}$が与えられた場合、提案分布$q(\bm{z}^{*}|\bm{z}^{(\tau)})$に基づいて$\bm{z}^{(\tau+1)}$の候補を提案。

\item 確率
\begin{eqnarray}
A(\bm{z}^{*},\bm{z}^{(\tau)}) = \min \left(1, \frac{\tilde{p}(\bm{z}^{*})q(\bm{z}^{(\tau)}|\bm{z}^{*})}{\tilde{p}(\bm{z}^{(\tau)})q(\bm{z}^{*}|\bm{z}^{\tau})} \right)
\end{eqnarray}
で$\bm{z}^{*}$を$\bm{z}^{(\tau+1)}$として受理し、棄却された場合は$\bm{z}^{(\tau+1)} = \bm{z}^{(\tau)}$とする。
\end{enumerate}
提案分布が対称な場合はMetropolis法と同一である。
遷移確率は
\begin{eqnarray}
T(\bm{z},\bm{z}') = q(\bm{z}'|\bm{z})A(\bm{z}',\bm{z})
\end{eqnarray}
であり、
\begin{eqnarray}
p(\bm{z})q(\bm{z}'|\bm{z})A(\bm{z}'|\bm{z}) &=& \min( p(\bm{z})q(\bm{z}'|\bm{z}) , p(\bm{z}')q(\bm{z}|\bm{z}') ) \notag \\
&=& \min( p(\bm{z}')q(\bm{z}|\bm{z}'), p(\bm{z})q(\bm{z}'|\bm{z}) ) \notag \\
&=& p(\bm{z}')q(\bm{z}|\bm{z}')A(\bm{z},\bm{z}')
\end{eqnarray}
となるため、詳細つり合いの条件を満たす。したがってMetropolis-Hastingsアルゴリズムで行うサンプリングは、分布$p(\bm{z})$に従う。

\section{ギブスサンプリング}
サンプリングしたい確率分布$p(\bm{z})=p(z_{1},\cdots,z_{M})$に対して、以下のアルゴリズムで行うサンプリングをギブスサンプリングという
\begin{enumerate}
\item $\bm{z}^{(0)}$を適当に選ぶ。
\item $\tau = 1,\cdots,T$に対して以下を行う。
\begin{itemize}
\item $z_{1}^{(\tau+1)}$ を$p(z_{1}|z_{2}^{(\tau)},\cdots,z_{M}^{(\tau)})$に従ってサンプルする。
\item[ ] $\vdots$
\item $z_{j}^{(\tau+1)}$ を$p(z_{j}|z_{1}^{(\tau+1)},\cdots,z_{j-1}^{(\tau+1)},\cdots,z_{j+1}^{(\tau)},\cdots,z_{M}^{(\tau)})$に従ってサンプルする。
\item[ ] $\vdots$
\item $z_{M}^{(\tau+1)}$ を$p(z_{M}|z_{1}^{(\tau+1)},\cdots,z_{M-1}^{(\tau+1)})$に従ってサンプルする。
\end{itemize}
\end{enumerate}
これは実は提案確率を$q(\bm{z}^{*}|\bm{z}) = p(z_{k}^{*}|\bm{z}_{\backslash k})$とするMetropolis-Hastingsサンプリングになっている。
ここで$\bm{z}_{\backslash k}$は$\bm{z}$から$z_{k}$を取り除いたものである。
確率による棄却判定が行われないのは
\begin{eqnarray}
\frac{p(\bm{z^{*}})q(\bm{z}|\bm{z}^{*})}{p(\bm{z})q(\bm{z}^{*}|\bm{z})} =
\frac{p(z_{k}^{*}|\bm{z}_{\backslash k}^{*}) p(\bm{z}_{\backslash k}^{*})p(z_{k}|\bm{z}_{\backslash k}^{*}) }{ p(z_{k}|\bm{z}_{\backslash k}) p(\bm{z}_{\backslash k}^{*})p(z_{k}^{*}|\bm{z}_{\backslash k}) } = 1
\end{eqnarray}
より$A(\bm{z}^{*},\bm{z})=1$となるためである。ここで$\bm{z}_{\backslash k}^{*} = \bm{z}_{\backslash k}$を用いた。

\section{スライスサンプリング}
正規化されていない分布$\tilde{p}(\bm{z})$に従ってサンプリングを行う方法として、付加的な変数$u$を用いた分布
\begin{eqnarray}
\hat{p}(z,u) =
\begin{cases}
1/Z_{p} & 0 \leq u \leq \tilde{p}(z)のとき \\
0 & それ以外
\end{cases}
\end{eqnarray}
からサンプリングを行う方法がある。これは$0 \leq u \leq \tilde{p}(z)$を満たす領域で定義された一様分布であり、このサンプリングを実現するのがスライスサンプリングである。
具体的には$z$の値が決められたときに範囲$0\leq u \leq \tilde{p}(z)$の領域から一様に$u$をサンプリングし、次に$u$を固定し$\{z:\tilde{p}(z)>u\}$で定義される領域から$z$を一様にサンプリングすることを繰り返す。
\textcolor{blue}{
本文に明示的に書かれていないが、これはギブスサンプリングの特別な場合に他ならない。
}

\section{ハイブリッドモンテカルロアルゴリズム}
省略
\subsection{力学系}
省略
\subsection{ハイブリッドモンテカルロアルゴリズム}
省略
\section{分配関数の推定}
省略

\chapter{連続潜在変数}
\section{主成分分析}

\subsection{分散最大化による定式化}
ここでは$D$次元の観測値のデータ集合$\{\bm{x}_{n}\}$を$M<D$次元の空間の上に射影することを考える。
まず１次元空間への射影を考える。
この空間の方向を$D$次元ベクトルを用いて、$\bm{u}_{1}$と表すことにすると、各データ点は$\bm{u}_{1}^{T}\bm{x}_{n}$に射影される。
その平均値は
\begin{eqnarray}
\bar{\bm{x}} = \frac{1}{N}\sum_{n=1}^{N}\bm{x}_{n}
\end{eqnarray}
を用いて$\bm{u}_{1}^{T}\bar{\bm{x}}_{n}$と表すことができる。また、その分散は
\begin{eqnarray}
\frac{1}{N} \sum_{n=1}^{N}\{ \bm{u}_{1}^{T}\bm{x}_{n}-\bm{u}_{1}^{T}\bar{\bm{x}} \}^{2} = \bm{u}_{1}^{T}\bm{S}\bm{u}_{1}
\end{eqnarray}
であらわされる。ここで
\begin{eqnarray}
\bm{S} = \frac{1}{N}\sum_{n=1}^{N}(\bm{x}_{n}-\bar{\bm{x}})(\bm{x}_{n}-\bar{\bm{x}})^{T}
\end{eqnarray}
である。これを最大化する$\bm{u}$を求めるためには、ラグランジュ未定乗数法を用いて
\begin{eqnarray}
\bm{u}_{1}^{T}\bm{S}\bm{u}_{1} + \lambda_{1}(1-\bm{u}_{1}^{T}\bm{u}_{1})
\end{eqnarray}
を微分し
\begin{eqnarray}
\bm{S}\bm{u}_{1} = \lambda_{1}\bm{u}_{1}
\end{eqnarray}
を得る。これに左から$\bm{u}_{1}^{T}$をかけると
\begin{eqnarray}
\bm{u}_{1}^{T}\bm{S}\bm{u}_{1} = \lambda_{1}
\end{eqnarray}
を得る。これらのことから、分散を最大にするには$\bm{u}_{1}$を$\bm{S}$の最大固有値に対応する固有ベクトルにすればよいことがわかる。

\subsection{誤差最小化による定式化}
今度は射影誤差の最小化に基づいた主成分分析の定式化を考える。
まず$D$次元の完全正規直交系
\begin{eqnarray}
\bm{u}_{i}^{T}\bm{u}_{j} = \delta_{ij}
\end{eqnarray}
を導入する。すると各データ点は一意的に
\begin{eqnarray}
\bm{x}_{n} = \sum_{i=1}^{D}\alpha_{ni}\bm{u}_{i}
\end{eqnarray}
と表すことができ、正規直交性より
\begin{eqnarray}
\bm{x}_{n} = \sum_{i=1}^{D}(\bm{x}_{n}^{T}\bm{u}_{i})\bm{u}_{i}
\end{eqnarray}
と書くことができる。
しかしここでは$M$個の変数で各データ点を近似することにあるのであって、各データ点を
\begin{eqnarray}
\tilde{\bm{x}}_{n} = \sum_{i=1}^{M}z_{ni}\bm{u}_{i} + \sum_{i=M+1}^{D}b_{i}\bm{u}_{i}
\end{eqnarray}
と近似する。ここで$\{z_{ni}\}$はデータ点に依存しているが、$\{b_{i}\}$はすべてのデータ点に共通である。
近似は、誤差関数
\begin{eqnarray}
J = \frac{1}{N}\sum_{n=1}^{N} || \bm{x}_{n}-\tilde{\bm{x}}_{n} ||^{2}
\end{eqnarray}
を最小化するように$\{\bm{u}_{i}\},\{z_{ni}\},\{b_{i}\}$を決めることによって行う。
$\{z_{ni}\}$と$b_{i}$については、$J$に$\tilde{\bm{x}}_{n}$の表式を代入して微分をすることで
\begin{eqnarray}
z_{nj} &=& \bm{x}_{n}^{T}\bm{u}_{j} \notag \\
b_{j} &=& \bar{\bm{x}}^{T}\bm{u}_{j}
\end{eqnarray}
を得る。これより
\begin{eqnarray}
\bm{x}_{n}-\tilde{\bm{x}}_{n} = \sum_{i=M+1}^{D}\{(\bm{x}_{n}-\bar{\bm{x}})^{T}\bm{u}_{i} \} \bm{u}_{i}
\end{eqnarray}
が従い、
\begin{eqnarray}
J = \frac{1}{N}\sum_{n=1}^{N}\sum_{i=M+1}^{D}(\bm{x}_{n}^{T}\bm{u}_{i}-\bar{\bm{x}}^{T}\bm{u}_{i})^{2} = \sum_{i=M+1}^{D}\bm{u}_{i}^{T}\bm{S}\bm{u}_{i}
\end{eqnarray}
となる。これを最小化するには$\bm{u}_{i}(i>M)$を$\bm{S}$の小さい固有値に対応する固有ベクトルに選べばよく、
\begin{eqnarray}
J = \sum_{i=M+1}^{D}\lambda_{i}
\end{eqnarray}
となる。

\subsection{主成分分析の応用}
主成分分析はデータベクトル$\bm{x}_{n}$に対する圧縮方法として利用することができる。これは
\begin{eqnarray}
\bar{\bm{x}} = \sum_{i=1}^{D}(\bar{\bm{x}}^{T}\bm{u}_{i})\bm{u}_{i}
\end{eqnarray}
より
\begin{eqnarray}
\tilde{\bm{x}}_{n} &=& \sum_{i=1}^{M} (\bm{x}_{n}^{T}\bm{u}_{i})\bm{u}_{i} + \sum_{i=M+1}^{D} (\bar{\bm{x}}_{n}^{T}\bm{u}_{i})\bm{u}_{i} \notag \\
&=& \bar{\bm{x}} + \sum_{i=1}^{M} (\bm{x}_{n}^{T}\bm{u}_{i} - \bar{\bm{x}}_{n}^{T}\bm{u}_{i})\bm{u}_{i}
\end{eqnarray}
となるため、$D$次元ベクトルを$M$次元ベクトルで置き換えたことになるためである。

また主成分分析はデータの前処理にも応用できる。例えば、データ集合に対して標準化された共分散行列は
\begin{eqnarray}
\rho_{ij} = \frac{1}{N}\sum_{n=1}^{N}\frac{(x_{ni}-\bar{x}_{i})}{\rho_{i}} \frac{(x_{nj}-\bar{x}_{j})}{\rho_{j}}
\end{eqnarray}
であらわすことができるが、$\bm{S}\bm{U} = \bm{U}\bm{L}$を満たす固有ベクトルの行列$\bm{U}$および、対角成分が固有値の行列$\bm{L}$を用いて、
\begin{eqnarray}
\bm{y}_{n} = \bm{L}^{-1/2}\bm{U}^{T}(\bm{x}_{n}-\bar{\bm{x}})
\end{eqnarray}
を定義すると、
\begin{eqnarray}
\frac{1}{N}\sum_{n=1}^{N}\bm{y}_{n} \bm{y}_{n}^{T} = \bm{I}
\end{eqnarray}
となるため、平均が$0$で標準化された共分散行列が単位行列となる。

\subsection{高次元データに対する主成分分析}
ここでは、データ点の数がベクトル空間の次元$D$よりも小さい場合を考える。
まず、$\bm{X}$を$n$番目の行が$(\bm{x}_{n}-\bar{\bm{x}})^{T}$である$N\times D$次元の行列とする。
すると$\bm{S} = N^{-1}\bm{X}^{T}\bm{X}$と書くことができて、固有ベクトルの方程式は
\begin{eqnarray}
\frac{1}{N}\bm{X}^{T}\bm{X}\bm{u}_{i} = \lambda_{i}\bm{u}_{i}
\end{eqnarray}
となる。これは、$D$次元の固有値方程式であって、$D$次元空間の$N$点の集合は高々$N-1$次元の部分空間を定義するので、$D-N+1$個の固有値は$0$になる。
$0$でない固有値は上の指揮に$\bm{X}$をかけて、
\begin{eqnarray}
\frac{1}{N}\bm{X}\bm{X}^{T}(\bm{X}\bm{u}_{i}) = \lambda (\bm{X}\bm{u}_{i})
\end{eqnarray}
とすることで$N$次元の固有値方程式の解として得ることができる。

\section{確率的主成分分析}
ここでは主成分分析が、ある確率的潜在変数モデルの最尤解として表現されることを示す。
まず、主部分空間に対応する$M$次元の潜在変数$\bm{z}$を明示的に導入し、ガウス事前分布$p(\bm{z})$に従うものとする。
また、観測変数$\bm{x}$については、潜在変数を条件とするガウス分布$p(\bm{x}|\bm{z})$に従うものとする。
具体的には
\begin{eqnarray}
p(\bm{z}) &=& \mathcal{N}(\bm{z}|\bm{0},\bm{I}) \notag \\
p(\bm{x}|\bm{z}) &=& \mathcal{N}(\bm{x}|\bm{W}\bm{z} + \bm{\mu}, \sigma^{2}\bm{I})
\end{eqnarray}
とする。
これより、周辺分布は
\begin{eqnarray}
p(\bm{x}) &=& \int p(\bm{x}|\bm{z})p(\bm{z})d\bm{z} = \mathcal{N}(\bm{x}|\bm{\mu},\bm{C}) \notag \\
\bm{C} &=& \bm{W}\bm{W}^{T} + \sigma^{2}\bm{I}
\end{eqnarray}
となる。予測分布$p(\bm{x})$はパラメータ$\bm{\mu},\bm{W},\sigma^{2}$によって規定されるが、$\bm{W}$については、$\bm{W}\bm{W}^{T}$だけが周辺分布に影響を与えるため、
任意の直交行列$\bm{R}$に対して、$\tilde{\bm{W}} = \bm{W}\bm{R}$は同じ周辺分布を与える。
事後分布については
\begin{eqnarray}
p(\bm{z}|\bm{x}) &=& \mathcal{N}(\bm{z}|\bm{M}^{-1}\bm{W}^{T}(\bm{x}-\bm{\mu}),\sigma^{2}\bm{M}^{-1}) \notag \\
\bm{M} &=& \bm{W}^{T}\bm{W} + \sigma^{2}\bm{I}
\end{eqnarray}
を得る。

\subsection{最尤法による主成分分析}
データ点の集合$\bm{X}=\{\bm{x}_{n}\}$が与えられた下で、対数尤度関数は
\begin{eqnarray}
\ln p(\bm{X}|\bm{\mu},\bm{W},\sigma^{2}) &=& \sum_{n=1}^{N}\ln p(\bm{x}_{n}|\bm{W},\bm{\mu},\sigma^{2}) \notag \\
&=& -\frac{ND}{2}\ln(2\pi) - \frac{N}{2}\ln|\bm{C}| - \frac{1}{2}\sum_{n=1}^{N}(\bm{x}_{n}-\bm{\mu})^{T}\bm{C}^{-1}(\bm{x}_{n}-\bm{\mu}) \notag \\
\end{eqnarray}
で与えられ、$\bm{\mu}$については
\begin{eqnarray}
\bm{\mu}_{\mathrm{ML}} = \frac{1}{N}\sum_{n=1}^{N}\bm{x}_{n}
\end{eqnarray}
で与えられる。また$\bm{W}$については、対数尤度関数の全ての停留点は
\begin{eqnarray}
\bm{W}_{\mathrm{ML}} = \bm{U}_{M}(\bm{L}_{M}-\sigma^{2}\bm{I})^{1/2}\bm{R}
\end{eqnarray}
で与えられる。ここで$\bm{U}_{M}$は$D\times M$行列で、その列ベクトルはデータ共分散行列$\bm{S}$の固有ベクトルのサイズ$M$の任意の部分集合であり、
$M$次元行列$\bm{L}_{M}$は対角成分に固有ベクトルに対応する固有値を持つ対角行列である。
また、尤度関数の最大値は、$M$個の固有ベクトルを固有値の大きい$M$個に選んだ場合に得られることが知られている。
この場合、$\bm{W}$の列ベクトルは、主成分分析の主部分空間を成す。
そして$\sigma^{2}$の最尤解は
\begin{eqnarray}
\sigma_{\mathrm{ML}}^{2} = \frac{1}{D-M}\sum_{i=M+1}^{D}\lambda_{i}
\end{eqnarray}
で与えられる。すなわち、切り捨てられた次元に関連した分散の平均である。

また、事後分布についての平均は
\begin{eqnarray}
\mathbb{E}[\bm{z}|\bm{x}] = \bm{M}^{-1}\bm{W}_{\mathrm{ML}}^{T}(\bm{x}-\bar{\bm{x}})
\end{eqnarray}
で与えられる。$\sigma^{2}\rightarrow 0$の極限では
\begin{eqnarray}
\mathbb{E}[\bm{z}|\bm{x}] &\rightarrow& (\bm{W}_{\mathrm{ML}}^{T}\bm{W}_{\mathrm{ML}})^{-1}\bm{W}_{\mathrm{ML}}^{T}(\bm{x}-\bar{\bm{x}}) \notag \\
&=& \bm{L}_{M}\bm{U}_{M}^{T}(\bm{x}-\bar{\bm{x}})
\end{eqnarray}
となって、本文(12.24)の結果と一致する。

\textcolor{blue}{ここに直交射影とあるが、いまひとつよくわからない。通常、直交射影は射影の値域と核が直交することをいう。
ここで核とは、写像$\bm{A}$に対して、$\{\bm{x}|\bm{A}\bm{x}=\bm{0}\}$のことである。
この場合、地域は$\mathbb{R}^{M}$内の部分空間で、核は$\mathbb{R}^{D}$内の部分空間のため、どのように直交を定義するのだろうか。
}

\subsection{EMアルゴリズムによる主成分分析}
ここではEMアルゴリズムによってパラメータを求める方法を考える。
すでに厳密な最尤パラメータは得られているが、高次元空間においては、サンプル共分散行列を扱うよりも計算量的に有利という場合もある。
EMアルゴリズムに必要な完全データの対数尤度関数は
\begin{eqnarray}
\ln p(\bm{X},\bm{Z}|\bm{\mu},\bm{W},\sigma^{2}) = \sum_{n=1}^{N}\{ \ln p(\bm{x}_{n}|\bm{z}_{n}) + \ln p(\bm{z}_{n}) \}
\end{eqnarray}
で与えられる。そして、その潜在変数の事後分布の期待値は
\begin{eqnarray}
\mathcal{Q}(\bm{\theta},\bm{\theta}^{\mathrm{old}}) &=& \mathbb{E}[\ln p(\bm{X},\bm{Z}|\bm{\mu},\bm{W},\sigma^{2}) ] \notag \\
&=& -\sum_{n=1}^{N}\left\{ \frac{D}{2}\ln(2\pi \sigma^{2}) + \frac{1}{2}\mathrm{Tr}(\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}]) \right. \notag \\
&+& \frac{1}{2\sigma^{2}}||\bm{x}_{n}-\bm{\mu}||^{2} - \frac{1}{\sigma^{2}}\mathbb{E}[\bm{z}]^{T}\bm{W}^{T}(\bm{x}_{n}-\bm{\mu}) \notag \\
&+& \left. \frac{1}{2\sigma^{2}}\mathrm{Tr}(\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}]\bm{W}^{T}\bm{W}) + \frac{M}{2}\ln(2\pi) \right\}
\end{eqnarray}
となる。$\bm{\mu}$に対しては最尤解がサンプル平均になることをすでに知っているので、以降では置き換えることにする。
古いパラメータに依存する潜在変数の期待値は
\begin{eqnarray}
\mathbb{E}[\bm{z}_{n}] &=& \bm{M}_{\mathrm{old}}^{-1}\bm{W}_{\mathrm{old}}^{T}(\bm{x}_{n}-\bar{\bm{x}}) \notag \\
\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}] &=& \sigma_{\mathrm{old}}^{2}\bm{M}_{\mathrm{old}}^{-1} + \mathbb{E}[\bm{z}_{n}]\mathbb{E}[\bm{z}_{n}]^{T}
\end{eqnarray}
である。これより$\mathcal{Q}$を最大化するように新しいパラメータを選ぶと
\begin{eqnarray}
\bm{W}_{\mathrm{new}} &=& \left[ \sum_{n=1}^{N}(\bm{x}_{n}-\bar{\bm{x}})\mathbb{E}[\bm{z}_{n}]^{T}\right] \left[ \sum_{n=1}^{N}\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}] \right]^{-1} \notag \\
\sigma^{\mathrm{new}} &=& \frac{1}{ND} \sum_{n=1}^{N} \left\{ ||\bm{x}_{n}-\bar{\bm{x}}||^{2} - 2\mathbb{E}[\bm{z}_{n}]^{T} \bm{W}_{\mathrm{new}}^{T}(\bm{x}_{n}-\bar{\bm{x}}) \right. \notag \\
&+& \left. \mathrm{Tr}(\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}]\bm{W}_{\mathrm{new}}^{T}\bm{W}_{\mathrm{new}}) \right\}
\end{eqnarray}
を得る。

\subsection{ベイズ的主成分分析}
省略

\subsection{因子分析}
因子分析のモデルは、潜在変数による線形ガウス分布で、分散が対角的ではあるが等方的でないものである。
すなわち$D$次元対角行列$\bm{\Psi}$を用いて
\begin{eqnarray}
p(\bm{x}|\bm{z}) = \mathcal{N}(\bm{x}|\bm{W}\bm{z}+\bm{\mu},\bm{\Psi})
\end{eqnarray}
と表されるものである。この場合周辺分布は
\begin{eqnarray}
p(\bm{x}) &=& \mathcal{N}(\bm{x}|\bm{\mu},\bm{C}) \notag \\
\bm{C} &=& \bm{W}\bm{W}^{T} + \bm{\Psi}
\end{eqnarray}
で与えられる。
このモデルの最尤解は$\bm{\mu}$についてはサンプル平均$\bar{\bm{x}}$であるが、$\bm{W},\bm{\Psi}$については厳密に求めることはできず、EMアルゴリズムで解くことになる。
Eステップについては、
\begin{eqnarray}
\mathbb{E}[\bm{z}_{n}] &=& \bm{G}_{\mathrm{old}}\bm{W}_{\mathrm{old}}^{T}\bm{\Psi}_{\mathrm{old}}^{-1}(\bm{x}_{n}-\bar{\bm{x}}) \notag \\
\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}] &=& \bm{G}_{\mathrm{old}} + \mathbb{E}[\bm{z}_{n}]\mathbb{E}[\bm{z}_{n}]^{T} \notag \\
\bm{G} &=& (\bm{I} + \bm{W}^{T}\bm{\Psi}^{-1}\bm{W})^{-1}
\end{eqnarray}
を得る。これにより、更新式として
\begin{eqnarray}
\bm{W}_{\mathrm{new}} &=& \left[ \sum_{n=1}^{N}(\bm{x}_{n}-\bar{\bm{x}})\mathbb{E}[\bm{z}_{n}]^{T}\right] \left[ \sum_{n=1}^{N}\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}] \right]^{-1} \notag \\
\bm{\Psi}_{\mathrm{new}} &=& \mathrm{diag}\left\{ \bm{S}-\bm{W}_{\mathrm{new}}\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}[\bm{z}_{n}](\bm{x}_{n}-\bar{\bm{x}})^{T} \right\}
\end{eqnarray}
を得る。

\section{カーネル主成分分析}
ここでは観測変数のデータ集合$\{\bm{x}_{n}\}$を考え、サンプル平均をあらかじめ引き去っているものとして$\sum_{n}\bm{x}_{n}=0$とする。
そして、$M$次元の特徴空間への非線形変換$\bm{\phi}(\bm{x})$を考える。
特徴空間での$M\times M$サンプル共分散行列は$\sum_{n}\bm{\phi}(\bm{x}_{n})=0$を仮定すると
\begin{eqnarray}
\bm{C} = \frac{1}{N}\sum_{n=1}^{N} \bm{\phi}(\bm{x}_{n}) \bm{\phi}(\bm{x}_{n})^{T}
\end{eqnarray}
であり、その固有ベクトル展開は$i=1,\cdots,M$に対して
\begin{eqnarray}
\bm{C}\bm{v}_{i} = \lambda_{i}\bm{v}_{i}
\end{eqnarray}
と定義される。これは$\bm{C}$の定義から
\begin{eqnarray}
\frac{1}{N}\sum_{n=1}^{N}\bm{\phi}(\bm{x}_{n})\{\bm{\phi}(\bm{x}_{n})^{T}\bm{v}_{i}\} = \lambda_{i}\bm{v}_{i}
\end{eqnarray}
で与えられる。上の式からベクトル$\bm{v}_{i}$は$\bm{\phi}(\bm{x}_{n})$の線形結合で与えられるので
\begin{eqnarray}
\bm{v}_{i} = \sum_{n=1}^{N}a_{in}\bm{\phi}(\bm{x}_{n})
\end{eqnarray}
と書くことができる。
\textcolor{blue}{($\{\bm{\phi}(\bm{x}_{n})\}$がはる部分空間の次元が$N$未満の場合、$a_{in}$の値は一意的には決まらない。)}
これを元の方程式に代入すると
\begin{eqnarray}
\frac{1}{N}\sum_{n=1}^{N}\bm{\phi}(\bm{x}_{n})\bm{\phi}(\bm{x}_{n})^{T}\sum_{m=1}^{N}a_{im}\bm{\phi}(\bm{x}_{m}) = \lambda_{i}\sum_{n=1}^{N}a_{in}\bm{\phi}(\bm{x}_{n})
\end{eqnarray}
となり、両辺に$\bm{\phi}(\bm{x}_{l})^{T}$をかけて$\bm{\phi}(\bm{x}_{n})^{T}\bm{\phi}(\bm{x}_{n})=k(\bm{x}_{n},\bm{x}_{m})$と表現することによって
\begin{eqnarray}
\frac{1}{N}\sum_{n=1}^{N}k(\bm{x}_{l},\bm{x}_{n})\sum_{m=1}^{N}a_{im}k(\bm{x}_{n},\bm{x}_{m}) = \lambda_{i}\sum_{n=1}^{N}a_{in}k(\bm{x}_{l},\bm{x}_{n})
\end{eqnarray}
を得る。これは行列記法では
\begin{eqnarray}
\bm{K}^{2}\bm{a}_{i} = \lambda_{i}N\bm{K}\bm{a}_{i}
\end{eqnarray}
と書くことができて、$0$でない固有値および固有ベクトルは
\begin{eqnarray}
\bm{K}\bm{a}_{i} = \lambda_{i}N\bm{a}_{i}
\end{eqnarray}
を解くことによって得られる。また規格化条件は
\begin{eqnarray}
1 = \bm{v}_{i}^{T}\bm{v}_{i} = \sum_{n=1}^{N}\sum_{m=1}^{N}a_{in}a_{im}\bm{\phi}(\bm{x}_{n})^{T}\bm{\phi}(\bm{x}_{m}) = \bm{a}_{i}^{T}\bm{K}\bm{a}_{i} = \lambda_{i}N\bm{a}_{i}^{T}\bm{a}_{i}
\end{eqnarray}
で与えられる。

固有値問題を解いた場合、点$\bm{x}$の固有ベクトル$i$の上への射影は
\begin{eqnarray}
y_{i}(\bm{x}) = \bm{\phi}(\bm{x})^{T}\bm{v}_{i} = \sum_{n=1}^{N}a_{in}\bm{\phi}(\bm{x})^{T}\bm{\phi}(\bm{x}) = \sum_{n=1}^{N}a_{in}k(\bm{x},\bm{x}_{n})
\end{eqnarray}
で与えられ、カーネル関数だけを通して表される。

ところで、一般には$\bm{\phi}(\bm{x}_{n})$の平均は$0$にならない。
その場合、共分散行列は
\begin{eqnarray}
\tilde{\bm{\phi}}(\bm{x}_{n}) = \bm{\phi}(\bm{x}_{n}) - \frac{1}{N}\sum_{l=1}^{N}\bm{\phi}(\bm{x}_{l})
\end{eqnarray}
を用いて
\begin{eqnarray}
\tilde{K}_{nm} &=& \tilde{\bm{\phi}}(\bm{x}_{n})^{T}\tilde{\bm{\phi}}(\bm{x}_{m}) \notag \\
&=& \bm(\bm{x}_{n})^{T}\bm{\phi}(\bm{x}_{m}) - \frac{1}{N}\sum_{l=1}^{N} \bm(\bm{x}_{n})^{T}\bm{\phi}(\bm{x}_{l}) \notag \\
&-& \frac{1}{N}\sum_{l=1}^{N} \bm(\bm{x}_{l})^{T}\bm{\phi}(\bm{x}_{m}) + \frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{l=1}^{N} \bm(\bm{x}_{j})^{T}\bm{\phi}(\bm{x}_{l}) \notag \\
&=& k(\bm{x}_{n},\bm{x}_{m}) - \frac{1}{N}\sum_{l=1}^{N} k(\bm{x}_{l},\bm{x}_{m}) \notag \\
&-& \frac{1}{N} \sum_{l=1}^{N} k(\bm{x}_{n},\bm{x}_{l}) + \frac{1}{N^2}\sum_{j=1}^{N}\sum_{l=1}^{N}k(\bm{x}_{j},\bm{x}_{l})
\end{eqnarray}
と表すことができる。これは行列表記だと
\begin{eqnarray}
\tilde{\bm{K}} = \bm{K} - \bm{I}_{N}\bm{K} - \bm{K}\bm{1}_{N} + \bm{1}_{N}\bm{K}\bm{1}_{N}
\end{eqnarray}
と書くことができる。ここで$\bm{1}_{N}$は全ての要素が$1/N$の行列である。

\section{非線形潜在変数モデル}

\subsection{独立成分分析}
省略

\subsection{自己連想ニューラルネットワーク}
省略

\subsection{非線形多様体のモデル化}
省略


\chapter{系列データ}

\section{マルコフモデル}
一般に観測系列の同時分布は、確率の積の規則を用いることで
\begin{eqnarray}
p(\bm{x}_{1},\cdots,\bm{x}_{N}) = p(\bm{x}_{1})\prod_{n=2}^{N}p(\bm{x}_{N}|\bm{x}_{1},\cdots,\bm{x}_{n-1})
\end{eqnarray}
と書くことができるが、
\begin{eqnarray}
p(\bm{x}_{1},\cdots,\bm{x}_{N}) = p(\bm{x}_{1})\prod_{n=2}^{N}p(\bm{x}_{n}|\bm{x}_{n-1})
\end{eqnarray}
のように、$n$番目のデータが$n-1$番目のデータのみに依存するものを一次マルコフ連鎖という。この場合
\begin{eqnarray}
p(\bm{x}_{n}|\bm{x}_{1},\cdots,\bm{x}_{n-1}) = p(\bm{x}_{n}|\bm{x}_{n-1})
\end{eqnarray}
が成り立つ。同様にして二次マルコフ連鎖は
\begin{eqnarray}
p(\bm{x}_{1},\cdots,\bm{x}_{N}) = p(\bm{x}_{1})p(\bm{x}_{2}|\bm{x}_{1})\prod_{n=3}^{N}p(\bm{x}_{n}|\bm{x}_{n-1},\bm{x}_{n-2})
\end{eqnarray}
で与えられるモデルである。

また、各観測値$\bm{x}_{n}$に対し、対応する潜在変数$\bm{z}_{n}$を導入し、これがマルコフ連鎖を構成するとするモデルを状態空間モデルといい、その同時分布は
\begin{eqnarray}
p(\bm{x}_{1},\cdots,\bm{x}_{N},\bm{z}_{1},\cdots,\bm{z}_{N}) = p(\bm{z}_{1})\left[ \prod_{n=2}^{N}p(\bm{z}_{n}|\bm{z}_{n-1}) \right] \prod_{n=1}^{N}p(\bm{x}_{n}|\bm{z}_{n})
\end{eqnarray}
で与えられる。これは$\bm{z}_{n}$が与えられたときに$\bm{z}_{n-1}$と$\bm{z}_{n+1}$が独立、
\begin{eqnarray}
\bm{z}_{n+1} \Perp \bm{z}_{n-1} | \bm{z}_{n}
\end{eqnarray}
という性質を満たす。
潜在変数が離散変数の場合隠れマルコフモデル、潜在変数と観測変数の両方がガウス分布に従うとき、線形動的システムという。

\section{隠れマルコフモデル}
隠れマルコフモデルでは$\bm{z}_{n}$は一対K符号化を用いると便利であり、その遷移確率は$A_{jk}\equiv p(z_{nk}=1|z_{n-1,j}=1)$によってあらわされる。すなわち
\begin{eqnarray}
p(\bm{z}_{n}|\bm{z}_{n-1},\bm{A}) = \prod_{k=1}^{K}\prod_{j=1}^{K}A_{jk}^{z_{n-1,j}z_{nk}}
\end{eqnarray}
が成り立つ。最初のノードは確率のベクトル$\bm{\pi}$で表される周辺分布
\begin{eqnarray}
p(\bm{z}_{1}|\bm{\pi}) = \prod_{k=1}^{K}\pi_{k}^{z_{1k}}
\end{eqnarray}
を持つ。また、$\bm{z}_{n}$が与えられた場合の$\bm{x}_{n}$の分布を支配するパラメータを$\bm{\phi}$とする。具体的には
\begin{eqnarray}
p(\bm{x}_{n}|\bm{z}_{n},\bm{\phi}) = \prod_{k=1}^{K}p(\bm{x}_{n}|\bm{\phi}_{k})^{z_{nk}}
\end{eqnarray}
とあらわされる。

均一なモデルとは$\bm{A}$および$\bm{\phi}$が$n$に依存しないモデルのことであり、
\begin{eqnarray}
p(\bm{X},\bm{Z}|\bm{\theta}) = p(\bm{z}_{1}|\bm{\pi})\left[ \prod_{n=2}^{N}p(\bm{z}_{n}|\bm{z}_{n-1},\bm{A}) \right] \prod_{m=1}^{N}p(\bm{x}_{m}|\bm{z}_{m},\bm{\phi})
\end{eqnarray}
と書くことができる。ここで$\bm{\theta}=\{ \bm{\pi},\bm{A},\bm{\phi} \}$はモデルを支配するパラメータである。

また、$\bm{A}$の$k\leq j$となる$A_{jk}$の成分をゼロとして得られるモデルをleft-to-rightHMMという。

\subsection{HMMの最尤推定}
データ集合$\bm{X}=\{\bm{x}_{1},\cdots,\bm{x}_{N}\}$が観測された場合のHMMのパラメータを最尤推定で決定することを考える。
そのためにEM法を用いることにする。この場合
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}^{old}) &=& \sum_{\bm{Z}}p(\bm{Z}|\bm{X},\bm{\theta}^{old})\ln p(\bm{X},\bm{Z}|\bm{\theta}) \notag \\
&=& \sum_{\bm{Z}}p(\bm{Z}|\bm{X},\bm{\theta}^{old})\left[ \ln p(\bm{z}_{1}|\bm{\pi}) + \sum_{n=2}^{N}\ln p(\bm{z}_{n}|\bm{z}_{n-1},\bm{A})
+ \sum_{n=1}^{N}\ln p(\bm{x}_{n}|\bm{z}_{n},\bm{\phi}) \right] \notag \\
\end{eqnarray}
となるが、
\begin{eqnarray}
\gamma(\bm{z}_{n}) &=& p(\bm{z}_{n}|\bm{X},\bm{\theta}^{old}) \notag \\
\xi(\bm{z}_{n-1},\bm{z}_{n}) &=& p(\bm{z}_{n-1},\bm{z}_{n}|\bm{X},\bm{\theta}^{old})
\end{eqnarray}
と表記し、さらに
\begin{eqnarray}
\gamma(z_{nk}) &=& p(z_{nk}=1|\bm{X},\bm{\theta}^{old}) \notag \\
\xi(z_{n-1,j},z_{nk}) &=& p(z_{n-1,j}=z_{nk}=1|\bm{X},\bm{\theta}^{old})
\end{eqnarray}
と書くことにすると、
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}^{old}) &=& \sum_{k=1}^{K}\gamma(z_{1k})\ln \pi_{k} + \sum_{n=2}^{N}\sum_{j=1}^{K}\sum_{k=1}^{K}\xi(z_{n-1,j},z_{nk})\ln A_{jk} \notag \\
&+& \sum_{n=1}^{N}\sum_{k=1}^{K}\gamma(z_{nk}) \ln p(\bm{x}_{n}|\bm{\phi}_{k})
\end{eqnarray}
を得る。
Mステップでは$\gamma(\bm{z}_{n})$と$\xi(\bm{z}_{n-1},\bm{z}_{n})$を定数とみなし、パラメータ$\bm{\theta}=\{ \bm{\pi},\bm{A},\bm{\phi} \}$に関して$Q(\bm{\theta},\bm{\theta}^{old})$を最大化するが、これはラグランジュ未定乗数法を使って
\begin{eqnarray}
\pi_{k} &=& \frac{\gamma(z_{1k})}{\sum_{j=1}^{K}\gamma(z_{1j})} \notag \\
A_{jk} &=& \frac{\sum_{n=2}^{N}\xi(z_{n-1,j}z_{nk})}{\sum_{l=1}^{K}\sum_{n=2}^{N}\xi(z_{n-1,j}z_{nl})}
\end{eqnarray}
を得る。

\subsection{フォワードバックワードアルゴリズム}
次にEMアルゴリズムのＥステップに対応する$\gamma(z_{nk})$と$\xi(z_{n-1,j},z_{nk})$を求める方法について検討する。
そのために、条件付き独立性を以下に書き下すと
\begin{eqnarray}
p(\bm{X}|\bm{z}_{n}) &=& p(\bm{x}_{1},\cdots,\bm{x}_{n}|\bm{z}_{n})p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n}) \notag \\
p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{x}_{n},\bm{z}_{n}) &=& p(\bm{x}_{1},\cdots, \bm{x}_{n-1}|\bm{z}_{n}) \notag \\
p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{z}_{n-1},\bm{z}_{n}) &=& p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{z}_{n-1}) \notag \\
p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n}\bm{z}_{n+1}) &=& p(\bm{z}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n+1}) \notag \\
p(\bm{x}_{n+2},\cdots,\bm{x}_{N}|\bm{z}_{n+1}\bm{x}_{n+1}) &=& p(\bm{z}_{n+2},\cdots,\bm{x}_{N}|\bm{z}_{n+1}) \notag \\
p(\bm{X}|\bm{z}_{n-1}\bm{x}_{n}) &=& p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{z}_{n-1})p(\bm{x}_{n}|\bm{z}_{n})p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n}) \notag \\
p(\bm{x}_{N+1}|\bm{X},\bm{z}_{N+1} &=& p(\bm{x}_{N+1}|\bm{z}_{N+1}) \notag \\
p(\bm{z}_{N+1}|\bm{z}_{N},\bm{X}) &=& p(\bm{z}_{N+1}|\bm{z}_{N})
\end{eqnarray}
となる。
そして$\gamma(\bm{z}_{n})$については、ベイズの定理と条件付き独立性より
\begin{eqnarray}
\gamma(\bm{z}_{n}) &=& p(\bm{z}_{n}|\bm{X}) = \frac{p(\bm{X}|\bm{z}_{n})p(\bm{z}_{n})}{p(\bm{X})} \notag \\
&=& \frac{p(\bm{x}_{1},\cdots,\bm{x}_{n},\bm{z}_{n})p(\bm{x}_{n+1},\cdots,\bm{x}_{N})}{p(\bm{X})}
= \frac{\alpha(\bm{z}_{n})\beta(\bm{z}_{n}) }{ p(\bm{X}) }
\end{eqnarray}
を得る。ただし
\begin{eqnarray}
\alpha(\bm{z}_{n}) \equiv p(\bm{x}_{1},\cdots,\bm{x}_{n},\bm{z}_{n}) \notag \\
\beta(\bm{z}_{n}) \equiv p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n})
\end{eqnarray}
である。そして、$\alpha,\beta$は再帰的に求めることができて、
\begin{eqnarray}
\alpha(\bm{z}_{n}) &=& p(\bm{x}_{1},\cdots,\bm{x}_{n},\bm{z}_{n}) \notag \\
&=& p(\bm{x}_{1},\cdots,\bm{x}_{n}|\bm{z}_{n})p(\bm{z}_{n}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n})  p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{z}_{n})p(\bm{z}_{n}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n}) \sum_{\bm{z}_{n-1}} p(\bm{x}_{1},\cdots,\bm{x}_{n-1},\bm{z}_{n-1},\bm{z}_{n}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n}) \sum_{\bm{z}_{n-1}} p(\bm{x}_{1},\cdots,\bm{x}_{n-1},\bm{z}_{n}|\bm{z}_{n-1}) p(\bm{z}_{n-1}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n}) \sum_{\bm{z}_{n-1}} p(\bm{x}_{1},\cdots,\bm{x}_{n-1}|\bm{z}_{n-1})p(\bm{z}_{n}|\bm{z}_{n-1})p(\bm{z}_{n-1}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n}) \sum_{\bm{z}_{n-1}} p(\bm{x}_{1},\cdots,\bm{x}_{n-1},\bm{z}_{n-1})p(\bm{z}_{n}|\bm{z}_{n-1}) \notag \\
&=& p(\bm{x}_{n}|\bm{z}_{n}) \sum_{\bm{z}_{n-1}} \alpha(\bm{z}_{n-1})p(\bm{z}_{n}|\bm{z}_{n-1})
\end{eqnarray}
この初期条件は
\begin{eqnarray}
\alpha(\bm{z}_{1}) = p(\bm{x}_{1},\bm{z}_{1}) = p(\bm{z}_{1})p(\bm{x}_{1}|\bm{z}_{1}) = \prod_{k=1}^{K}\{ \pi_{k}p(\bm{x}_{1}|\bm{\phi}_{k})\}^{z_{1k}}
\end{eqnarray}
で与えられる。

同様に$\beta(\bm{z}_{n})$についても
\begin{eqnarray}
\beta(\bm{z}_{n}) &=& p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{n+1}}p(\bm{x}_{n+1},\cdots,\bm{x}_{N},\bm{z}_{n+1}|\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{n+1}}p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n},\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{n+1}}p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{n+1}}p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n+1})p(\bm{x}_{n+1}|\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{n+1}}\beta(\bm{z}_{n+1})p(\bm{x}_{n+1}|\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n})
\end{eqnarray}
を得る。初期値については本文(13.33)において$n=N$とおき、$\alpha$の定義で置き換えると、
\begin{eqnarray}
p(\bm{z}_{N}|\bm{X}) = \frac{p(\bm{X},\bm{z}_{N})\beta(\bm{z}_{N})}{p(\bm{X})}
\end{eqnarray}
となることから$\beta(\bm{z}_{N})=1$とすればよいことがわかる。

また本文(13..33)の両辺において$\bm{z}_{n}$について和を取ると
\begin{eqnarray}
p(\bm{X}) &=& \sum_{\bm{z}_{n}} \alpha(\bm{z}_{n}) \beta(\bm{z}_{n}) \notag \\
&=& \sum_{\bm{z}_{N}}\alpha(\bm{z}_{N})
\end{eqnarray}
を得る。

次に$\xi(\bm{z}_{n-1},\bm{z}_{n})$については
\begin{eqnarray}
\xi(\bm{z}_{n-1},\bm{z}_{n}) &=& p(\bm{z}_{n-1},\bm{z}|\bm{X}) \notag \\
&=& \frac{p(\bm{X}|\bm{z}_{n-1},\bm{z}_{n})p(\bm{z}_{n-1},\bm{z}_{n})}{p(\bm{X})} \notag \\
&=& \frac{p(\bm{x},\cdots,\bm{x}_{n-1}|\bm{z}_{n-1})p(\bm{x}_{n}|\bm{z}_{n})p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n})p(\bm{z}_{n}|\bm{z}_{n-1})p(\bm{z}_{n-1}) }{p(\bm{X})} \notag \\
&=& \frac{\alpha(\bm{z}_{n-1}p(\bm{x}_{n}|\bm{z}_{n})p(\bm{z}_{n}|\bm{z}_{n-1})\beta(\bm{z}_{n})}{p(\bm{X})}
\end{eqnarray}
を得る。

最後に予測分布については
\begin{eqnarray}
p(\bm{x}_{N+1}|\bm{X}) &=& \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1},\bm{z}_{N+1}|\bm{X}) \notag \\
&=& \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1}|\bm{z}_{N+1})p(\bm{z}_{N+1}|\bm{X}) \notag \\
&=& \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1}|\bm{z}_{N+1})\sum_{\bm{z}_{N}}p(\bm{z}_{N+1},\bm{z}_{N}|\bm{X}) \notag \\
&=& \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1}|\bm{z}_{N+1})\sum_{\bm{z}_{N}}p(\bm{z}_{N+1}|\bm{z}_{N})p(\bm{z}_{N}|\bm{X}) \notag \\
&=& \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1}|\bm{z}_{N+1})\sum_{\bm{z}_{N}}p(\bm{z}_{N+1}|\bm{z}_{N}) \frac{p(\bm{z}_{N},\bm{X})}{p(\bm{X})} \notag \\
&=& \frac{1}{p(\bm{X})} \sum_{\bm{z}_{N+1}}p(\bm{x}_{N+1}|\bm{z}_{N+1})\sum_{\bm{z}_{N}}p(\bm{z}_{N+1}|\bm{z}_{N}) \alpha(\bm{z}_{N})
\end{eqnarray}

\subsection{HMMの積和アルゴリズム}
省略

\subsection{スケーリング係数}
実際にフォワードバックワードアルゴリズムを利用する場合、値が指数関数的に小さくなってしまう場合がある。
そこで$\alpha(\bm{z}_{n})$の規格化された表式
\begin{eqnarray}
\hat\alpha(\bm{z}_{n}) = p(\bm{z}_{n}|\bm{x}_{1},\cdots,\bm{x}_{n}) = \frac{\alpha(\bm{z}_{n})}{p(\bm{x}_{1},\cdots,\bm{x}_{n})}
\end{eqnarray}
を導入する。さらに
\begin{eqnarray}
c_{n} = p(\bm{x}_{n}|\bm{x}_{1},\cdots,\bm{x}_{n-1})
\end{eqnarray}
を定義すると、乗法定理により
\begin{eqnarray}
p(\bm{x}_{1},\cdots,\bm{x}_{n}) = \prod_{m=1}^{n}c_{m}
\end{eqnarray}
を得る。これより
\begin{eqnarray}
\alpha(\bm{z}_{n}) = p(\bm{z}_{n}|\bm{x}_{1},\cdots,\bm{x}_{n})p(\bm{x}_{1},\cdots,\bm{x}_{n}) = \left( \prod_{m=1}^{n}c_{m} \right) \hat{\alpha}(\bm{z}_{n})
\end{eqnarray}
が得られるため、$\alpha$の再帰式に代入することで
\begin{eqnarray}
c_{n}\hat{\alpha}(\bm{z}_{n}) = p(\bm{x}_{n}|\bm{z}_{n})\sum_{\bm{z}_{n-1}}\hat{\alpha}(\bm{z}_{n-1})p(\bm{z}_{n}|\bm{z}_{n-1})
\end{eqnarray}
を得る。同様にして$\beta$についても
\begin{eqnarray}
\hat{\beta}(\bm{z}_{n}) = \frac{\beta(\bm{z}_{n})}{\prod_{m=n+1}^{N}c_{m}} = \frac{p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{z}_{n})}{p(\bm{x}_{n+1},\cdots,\bm{x}_{N}|\bm{x}_{1},\cdots,\bm{x}_{n})}
\end{eqnarray}
と定義すると再帰式は
\begin{eqnarray}
c_{n+1}\hat{\beta}(\bm{z}_{n}) = \sum_{\bm{z}_{n+1}}\hat{\beta}(\bm{z}_{n+1})p(\bm{x}_{n+1}|\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n})
\end{eqnarray}
となり、尤度関数と周辺確率は
\begin{eqnarray}
p(\bm{X}) &=& \prod_{n=1}^{N}c_{n} \notag \\
\gamma(\bm{z}_{n}) &=& \hat{\alpha}(\bm{z}_{n})\hat{\beta}(\bm{z}_{n}) \notag \\
\xi(\bm{z}_{n-1},\bm{z}_{n}) &=& (c_{n})^{-1} \hat{\alpha}(\bm{z}_{n-1})p(\bm{x}_{n}|\bm{z}_{n})p(\bm{z}_{n}|\bm{z}_{n-1})\hat{\beta}(\bm{z}_{n})
\end{eqnarray}
となる。

\subsection{Viterbiアルゴリズム}
ここでは観測データ$\{\bm{x}_{1},\cdots,\bm{x}_{N}\}$が与えられた場合に、最も確からしい$\bm{z}_{n}$の系列を求めることを考える。
\textcolor{blue}{
そこで
\begin{eqnarray}
w(\bm{z}_{n}) = \max_{\bm{z}_{1},\cdots,\bm{z}_{n-1}} \ln p(\bm{x}_{1},\cdots,\bm{x}_{n},\bm{z}_{1},\cdots,\bm{z}_{n})
\end{eqnarray}
と定義すると
\begin{eqnarray}
w(\bm{z}_{n+1}) &=& \max_{\bm{z}_{1},\cdots,\bm{z}_{n}} \ln p(\bm{x}_{1},\cdots,\bm{x}_{n+1},\bm{z}_{1},\cdots,\bm{z}_{n+1}) \notag \\
&=&\max_{\bm{z}_{1},\cdots,\bm{z}_{n}} \left[ \ln p(\bm{x}_{1},\cdots,\bm{x}_{n},\bm{z}_{1},\cdots,\bm{z}_{n}) + \ln p(\bm{z}_{n+1}|\bm{z}_{n}) + \ln p(\bm{x}_{n+1}|\bm{z}_{n+1}) \right] \notag \\
&=& \ln p(\bm{x}_{n+1}|\bm{z}_{n+1}) + \max_{\bm{z}_{n}} \{ \ln p(\bm{z}_{n+1}|\bm{z}_{n}) + w(\bm{z}_{n}) \}
\end{eqnarray}
を得る。また、
\begin{eqnarray}
w(\bm{z}_{1}) = \ln p(\bm{z}_{1}) + \ln p(\bm{x}_{1}|\bm{z}_{1})
\end{eqnarray}
であるため、$n=1$から順番に求めていけば最終的に求めたい量である
\begin{eqnarray}
\max_{\bm{z}_{n}}w(\bm{z}_{n}) = \max_{\bm{Z}}p(\bm{X},\bm{Z})
\end{eqnarray}
を求めることが可能になる。
}

\subsection{隠れマルコフモデルの拡張}
省略

\section{線形動的システム}
この節では以下の確率分布を持つ線形動的システムを考えることにする。
\begin{eqnarray}
p(\bm{z}_{n}|\bm{z}_{n-1}) &=& \mathcal{N}(\bm{z}_{n}|\bm{A}\bm{z}_{n-1},\bm{\Gamma}) \notag \\
p(\bm{x}_{n}|\bm{z}_{n}) &=& \mathcal{N}(\bm{x}_{n}|\bm{C}\bm{z}_{n},\bm{\Sigma}) \notag \\
p(\bm{z}_{1}) &=& \mathcal{N}(\bm{z}_{1}|\bm{\mu}_{0},\bm{P}_{0})
\end{eqnarray}
したがってパラメータ集合は$\bm{\theta}=\{ \bm{A},\bm{\Gamma},\bm{C},\bm{\Sigma},\bm{\mu}_{0},\bm{P}_{0}\}$と表すことができる。

\subsection{LDSにおける推論}
\textcolor{blue}{（LDSはLinear Dynamical Systemのこと。）}
ここでは観測系列で条件付けられた潜在変数の周辺分布を求める問題を考える。HMMの時と同様に
\begin{eqnarray}
\hat\alpha(\bm{z}_{n}) = p(\bm{z}_{n}|\bm{x}_{1},\cdots,\bm{x}_{n})
\end{eqnarray}
を定義しこれを
\begin{eqnarray}
\hat\alpha(\bm{z}_{n}) = \mathcal{N}(\bm{z}_{n}|\bm{\mu}_{n},\bm{V}_{n})
\end{eqnarray}
とおくことにする。その再帰式はHMMの時と同じようにして
\begin{eqnarray}
c_{n}\hat{\alpha}(\bm{z}_{n}) = p(\bm{x}_{n}|\bm{z}_{n}) \int \hat{\alpha}(\bm{z}_{n-1})p(\bm{z}_{n}|\bm{z}_{n-1})d\bm{z}_{n-1}
\end{eqnarray}
となり、これを計算することで
\begin{eqnarray}
\bm{\mu}_{n} &=& \bm{A}\bm{\mu}_{n-1} + \bm{K}_{n}(\bm{x}_{n}-\bm{C}\bm{A}\bm{\mu}_{n-1}) \notag \\
\bm{V}_{n} &=& (\bm{I}-\bm{K}_{n}\bm{C})\bm{P}_{n-1} \notag \\
c_{n} &=& \mathcal{N}(\bm{x}_{n}|\bm{C}\bm{A}\bm{\mu}_{n-1},\bm{C}\bm{P}_{n-1}\bm{C}^{T}+\bm{\Sigma})
\end{eqnarray}
を得る。ただし
\begin{eqnarray}
\bm{P}_{n-1} &=& \bm{A}\bm{V}_{n-1}\bm{A}^{T} + \bm{\Gamma} \notag \\
\bm{K}_{n} &=& \bm{P}_{n-1}\bm{C}^{T}(\bm{C}\bm{P}_{n-1}\bm{C}^{T}+\bm{\Sigma})^{-1}
\end{eqnarray}
である。
初期条件については
\begin{eqnarray}
c_{1}\hat{\alpha}(\bm{z}_{1}) = p(\bm{z}_{1})p(\bm{x}_{1}|\bm{z}_{1})
\end{eqnarray}
であることから
\begin{eqnarray}
\bm{\mu}_{1} &=& \bm{\mu}_{0} + \bm{K}_{1}(\bm{x}_{1}-\bm{C}\bm{\mu}_{0}) \notag \\
\bm{V}_{1} &=& (\bm{I}-\bm{K}_{1}\bm{C})\bm{P}_{0} \notag \\
c_{1} &=& \mathcal{N}(\bm{x}_{1}|\bm{C}\bm{\mu}_{0},\bm{C}\bm{P}_{n-1}\bm{C}^{T}+\bm{\Sigma}) \notag \\
\bm{K}_{1} &=& \bm{P}_{0}\bm{C}^{T}(\bm{C}\bm{P}_{0}\bm{C}^{T}+\bm{\Sigma})^{-1}
\end{eqnarray}
である。

同様に
\begin{eqnarray}
\gamma(\bm{z}_{n}) = p(\bm{z}_{n}|\bm{X})
\end{eqnarray}
を定義し、
\begin{eqnarray}
\gamma(\bm{z}_{n}) = \hat{\alpha}(\bm{z}_{n})\hat{\beta}(\bm{z}_{n}) = \mathcal(\bm{z}_{n}|\hat{\bm{\mu}}_{n},\hat{\bm{V}}_{n})
\end{eqnarray}
によって$\hat{\beta}(\bm{z}_{n})$および$\hat{\bm{\mu}}_{n},\hat{\bm{V}_{n}}$を定義すると再帰式は
\begin{eqnarray}
c_{n+1}\hat{\beta}(\bm{z}_{n}) = \int \hat{\beta}(\bm{z}_{n+1})p(\bm{x}_{n+1}|\bm{z}_{n+1})p(\bm{z}_{n+1}|\bm{z}_{n})d\bm{z}_{n+1}
\end{eqnarray}
となる。これを計算すると
\begin{eqnarray}
\hat{\bm{\mu}} &=& \bm{\mu}_{n} + \bm{J}_{n}(\hat{\bm{\mu}}_{n+1} - \bm{A}\bm{\mu}_{n}) \notag \\
\hat{\bm{V}}_{n} &=& \bm{V}_{n} + \bm{J}_{n}(\hat{V}_{n+1}- \bm{P}_{n}) \bm{J}_{n}^{T}
\end{eqnarray}
を得る。
ここで
\begin{eqnarray}
\bm{J}_{n} = \bm{V}_{n}\bm{A}^{T}(\bm{P}_{n})^{-1}
\end{eqnarray}
である。また２つ組の事後周辺分布については
\begin{eqnarray}
\xi(\bm{z}_{n-1},\bm{z}_{n}) &=& (c_{n})^{-1}\hat{\alpha}(\bm{z}_{n-1})p(\bm{x}_{n}|\bm{z}_{n})p(\bm{z}_{n}|\bm{z}_{n-1})\hat{\beta}(\bm{z}_{n}) \notag \\
&=& \frac{ \mathcal{N}(\bm{z}_{n-1}|\bm{\mu}_{n-1},\bm{V}_{n-1}) \mathcal{N}(\bm{z}_{n}|\bm{A}\bm{z}_{n-1},\bm{\Gamma})
\mathcal{N}(\bm{x}_{n}|\bm{Cz}_{n},\bm{\Sigma})\mathcal{N}(\bm{z}_{n}|\hat{\bm{\mu}}_{n},\hat{\bm{V}}_{n}) }{c_{n}\hat{\alpha}(\bm{z}_{n})} \notag \\
\end{eqnarray}
が成り立つ。これは平均が$[\hat{\bm{\mu}}_{n-1},\hat{\bm{\mu}}_{n}]^{T}$で与えられるガウス分布であり、さらに
\begin{eqnarray}
\mathrm{cov}[\bm{z}_{n-1},\bm{z}_{n}] = \bm{J}_{n-1}\hat{\bm{V}}_{n}
\end{eqnarray}
が成り立つ。

\subsection{LDSの学習}
ここではEMアルゴリズムを用いてパラメータ$\bm{\theta}=\{\bm{A},\bm{\Gamma},\bm{C},\bm{\Sigma},\bm{\mu}_{0},\bm{P}_{0}\}$を推定することを考える。
そのためには確率分布$p(\bm{Z}|\bm{X},\bm{\theta}^{\mathrm{old}})$に基づいた潜在変数の期待値が必要であり、それらは今までの結果から
\begin{eqnarray}
\mathbb{E}[\bm{z}_{n}] &=& \hat{\bm{\mu}_{n}} \notag \\
\mathbb{E}[\bm{z}_{n}\bm{z}_{n-1}^{T}] &=& \hat{\bm{V}}_{n}\bm{J}_{n-1}^{T} + \hat{\bm{\mu}}_{n}\hat{\bm{\mu}}_{n-1}^{T} \notag \\
\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}] &=& \hat{\bm{V}}_{n} + \hat{\bm{\mu}}_{n}\hat{\bm{\mu}}_{n}^{T}
\end{eqnarray}
となる。また、完全データの対数尤度関数は
\begin{eqnarray}
\ln p(\bm{X},\bm{Z}|\bm{\theta}) &=& \ln p(\bm{z}_{1}|\bm{\mu}_{0},\bm{P}_{0}) + \sum_{n=2}^{N}\ln p(\bm{z}_{n}|\bm{z}_{n-1},\bm{A},\bm{\Gamma}) \notag \\
&+& \sum_{n=1}^{N}\ln p(\bm{x}_{n}|\bm{z}_{n},\bm{C},\bm{\Sigma})
\end{eqnarray}
で与えられる。事後分布での期待値は
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}^{\mathrm{old}}) &=& \int \ln p(\bm{X},\bm{Z}|\bm{\theta}) p(\bm{Z}|\bm{X},\bm{\theta}^{\mathrm{old}})d\bm{Z} \notag \\
&=& \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} [\ln p(\bm{X},\bm{Z}|\bm{\theta}) ]
\end{eqnarray}
である。最初に$\bm{\mu}_{0}$と$\bm{P}_{0}$について考える。これらに依存しない項を定数項とすると
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}_{old}) &=& -\frac{1}{2}\ln |\bm{P}_{0}| - \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[ \frac{1}{2}(\bm{z}_{1}-\bm{\mu}_{0})^{T}\bm{P}_{0}^{-1}(\bm{z}_{1}-\bm{\mu}_{0}) \right] + \mathrm{const} \notag \\
\end{eqnarray}
となる。これを最大化するには2.3.4節の結果から
\begin{eqnarray}
\bm{\mu}_{0}^{\mathrm{new}} &=& \mathbb{E}[\bm{z}_{1}] \notag \\
\bm{P}_{0}^{\mathrm{new}} &=& \mathbb{E}[\bm{z}_{1}\bm{z}_{1}^{T}] - \mathbb{E}[\bm{z}_{1}] \mathbb{E}[\bm{z}_{1}^{T}]
\end{eqnarray}
同様に$\bm{A}$と$\bm{\Gamma}$については
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}_{old}) = -\frac{N-1}{2}\ln |\bm{\Gamma}| - \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[\frac{1}{2}\sum_{n=2}^{N}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1})^{T}\bm{\Gamma}^{-1}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1}) \right] + \mathrm{const}
\notag \\
\end{eqnarray}
であり、$\bm{A}$については
\begin{eqnarray}
\bm{A}^{\mathrm{new}} = \left( \sum_{n=2}^{N}\mathrm{E}[\bm{z}_{n}\bm{z}_{n-1}^{T}] \right)  \left( \sum_{n=2}^{N}\mathrm{E}[\bm{z}_{n-1}\bm{z}_{n-1}^{T}] \right)^{-1}
\end{eqnarray}
となる。
\textcolor{blue}{
これは
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}_{old}) &=& -\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[\frac{1}{2}\sum_{n=2}^{N} [z(n)_{i}-z(n-1)_{h}A_{ih} ]\Gamma^{-1}_{ij}[ z(n)_{j}-A_{jk}z(n-1)_{k} ] \right]  + \mathrm{const} \notag \\
\frac{\partial}{\partial A_{\alpha\beta}} Q(\bm{\theta},\bm{\theta}_{old}) &=& -\frac{1}{2}\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[ \sum_{n=2}^{N} z(n-1)_{\beta}\Gamma^{-1}_{\alpha j} [ z(n)_{j}-A_{jk}z(n-1)_{k} ] \right. \notag \\
&+& \left. [ z(n)_{i}-z(n-1)_{h}A_{ih} ]\Gamma^{-1}_{i\alpha} z(n-1)_{\beta} \right] \notag \\
&=& -\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} \left[ \sum_{n=2}^{N} z(n-1)_{\beta}\Gamma^{-1}_{\alpha j}
[z(n)_{j}-A_{jk}z(n-1)_{k}] \right] \notag \\
\frac{\partial}{\partial \bm{A}}Q(\bm{\theta},\bm{\theta}_{old}) &=& -\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} \left[ \sum_{n=2}^{N}\bm{z}_{n-1}[\bm{\Gamma}^{-1}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1})]^{T} \right] \notag \\
&=& -\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} \left[ \sum_{n=2}^{N}\bm{z}_{n-1}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1})^{T}(\bm{\Gamma}^{-1})^{T} \right] = 0
\end{eqnarray}
より、
\begin{eqnarray}
\mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}  \left[ \sum_{n=2}^{N}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1})\bm{z}_{n-1}^{T} \right] = 0
\end{eqnarray}
が成り立つことから得られる。ここで成分表記するときはベクトルの足と区別するため$\bm{z}_{n}$を$z(n)$と書いた。
}
また$\bm{\Gamma}$については
\begin{eqnarray}
\bm{\Gamma}^{\mathrm{new}} &=& \frac{1}{N-1} \sum_{n=2}^{N} \left\{
\mathbb{E} [\bm{z}_{n}\bm{z}_{n}^{T}] - \bm{A}^{\mathrm{new}}\mathbb{E}[\bm{z}_{n-1}\bm{z}_{n}^{T}] \right. \notag \\
&-& \left. \mathbb{E} [\bm{z}_{n}\bm{z}_{n-1}^{T}](\bm{A}^{\mathrm{new}})^{T} +
\bm{A}^{\mathrm{new}}\mathbb{E}[\bm{z}_{n-1}\bm{z}_{n-1}^{T}](\bm{A}^{\mathrm{new}})^{T}
\right\}
\end{eqnarray}
これは
\textcolor{blue}{
\begin{eqnarray}
\frac{\partial }{\partial \Gamma_{\alpha\beta}}\Gamma^{-1}_{ij} = - \Gamma^{-1}_{i\alpha}\Gamma^{-1}_{\beta j} \quad
\frac{\partial }{\partial \Gamma_{\alpha\beta}} \ln |\bm{\Gamma}| = \Gamma^{-1}_{\beta\alpha}
\end{eqnarray}
を用いて
\begin{eqnarray}
\frac{\partial}{\partial \Gamma_{\alpha\beta}}Q(\bm{\theta},\bm{\theta}_{old}) &=&
-\frac{N-1}{2} \Gamma^{-1}_{\beta\alpha} + \frac{1}{2} \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[ \sum_{n=2}^{N}[z(n)_{i}-z(n-1)_{h}A^{T}_{hi}]\Gamma^{-1}_{i\alpha}\Gamma^{-1}_{\beta j}[z(n)_{j}-A_{jk}z(n-1)_{k}] \right] \notag \\
\frac{\partial}{\partial \bm{\Gamma}}Q(\bm{\theta},\bm{\theta}_{old}) &=&
-\frac{N-1}{2} \bm{\Gamma}^{-1} + \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} \left[ \sum_{n=2}^{N} \bm{\Gamma}^{-1}(\bm{z}_{n}-\bm{A}\bm{z}_{n-1})(\bm{z}_{n}^{T}-\bm{z}_{n-1}^{T}\bm{A}^{T})\bm{\Gamma}^{-1}\right] = 0
\end{eqnarray}
となり
\begin{eqnarray}
\bm{\Gamma} = \frac{1}{N-1} \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}} \left[ (\bm{z}_{n}-\bm{A}\bm{z}_{n-1})(\bm{z}_{n}^{T}-\bm{z}_{n-1}^{T}\bm{A}^{T}) \right]
\end{eqnarray}
となることから得られる。
}
最後に$\bm{C}$と$\bm{\Sigma}$については
\begin{eqnarray}
Q(\bm{\theta},\bm{\theta}_{old}) = -\frac{N}{2} \ln |\bm{\Sigma}| - \mathbb{E}_{\bm{Z}|\bm{\theta}^{\mathrm{old}}}
\left[ \frac{1}{2}\sum_{n=2}^{N}(\bm{x}_{n}-\bm{C}\bm{z}_{n})^{T}\bm{\Sigma}^{-1}(\bm{x}_{n}-\bm{C}\bm{z}_{n}) \right] + \mathrm{const} \notag \\
\end{eqnarray}
を微分して
\begin{eqnarray}
\bm{C}^{\mathrm{new}} &=& \left( \sum_{n=1}^{N}\bm{x}_{n}\mathrm{E}[\bm{z}_{n}^{T}] \right)  \left( \sum_{n=1}^{N}\mathrm{E}[\bm{z}_{n}\bm{z}_{n}^{T}] \right)^{-1} \notag \\
\bm{\Sigma}^{\mathrm{new}} &=& \frac{1}{N}\sum_{n=1}^{N} \left\{
\mathbb{E} [\bm{x}_{n}\bm{x}_{n}^{T}] - \bm{C}^{\mathrm{new}}\mathbb{E}[\bm{z}_{n}]\bm{x}_{n}^{T} \right. \notag \\
&-& \left. \bm{x}_{n}\mathbb{E} [\bm{z}_{n}^{T}](\bm{C}^{\mathrm{new}})^{T} +
\bm{C}^{\mathrm{new}}\mathbb{E}[\bm{z}_{n}\bm{z}_{n}^{T}](\bm{C}^{\mathrm{new}})^{T}
\right\}
\end{eqnarray}
を得る。

\subsection{LDSの拡張}
省略

\subsection{粒子フィルタ}
ここでは、ガウス分布でない出力密度を用いる動的システムでも利用可能な、サンプリング手法を用いた推論アルゴリズムを考える。
観測値$\bm{X}_{n}$が与えられたときに、事後分布$p(\bm{z}_{n}|\bm{X}_{n})$から$L$個のサンプルを生成すると、
$\bm{z}_{n}$の関数の期待値は
\begin{eqnarray}
\mathbb{E}[f(\bm{z}_{n})] &=& \int f(\bm{z}_{n})p(\bm{z}_{n}|\bm{X}_{n})d\bm{z}_{n} \notag \\
&=& \int f(\bm{z}_{n})p(\bm{z}_{n}|\bm{x}_{n},\bm{X}_{n-1})d\bm{z}_{n} \notag \\
&=& \int f(\bm{z}_{n})\frac{p(\bm{x}_{n},\bm{z}_{n}|\bm{X}_{n-1})}{p(\bm{x}_{n}|\bm{X}_{n-1})} d\bm{z}_{n} \notag \\
&=& \frac{\int f(\bm{z}_{n})p(\bm{x}|\bm{z}_{n})p(\bm{z}_{n}|\bm{X}_{n-1})d\bm{z}_{n}}{ \int p(\bm{x}_{n}|\bm{z}_{n})p(\bm{z}_{n}|\bm{X}_{n-1})d\bm{z}_{n} } \notag \\
&\approx& \frac{\sum_{l=1}^{l}f(\bm{z}_{n}^{(l)})p(\bm{x}_{n}|\bm{z}_{n}^{(l)}) }{\sum_{l=1}^{l} p(\bm{x}_{n}|\bm{z}_{n}^{(l)}) } \notag \\
&=& \sum_{l=1}^{l}w_{n}^{(l)}f(\bm{z}_{n}^{(l)})
\end{eqnarray}
となる。ここで
\begin{eqnarray}
w_{n}^{(l)} = \frac{p(\bm{x}_{n}|\bm{z}_{n}^{(l)})}{\sum_{m=1}^{L}p(\bm{x}_{n}|\bm{z}_{n}^{(m)})}
\end{eqnarray}
である。
以下保留。

\chapter{モデルの結合}
\section{ベイズモデル平均化}
混合ガウス分布の場合、周辺分布は
\begin{eqnarray}
p(\bm{x}) &=& \sum_{\bm{z}}p(\bm{x},\bm{z}) \notag \\
&=& \sum_{k=1}^{K}\pi_{k}\mathcal{N}(\bm{x}|\bm{\mu}_{k},\bm{\Sigma}_{k})
\end{eqnarray}
と表され、データ集合$\bm{X}=\{\bm{x}_{1},\cdots,\bm{x}_{N}\}$の周辺確率は
\begin{eqnarray}
p(\bm{X}) = \prod_{n=1}^{N}p(\bm{x}_{n}) = \prod_{n=1}^{N}\left[ \sum_{\bm{z}_{n}}p(\bm{x}_{n},\bm{z}_{n}) \right]
\end{eqnarray}
で与えられ、観測されたデータ点$\bm{x}_{n}$ごとに対応する潜在変数$\bm{z}_{n}$が存在する。

一方で、ベイズモデル平均化では$h=1,\cdots,H$で番号付けされたいくつかの異なるモデルが存在し、その事前確立$p(h)$が存在するとする。
この場合周辺分布は
\begin{eqnarray}
p(\bm{X})=\sum_{k=1}^{H}p(\bm{X}|h)p(h)
\end{eqnarray}
で与えれられる。データ集合のサイズが大きくなれば、不確実性は減少し事後確率$p(h|\bm{X})$は漸近的に一つのモデルに集中する。

\section{コミッティ}
あるデータ集合$\bm{X}=\{\bm{x}_{1},\cdots,\bm{x}_{N}\}$に対して、$M$個のブートストラップ集合（元の集合からランダムに$N$個のデータを選んだもの。元のデータ数が$N$なので、一つの集合に重複したデータが含まれたりもする。）を生成し、$M$個の予測モデル$y_{m}(\bm{x})$を訓練したとすると、コミッティの予測は
\begin{eqnarray}
y_{\mathrm{COM}}(\bm{x}) = \frac{1}{M}\sum_{m=1}y_{m}(\bm{x})
\end{eqnarray}
で与えられる。予測しようとする本当の回帰関数が$h(\bm{x})$で与えられる場合、それぞれのモデルの誤差は
\begin{eqnarray}
y_{m}(\bm{x}) = h(\bm{x}) + \epsilon_{m}(\bm{x})
\end{eqnarray}
で与えられる。すると平均二乗誤差は
\begin{eqnarray}
\mathbb{E}_{\bm{x}}\left[\{y_{m}(\bm{x})-h(\bm{x})\}^{2}\right] = \mathbb{E}_{\bm{x}}[\epsilon_{m}(\bm{x})^{2}]
\end{eqnarray}
で与えられ、各モデルの平均二乗誤差の平均は
\begin{eqnarray}
E_{\mathrm{AV}} = \frac{1}{M}\sum_{m=1}^{M} \mathbb{E}_{\bm{x}}[\epsilon_{m}(\bm{x})^{2}]
\end{eqnarray}
となる。一方で、コミッティについての誤差の期待値は
\begin{eqnarray}
E_{\mathrm{COM}} &=& \mathbb{E}_{\bm{x}} \left[ \left\{ \frac{1}{M}\sum_{m=1}^{M}y_{m}(\bm{x})-h(\bm{x}) \right\}^{2} \right] \notag \\
&=& \mathbb{E}_{\bm{x}} \left[ \left\{ \frac{1}{M}\sum_{m=1}^{M}\epsilon_{m}(\bm{x}) \right\}^{2} \right]
\end{eqnarray}
となる。もし、誤差の平均が$0$で無相関、すなわち
\begin{eqnarray}
\mathbb{E}_{\bm{x}}[\epsilon_{m}(\bm{x})] &=& 0 \notag \\
\mathbb{E}_{\bm{x}}[\epsilon_{m}(\bm{x})\epsilon_{l}(\bm{x})] &=& 0
\end{eqnarray}
である場合
\begin{eqnarray}
E_{\mathrm{COM}} = \frac{1}{M} E_{\mathrm{AV}}
\end{eqnarray}

\section{ブースティング}
ここでは訓練データが二値の目標変数$t_{1},\cdots,t_{N}$を伴う入力ベクトル$\bm{x}_{1},\cdots,\bm{x}_{N}$で、分類器$y(\bm{x})\in \{-1,1\}$を構成するAdaBoostアルゴリズムを扱う。具体的には以下のようにする。
\begin{enumerate}
\item $n=1,\cdots,N$のデータ重み係数$\{w_{n}\}$を$w_{n}^{(1)}=1/N$に初期化
\item $m=1,\cdots,M$について以下を繰り返す。
\begin{itemize}
\item[(a)] 分類器$y_{m}(\bm{x})$を
\begin{eqnarray}
J_{m} = \sum_{n=1}^{N}w_{n}^{(m)}I(y_{m}(\bm{x}_{n})\neq t_{n})
\end{eqnarray}
を最小化するように最適化。（これはできると仮定する。）ここで$I(y_{m}(\bm{x}_{n})\neq t_{n})$は$y_{m}(\bm{x}_{n})\neq t_{n}$のときには$1$でそうでないときは$0$である。
\item[(b)] 次の値を計算する
\begin{eqnarray}
\epsilon_{m} &=& \frac{\sum_{n=1}^{N}w_{n}^{(m)}I(y_{m}(\bm{x}_{n})\neq t_{n})}{\sum_{n=1}^{N}w_{n}^{(m)}} \notag \\
\alpha_{m} &=& \ln\left\{\frac{1-\epsilon_{m}}{\epsilon_{m}}\right\}
\end{eqnarray}
\item[(c)] データ点の重み係数を以下の式で更新
\begin{eqnarray}
w_{n}^{(m+1)} = w_{n}^{(m)}\exp\{\alpha_{n}I(y_{m}(\bm{x}_{n})\neq t_{n})\}
\end{eqnarray}
\end{itemize}
\item 以下の式で最終モデルの予測を構成する
\begin{eqnarray}
Y_{M}(\bm{x}) = \mathrm{sign}\left(\sum_{m=1}^{M}\alpha_{m}y_{m}(\bm{x})\right)
\end{eqnarray}
\end{enumerate}

\subsection{指数誤差の最小化}
ここでは前節のAdaBoostが指数誤差関数の逐次的最適化から解釈できることを見る。
指数誤差関数は
\begin{eqnarray}
E = \sum_{n=1}^{N}\exp\{-t_{n}f_{m}(\bm{x}_{n})\}
\end{eqnarray}
で定義され、
\begin{eqnarray}
f_{m}(\bm{x}) &=& \frac{1}{2}\sum_{l=1}^{m}\alpha_{l}y_{l}(\bm{x})
\end{eqnarray}
である。ここで、
ベース分類器$y_{1}(\bm{x}),\cdots,y_{m-1}(\bm{x})$とそれらの係数$\alpha_{1},\cdots,\alpha_{m-1}$が固定されていると仮定し、$\alpha_{m}$と$y_{m}(\bm{x})$について最小化を行う。
誤差関数から$y_{m}(\bm{x})$の寄与を分離すると
\begin{eqnarray}
E &=& \sum_{n=1}^{N} \exp \left\{ -t_{n}f_{m-1}(\bm{x}_{n}) - \frac{1}{2}t_{n}\alpha_{m}y_{m}(\bm{x}_{n}) \right\} \notag \\
&=& \sum_{n=1}^{N} w_{n}^{(m)} \exp \left\{  - \frac{1}{2}t_{n}\alpha_{m}y_{m}(\bm{x}_{n}) \right\} \notag \\
&=& (e^{\alpha_{m}/2} - e^{-\alpha_{m}/2})\sum_{n=1}^{N}w_{n}^{(m)}I(y_{m}(\bm{x}_{n})\neq t_{n}) + e^{-\alpha_{m}/2}\sum_{n=1}^{N}w_{n}^{(m)}
\end{eqnarray}
となる。これを$y_{m}(\bm{x})$について最小化することはAdaBoostの2.(a)と等価であり、$\alpha_{m}$の最小化は2.(b)に帰着する。
そして新しい重みについては
\begin{eqnarray}
w_{n}^{(m+1)} &=& w_{n}^{(m)}\exp\left\{ -\frac{1}{2}t_{n}\alpha_{m}y_{m}(\bm{x}_{n}) \right\} \notag \\
&=& w_{n}^{(m)}\exp(-\alpha_{m}/2)\exp\{\alpha_{m}I(y_{m}(\bm{x}_{n})\neq t_{n})\}
\end{eqnarray}
となり、$n$に独立な$\exp(-\alpha_{m}/2)$が無視できるため、やはりAdaBoostのものに一致する。

\subsection{ブースティングのための誤差関数}
省略

\section{木構造モデル}
省略

\section{条件付き混合モデル}
省略





\end{document}