ベイズ線形回帰 - 機械学習基礎理論独習

やりたい事

データが与えられたときに予測をしたいのですが、
そのためには、データにフィットする曲線を求めればよさそうです。
出力値が入力値の関数で表せるならば、新たな入力に対して予測ができます。
単なる点推定ではなく、不確実性を表現できると予測としては最高です。
これがベイズ線形回帰で実現できます。

概要

本記事では説明することが多いので、流れを説明します。
実際の項目とは異なります。

「モデルの決定」どんなモデルを使うか決めます。
「パラメータの導出」モデルのパラメータを訓練データから求めます。
「モデル選択」ハイパーパラメータ $M$ の値についてざっくり検討します。
「過学習と正則化」過学習と正則化について簡単に触れます。
「最尤推定」確率モデルについて最尤推定します。
「事後分布とMAP推定」事後分布を求め、MAP推定します。事前分布を2パターン設けて説明します。
「予測分布」予測分布を導出します。
「周辺尤度最大化によるモデル選択」周辺尤度最大化により、ハイパーパラメータを推定します。

データ

データは $N$ の組 $(x_n,t_n)$ からなり、
$x_n\in[0,1]$ に対して、目標値 $t_n$ は $y=\sin 2\pi x$ に対応する値に、標準偏差が $0.3$ のガウス分布に従うノイズを加算して作られます。

変数

変数でよく登場するや重要なものをまとめておきます。
$N\in{\mathbb N}:$ データの個数
$D\in{\mathbb N}:$ 入力空間の次元
$M\in{\mathbb N}:$ 特徴空間の次元
${\bf x}_n\in{\mathbb R}^D,\ n\in{\mathbb N}:$ 訓練データの $n$ 番目の入力値
$t_n\in{\mathbb R},\ n\in{\mathbb N}:$ 訓練データの $n$ 番目の目標値
${\bf X}=\{{\bf x}_1,\cdots{\bf x}_N\}$ : 訓練データの入力値全体の集合
${\bf t}=(t_1,\cdots,t_N)^\top\in{\mathbb R}^N$ : 訓練データの目標値からなるベクトル
${\bf x}=(x_1,\cdots,x_D)^\top\in{\mathbb R}^D:$ 入力値
$t\in{\mathbb R}:$ 目標値
${\bf w}=(w_0,\dots,w_{M-1})^\top\in{\mathbb R}^M:$ パラメータ
$\phi_j:{\mathbb R}^D \rightarrow\ {\mathbb R}, \ j\in[0,M-1]\:$ 基底関数( $j$ は基底関数のパラメータ)
${\boldsymbol\phi}: {\mathbb R}^D\ \rightarrow\ {\mathbb R}^M, {\boldsymbol\phi}({\bf x})=(\phi_0({\bf x}),\dots,\phi_{M-1}({\bf x}))^\top:$ 特徴量関数
${\boldsymbol\Phi}=({\boldsymbol\phi}({\bf x}_1),\cdots,{\boldsymbol\phi}({\bf x}_N))^\top\in{\mathbb R}^{N\times M}$ : $n$ 行目が ${\boldsymbol\phi}({\bf x}_n)$ からなる行列(訓練データの入力値からなる行列)
$\lambda\in{\mathbb R}:$ 正則化の係数
$\beta\in{\mathbb R}:$ $t$ の予測値の精度(分散の逆数)
※ ${\bf w}$ の事前分布(ガウス分布)の平均は $({\bf 0},\alpha{\bf I})$ を使う場合と、 $({\bf m}_0,{\bf S}_0)$ を使う場合があります。
$\alpha\in{\mathbb R}:$ ${\bf w}$ の事前分布(多次元ガウス分布)の分散共分散が等方的な行列(スカラー行列)の時のスカラー
${\bf m}_0\in{\mathbb R}^M:$ ${\bf w}$ の事前分布(多次元ガウス分布)の平均
${\bf S}_0\in{\mathbb R}^{M\times M}:$ ${\bf w}$ の事前分布(多次元ガウス分布)の分散共分散行列
${\bf m}_N\in{\mathbb R}^M:$ ${\bf w}$ の事後分布(多次元ガウス分布)の平均
${\bf S}_N\in{\mathbb R}^{M\times M}:$ ${\bf w}$ の事後分布(多次元ガウス分布)の分散共分散行列
$m({\bf x})\in{\mathbb R}:$ ${\bf x}$ が与えられたときの $t$ の予測分布(1次元ガウス分布)の平均
$s^2({\bf x})\in{\mathbb R}:$ ${\bf x}$ が与えられたときの $t$ の予測分布(1次元ガウス分布)の分散

注意

入力値 ${\bf x}$ は ${\bf x}=(x_1,\dots,x_D)\in\mathbb{R}^D$ とします。
入力値はベクトルなのですが、図は細字でスカラーで書かれていたりして、細字と太字が混在しています。
すみませんが適宜読み分けてください。

線形基底関数モデル

最も単純な線形回帰モデルは入力変数の線形結合

$\begin{eqnarray} y({\bf x}, {\bf w})=w_0+w_1x_1+\cdots+w_D x_D\tag{1} \end{eqnarray}$

で表されますが、入力変数 $x_i$ に関して線形関数になっているため表現能力に乏しいです。
例えば、 $D=1$ のとき $y({\bf x},{\bf w})=w_0+w_1x$ となり、直線しか表現できません。
そこで、以下のようにこのモデルを拡張します。

$\begin{eqnarray} y({\bf x}, {\bf w})&=&w_0+w_1\phi_1({\bf x})+\cdots+w_{M-1}{\phi}_{M-1}({\bf x})\\ &=&w_0+\sum_{j=1}^{M-1}w_j\phi_j({\bf x})\tag{2} \end{eqnarray}$

$\phi_i$ を基底関数といいます。
ここで、式 $(2)$ のように拡張することによりパラメータ数は $M$ になったことに注意しましょう。

ダミーの基底関数 $\phi_0({\bf x})=1$ を追加すると、式 $(2)$ は以下のように簡単に書くことができます。

$\begin{eqnarray} y({\bf x}, {\bf w})&=&w_0\phi_0({\bf x})+w_1\phi_1({\bf x})+\cdots+w_{M-1}{\phi}_{M-1}({\bf x})\\ &=&\sum_{j=0}^{M-1}w_j\phi_j({\bf x})\\ &=&{\bf w}^\top{\boldsymbol\phi}({\bf x})\tag{3} \end{eqnarray}$

ただし、 ${\bf w}=(w_0,\dots,w_{M-1})^\top,\ {\boldsymbol\phi}({\bf x})=(\phi_0({\bf x}),\dots,\phi_{M-1}({\bf x}))^\top$ とします。
${\boldsymbol\phi}$ という写像で $\bf x$ を変換しているわけですが、これはパターン認識の実用場面でよくみられる前処理(特徴抽出)に相当します。

基底関数には非線形関数を用います。
例えば、 $D=1$ のとき、 $\phi_j(x)=x^j$ は $x$ について非線形です。
$D=1,\ M=3, \phi_j(x)=x^j$ のとき、式 $(3)$ は $y({\bf x},{\bf w})=w_0+w_1x+w_2x^2+w_3x^3$ と $x$ の $3$ 次関数となっており、式 $(1)$ と比べて表現能力が向上していることが分かります。

本記事では基底関数の選択に依存しない内容なので、基底関数を限定しませんが
基底関数として $\phi_j(x)=x^j$ をイメージするとよいと思います。

曲線フィッティング

$N$ 個のデータ $({\bf x}_1,t_1),\ldots,({\bf x}_N,t_N)$ が与えられていて、 ${\bf x}_n\in\mathbb{R},\ t_n\in\mathbb{R}$ とします。
${\bf x}_n$ のように入力値をベクトルで表記しますが、実際は $D=1$ のスカラーとします。

このデータに当てはまる曲線を求めてみます。

式 $(3)$ より、モデルを $y({\bf x}, {\bf w})={\bf w}^\top{\boldsymbol\phi}({\bf x})$ とします。

目的関数の定義

パラメータは最小二乗法によって求めます。

あるデータ $({\bf x}_n,t_n)$ について着目します。
${\bf x}_n$ の予測値は ${\bf w}^\top{\boldsymbol\phi}({\bf x}_n)$ なので、
$t_n$ との差の $2$ 乗は $({\bf w}^\top{\boldsymbol\phi}({\bf x}_n)-t_n)^2$ となります。
これを全データについて足し合わせたものに $1/2$ を乗じたものを目的関数とします。

$\begin{eqnarray} E({\bf w})=\frac{1}{2}\sum_{n=1}^N\left({\bf w}^\top{\boldsymbol\phi}({\bf x}_n)-t_n\right)^2\tag{4} \end{eqnarray}$

ここで、

$\begin{eqnarray} {\boldsymbol\Phi}=\begin{pmatrix} \boldsymbol\phi({\bf x}_1)^\top \\\ \vdots \\\ \boldsymbol\phi({\bf x}_N)^\top\ \end{pmatrix},\ {\bf t}=(t_1,\cdots,t_N)^\top\tag{5} \end{eqnarray}$

とおくと、目的関数は

$\begin{eqnarray} E({\bf w})=\frac{1}{2}||\boldsymbol\Phi{\bf w}-{\bf t}||^2\tag{6} \end{eqnarray}$

と書けます。

目的関数の最小化

目的関数を最小化するために ${\bf w}$ で偏微分し、 $={\bf 0}$ とし、
これについて解くと ${\bf w}$ が定まり、予測曲線が決定します。
誤差の二乗和からなる目的関数を最小化するので、これを最小二乗法と言います。

以下、導出です。

