TL;DR

ハイパーパラメータをどうやって決める？ ⇒ 周辺対数尤度（＝モデルエビデンス）を最大にするようなハイパーパラメータの値を推定値とする方法があり、エビデンス近似という

3.5 エビデンス近似

これまでに重みパラメータ$w$以外に、ハイパーパラメータ$ \alpha, \beta $を導入していた。

\begin{align}

p( t | x, w, \beta ) = N( t | y(x, w), \beta^{-1} ) \tag{3.8} \\
p(w | \alpha ) = N(w | 0, \alpha^{-1} I ) \tag{3.52}

\end{align}

つまり、$ \alpha $は重みパラメータ$w$の事前分布がどの程度ばらつくかの仮定を示すガウス分布の精度、$ \beta $は観測値$t$にどの程度のガウスノイズが乗っているかの仮定を示すガウス分布の精度である。

本当は$w$だけでなく、これらのハイパーパラメータも周辺化した下式の予測分布を使いたい。が、積分3つもあることから予想できるように、この計算は非常に難しい。

$$
p(t | \mathbf{t} ) = \int \int \int p(t | w, \beta) p(w | \mathbf{t}, \alpha, \beta) p(\alpha, \beta | \mathbf{t} ) dw d\alpha d\beta \tag{3.74}
$$

なので、$w$に関して得られた周辺尤度（＝前節のモデルエビデンス）を最大にするようなハイパーパラメータを求めるという近似法について議論する。この近似法はエビデンス近似と呼ばれる。

さて、ハイパーパラメータ$ \alpha, \beta $が分布でなく、ある1つの値$ \hat{\alpha}, \hat{\beta} $に決まったと仮定すれば式(3.75)の予測分布は以下のように近似できる。

$$
p(t | \mathbf{t} ) \simeq p( t | \mathbf{t}, \hat{\alpha}, \hat{\beta} )
= \int p(t | w, \hat{\beta} ) p(w | \mathbf{t}, \hat{\alpha}, \hat{\beta} ) dw \tag{3.75}
$$

ハイパーパラメータに関する事後分布$ p(\alpha, \beta | \mathbf{t} ) $はベイズの定理により

$$
p( \alpha, \beta | \mathbf{t} ) \propto p( \mathbf{t} | \alpha, \beta ) p( \alpha, \beta ) \tag{3.76}
$$

だが、データを見る前に$ \alpha, \beta $がどんな値をとりそうか、といった知見は普通ない。すなわち、事前分布$ p( \alpha, \beta ) $は未知なので平坦な分布であると仮定して、事後分布を最大にする代わりに尤度$ p( \mathbf{t} | \alpha, \beta ) $を最大にするような$ \alpha, \beta $を推定するアプローチをとる。

以上の関係性を図に整理した。

3.5.1 エビデンス関数の評価

一旦$ \alpha, \beta $の推定から離れて、周辺尤度（＝前節のモデルエビデンス）の導出および周辺尤度によってどのようなモデルが選択されるか（3.4節の続き）をみる。

$ \alpha $が決まったとすると$w$の分布が決まり、$w$と$ \beta $が決まったとするとその重み$w$での予測分布が決まる。これらの同時分布の$w$についての平均をとる＝$w$について周辺化することで周辺尤度が得られる。

\begin{align}

p( \mathbf{t} | \alpha, \beta ) &= \int p( \mathbf{t} | w, \beta ) p(w | \alpha ) dw \tag{3.77} \\

\end{align}

さて、式(3.77)の積分を解析的に求めるために違う形式に書き直そう（演習3.17）。

\begin{align}

\ln p( \mathbf{t} | w, \beta ) &= \frac{N}{2} \ln \beta - \frac{N}{2} \ln ( 2 \pi ) - \frac{\beta}{2} \sum_n \{ t_n - w^T \phi (x_n) \}^2 \tag{3.11} \\

p( \mathbf{t} | w, \beta ) &= \exp \left[
 \frac{N}{2} \ln \beta - \frac{N}{2} \ln ( 2 \pi ) - \frac{\beta}{2} \sum_n \{ t_n - w^T \phi (x_n) \}^2
\right] \\

&= \exp \left[
\ln \left( \frac{\beta}{2 \pi} \right)^{N/2} - \frac{\beta}{2} \sum_n \{ t_n - w^T \phi (x_n) \}^2
\right] \\

&= \left( \frac{\beta}{2 \pi} \right)^{N/2} \exp \left[
- \frac{\beta}{2} \sum_n \{ t_n - w^T \phi (x_n) \}^2
\right] \\

&= \left( \frac{\beta}{2 \pi} \right)^{N/2} \exp \left[
- \frac{\beta}{2} || \mathbf{t} - \Phi w ||^2
\right] \\ \tag{A} 

\end{align}

以下の事前分布のガウス分布の書き下しにおいて、行列式部分の書き下し方はpart2参照。
参照：part2, MAP推定のところ

\begin{align}

