TL;DR
- ガウス的なアプローチで、予測分布を出したいときは事後確率分布をパラメータ$w$に関して積分する計算が必要。だが、事後確率分布がシグモイド関数/ソフトマックス関数がごちゃごちゃするような形だと積分できない... ⇒ややこしい形の分布を積分可能なガウス分布に近似することで計算を簡単にしよう！
4.4 ラプラス近似
入力$x$を与えたときのあるクラス$C_k$に対する予測分布は、事後確率分布$p(w|\mathbf{t})$を$w$について積分する計算が必要である（参照⇒part10。これは3章回帰問題の予測分布だが、4.5章で識別モデルの予測分布も同様であることをみる）。
しかし、事後確率がシグモイド関数やソフトマックス関数で与えられるとモデル化しており（参照⇒part15）、ガウス分布ではないので、$w$に関して正確に積分することができない。
なので、事後確率分布をガウス分布で近似することを考える。このように確率密度分布をガウス分布で近似することをラプラス近似という。
まず、連続な1変数zの場合
以下の式で定義される分布$p(z)$を仮定する。
$$
p(z) = \frac{1}{Z} f(z) \tag{4.125}
$$
ここで、$Z$は$f(z)$を$z$について積分した値で、$ \int p(z) dz = 1$となるように存在する正規化係数である。
ラプラス近似でやりたいことは、（任意の）分布$p(z)$の最頻値を中心とするガウス分布による近似$q(z)$を見つけること。なので$p(z)$の最頻値を見つけよう。
最頻値より$z$の値を少しでも大きく/小さくする（グラフでいえば少しでも右に/左に動く）と$z$の値は小さくなるのだから、最頻値$z_0$での$p(z_0)$の傾きは0であり、微分は0になることを頭に入れておく。すなわち、
$$
p(z_0)' = f(z_0)' = 0 \tag{A}
$$
ガウス分布は、その対数が変数の2次関数で書ける1。なので、任意の分布$f(z)$の対数をテイラー展開を使って、2次関数で表そう。こうすると、近似分布$q(z)$がガウス分布の形になる。
最頻値$z_0$を中心とした$ \ln f(z) $のテイラー展開は
\begin{align}
\ln f(z) =& \ln f(z_0) + \frac{1}{1!} \{ \ln f(z_0) \}' (z-z_0) +
\frac{1}{2!} \{ \ln f(z_0) \}'' (z-z_0)^2 + (これ以降は微小のため無視)
\end{align}
ここで、式(A)より
$$
\{ \ln f(z_0) \}' = \frac{ f(z_0)' }{ f(z_0) } = 0
$$
となるため、
\begin{align}
\ln f(z) \simeq \ln f(z_0) - \frac{1}{2} A (z - z_0)^2 \tag{4.127} \\
A = - \{ \ln f(z_0) \}'' \tag{4.128}
\end{align}
このテイラー展開を指数をとると、
\begin{align}
f(z) &\simeq \exp \left\{
\ln f(z_0) - \frac{A}{2} (z - z_0)^2
\right\} \\
&= \exp \{ \ln f(z_0) \} \cdot \exp \{ - \frac{A}{2} (z - z_0)^2 \} \\
&= f(z_0) \exp \{ - \frac{A}{2} (z - z_0)^2 \} \tag{4.129}
\end{align}
任意の分布$f(z)$がガウス分布と同じexpの2次関数の形で書くことができた。なので、近似分布$q(z)$は式
(4.129)に正規化係数を足した下式となることが導出できた。
$$
q(z) = \left( \frac{A}{2 \pi} \right)^{1/2} \exp
\left\{ - \frac{A}{2} (z - z_0)^2 \right\} \tag{4.130}
$$
ただし、式(4.130)のルート（1/2乗）の計算ができなくなるので$A>0$、つまり$f(z_0)$での2階微分が負の場合しかこの近似式は使えない。
この式を使って（事後確率の）シグモイド関数とexp関数の積をガウス分布に近似したものが下図。オレンジがシグモイド関数で、赤が近似ガウス分布。
多次元（M次元）の場合
1変数の場合と同様、任意の分布$p(z) = f(z)/Z$（$z$は多次元）の最頻値$z_0$を求めよう。点$z_0$では、勾配$\nabla f(z)=0$となる。$z_0$の周りで$ \ln f(z) $をテイラー展開すると
\begin{align}
\ln f(z) \simeq \ln f(z_0) - \frac{1}{2} (z - z_0)^T A (z - z_0) \tag{4.131} \\
\end{align}
ここで、1変数のときの議論より、行列$A$は$ \ln f(z_0) $の2階微分に相当するものでヘッセ行列となる。つまり、
$$
A = - \nabla \nabla \ln f(z_0) \tag{4.132}
$$
これまた1変数のときと同様に、式(4.131)の対数をとり、正規化係数を求めることで以下の近似分布が得られる。
$$
q(z) = \frac{|A|^{1/2}}{(2 \pi)^{M/2}} \exp \left\{
- \frac{1}{2} (z - z_0)^T A (z - z_0)
\right\} = N(z | z_0, A^{-1} ) \tag{4.134}
$$
行列式$|A|$が正となるときに、この近似式が使える。行列式が正となるのは、$z_0$が局所最大であって、局所最小または鞍点でない場合らしい。
ラプラス近似の利点と欠点
利点
- 観測データが増えるほど、ガウス分布による近似が良くなることが期待できる（データが多い状況で有用）
欠点
- 現実の分布は多峰的であることがおおく、単峰のガウス分布による近似では分布を局所的にしか捉えられない
- ガウス分布による近似のため、実数変数の場合のみ適用可能
4.4.1 モデルの比較とBIC
上節で導出した近似を使い、関数$f(z)$の積分値である正規化係数$Z$を近似することができる。
expの部分の積分はガウス分布の正規化係数の逆数が出てくることを利用している（⇒part4）。
つまり、ガウス分布の積分は1になることから
～～～
$\int$ (正規化係数) exp(ほにゃほにゃ) = 1
$ \int $ exp(ほにゃほにゃ) = (正規化係数)^{-1}
～～～
ことを利用して：
\begin{align}
Z = \int f(z) dz &\simeq f(z_0) \int \exp \left\{
- \frac{1}{2} (z - z_0)^T A (z - z_0)
\right\} dz \\
&= f(z_0) \frac{(2 \pi)^{M/2}}{|A|^{1/2}} \tag{4.135}
\end{align}
つづいて、データ集合$D$と（複数の）パラメータ$ \{ \theta_i \}$を持つモデル集合$ \{ M_i \} $を考える。各モデルに対し、尤度関数$ p(D | \theta_i, M_i)$を定義する。
パラメータ$ \{ \theta_i \}$に関する事前確率$p(\theta_i | M_i) $を導入するなら、モデルエビデンスが計算できる：
$$
p(D | M_i) = \int p(D | \theta, M_i) p(\theta | M_i) d\theta \tag{4.136}
$$
（参照⇒part11。式(3.68)の$w$を$\theta$に置き換えれば同じ式！）
いろんなモデルからどのモデルがいいかを選択するとき、このモデルエビデンスが重要な働きを示すということを3.4節「ベイズモデル比較」で学んだので、もう少し見てみよう。
演習4.22 ラプラス近似の下での対数モデルエビデンスの式を導出
$ p(D | M_i) $を$Z$, $ p(D | \theta, M_i) p(\theta | M_i) $を$ f(\theta) $として、式(4.135)を使うと：
\begin{align}
p(D) & \simeq f( \theta_{MAP} ) \frac{(2 \pi)^{M/2}}{|A|^{1/2}} =
p(D | \theta_{MAP}, M_i) p( \theta_{MAP} | M_i) \frac{(2 \pi)^{M/2}}{|A|^{1/2}} \\
\\
\ln p(D) & \simeq \ln p(D | \theta_{MAP}, M_i) + p( \theta_{MAP} | M_i) + \frac{M}{2} \ln (2 \pi) - \frac{1}{2} \ln |A| \tag{4.137}
\end{align}
ここで、$ \theta_{MAP} $は事後確率分布の最頻値での$ \theta $の値である。
$A$は式(4.128)の$A = - \{ \ln f(z_0) \}''$という定義から
$$
A = - \nabla \nabla p(D | \theta_{MAP} ) p(\theta_{MAP})
= - \nabla \nabla \ln p( \theta_{MAP} | D) \tag{4.138}
$$
中辺までの式変形は定義そのままなのでわかる。中辺から右辺への式変形が分からない。右辺には$p(D)$がかからないだろうか...？
ひとまず、この式を認めて行列$A$はパラメータに関する事後確率$p (\theta_{MAP} | D)$の負の対数の2階微分で求まると理解しておこう。
式(4.137)の右辺第1項は、最適なパラメータを使用して評価した対数尤度であり、残り3つの項はモデルの複雑さ対するペナルティであり、オッカム係数と呼ばれる。
この式に荒い近似を使うと、モデル選択の指標としてよく知られる以下のBICが導出できるらしい（演習問題だがわからないので略）。
\begin{align}
\text{BIC} &= \ln p(D) \simeq \ln p(D | \theta_{MAP} ) - \frac{1}{2} M \ln N \tag{4.139} \\
\text{AIC} &= \ln p(D | \theta_{MAP} ) - M \tag{1.73}
\end{align}
AIC, BICとも最適なモデルを選択するための指標としてよく用いられる。AICと比較し、BICの方がモデルの複雑さにより重いペナルティを科している。
4.5 ベイズロジスティック回帰
ロジスティック回帰モデルには、$w$に関する積分が大変なシグモイド関数/ソフトマックス関数がいっぱいでてくるので、厳密にベイズ推論を適用するのは難しい。そこで、ラプラス近似を使っていこう。
4.5.1 ラプラス近似
計算を簡単にするため、パラメータ$w$に関する事後分布をガウス分布で近似したい。なので、事前分布もガウス分布と仮定しよう。
事前ガウス分布は一般形：
$$
p(w) = N(w | m_0, S_0) \tag{4.140}
$$
ここで、$m_0, S_0$はある固定のハイパーパラメータである。$w$の事後分布は、
$$
p(w | \mathbf{t} ) \propto p(w) p(\mathbf{t} | w) \tag{4.141}
$$
で与えられる。両辺の対数をとって、対数事後確率分布を導出する。
\begin{align}
尤度関数p( \mathbf{t} | w) &= \prod_{n=1}^N \sigma(w^T \phi_n)^{t_n} \{ 1 - \sigma(w^T \phi_n) \}^{1-t_n} \tag{4.89} \\
\\
事後確率分布p(w | \mathbf{t} ) &= - \frac{1}{2} (w - m_0)^T S_0^{-1} (w - m_0) +
\sum_{n=1}^N \{ t_n \ln y_n + (1 - t_n) \ln (1 - y_n) \} + const. \tag{4.142}
\end{align}
ここで、$y_n = \sigma(w^T \phi_n) $である。
式(4.134)より、
- ラプラス近似したガウス分布の平均は式(4.134)でいう$z_0$であり、$z_0$は最頻値のことである。なので、ガウス分布の形を仮定している事後分布を最大化するMAP（最大事後確率）解$w_{MAP}$を求め、それがガウス分布の平均を定義する。
- ラプラス近似したガウス分布の精度行列は行列$A$であり、式(4.138)より負の対数事後確率の2階微分で与えられる。つまり、 $$ S_N^{-1} = - \nabla \nabla \ln p(w| \mathbf{t} ) = S_0^{-1} + \sum_{n=1}^N y_n ( 1 - y_n) \phi_n \phi_n^T \tag{4.143} $$
よって、事後確率分布のガウス分布による近似は：
$$
q(w) = N(w | w_{MAP}, S_N) \tag{4.144}
$$
あとは、予測分布を得るために、この近似ガウス分布を周辺化（積分）する作業がある ⇒ 次節
4.5.2 予測分布
＜おことわり＞式変形断念しました。結果だけ簡単に記すのみになってます
2クラス問題について考える。
特徴ベクトル$\phi(x)$が入力された際のクラス$C_1$に対する予測分布は、事後確率分布$p(w|\mathbf{t})$に関して周辺化して得られる。つまり、
$$
p(C_1 | \phi, \mathbf{t} ) = \int p(C_1 | \phi, w)p(w|\mathbf{t})dw
\simeq \int \sigma(w^T \phi) q(w) dw \tag{4.145}
$$
クラス$C_2$に対する予測分布は$ p(C_2 | \phi, \mathbf{t} ) = 1 - p(C_1 | \phi, \mathbf{t} ) $で与えられる。
事後分布をガウス分布に近似したが、式(4.145)はシグモイド関数×ガウス関数の積分となっており、まだややこしい。なので、いろいろ数学テクニックを駆使して導出した予測分布の近似結果のみ以下に記す2。
\begin{align}
p(C_1 | \phi, \mathbf{t} ) &= \int \sigma(w^T \phi) N(w^T \phi | \mu_a, \sigma_a^2 ) \tag{4.151} \\
\mu_a &= w_{MAP}^T \phi \tag{4.149} \\
\sigma_a^2 &= \phi^T S_N \phi \tag{4.150}
\end{align}
ただ、この予測分布の近似結果式も解析的に評価できないので、プロビット関数の逆関数というもの3を使って近似式を用いたりするらしい。近似のオンパレードになるのでした。