$\begin{eqnarray} &&\nabla E({\bf w})={\bf 0}\\ &&\Leftrightarrow\frac{\partial}{\partial {\bf w}}E({\bf w})={\bf 0}\\ &&\Leftrightarrow\frac{\partial}{\partial {\bf w}}\frac{1}{2}||\boldsymbol\Phi{\bf w}-{\bf t}||^2={\bf 0}\\ &&\Leftrightarrow\frac{1}{2}\frac{\partial}{\partial {\bf w}}(\boldsymbol\Phi{\bf w}-{\bf t})^\top(\boldsymbol\Phi{\bf w}-{\bf t})={\bf 0}\\ &&\Leftrightarrow\frac{1}{2}\frac{\partial}{\partial {\bf w}}({\bf w}^\top\boldsymbol\Phi^\top-{\bf t}^\top)(\boldsymbol\Phi{\bf w}-{\bf t})={\bf 0}\\ &&\Leftrightarrow\frac{1}{2}\frac{\partial}{\partial {\bf w}}({\bf w}^\top\boldsymbol\Phi^\top\boldsymbol\Phi{\bf w}-{\bf w}^\top\boldsymbol\Phi^\top{\bf t}-{\bf t}^\top\boldsymbol\Phi{\bf w}+{\bf t}^\top{\bf t})={\bf 0}\\ &&\Leftrightarrow\frac{1}{2}\frac{\partial}{\partial {\bf w}}({\bf w}^\top\boldsymbol\Phi^\top\boldsymbol\Phi{\bf w}-2{\bf w}^\top\boldsymbol\Phi^\top{\bf t}+||{\bf t}||^2)={\bf 0}\\ &&\Leftrightarrow\frac{1}{2}(2\boldsymbol\Phi^\top\boldsymbol\Phi{\bf w}-2\boldsymbol\Phi^\top{\bf t})={\bf 0}\\ &&\Leftrightarrow \boldsymbol\Phi^\top\boldsymbol\Phi{\bf w}=\boldsymbol\Phi^\top{\bf t}\\ &&\Leftrightarrow {\bf w}=(\boldsymbol\Phi^\top\boldsymbol\Phi)^{-1}\boldsymbol\Phi^\top{\bf t}\tag{7}\\ \end{eqnarray}$

$\boldsymbol\Phi^\top\boldsymbol\Phi$ は正則であるとします。

目的関数は凸関数である

目的関数が凸関数であること示すには、目的関数のヘッセ行列が半正定値行列であることを示せばよいです。
以下、証明です。

まず、目的関数のヘッセ行列を求めます。

$\begin{eqnarray} \nabla\nabla E({\bf w})&=&\frac{\partial}{\partial {\bf w}}\left(\frac{\partial}{\partial {\bf w}^\top}E({\bf w})\right)\\ &=&\frac{\partial}{\partial {\bf w}}\left(\frac{\partial}{\partial {\bf w}}E({\bf w})\right)^\top\\ &=&\frac{\partial}{\partial {\bf w}}\left(\boldsymbol\Phi^\top\boldsymbol\Phi{\bf w}-\boldsymbol\Phi^\top{\bf t}\right)^\top\\ &=&\frac{\partial}{\partial {\bf w}}({\bf w}^\top\boldsymbol\Phi^\top\boldsymbol\Phi-{\bf t}^\top\boldsymbol\Phi)\\ &=&(\boldsymbol\Phi^\top\boldsymbol\Phi)^\top\\ &=&\boldsymbol\Phi^\top\boldsymbol\Phi\tag{8}\\ \end{eqnarray}$

次に、ヘッセ行列が半正定値行列であることを示します。
${\bf x}\in\mathbb{R}^D$ とします。

$\begin{eqnarray} {\bf x}^\top\boldsymbol\Phi^\top\boldsymbol\Phi{\bf x}&=&(\boldsymbol\Phi{\bf x})^\top(\boldsymbol\Phi{\bf x})\\ &=&||\boldsymbol\Phi{\bf x}||^2\geq0\tag{9}\\ \end{eqnarray}$

目的関数のヘッセ行列は半正定値行列であることが示せました。
よって、目的関数は凸関数です。

モデル比較

与えられたデータに対して、 ${\bf w}$ が計算できることは分かりましたが、パラメータ $M$ を選ぶ問題が残っています。
基底関数を $\phi_j({\bf x})=x^j$ としたときの $M=0,1,3,9$ の場合の予測曲線は以下のようになります。

$M=0,1$ の場合は、データへの当てはまりがあまり良くないことが分かります。
$M=3$ の場合が、もっともよく当てはまっているように見えます。
$M=9$ の場合、訓練データには非常によく当てはまっています。データの数 $N$ が $10$ であり、予測曲線が $9$ 次関数であるため $E({\bf w}^\star)=0$ になっています。しかし、当てはめた曲線は無茶苦茶に発振したようになっています。このようなふるまいを過学習といいます。

私たちの目標は、新たなデータに対して正確な予測を行える汎化性能を達成することです。
訓練データと訓練データとは独立したテストデータを用いて平均二乗平方根誤差 $E_{RMS}=\sqrt{2E({\bf w}^\star)/N}$ を計算した結果が以下です。

$M=3,4,5,6,7,8$ の場合は訓練データ、テストデータの平均二乗平方根誤差が小さいので、良さそうなMであることが分かります。
$M=9$ の場合は、訓練データの平均二乗平方根誤差は $0$ ですが、テストデータの時の値が非常に大きく、過学習していることが見て取れます。

次に、モデルの次数 $M$ は固定し、データ集合のサイズを変えてみたときの振る舞いを以下の図 $(M=9)$ に示します。

モデルの複雑さを固定したとき、データ集合のサイズが大きくなるにつれて過学習の問題は深刻でなくなっていることが分かります。

リッジ回帰(L2正則化)

過学習と正則化

過学習とはモデルがデータに適合しすぎてしまうことを言います。

過学習を防ぐために目的関数に正則化項(罰則項)を加えて最適化します。
このことを正則化と言います。

この記事では加える正則化項がL2ノルムなのでL2正則化といいます。リッジ回帰とも言います。

$\begin{eqnarray} E({\bf w})=\frac{1}{2}||{\boldsymbol\Phi}{\bf w}-{\bf t}||^2+\frac{\lambda}{2}||{\bf w}||^2\tag{10} \end{eqnarray}$

正則化の効果

正則化の効果をグラフで見てみます。
$\lambda=e^{-18}$ としました。
グラフより正則化有の方がいい感じに曲線にフィットしているのが分かります。
また、 ${\bf w}$ の要素の絶対値も正則化有の方が小さいことが分かります。

解析解

目的関数を最小化するような ${\bf w}$ を求めてみます。

$\begin{eqnarray} &&\frac{\partial}{\partial{\bf w}}E({\bf w})={\bf 0}\\ &&\Leftrightarrow\frac{\partial}{\partial{\bf w}}\left(\frac{1}{2}||{\boldsymbol\Phi}{\bf w}-{\bf t}||^2+\frac{\lambda}{2}||{\bf w}||^2\right)={\bf 0}\\ &&\Leftrightarrow\frac{\partial}{\partial{\bf w}}\left(\frac{1}{2}{\bf w}^T\boldsymbol\Phi^T\boldsymbol\Phi{\bf w}-{\bf w}^T\boldsymbol\Phi^T{\bf t}+\frac{1}{2}||{\bf t}||^2+\frac{\lambda}{2}||{\bf w}||^2\right)={\bf 0}\\ &&\Leftrightarrow\boldsymbol\Phi^T\boldsymbol\Phi{\bf w}-\boldsymbol\Phi^T{\bf t}+\lambda{\bf w}={\bf 0}\\ &&\Leftrightarrow(\boldsymbol\Phi^T\boldsymbol\Phi+\lambda{\bf I}){\bf w}=\boldsymbol\Phi^T{\bf t}\\ &&\Leftrightarrow{\bf w}=(\boldsymbol\Phi^T\boldsymbol\Phi+\lambda{\bf I})^{-1}\boldsymbol\Phi^T{\bf t}\tag{11}\\ \end{eqnarray}$

$\boldsymbol\Phi^T\boldsymbol\Phi+\lambda{\bf I}$ は正則であるとします。
ただし、 $w_0$ は正則化項から外すことも多いので注意してください。
$w_0$ を正則化項から外す場合は、式 $(11)$ の ${\bf I}$ の1行1列成分を0にしてください。

最尤推定

目標変数 $t$ が $y({\bf x},{\bf w})$ と加法性のガウスノイズの和で表される場合を考えます。

$\begin{eqnarray} t=y({\bf x},{\bf w})+\epsilon\tag{12}\\ \end{eqnarray}$

$\epsilon$ は期待値が $0$ で精度(分散の逆数)が $\beta$ のガウス分布に従う確率変数です。すなわち、

$\begin{eqnarray} p(t|{\bf x},{\bf w},\beta)=\mathcal{N}(t|y({\bf x},{\bf w}),\beta^{-1})\tag{13}\\ \end{eqnarray}$

と表せます。
式 $(13)$ のイメージ図を記載します。

ここで尤度関数を求めます。

$\begin{eqnarray} p({\bf t}|{\bf X},{\bf w},\beta)=\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\tag{14}\\ \end{eqnarray}$

式 $(14)$ で ${\bf X}=\{{\bf x}_1,\ldots,{\bf x}_N\},\ {\bf t}=\{t_1,\ldots,t_N\}$ としました。
式 $(14)$ に対数を取って、対数尤度関数を求めます。

$\begin{eqnarray} \ln p({\bf t}|{\bf X},{\bf w},\beta)&=&\sum_{n=1}^N\ln\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\\ &=&\frac{N}{2}\ln\beta-\frac{N}{2}\ln(2\pi)-\frac{\beta}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2\tag{15}\\ \end{eqnarray}$

なんと、 $(15)$ の第三項に二乗和誤差関数が現れました。
すなわち、条件付ガウスノイズ分布の下で線形モデルに対する尤度関数の最大化は、二乗和誤差関数を最小化するのと等価であることが分かります。

式 $(15)$ を ${\bf w}$ で偏微分し、 $={\bf 0}$ とおいて解くと(導出については式 $(7)$ をご覧ください。全く同じです。)

$\begin{eqnarray} {\bf w}_{\rm ML}=(\boldsymbol\Phi^\top\boldsymbol\Phi)^{-1}\boldsymbol\Phi^\top{\bf t}\tag{16}\\ \end{eqnarray}$