p(w | \alpha) &= N(w | 0, \alpha^{-1} I ) \\
&= \left( \frac{1}{ 2 \pi} \right)^{M/2} \frac{1}{ \alpha^{- M/2 } } \exp \left( 
- \frac{\alpha}{2} w^T w
\right) \\

&= \left( \frac{ \alpha }{ 2 \pi} \right)^{M/2} \exp \left( 
- \frac{\alpha}{2} w^T w \tag{B}
\right) \\

& \text{周辺尤度の積分の中身は式(A), (B)を掛けたもの。wに依存しない項を積分の外に出すと} \\

p( \mathbf{t} | \alpha, \beta ) &= \left( \frac{\beta}{2 \pi} \right)^{N/2}
\left( \frac{ \alpha }{ 2 \pi} \right)^{M/2}
\int \exp \{ -E(w) \} dw \tag{3.78} \\

E(w) &= \beta E_D (w) + \alpha E_w (w) \\
&= \frac{ \beta}{2} || \mathbf{t} - \Phi w ||^2 + \frac{ \alpha}{2} w^T w \tag{3.79}

\end{align}

つづいて式(3.79)を$w$について平方完成することで以下の式に変形する（平方完成はけっこうややこしいが、以下の変形後の式と変形前の式(3.79)が同じ式であることは比較的簡単に確かめられる）。

\begin{align}

E(w) &= E(m_N) + \frac{1}{2} (w - m_N)^T A (w - m_N) \tag{3.80} \\
A &= \alpha I + \beta \Phi^T \Phi \tag{3.81} \\
E(m_N) &= \frac{\beta}{2} || \mathbf{t} - \Phi m_N ||^2 + \frac{\alpha}{2} m_N^T m_N \tag{3.82} \\
m_N &= \beta A^{-1} \Phi^T \mathbf{t} \tag{3.84}

\end{align}

式(3.84)は3.3節「ベイズ線形回帰」で行ったパラメータ$w$に関する事後分布の平均と同じ式になっているので、事後分布の平均を表している。

以上の変形により式(3.78)の$w$に関する積分の中身の$E(w)$が$w$によらない項$E(m_N)$と$w$による項に分けることが出来た。
$E(m_N)$は積分の外に出すことができる。$w$による項はガウス分布の正規化係数がない形をしているので、正規化係数つきガウス分布の積分が1になることを用いて、$w$に関する積分を次のように評価できる（演習3.19）。

\begin{align}

ガウス分布の正規化係数をZとして \\
\int Z \exp \{ \cdots \} dw &= 1 \\
\int \exp \{ \cdots \} dw &= 1 / Z \\

さて \exp \left\{ - \frac{1}{2} (w - m_N)^T A (w - m_N) \right\} を指数部分にもつ \\
ガウス分布の正規化係数Zは\frac{1}{ (2 \pi)^{M/2} } \frac{1}{ |A|^{1/2} }であるから \\
\int \exp \left\{ - \frac{1}{2} (w - m_N)^T A (w - m_N) \right\} dw &= 
( 2 \pi )^{M/2} |A|^{-1/2} \\

周辺尤度は p( \mathbf{t} | \alpha, \beta ) &= \left( \frac{\beta}{2 \pi} \right)^{N/2}
\left( \frac{ \alpha }{ 2 \pi} \right)^{M/2}
\int \exp \{ -E(w) \} dw \tag{3.78} \\

&= \left( \frac{\beta}{2 \pi} \right)^{N/2}
\left( \frac{ \alpha }{ 2 \pi} \right)^{M/2}
\exp \{ -E(m_N) \} \int \exp \left\{ - \frac{1}{2} (w - m_N)^T A (w - m_N) \right\} \\

&= \left( \frac{\beta}{2 \pi} \right)^{N/2}
\left( \frac{ \alpha }{ 2 \pi} \right)^{M/2}
\exp \{ -E(m_N) \} ( 2 \pi )^{M/2} |A|^{-1/2} \\

対数をとると \\

\ln p( \mathbf{t} | \alpha, \beta ) &= \frac{M}{2} \ln \alpha + \frac{N}{2} \ln \beta
- E(m_N) - \frac{1}{2} \ln |A| - \frac{N}{2} \ln (2 \pi ) \tag{3.86}

\end{align}

ようやく対数周辺尤度（＝対数モデルエビデンス）を積分表現を含まない具体的な式の形に変形できたので、3.4節「ベイズモデル比較」で議論したようにこの値を基にどの複雑さのモデルが選択されるかを見てみよう。

サイン関数へのフィッティングを何次の多項式モデルですればいいかという問題を考える。ここでは、ハイパーパラメータ$\alpha$（とたぶん$\beta$も）がある値であると決め打ちする。

横軸にモデルの次元数、縦軸に対数モデルエビデンスをプロットすると下図左になる。

3次の項まで使ったとき、真のサイン関数（下図右の緑線）のような「値が上がって下がってまた上がる」が表現できるようになり、フィッティング度合いが大きく向上するため、対数モデルエビデンスも大きく増加している。4次以上を使うと、フィッティング度合いは若干改善するものの、モデルの複雑度に対するペナルティ項の効果の方が強く働き、対数モデルエビデンスは減少している。こうして、中程度の複雑度のモデルが選択される。