式 $(16)$ は当然、最小二乗法の解と同じです。
式 $(15)$ を $\beta$ で偏微分し、 $=0$ とおいて解くと $\beta$ が求まります。

$\begin{eqnarray} {\beta_{\rm ML}}^{-1}=\frac{1}{N}\sum_{n=1}^N(t_n-{\bf w}_{\rm ML}^\top{\boldsymbol\phi}({\bf x}_n))^2\tag{17}\\ \end{eqnarray}$

事後分布とMAP推定

wの事前分布の平均が0で共分散行列が等方的な場合

MAP推定とは事後分布を最大するようなパラメータ ${\bf w}$ を点推定することです。
事後分布を求めるために事前分布を定めます。
尤度は1次元ガウス分布の積であり、多次元ガウス分布であるため、共役事前分布として多次元ガウス分布とします。
恣意的ですが、平均は ${\bf 0}$ 、共分散行列は等方的な分布 $\alpha^{-1}{\bf I}$ とします。 $\alpha$ は定数です。

$\begin{eqnarray} p({\bf w}|\alpha)=\mathcal{N}({\bf w}|{\bf 0},\alpha^{-1}{\bf I})\tag{18} \end{eqnarray}$

尤度関数は式 $(14)$ より以下のようになります。 $\beta$ は定数です。

$\begin{eqnarray} p({\bf t}|{\bf X},{\bf w},\beta)=\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\tag{19}\\ \end{eqnarray}$

あとは、ベイズの定理を用いれば、事後分布 $p({\bf w}|{\bf X},{\bf t},\alpha,\beta)$ は求まるのですが、
ベイズの定理というのは、そもそも同時分布を異なる条件付き確率の積にして、変形しただけのものですから、
それにならって同時分布 $p({\bf X},{\bf t},{\bf w}|\alpha,\beta)$ から考えていきたいと思います。
その同時分布を考えるために、グラフィカルモデルという図がとても便利なので、今回のグラフィカルモデルを以下に示します。

まず、同時分布 $p({\bf X},{\bf t},{\bf w}|\alpha,\beta)$ を以下のように変形します。

$\begin{eqnarray} p({\bf X},{\bf t},{\bf w}|\alpha,\beta)&=&p({\bf t},{\bf w}|{\bf X},\alpha,\beta)p({\bf X}|\alpha,\beta)\\ &=&p({\bf w}|{\bf X},{\bf t},\alpha,\beta)p({\bf t}|{\bf X},\alpha,\beta)p({\bf X}|\alpha,\beta)\\ &=&p({\bf w}|{\bf X},{\bf t},\alpha,\beta)p({\bf t}|{\bf X},\beta)p({\bf X})\tag{20} \end{eqnarray}$

次に、同時分布 $p({\bf X},{\bf t},{\bf w}|\alpha,\beta)$ を以下のように式 $(20)$ とは異なる変形します。

$\begin{eqnarray} p({\bf X},{\bf t},{\bf w}|\alpha,\beta)&=&p({\bf t},{\bf w}|{\bf X},\alpha,\beta)p({\bf X}|\alpha,\beta)\\ &=&p({\bf t}|{\bf X},{\bf w},\alpha,\beta)p({\bf w}|{\bf X},\alpha,\beta)p({\bf X}|\alpha,\beta)\\ &=&p({\bf t}|{\bf X},{\bf w},\beta)p({\bf w}|\alpha)p({\bf X})\tag{21} \end{eqnarray}$

式 $(20),(21)$ を等号で結んで、事後分布 $p({\bf w}|{\bf X},{\bf t},\alpha,\beta)$ について解きます。

$\begin{eqnarray} p({\bf w}|{\bf X},{\bf t},\alpha,\beta)&=&\dfrac{p({\bf t}|{\bf X},{\bf w},\beta)p({\bf w}|\alpha)p({\bf X})}{p({\bf t}|{\bf X},\beta)p({\bf X})}\\ &=&\dfrac{p({\bf t}|{\bf X},{\bf w},\beta)p({\bf w}|\alpha)}{p({\bf t}|{\bf X},\beta)}\\ &\propto& p({\bf t}|{\bf X},{\bf w},\beta)p({\bf w}|\alpha)\\ &\propto&\left(\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\right)\mathcal{N}({\bf w}|{\bf 0},\alpha^{-1}{\bf I})\tag{22} \end{eqnarray}$

式 $(22)$ に対数を取ります。

$\begin{eqnarray} \ln p({\bf w}|{\bf X},{\bf t},\alpha,\beta)&=&\left(\sum_{n=1}^N\ln\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\right)+\ln\mathcal{N}({\bf w}|{\bf 0},\alpha^{-1}{\bf I})+{\rm const.}\\ &=&-\frac{N}{2}\ln2\pi+\frac{N}{2}\ln\beta-\frac{\beta}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2-\frac{1}{2}\ln2\pi+\frac{1}{2}\ln\alpha-\frac{\alpha}{2}{\bf w}^\top{\bf w}+{\rm const.}\\ &=&-\frac{\beta}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2-\frac{\alpha}{2}{\bf w}^\top{\bf w}+{\rm const.}\\ &=&-\beta\left(\frac{1}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2+\frac{\alpha}{2\beta}{\bf w}^\top{\bf w}\right)+{\rm const.}\\ &=&-\beta\left(\frac{1}{2}||{\boldsymbol\Phi}{\bf w}-{\bf t}||^2+\frac{1}{2}\cdot\frac{\alpha}{\beta}||{\bf w}||^2\right)+{\rm const.}\tag{23}\\ \end{eqnarray}$

なんと、式 $(23)$ にリッジ回帰(L2正則化)と同じ目的関数が出てきました。

すなわち、この事後分布を ${\bf w}$ について最大化することは、二乗和誤差関数と二次正則化項の和を最小化することと等価であることが分かります。
このような結果になったのは、事前分布の平均を ${\bf 0}$ 、共分散行列を等方的な分布 $\alpha^{-1}{\bf I}$ としたためです。
式 $(23)$ がガウス分布であることは、以下の「wの事前分布に特別な仮定を置かない場合」で示します。

wの事前分布に特別な仮定を置かない場合

今度は、 ${\bf w}$ の事前分布に仮定を設けずに、 ${\bf w}$ の事後分布を求めていきます。
${\bf w}$ の事前分布の平均を ${\bf m}_0$ 、共分散を ${\bf S}_0$ とします。

$\begin{eqnarray} &&p({\bf w}|{\bf m}_0,{\bf S}_0)=\mathcal{N}({\bf w}|{\bf m}_0,{\bf S}_0)\tag{24}\\ \end{eqnarray}$

以下にこの時のグラフィカルモデルを示します。

先程と同じ要領で、事後分布を求めます。

$\begin{eqnarray} p({\bf w}|{\bf X},{\bf t},{\bf m}_0,{\bf S}_0,\beta)&\propto&p({\bf t}|{\bf X},{\bf w},\beta)p({\bf w}|{\bf m}_0,{\bf S}_0)\\ &=&\left(\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\right)\mathcal{N}({\bf w}|{\bf m}_0,{\bf S}_0)\tag{25}\\ \end{eqnarray}$

対数を取って、 $\bf w$ についてまとめます。

$\begin{eqnarray} &&\ln p({\bf w}|{\bf X},{\bf t},{\bf m}_0,{\bf S}_0,\beta)\\ &=&\sum_{n=1}^N\ln\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})+\ln\mathcal{N}({\bf w}|{\bf m}_0,{\bf S}_0)+{\rm const.}\\ &=&-\frac{N}{2}\ln2\pi+\frac{N}{2}\ln\beta-\frac{\beta}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2-\frac{D}{2}\ln2\pi-\frac{D}{2}\ln|{\bf S}_0|-\frac{1}{2}({\bf w}-{\bf m}_0)^\top{\bf S}_0^{-1}({\bf w}-{\bf m}_0)+{\rm const.}\\ &=&-\frac{\beta}{2}\sum_{n=1}^N(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2-\frac{1}{2}({\bf w}-{\bf m}_0)^\top{\bf S}_0^{-1}({\bf w}-{\bf m}_0)+{\rm const.}\\ &=&-\frac{\beta}{2}||{\bf t}-{\boldsymbol\Phi}{\bf w}||^2-\frac{1}{2}({\bf w}-{\bf m}_0)^\top{\bf S}_0^{-1}({\bf w}-{\bf m}_0)+{\rm const.}\\ &=&-\frac{\beta}{2}({\bf w}^\top{\boldsymbol\Phi}^\top{\boldsymbol\Phi}{\bf w}-2{\bf w}^\top{\boldsymbol\Phi}^\top{\bf t})-\frac{1}{2}({\bf w}^\top{\bf S}_0^{-1}{\bf w}-2{\bf w}{\bf S}_0^{-1}{\bf m}_0)+{\rm const.}\\ &=&-\frac{1}{2}({\bf w}^\top({\boldsymbol\Phi}^\top\beta{\boldsymbol\Phi}+{\bf S}_0^{-1}){\bf w}-2{\bf w}^\top(\beta{\boldsymbol\Phi}^\top{\bf t}-2{\bf S}_0^{-1}{\bf m}_0))+{\rm const.}\tag{26}\\ \end{eqnarray}$

$p({\bf w}|{\bf X},{\bf t},{\bf m}_0,{\bf S}_0,\beta)$ は多次元ガウス分布なので、平均を ${\bf m}_N$ 、共分散行列を ${\bf S}_N$ とおき、 $\bf w$ についてまとめます。

$\begin{eqnarray} \ln \mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N)&=&-\frac{1}{2}({\bf w}^\top{\bf S}_N^{-1}{\bf w}-2{\bf w}^\top{\bf S}_N^{-1}{\bf m}_N)+{\rm const.}\tag{27}\\ \end{eqnarray}$

式 $(4),(5)$ を係数比較します。

$\begin{eqnarray} &&{\bf m}_N={\bf S}_N({\bf S}_0^{-1}{\bf m}_0+\beta{\boldsymbol\Phi}^\top{\bf t})\tag{28}\\ &&{\bf S}_N^{-1}={\bf S}_0^{-1}+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}\tag{29}\\ \end{eqnarray}$

MAP推定は事後分布を最大にするようなパラメータ $\bf w$ を点推定することであり、
多次元ガウス分布は平均が最大値をとるので、その値を ${\bf w}_{\rm MAP}$ と書くと、

$\newcommand{\argmax}{\mathop{\rm arg~max}\limits}\begin{eqnarray} {\bf w}_{\rm MAP}&=&\argmax_{\bf w}\ln \mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N)\\\ &=&{\bf m}_N\\ &=&{\bf S}_N({\bf S}_0^{-1}{\bf m}_0+\beta{\boldsymbol\Phi}^\top{\bf t})\tag{30}\\ \end{eqnarray}$

予測分布

この「予測分布」の説明では、条件付き確率の条件 ${\beta},{\bf m}_0,{\bf S}_0$ を省略します。

予測分布とは、新たな入力値 $\bf x$ に対する $t$ の分布のことです。

$\begin{eqnarray} &&p(t|{\bf x},{\bf X},{\bf t})\tag{31} \end{eqnarray}$

この時のグラフィカルモデルを以下に示します。

予測分布の導出について2通り示そうと思います。

導出その1 - ガウス分布の周辺化に関する公式を使う方法

予測分布は $t,{\bf w}$ の条件付き同時分布を $\bf w$ について周辺化することにより求まります。

$\begin{eqnarray} &&p(t|{\bf x},{\bf X},{\bf t})=\int p(t,{\bf w}|{\bf x},{\bf X},{\bf t}){\rm d}{\bf w}\tag{32} \end{eqnarray}$

$p(t,{\bf w}|{\bf x},{\bf X},{\bf t})$ を変形します。

$\begin{eqnarray} p(t,{\bf w}|{\bf x},{\bf X},{\bf t})&=&p(t|{\bf w},{\bf x},{\bf X},{\bf t})p({\bf w}|{\bf x},{\bf X},{\bf t})\\ &=&p(t|{\bf w},{\bf x})p({\bf w}|{\bf X},{\bf t})\tag{33} \end{eqnarray}$

$(33)$ を $(32)$ に代入します。

$\begin{eqnarray} &&p(t|{\bf x},{\bf X},{\bf t})=\int p(t|{\bf w},{\bf x})p({\bf w}|{\bf X},{\bf t}){\rm d}{\bf w}\tag{34} \end{eqnarray}$

ここで、

$\begin{eqnarray} &&p(t|{\bf w},{\bf x})=\mathcal{N}(t|{\bf w}^\top{\boldsymbol\phi}({\bf x}),\beta^{-1})\tag{35}\\ &&p({\bf w}|{\bf X},{\bf t})=\mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N)\tag{36}\\ \end{eqnarray}$

$(35)$ は仮定であり、 $(36)$ は事後分布で式 $(28),(29)$ で求めています。
$(35),(36)$ を $(34)$ に代入します。

$\begin{eqnarray} &&p(t|{\bf x},{\bf X},{\bf t})=\int \mathcal{N}(t|{\bf w}^\top{\boldsymbol\phi}({\bf x}),\beta^{-1})\mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N){\rm d}{\bf w}\tag{37} \end{eqnarray}$

$(37)$ は以下の公式を使えば、積分計算せずに求まります。

$\begin{eqnarray} &&p(\bf x) &=& \mathcal{N}(\mathbf x | \boldsymbol\mu, \mathbf\Lambda^{-1})\tag{38}\\ &&p(\bf y | \bf x) &=& \mathcal{N}(\mathbf y | \mathbf A \mathbf x + \mathbf b, \mathbf{L}^{-1}) \tag{39}\\ \end{eqnarray}$

のとき、

$\begin{eqnarray} p(\mathbf y) = \mathcal{N}(\mathbf y | \mathbf A {\boldsymbol\mu} + \mathbf b , \mathbf{L}^{-1} + \mathbf A \mathbf \Lambda^{-1} \mathbf A^{\top}) \tag{40} \end{eqnarray}$

が成り立ちます。

この公式に ${\bf x}={\bf w},\ {\boldsymbol\mu}={\bf m}_N,\ {\boldsymbol\Lambda}^{-1}={\bf S}_N,\ {\bf y}=t,\ {\bf A}{\bf x}={\boldsymbol\phi}({\bf x})^\top{\bf w},\ {\bf L}^{-1}=\beta^{-1},\ {\bf b}=0$ のように当てはめると
$(37)$ は以下のようになります。

$\begin{eqnarray} &&p(t|{\bf x},{\bf X},{\bf t})=\mathcal{N}(t|m({\bf x}),s^2({\bf x}))\tag{41}\\ \end{eqnarray}$

ただし、

$\begin{eqnarray} &&m({\bf x})={\boldsymbol\phi}({\bf x})^\top{\bf m}_N\tag{42}\\ &&s^2({\bf x})=\beta^{-1}+{\boldsymbol\phi}({\bf x})^\top{\bf S}_N{\boldsymbol\phi}({\bf x})\tag{43}\\ \end{eqnarray}$

です。

以下に、予測分布のグラフを3つ示します。

青色の丸は訓練データを、点線は真の値を、赤色の線は予測分布の平均を、青色の線は予測分布±予測分布の標準偏差を表します。
データがある箇所では、予測分布の分散が小さくなっていることが分かります。

導出その2

同時分布 $p(t,{\bf w}|{\bf x},{\bf X},{\bf t})$ を変形します。

$\begin{eqnarray} p(t,{\bf w}|{\bf x},{\bf X},{\bf t})=p({\bf w}|{\bf x},t,{\bf X},{\bf t})p(t|{\bf x},{\bf X},{\bf t})\tag{44} \end{eqnarray}$

式 $(33),(44)$ より、 $p(t|{\bf x},{\bf X},{\bf t})$ は以下のようになります。

$\begin{eqnarray} p(t|{\bf x},{\bf X},{\bf t})&=&\dfrac{p(t|{\bf w},{\bf x})p({\bf w}|{\bf X},{\bf t})}{p({\bf w}|{\bf x},t,{\bf X},{\bf t})}\\ &\propto&\dfrac{p(t|{\bf w},{\bf x})}{p({\bf w}|{\bf x},t,{\bf X},{\bf t})}\tag{45} \end{eqnarray}$

式 $(45)$ に対数を取ります。

$\begin{eqnarray} \ln p(t|{\bf x},{\bf X},{\bf t})&=&\ln p(t|{\bf x},{\bf w})-\ln p({\bf w}|{\bf x},t,{\bf X},{\bf t})+{\rm const.}\tag{46} \end{eqnarray}$

式 $(46)$ の第 $1$ 項の $p(t|{\bf x},{\bf w})$ は式 $(35)$ と同じです。
式 $(46)$ の第 $2$ 項の $p({\bf w}|{\bf x},t,{\bf X},{\bf t})$ は ${\bf w}$ の事後分布のパラメータを事前分布のパラメータとして
1件からなるデータ $\{({\bf x},t)\}$ で事後分布がさらに更新されると考えればよいので、

$\begin{eqnarray} &&p({\bf w}|{\bf x},t,{\bf X},{\bf t})={\mathcal N}({\bf w}|{\bf m}_*,{\bf S}_*)\tag{47}\\ &&{\bf m}_*={\bf S}_*({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))\tag{48}\\ &&{\bf S}_*^{-1}={\bf S}_N^{-1}+\beta{\boldsymbol\phi}({\bf x}){\boldsymbol\phi}({\bf x})^\top\tag{49}\\ \end{eqnarray}$

となります。(式 $(28),(29)$ 参照。)
式 $(35),(47)$ を式 $(46)$ に代入します。

$\begin{eqnarray} &&\ln p(t|{\bf x},{\bf X},{\bf t})\\ &=&\ln \mathcal{N}(t|{\bf w}^\top{\boldsymbol\phi}({\bf x}),\beta^{-1})-\ln {\mathcal N}({\bf w}|{\bf m}_*,{\bf S}_*)+{\rm const.}\\ &=&-\frac{1}{2}\beta(t-{\bf w}^\top{\boldsymbol\phi}({\bf x}))^2+\frac{1}{2}({\bf w}-{\bf m}_*)^\top{\bf S}^{-1}({\bf w}-{\bf m}_*)+{\rm const.}\\ &=&-\frac{1}{2}\beta(t^2-2t{\bf w}^\top{\boldsymbol\phi}({\bf x})+({\bf w}^\top{\boldsymbol\phi}({\bf x}))^2)+\frac{1}{2}({\bf w}^\top{\bf S}^{-1}{\bf w}-2{\bf w}^\top{\bf S}_*^{-1}{\bf m}_*+{\bf m}_*^\top{\bf S}_*^{-1}{\bf m}_*)+{\rm const.}\\ &=&-\frac{1}{2}\beta(t^2-2t{\bf w}^\top{\boldsymbol\phi}({\bf x}))+\frac{1}{2}\underbrace{(-2{\bf w}^\top{\bf S}_*^{-1}{\bf m}_*+{\bf m}_*^\top{\bf S}_*^{-1}{\bf m}_*)}_{:=A}+{\rm const.}\tag{50} \end{eqnarray}$

$A:=-2{\bf w}^\top{\bf S}_*^{-1}{\bf m}_*+{\bf m}_*^\top{\bf S}_*^{-1}{\bf m}_*$ を計算します。