3.5.2 エビデンス関数の最大化

式(3.86)でようやく積分表現を含まない具体的な形の対数周辺尤度の式が導出できたので、これが最大となるような$ \alpha, \beta $の推定式を導出しよう。

まずは$ \alpha $の推定から。
流れは最尤推定と同じように $\alpha$について偏微分＝０の方程式を解く。
まず、対数周辺尤度の$ -\frac{1}{2} \ln |A| $項の偏微分について考えるが、以下の数学テクニックを使う。

行列式は固有ベクトルの積で書ける¹
対称行列$A, B$の和$ C = A + B$である行列$C$の固有値の和は、$A$の固有値の和と$B$の固有値の和を足したものと等しい²

式(3.81)右辺に現れる行列$ \beta \Phi^T \Phi $の固有値が$ \lambda_i $としよう。この行列は対称行列であることに留意。なので、対称行列の和である行列Aの固有値は$ \prod_i ( \lambda_i + \alpha) $と書ける。

よって

$$
\frac{d}{d \alpha} \ln |A| = \frac{d}{d \alpha} \ln \prod_i ( \lambda_i + \alpha)
= \frac{d}{d \alpha} \ln \sum_i ( \lambda_i + \alpha)
= \sum_i \frac{1}{ \lambda_i + \alpha } \tag{3.88}
$$

この結果を用いて、対数周辺尤度の$\alpha$に関する偏微分＝０　を解くと

\begin{align}

0 &= \frac{M}{2 \alpha} - \frac{1}{2} m_N^T m_N - \frac{1}{2} \sum_i \frac{1}{ \lambda_i + \alpha} \tag {3.89} \\
2 \alphaを掛け、整理すると \\
\alpha m_N^T m_N &= M - \alpha \sum_i \frac{1}{ \lambda_i + \alpha}\\
&= \sum_i \left( \frac{ \lambda_i + \alpha }{ \lambda_i + \alpha } - \frac{ \alpha }{ \lambda_i + \alpha } \right)\\

&= \sum_i \frac{ \lambda_i }{ \alpha + \lambda_i } = \gamma \tag{3.90, 3.91} \\

\therefore \alpha &= \frac{ \gamma }{ m_N^T m_N } \tag{3.92}


\end{align}

式(3.92)を満たす$ \alpha $が周辺尤度を最大にする$ \alpha $である。しかし、右辺の$ \gamma, m_N $も$ \alpha $に依存するため、1発でこの式を解くことができない。なので、次のような繰り返し手順により$ \alpha $を求める。

$ \alpha $の初期値を適当に決める
$ m_N = \beta S_N \Phi^T \mathbf{t} $ (3.54) から$m_N$を求める
式(3.90)から$\gamma$を求める
$ m_N, \gamma $から$ \alpha $を推定
推定された$ \alpha $から$m_N, \gamma, \alpha$を再推定する過程を収束するまで繰り返す

続いて$ \beta $の推定だが、こちらは$ \alpha $の推定手順と同じ。

行列$ \beta \Phi^T \Phi $の固有値が$ \lambda_i $としたことから、$ \lambda_i $が$ \beta $に比例することに注意すれば、$ d \lambda_i / d \beta = \lambda_i / \beta $より

$$
\frac{d}{d \beta} \ln |A| = \frac{d}{d \beta} \sum_i \ln ( \lambda_i + \alpha )
= \frac{1}{\beta} \sum_i \frac{ \lambda_i }{ \lambda_i + \alpha } = \frac{ \gamma}{ \beta } \tag{3.93}
$$

この結果を用いて対数尤度の$\beta$に関する偏微分＝０、の式を整理すれば

$$
\frac{1}{\beta} = \frac{1}{ N - \gamma } \sum_{n=1}^N
\{ t_n - m_N^T \phi (x_n) \}^2 \tag{3.95}
$$

3.5.3 有効パラメータ数

（よく分からないので省略。今後の章の内容とは独立しているのできっと大丈夫）

3.6 固定された基底関数の限界

本章の線形回帰モデルでは、固定された非線形基底関数を線形結合したモデルを考えてきた。
これにより、"パラメータに関する微分＝０"の方程式を解くとパラメータの最尤推定解が得られるなど、メリットがあった一方、問題点もある。

問題点の１つは、次元の呪いである。入力空間の次元数$D$に応じて指数的に基底関数の数を増やしていく必要がある。

この問題に対するアプローチとして、後の章ではニューラルネットなどのより複雑な問題なモデルを見ていきますよ～、という予告。

テキストの付録式(C.47)やp.79 @ 2.3節ガウス分布に記述あり ↩
テキストにはこの説明が見つけられなかったが、ググるとこんなページがありました。https://mathoverflow.net/questions/4224/eigenvalues-of-matrix-sums ↩