$\begin{eqnarray} A&=&-2{\bf w}^\top{\bf S}_*^{-1}{\bf m}_*+{\bf m}_*^\top{\bf S}_*^{-1}{\bf m}_*\\ &=&-2{\bf w}^\top{\bf S}_*^{-1}{\bf S}_*({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))+({\bf S}_*({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x})))^\top{\bf S}_*^{-1}{\bf S}_*({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))\\ &=&-2{\bf w}^\top({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))+({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))^\top{\bf S}_*({\bf S}_N^{-1}{\bf m}_N+\beta t{\boldsymbol\phi}({\bf x}))\\ &=&-2{\bf w}^\top{\bf S}_N^{-1}{\bf m}_N-2\beta t{\bf w}^\top{\boldsymbol\phi}({\bf x})+\beta^2 t^2{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\boldsymbol\phi}({\bf x})+2\beta t{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N+({\bf S}_N^{-1}{\bf m}_N)^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N\\ &=&-2\beta t{\bf w}^\top{\boldsymbol\phi}({\bf x})+\beta^2 t^2{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\boldsymbol\phi}({\bf x})+2\beta t{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N+{\rm const.}\tag{51} \end{eqnarray}$

式 $(51)$ を式 $(50)$ に代入します。

$\begin{eqnarray} \ln p(t|{\bf x},{\bf X},{\bf t})&=&-\frac{1}{2}\beta(t^2-2t{\bf w}^\top{\boldsymbol\phi}({\bf x}))+\frac{1}{2}(-2\beta t{\bf w}^\top{\boldsymbol\phi}({\bf x})+\beta^2 t^2{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\boldsymbol\phi}({\bf x})+2\beta t{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N)+{\rm const.}\\ &=&-\frac{1}{2}(\beta t^2-\beta^2 t^2{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\boldsymbol\phi}({\bf x})-2\beta t{\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N)+{\rm const.}\\ &=&-\frac{1}{2}((\underbrace{\beta -\beta^2 {\boldsymbol\phi}({\bf x})^\top {\bf S}_* {{\boldsymbol\phi}({\bf x})}}_{:=B}) t^2\underbrace{-2\beta {\boldsymbol\phi}({\bf x})^\top{\bf S}_*{\bf S}_N^{-1}{\bf m}_N}_{:=C} t)+{\rm const.}\tag{52}\\ \end{eqnarray}$

${\bf S}_{*}$ を計算するため、以下のウッドベリーの公式を用います。

${\bf A}\in{\mathbb R}^{n\times n},\ {\bf B}\in{\mathbb R}^{n\times k},\ {\bf C}\in{\mathbb R}^{k\times n},\ {\bf D}\in{\mathbb R}^{k\times k}$ で ${\bf A}^{-1},\ {\bf D}^{-1}$ が存在するとき

$\begin{eqnarray} ({\bf A}+{\bf BDC})^{-1}={\bf A}^{-1}-{\bf A}^{-1}{\bf B}({\bf D}^{-1}+{\bf C}{\bf A}^{-1}{\bf B})^{-1}{\bf C}{\bf A}^{-1}\tag{53} \end{eqnarray}$

が成り立つ。
式 $(49)$ より、 ${\bf S}_*=({\bf S}_N^{-1}+\beta{\boldsymbol\phi}({\bf x}){\boldsymbol\phi}({\bf x})^\top)^{-1}$ なので、
${\bf A}={\bf S}_N^{-1},\ {\bf B}={\boldsymbol\phi}({\bf x}),\ {\bf C}={\boldsymbol\phi}({\bf x})^\top,\ {\bf D}=\beta$ とウッドベリーの公式に対応づけると

$\begin{eqnarray} {\bf S}_{*} = {\bf S}_N - {\bf S}_N {\boldsymbol\phi}({\bf x})( {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) )^{-1} {\boldsymbol\phi} ({\bf x}) ^\top {\bf S}_N\tag{54} \end{eqnarray}$

となります。

式 $(52)$ の $B:=\beta -\beta^2 {\boldsymbol\phi}({\bf x})^\top {\bf S}_* {{\boldsymbol\phi}({\bf x})}$ を計算します。

$\begin{eqnarray} B&=&\beta -\beta^2 {\boldsymbol\phi}({\bf x})^\top {\bf S}_* {{\boldsymbol\phi}({\bf x})}\\ &=&\beta -\beta^2 {\boldsymbol\phi}({\bf x})^\top ({\bf S}_N - {\bf S}_N {\boldsymbol\phi}({\bf x})( {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) )^{-1} {\boldsymbol\phi} ({\bf x}) ^\top {\bf S}_N) {{\boldsymbol\phi}({\bf x})}\\ &=&\beta -\frac{\beta^2}{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } {\boldsymbol\phi}({\bf x})^\top (( {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) ){\bf S}_N - {\bf S}_N {\boldsymbol\phi}({\bf x}) {\boldsymbol\phi} ({\bf x}) ^\top {\bf S}_N) {{\boldsymbol\phi}({\bf x})}\\ &=&\beta -\frac{\beta^2}{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } {\boldsymbol\phi}({\bf x})^\top ( {\beta}^{-1}{\bf S}_N+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) {\bf S}_N - {\bf S}_N {\boldsymbol\phi}({\bf x}) {\boldsymbol\phi} ({\bf x}) ^\top {\bf S}_N) {{\boldsymbol\phi}({\bf x})}\\ &=&\beta -\frac{\beta^2}{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } ( {\beta}^{-1} {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) + ({\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) )^2 - ({{\boldsymbol\phi} ({\bf x})} ^\top {\bf S}_N {{\boldsymbol\phi}({\bf x})})^2 ) \\ &=&\beta -\frac{\beta^2}{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } {\beta}^{-1} {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) \\ &=&\beta -\frac{\beta {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } \\ &=&\frac{1 + \beta {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) - \beta {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } \\ &=&\frac{ 1 }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } \tag{55}\\ \end{eqnarray}$

式 $(52)$ の $C:= -2 \beta {\boldsymbol\phi}({\bf x})^\top {\bf S}_* {\bf S}_N^{-1} {\bf m}_N$ を計算します。

$\begin{eqnarray} C&=& -2 \beta {\boldsymbol\phi}({\bf x})^\top {\bf S}_* {\bf S}_N^{-1} {\bf m}_N\\ &=& -2 \beta {\boldsymbol\phi}({\bf x})^\top({\bf S}_N - {\bf S}_N {\boldsymbol\phi}({\bf x})( {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) )^{-1} {\boldsymbol\phi} ({\bf x}) ^\top {\bf S}_N) {\bf S}_N^{-1} {\bf m}_N\\ &=& -2 \beta {\boldsymbol\phi}({\bf x})^\top({\bf I} - {\bf S}_N {\boldsymbol\phi}({\bf x})( {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) )^{-1} {\boldsymbol\phi} ({\bf x}) ^\top ) {\bf m}_N\\ &=& \frac{ -2 \beta }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } { {\boldsymbol\phi}({\bf x}) }^\top( ({\beta}^{-1}+ { {\boldsymbol\phi}({\bf x}) } ^\top {\bf S}_N { {\boldsymbol\phi} ({\bf x}) } ) {\bf I} - {\bf S}_N { {\boldsymbol\phi}({\bf x}) } { {\boldsymbol\phi} ({\bf x}) } ^\top ) {\bf m}_N\\ &=& \frac{ -2 \beta }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } { {\boldsymbol\phi}({\bf x}) }^\top( {\beta}^{-1} {\bf I} + { {\boldsymbol\phi}({\bf x}) } ^\top {\bf S}_N { {\boldsymbol\phi} ({\bf x}) } - {\bf S}_N { {\boldsymbol\phi}({\bf x}) } { {\boldsymbol\phi} ({\bf x}) } ^\top ) {\bf m}_N\\ &=& \frac{ -2 \beta }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } ( {\beta}^{-1} { {\boldsymbol\phi}({\bf x}) }^\top {\bf m}_N + { {\boldsymbol\phi}({\bf x}) }^\top { {\boldsymbol\phi}({\bf x}) } ^\top {\bf S}_N { {\boldsymbol\phi} ({\bf x}) } {\bf m}_N - { {\boldsymbol\phi}({\bf x}) }^\top {\bf S}_N { {\boldsymbol\phi}({\bf x}) } { {\boldsymbol\phi} ({\bf x}) } ^\top {\bf m}_N ) \\ &=& \frac{ -2 }{ {\beta}^{-1}+ {\boldsymbol\phi}({\bf x})^\top {\bf S}_N {\boldsymbol\phi} ({\bf x}) } { {\boldsymbol\phi}({\bf x}) } ^\top {\bf m}_N \tag{56} \end{eqnarray}$

式 $(55),(56)$ を式 $(52)$ に代入します。

$\begin{eqnarray} \ln p(t|{\bf x},{\bf X},{\bf t})&=&-\frac{1}{2}\left(({\beta}^{-1}+ {\boldsymbol\phi}({\bf x}) ^\top {\bf S}_N { {\boldsymbol\phi} ({\bf x}) })^{-1} t^2 -2 ({\beta}^{-1}+ { {\boldsymbol\phi}({\bf x}) } ^\top {\bf S}_N { {\boldsymbol\phi} ({\bf x}) } )^{-1} { {\boldsymbol\phi}({\bf x}) } ^\ top {\bf m}_N t\right)+{\rm const.}\tag{57} \end{eqnarray}$

式 $(57)$ より、予測分布はガウス分布なので、平均を $m({\bf x})$ 、分散を $s^2({\bf x})$ とおきます。

$\begin{eqnarray} \ln {\mathcal N} (t | m({\bf x}), {s^2}({\bf x}) ) =-\frac{1}{2}\left(\frac{1}{s^2({\bf x})}t^2-2\frac{1}{s^2({\bf x})} m({\bf x}) t\right)+{\rm const.}\tag{58} \end{eqnarray}$

$(57)$ と $(58)$ を係数比較します。

$\begin{eqnarray} &&m({\bf x})={\boldsymbol\phi}({\bf x})^\top{\bf m}_N\tag{59}\\ &&s^2({\bf x})=\beta^{-1}+{\boldsymbol\phi}({\bf x})^\top{\bf S}_N{\boldsymbol\phi}({\bf x})\tag{60}\\ \end{eqnarray}$

解法1と一致しました。

ベイズモデル比較

ベイズの立場からモデル比較の問題を考えます。

周辺尤度

$L$ 個のモデル $\{\mathcal{M}_i\}(i=1,\ldots,L)$ を比較する場合を考えます。
観測されたデータを $\mathcal{D}$ とします。
モデルの不確かさは、事前確率分布 $p({\mathcal{M}_i})$ で表します。
この時、モデルの事後分布は

$\begin{eqnarray} p(\mathcal{M}_i|\mathcal{D})\propto p(\mathcal{M}_i)p(\mathcal{D}|\mathcal{M}_i)\tag{61} \end{eqnarray}$

となります。

簡単のため、すべてのモデルの事前分布 $p(\mathcal{M}_i)$ が等しい場合を考えます。
ベイズモデル比較では、モデルエビデンス $p(\mathcal{D}|\mathcal{M}_i)$ と呼ばれる項が重要な働きをします。
モデルエビデンスは周辺尤度とも呼ばれます。
2つのモデルに対するエビデンスの比 $p(\mathcal{D}|\mathcal{M}_i)/p(\mathcal{D}|\mathcal{M}_j)$ はベイズ因子と呼ばれます。

モデルの事後分布が分かれば、予測分布は以下のようになります。

$\begin{eqnarray} p(t|{\bf x},\mathcal{D})=\sum_{i=1}^Lp(t|{\bf x},\mathcal{M}_i,\mathcal{D})p(\mathcal{M}_i|\mathcal{D})\tag{62} \end{eqnarray}$

モデル平均の単純な近似は、一番尤もらしいモデルを1つ選ぶ方法です。
これをモデル選択と言います。

パラメータ $\bf w$ を持つモデルに対しては、周辺尤度は以下で与えられます。

$\begin{eqnarray} p(\mathcal{D}|\mathcal{M}_i)=\int p(\mathcal{D}|{\bf w},\mathcal{M}_i)p({\bf w}|\mathcal{M}_i){\rm d}{\bf w}\tag{63} \end{eqnarray}$

また、周辺尤度は事後分布を計算するときの分母に現れる正規化定数でもあります。

$\begin{eqnarray} p({\bf w}|\mathcal{D},\mathcal{M}_i)=\frac{p(\mathcal{D}|{\bf w},\mathcal{M}_i)p({\bf w}|\mathcal{M}_i)}{p(\mathcal{D}|\mathcal{M}_i)}\tag{64} \end{eqnarray}$

周辺尤度が大きい方が正しいモデルである理由

ベイズモデルの比較の枠組みでは、真の分布が含まれていると暗に仮定しています。
2つのモデル $\mathcal{M}_1,\mathcal{M}_2$ があり、 $\mathcal{M}_1$ が正しいモデルだと仮定します。
ベイズ因子をデータの集合の分布に対して、平均します。

$\begin{eqnarray} \int p(\mathcal{D}|\mathcal{M}_1)\ln\frac{p(\mathcal{D}|\mathcal{M}_1)}{p(\mathcal{D}|\mathcal{M}_2)}{\rm d}\mathcal{D}\geq 0\tag{65} \end{eqnarray}$

$(66)$ はKLダイバージェンスなので、0以上です。

つまり、正しいモデルに対して、平均的に

$\begin{eqnarray} \left\langle\frac{p(\mathcal{D}|\mathcal{M}_1)}{p(\mathcal{D}|\mathcal{M}_2)}\right\rangle_{p(\mathcal{D}|\mathcal{M}_1)}\geq 1\tag{66}\\ \end{eqnarray}$

なので、 $p(\mathcal{D}|\mathcal{M}_2)$ が大きいほど正しいモデルに近いことが分かります。

周辺尤度によるモデル比較 - M

一般の周辺尤度

一般に、観測データ $\mathcal{D}$ 、パラメータ $\theta$ のとき、 $\theta$ の事後分布は以下の式で表されます。

$\begin{eqnarray} p(\theta|{\mathcal D})=\frac{p({\mathcal D}|\theta)p(\theta)}{p(\mathcal D)}\tag{67} \end{eqnarray}$

$(67)$ の $p(\mathcal D)$ を周辺尤度といい、モデルからデータ $\mathcal D$ が出現する尤もらしさを表しています。
$(67)$ を $p(\mathcal D)$ について解きます。

$\begin{eqnarray} p(\mathcal D)=\frac{p({\mathcal D}|\theta)p(\theta)}{p(\theta|{\mathcal D})}\tag{68} \end{eqnarray}$

このモデルの周辺尤度

今回のモデルでは、 $\mathcal{D}=\{{\bf X},{\bf t}\},\ \theta=\{{\bf w}\}$ となるので、周辺尤度は

$\begin{eqnarray} p(\mathcal D)&=&p({\bf X},{\bf t})\\ &=&p({\bf t}|{\bf X})p({\bf X})\\ &\propto&p({\bf t}|{\bf X})\tag{69} \end{eqnarray}$

となります。
式 $(69)$ で $p({\bf X})$ を無視しているのは、 ${\bf X}$ は今回のモデルでは常に与えられているものだからで、
興味があるのは、周辺尤度の値ではないからです。(ベイズ推論による機械学習入門 p109 より引用)
そのため、周辺尤度として $p({\bf t}|{\bf X})$ を用います。

周辺尤度の計算

ベイズの定理より、周辺尤度は以下のようになります。

$\begin{eqnarray} p({\bf t}|{\bf X})=\frac{p({\bf t}|{\bf X},{\bf w})p({\bf w})}{p({\bf w}|{\bf X},{\bf t})}\tag{70} \end{eqnarray}$

対数を取ります。

$\begin{eqnarray} \ln p({\bf t}|{\bf X})=\ln p({\bf t}|{\bf X},{\bf w})+\ln p({\bf w})-\ln p({\bf w}|{\bf X},{\bf t})\tag{71} \end{eqnarray}$

$p({\bf t}|{\bf X},{\bf w})$ は尤度を、 $p({\bf w})$ は事前分布を、 $p({\bf w}|{\bf X},{\bf t})$ は事後分布を表します。

$\begin{eqnarray} &&p({\bf t}|{\bf X},{\bf w})=\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\tag{72}\\ &&p({\bf w})=\mathcal{N}({\bf w}|{\bf m}_0,{\bf S}_0)\tag{73}\\ &&p({\bf w}|{\bf X},{\bf t})=\mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N)\tag{74}\\ \end{eqnarray}$

$\ln p({\bf t}|{\bf X},{\bf w})$ を展開します。

$\begin{eqnarray} \ln p({\bf t}|{\bf X},{\bf w})&=&\ln\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\\ &=&\sum_{n=1}^N\ln\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\\ &=&-\frac{1}{2}\sum_{n=1}^N\left(\beta(t_n-{\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2+\ln\beta^{-1}+\ln2\pi\right)\\ &=&-\frac{1}{2}\sum_{n=1}^N\left(\beta t_n^2-2\beta t_n{\bf w}^\top{\boldsymbol\phi}({\bf x}_n)+\beta({\bf w}^\top{\boldsymbol\phi}({\bf x}_n))^2-\ln\beta+\ln2\pi\right)\tag{75}\\ \end{eqnarray}$

$\ln p({\bf w})$ を展開します。

$\begin{eqnarray} \ln p({\bf w})&=&\ln\mathcal{N}({\bf w}|{\bf m}_0,{\bf S}_0)\\ &=&-\frac{1}{2}\left(({\bf w}-{\bf m}_0)^\top{\bf S}_0^{-1}({\bf w}-{\bf m}_0)+\ln|{\bf S}_0|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top{\bf S}_0^{-1}{\bf w}-2{\bf w}^\top{\bf S}_0^{-1}{\bf m}_0+{\bf m}_0^\top{\bf S}_0^{-1}{\bf m}_0+\ln|{\bf S}_0|+D\ln2\pi\right)\tag{76}\\ \end{eqnarray}$

$\ln p({\bf w}|{\bf X},{\bf t})$ を展開します。

$\begin{eqnarray} \ln p({\bf w}|{\bf X},{\bf t})&=&\ln\mathcal{N}({\bf w}|{\bf m}_N,{\bf S}_N)\\ &=&-\frac{1}{2}\left(({\bf w}-{\bf m}_N)^\top{\bf S}_0^{-1}({\bf w}-{\bf m}_N)+\ln|{\bf S}_N|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top{\bf S}_N^{-1}{\bf w}-2{\bf w}^\top{\bf S}_N^{-1}{\bf m}_N+{\bf m}^\top{\bf S}_N^{-1}{\bf m}_N+\ln|{\bf S}_N|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top({\bf S}_0^{-1}+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}){\bf w}-2{\bf w}^\top({\bf S}_0^{-1}{\bf m}_0+\beta{\boldsymbol\Phi{\bf t}})+{\bf m}^\top{\bf S}_N^{-1}{\bf m}_N+\ln|{\bf S}_N|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top\left({\bf S}_0^{-1}+\beta\sum_{n=1}^N{\boldsymbol\phi}({\bf x}_n){\boldsymbol\phi}({\bf x}_n)^\top\right){\bf w}-2{\bf w}^\top\left({\bf S}_0^{-1}{\bf m}_0+\beta\sum_{n=1}^Nt_n{\boldsymbol\phi}({\bf x}_n)\right)+{\bf m}^\top{\bf S}_N^{-1}{\bf m}_N+\ln|{\bf S}_N|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top{\bf S}_0^{-1}{\bf w}+\beta\sum_{n=1}^N{\bf w}^\top{\boldsymbol\phi}({\bf x}_n){\boldsymbol\phi}({\bf x}_n)^\top{\bf w}-2{\bf w}^\top{\bf S}_0^{-1}{\bf m}_0-2\beta\sum_{n=1}^Nt_n{\bf w}^\top{\boldsymbol\phi}({\bf x}_n)+{\bf m}^\top{\bf S}_N^{-1}{\bf m}_N+\ln|{\bf S}_N|+D\ln2\pi\right)\\ &=&-\frac{1}{2}\left({\bf w}^\top{\bf S}_0^{-1}{\bf w}+\beta\sum_{n=1}^N\left({\bf w}^\top{\boldsymbol\phi}({\bf x}_n)\right)^2-2{\bf w}^\top{\bf S}_0^{-1}{\bf m}_0-2\beta\sum_{n=1}^Nt_n{\bf w}^\top{\boldsymbol\phi}({\bf x}_n)+{\bf m}^\top{\bf S}_N^{-1}{\bf m}_N+\ln|{\bf S}_N|+D\ln2\pi\right)\tag{77}\\ \end{eqnarray}$

$(75),(76),(77)$ を $(71)$ に代入します。

$\begin{eqnarray} \ln p({\bf t}|{\bf X})=-\frac{1}{2}\left(\sum_{n=1}^N(\beta t_n^2-\ln\beta+\ln2\pi)+{\bf m}_0^\top{\bf S}_0^{-1}{\bf m}_0+\ln|{\bf S}_0|-{\bf m}_N^\top{\bf S}_N^{-1}{\bf m}_N-\ln|{\bf S}_N|\right)\tag{78} \end{eqnarray}$

これで対数周辺尤度が計算できます。

このモデルのパラメータ $M$ を変えて、対数周辺尤度を計算した結果が以下です。

$M=3$ の時が最大になっています。

周辺尤度によるモデル比較 - α、β

周辺尤度

事前分布を平均は ${\bf 0}$ 、共分散行列は等方的な分布 $\alpha^{-1}{\bf I}_M$ として、 $\alpha,\beta$ を推定します。

事前分布は

$\begin{eqnarray} &&p({\bf w}|\alpha)=\mathcal{N}({\bf w}|{\bf 0},\alpha^{-1}{\bf I}_M)\tag{79}\\ \end{eqnarray}$

です。

尤度は

$\begin{eqnarray} p({\bf t}|{\bf w},\beta)&=&\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\boldsymbol\phi}({\bf x}_n),\beta^{-1})\\ &=&\prod_{n=1}^N\left(\frac{\beta}{2\pi}\right)^{\frac{1}{2}}\exp\left(-\frac{\beta}{2}\left(t_n-{\boldsymbol\phi}({\bf x}_n)^\top{\bf w}\right)^2\right)\\ &=&\frac{1}{(2\pi\beta^{-1})^{\frac{N}{2}}}\exp\left(-\frac{\beta}{2}\sum_{n=1}^N \left(t_n-{\boldsymbol\phi}({\bf x}_n)^\top{\bf w}\right)^2\right)\\ &=&\frac{1}{(2\pi)^{\frac{N}{2}}|\beta^{-1}{\bf I}_N|^{\frac{1}{2}}}\exp\left(-\frac{\beta}{2}({\bf t}-{\boldsymbol\Phi}{\bf w})^\top({\bf t}-{\boldsymbol\Phi}{\bf w})\right)\\ &=&\mathcal{N}({\bf t}|{\boldsymbol\Phi}{\bf w},\beta^{-1}{\bf I}_N)\tag{80}\\ \end{eqnarray}$

です。(多次元正規分布で書き換えました)

周辺尤度を計算します。(先ほど求めた方法とは異なる方法で求めます)

$\begin{eqnarray} p({\bf t}|\alpha,\beta)&=&\int p({\bf t},{\bf w}|\alpha,\beta){\rm d}{\bf w}\\ &=&\int p({\bf t}|{\bf w},\beta)p({\bf w}|\alpha){\rm d}{\bf w}\\ &=&\int \mathcal{N}({\bf t}|{\boldsymbol\Phi}{\bf w},\beta^{-1}{\bf I}_N)\mathcal{N}({\bf w}|{\bf 0},\alpha^{-1}{\bf I}_M){\rm d}{\bf w}\tag{81}\\ \end{eqnarray}$

$(81)$ は以下の公式を使えば、積分計算せずに求まります。

$\begin{eqnarray} &&p(\bf x) &=& \mathcal{N}(\mathbf x | \boldsymbol\mu, \mathbf\Lambda^{-1})\tag{82}\\ &&p(\bf y | \bf x) &=& \mathcal{N}(\mathbf y | \mathbf A \mathbf x + \mathbf b, \mathbf{L}^{-1}) \tag{83}\\ \end{eqnarray}$

のとき、

$\begin{eqnarray} p(\mathbf y) = \mathcal{N}(\mathbf y | \mathbf A {\boldsymbol\mu} + \mathbf b , \mathbf{L}^{-1} + \mathbf A \mathbf \Lambda^{-1} \mathbf A^{\top}) \tag{84} \end{eqnarray}$

が成り立ちます。

この公式に ${\bf x}={\bf w},{\boldsymbol\mu}={\bf 0},{\boldsymbol\Lambda}^{-1}=\alpha^{-1}{\bf I}_M,{\bf y}={\bf t},{\bf A}={\boldsymbol\Phi},{\bf L}^{-1}=\beta^{-1}{\bf I}_N,{\bf b}={\bf 0}$ のように当てはめると

$\begin{eqnarray} &&p({\bf t}|\alpha,\beta)=\mathcal{N}({\bf t}|{\bf 0},\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top)\tag{85}\\ \end{eqnarray}$

となります。
対数を取ります。

$\begin{eqnarray} \ln p({\bf t}|\alpha,\beta)&=&\ln \mathcal{N}({\bf t}|{\bf 0},\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top)\\ &=&-\frac{N}{2}\ln 2\pi-\frac{1}{2}\ln|\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top|-\frac{1}{2}{\bf t}^\top(\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top)^{-1}{\bf t}\tag{86}\\ \end{eqnarray}$

ここで公式 $|a{\bf A}|=a^N|{\bf A}|,|{\bf I}_N+{\bf AB}^\top|=|{\bf I}_M+{\bf A}^\top{\bf B}|$ より、 $(86)$ の第2項は

$\begin{eqnarray} |\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top|&=&\beta^{-N}|{\bf I}_N+\beta\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top|\\ &=&\beta^{-N}|{\bf I}_M+\beta\alpha^{-1}{\boldsymbol\Phi}^\top{\boldsymbol\Phi}|\\ &=&\beta^{-N}\alpha^{-M}|\alpha{\bf I}_M+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}|\\ &=&\beta^{-N}\alpha^{-M}|{\bf A}|\tag{87} \end{eqnarray}$

$(87)$ で ${\bf A}=\alpha{\bf I}_M+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}$ とおきました。

ウッドベリーの公式

$\begin{eqnarray} ({\bf A}+{\bf B}{\bf D}^{-1}{\bf C})^{-1}={\bf A}^{-1}-{\bf A}^{-1}{\bf B}({\bf D}+{\bf C}{\bf A}^{-1}{\bf B})^{-1}{\bf C}{\bf A}^{-1}\tag{88} \end{eqnarray}$

を ${\bf A}=\beta^{-1}{\bf I}_N,{\bf B}=\alpha^{-1}{\boldsymbol\Phi},{\bf C}={\boldsymbol\Phi}^\top,{\bf D}^{-1}={\bf I}_M$ として、 $(86)$ の第3項に適用します。

$\begin{eqnarray} -\frac{1}{2}{\bf t}^\top(\beta^{-1}{\bf I}_N+\alpha^{-1}{\boldsymbol\Phi}{\boldsymbol\Phi}^\top)^{-1}{\bf t}&=&-\frac{1}{2}{\bf t}^\top(\beta{\bf I}_N-\beta\alpha^{-1}{\boldsymbol\Phi}({\bf I}_M+{\boldsymbol\Phi}^\top\beta\alpha^{-1}{\boldsymbol\Phi})^{-1}{\boldsymbol\Phi}^\top\beta){\bf t}\\ &=&-\frac{1}{2}{\bf t}^\top(\beta{\bf I}_N-\beta^2\alpha^{-1}{\boldsymbol\Phi}(\alpha^{-1}(\alpha{\bf I}_M+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}))^{-1}{\boldsymbol\Phi}^\top){\bf t}\\ &=&-\frac{1}{2}{\bf t}^\top(\beta{\bf I}_N-\beta^2\alpha^{-1}{\boldsymbol\Phi}(\alpha(\alpha{\bf I}_M+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi})^{-1}){\boldsymbol\Phi}^\top){\bf t}\\ &=&-\frac{1}{2}{\bf t}^\top(\beta{\bf I}_N-\beta^2{\boldsymbol\Phi}(\alpha{\bf I}_M+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi})^{-1}{\boldsymbol\Phi}^\top){\bf t}\\ &=&-\frac{1}{2}{\bf t}^\top(\beta{\bf I}_N-\beta^2{\boldsymbol\Phi}{\bf A}^{-1}{\boldsymbol\Phi}^\top){\bf t}\\ &=&-\frac{1}{2}(\beta{\bf t}^\top{\bf t}-{\beta^2}{\bf t}^\top{\boldsymbol\Phi}{\bf A}^{-1}{\boldsymbol\Phi}^\top{\bf t})\\ &=&-\frac{1}{2}(\beta{\bf t}^\top{\bf t}-{\beta^2}{\bf t}^\top{\boldsymbol\Phi}{\bf A}^{-1}{\bf A}{\bf A}^{-1}{\boldsymbol\Phi}^\top{\bf t})\\ &=&-\frac{1}{2}(\beta{\bf t}^\top{\bf t}-({\beta}{\bf t}^\top{\boldsymbol\Phi}{\bf A}^{-1}){\bf A}(\beta{\bf A}^{-1}{\boldsymbol\Phi}^\top{\bf t}))\\ &=&-\frac{1}{2}(\beta{\bf t}^\top{\bf t}-(\beta{\bf A}^{-1}{\boldsymbol\Phi}^\top{\bf t})^\top{\bf A}(\beta{\bf A}^{-1}{\boldsymbol\Phi}^\top{\bf t}))\\ &=&-\frac{1}{2}(\beta{\bf t}^\top{\bf t}-{\bf m}_N^\top{\bf A}{\bf m}_N)\\ &=&-\frac{1}{2} \left ( \beta \mathbf t^\top \mathbf t-2 \mathbf m^\top_N \mathbf A \mathbf m_N + \mathbf m^\top_N \mathbf A \mathbf m_N \right ) \\ &=& - \frac { 1 } { 2 } \left ( \beta \mathbf { t } ^\top \mathbf { t } - 2 \mathbf { m } ^ \top_ { N } \mathbf { A } \left ( \beta \mathbf { A } ^ { - 1 } \mathbf { \Phi } ^\top \mathbf { t } \right ) + \mathbf { m } ^\top_ { N } \left ( \alpha \mathbf { I } _ { M } + \beta \mathbf { \Phi } ^\top \mathbf { \Phi } \right ) \mathbf { m } _ { N } \right ) \\ &=&-\frac{1}{2} \left ( \beta \mathbf t^\top \mathbf t-2 \mathbf m^\top_N \boldsymbol \Phi ^\top \mathbf t \beta + \beta \mathbf m^\top_N \boldsymbol \Phi ^\top\boldsymbol \Phi\mathbf m_N+\alpha \mathbf m^\top_N \mathbf m_N \right ) \\ &=&-\frac{1}{2} \left ( \beta \left ( \mathbf t-\boldsymbol \Phi \mathbf m_N \right ) ^\top \left ( \mathbf t-\boldsymbol \Phi \mathbf m_N \right ) +\alpha \mathbf m_N^\top \mathbf m_N \right ) \\ &=& -\frac{\beta}{2} \left | \left | \mathbf t-\boldsymbol \Phi \mathbf m_N \right | \right | ^ 2-\frac{\alpha}{2}\mathbf m_N^\top \mathbf m_N\tag{89} \end{eqnarray}$

$(87),(89)$ を $(86)$ に代入します。

$\begin{eqnarray} \ln p({\bf t}|\alpha,\beta)&=&-\frac{N}{2}\ln 2\pi-\frac{1}{2}\ln\beta^{-N}\alpha^{-M}|{\bf A}|-\frac{\beta}{2} \left | \left | \mathbf t-\boldsymbol \Phi \mathbf m_N \right | \right | ^ 2-\frac{\alpha}{2}\mathbf m_N^\top \mathbf m_N\\ &=&\frac{M}{2}\ln\alpha+\frac{N}{2}\ln\beta-E({\bf m}_N)-\frac{1}{2}\ln|{\bf A}|-\frac{N}{2}\ln(2\pi)\tag{90}\\ \end{eqnarray}$

$(90)$ で $E({\bf m}_N)$ は以下のようにおきました。

$\begin{eqnarray} E({\bf m}_N)=\frac{\beta}{2} \left | \left | \mathbf t-\boldsymbol \Phi \mathbf m_N \right | \right | ^ 2+\frac{\alpha}{2}\mathbf m_N^\top \mathbf m_N\tag{91} \end{eqnarray}$

周辺尤度の最大化

対数周辺尤度を最大化する $\alpha,\beta$ を求めるのですが、解析的に求めるのが難しい為、
$\alpha,\beta$ の更新式を導出し、繰り返し法で解を求めます。

行列 $\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}$ を行列 $\bf P$ で対角化します。

$\begin{eqnarray} {\bf P}^{-1}(\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}){\bf P}={\rm diag}(\lambda_1,\ldots,\lambda_M)\tag{92} \end{eqnarray}$

すると、行列 ${\bf A}$ にも対角化ができます。

$\begin{eqnarray} {\bf P}^{-1}{\bf A}{\bf P}={\rm diag}(\alpha+\lambda_1,\ldots,\alpha+\lambda_M)\tag{93} \end{eqnarray}$

式 $(93)$ の導出については、後で説明します。
よって、

$\begin{eqnarray} |{\bf A}|=\prod_{i=1}^M(\lambda_i+\alpha)\tag{94} \end{eqnarray}$

となります。

$\ln|{\bf A}|$ を $\alpha$ で微分します。

$\begin{eqnarray} \frac{d}{d\alpha}\ln|{\bf A}|&=&\frac{d}{d\alpha}\ln\prod_{i=1}^M(\lambda_i+\alpha)\\ &=&\frac{d}{d\alpha}\sum_{i=1}^M\ln(\lambda_i+\alpha)\\ &=&\sum_{i=1}^M\frac{1}{\lambda_i+\alpha}\tag{95}\\ \end{eqnarray}$

対数周辺尤度 $(90)$ を $\alpha$ で微分して $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial}{\partial\alpha}\ln p({\bf t}|\alpha,\beta)=0\\ &&\Leftrightarrow \frac{M}{2\alpha}-\frac{1}{2}{\bf m}_N^\top{\bf m}_N-\frac{1}{2}\sum_{i=1}^M\frac{1}{\lambda_i+\alpha}=0\\ &&\Leftrightarrow\alpha{\bf m}_N^\top{\bf m}_N=M-\sum_{i=1}^M\frac{\alpha}{\lambda_i+\alpha}\\ &&\Leftrightarrow\alpha{\bf m}_N^\top{\bf m}_N=\sum_{i=1}^M\frac{\lambda_i}{\lambda_i+\alpha}\\ &&\Leftrightarrow\alpha=\frac{\gamma}{{\bf m}_N^\top{\bf m}_N}\tag{96}\\ \end{eqnarray}$

$(96)$ で $\gamma$ を以下のようにおきました。

$\begin{eqnarray} \gamma=\sum_{i=1}^M\frac{\lambda_i}{\lambda_i+\alpha}\tag{97} \end{eqnarray}$

$(96)$ が $\alpha$ の更新式です。
$(96)$ の右辺の $\gamma$ は $\alpha$ を含むので、 $(96)$ は $\alpha$ に関する陰関数になっていることに注意してください。

固有値 $\lambda_i$ は $\beta$ に比例するので、 $\lambda_i=\beta a_i$ とおいて、 $\beta$ で微分すると、

$\begin{eqnarray} \frac{d}{d\beta}\lambda_i&=&\frac{d}{d\beta}\beta a_i\\ &=&a_i\\ &=&\frac{\lambda_i}{\beta}\tag{98} \end{eqnarray}$

となります。
$\ln|{\bf A}|$ を $\beta$ で微分します。

$\begin{eqnarray} \frac{d}{d\beta}\ln|{\bf A}|&=&\frac{d}{d\beta}\ln\prod_{i=1}^M(\lambda_i+\alpha)\\ &=&\frac{d}{d\beta}\sum_{i=1}^M\ln(\lambda_i+\alpha)\\ &=&\sum_{i=1}^M\frac{d}{d(\lambda_i+\alpha)}\ln(\lambda_i+\alpha)\frac{d}{d\beta}(\lambda_i+\alpha)\\ &=&\frac{1}{\beta}\sum_{i=1}^M\frac{\lambda_i}{\lambda_i+\alpha}\\ &=&\frac{\gamma}{\beta}\tag{99} \end{eqnarray}$

対数周辺尤度 $(90)$ を $\beta$ で微分して $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial}{\partial\beta}\ln p({\bf t}|\alpha,\beta)=0\\ &&\Leftrightarrow\frac{N}{2\beta}-\frac{1}{2}||{\bf t}-{\boldsymbol\Phi}{\bf m}_N||^2-\frac{\gamma}{2\beta}=0\\ &&\Leftrightarrow\frac{1}{\beta}=\frac{1}{N-\gamma}||{\bf t}-{\boldsymbol\Phi}{\bf m}_N||^2\tag{100}\\ \end{eqnarray}$

$(100)$ は $\beta$ の更新式です。
$(100)$ の右辺の $\gamma$ は $\lambda_i$ を含んでおり、 $\lambda_i$ は $\beta$ に依存するので、
$(100)$ は $\beta$ に関する陰関数になっています。

$\alpha,\beta$ に関する更新式が求まりましたので、収束するまで繰り返せば、近似解 $\hat{\alpha},\hat{\beta}$ が求まります。
MAP推定とリッジ回帰の関係の記事より、リッジ回帰の $\lambda$ は $\alpha / \beta$ なので、

$\begin{eqnarray} \hat{\lambda}=\frac{\hat{\alpha}}{\hat{\beta}}\tag{101} \end{eqnarray}$

となります。

式 $(93)$ の導出

$\alpha{\bf I}$ の固有値は全て $\alpha$ (重解)なので、固有ベクトルを $\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}$ の固有ベクトルに合わせることができます。
よって、対角化する直交行列は $\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}$ を対角化する直交行列 ${\bf P}$ と同じものとなります。

以下で、式 $(93)$ を導出します。

$\begin{eqnarray} &&{\bf P}^{-1}(\alpha{\bf I}){\bf P}+{\bf P}^{-1}(\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}){\bf P}={\rm diag}(\alpha,\ldots,\alpha)+{\rm diag}(\lambda_1,\ldots,\lambda_M)\\ &&\Leftrightarrow{\bf P}^{-1}(\alpha{\bf I}+\beta{\boldsymbol\Phi}^\top{\boldsymbol\Phi}){\bf P}={\rm diag}(\alpha+\lambda_1,\ldots,\alpha+\lambda_M)\\ &&\Leftrightarrow{\bf P}^{-1}{\bf A}{\bf P}={\rm diag}(\alpha+\lambda_1,\ldots,\alpha+\lambda_M)\tag{102} \end{eqnarray}$

最後に

本当は式 $(78)$ で ${\bf m}_0={\bf 0},{\bf S}_0=\alpha^{-1}{\bf I}$ を代入して式 $(90)$ を導きたかったんですが、計算が合わなかったので
昔の記事で計算したものを載せました。
計算が合ったら書き換えます。計算が合った方はコメントください。

参考文献

パターン認識と機械学習上巻
ベイズ推論による機械学習入門

やりたい事

概要

データ

変数

注意

線形基底関数モデル

曲線フィッティング

目的関数の定義

目的関数の最小化

目的関数は凸関数である

モデル比較

リッジ回帰(L2正則化)

過学習と正則化

正則化の効果

解析解

最尤推定

事後分布とMAP推定

wの事前分布の平均が0で共分散行列が等方的な場合

wの事前分布に特別な仮定を置かない場合

予測分布

導出その1 - ガウス分布の周辺化に関する公式を使う方法

導出その2

ベイズモデル比較

周辺尤度

周辺尤度が大きい方が正しいモデルである理由

周辺尤度によるモデル比較 - M

一般の周辺尤度

このモデルの周辺尤度

周辺尤度の計算

周辺尤度によるモデル比較 - α、β

周辺尤度

周辺尤度の最大化

式の導出

最後に

参考文献

式 $(93)$ の導出