PRML下巻を読んでいて、線形回帰を何度も例に出しながら各種の推定手法を説明するやり方がわかりやすいと感じたので、その部分をまとめてみます。
なので内容はPRML下巻からかいつまんだ焼き直しです（一部の記法や説明の仕方は書き換えています）。
EMアルゴリズムをまとめた次回の記事はこちらです
→線形回帰で理解するカーネル法・ガウス過程・RVM・EMアルゴリズム・変分ベイズ(2)
まとめる手法はタイトルにある通り
- カーネル法
- ガウス過程回帰
- RVM(関連ベクトルマシン)
- EMアルゴリズム
- 変分ベイズ
です。今回の記事はRVMまでです。
自分の理解が完璧ではなく誤りを含む可能性があるので、誤りを見つけた方はコメントで教えていただけると非常に助かります…。
線形回帰
まずは、すべての手法で共通して用いる設定を述べます（といいつつ、各手法では設定を少しずつ変更しています）。
$ n = 1,\ \ldots ,\ N $ に対して目標変数 $ t_n \in \mathbb{R} $ をパラメータ $
\boldsymbol{w} \in \mathbb{R}^d $ の線形結合 $ \boldsymbol{w}^{\top}\boldsymbol\phi(\boldsymbol x_n) $ で回帰することを考えます。
ただし $ \boldsymbol x_n $ は観測される入力変数、$ \phi(\boldsymbol x_n) \in \mathbb R^d $ は 何らかの写像 $ \boldsymbol \phi $ で入力変数を特徴空間に移したベクトルです。
例えば $ \boldsymbol \phi(x) = (1,\ x,\ x^2) $ とすれば2次の多項式回帰となります。
考える統計モデルは
\begin{align}
t_n &= \boldsymbol w^{\top} \boldsymbol \phi(x_n) + \varepsilon_n\ , \\
\varepsilon_n &\sim \mathcal{N}(0,\ \beta^{-1})\ , \\
\boldsymbol w &\sim \mathcal{N}(\boldsymbol 0,\ \alpha^{-1}I) \tag{1}
\end{align}
です。 $ \varepsilon_n\ (n = 1,\ \ldots ,\ N) $ は独立とします。
すなわち、誤差 $\varepsilon_n$ は独立に平均 $0$, 分散 $\beta^{-1}$ の正規分布に従い、
パラメータ $ \boldsymbol w $ の事前分布は平均 $\boldsymbol 0$, 共分散行列 $\alpha^{-1}I$ の $ d $ 次元正規分布とします。
後のために式 (1) をベクトル表記しておくと
\begin{align}
{\bf t} &= \Phi\boldsymbol w + \boldsymbol \varepsilon\ , \\
\boldsymbol \varepsilon &\sim \mathcal{N}(\boldsymbol 0,\ \beta^{-1}I)\ , \\
\boldsymbol w &\sim \mathcal{N}(\boldsymbol 0,\ \alpha^{-1}I) \tag{2}
\end{align}
となります。ただし
{\bf t} =
\begin{pmatrix}
t_1 \\
\vdots \\
t_N
\end{pmatrix}
\in \mathbb{R}^{N} , \quad
\boldsymbol \varepsilon =
\begin{pmatrix}
\varepsilon_1 \\
\vdots \\
\varepsilon_N
\end{pmatrix}
\in \mathbb{R}^{N} , \quad
\Phi =
\begin{pmatrix}
\boldsymbol \phi(\boldsymbol x_1)^{\top} \\
\vdots \\
\boldsymbol \phi(\boldsymbol x_N)^{\top}
\end{pmatrix}
\in \mathbb{R}^{N \times d}
です。
このとき、$\boldsymbol w$ で条件づけた ${\bf t}$ の条件付き分布と ${\bf t}$ の周辺分布は
\begin{align}
p({\bf t}| \boldsymbol w) &= \mathcal{N}(\Phi\boldsymbol w, \beta^{-1}I)\ , \\
p({\bf t}) &= \mathcal{N}(\boldsymbol 0, \alpha^{-1}\Phi\Phi^{\top} + \beta^{-1}I)
\end{align}
となります（2つ目の式は The Matrix Cookbook 式 (355) より）。
目標は、新たな入力 $ \boldsymbol x $ に対する $ t $ の予測値、あるいは予測分布 $p(t | {\bf t}) $ を求めることです。
MAP推定
カーネルを導入する前に、上で述べたモデルのパラメータ $ \boldsymbol w $ をMAP推定で求めてみます。
MAP推定は $ \boldsymbol w $ の事後分布 $ p(\boldsymbol w | {\bf t})$ を最大化することで $ \boldsymbol w $ を点推定します。
事後分布の対数をとると
\begin{align}
\log p(\boldsymbol w| {\bf t})
&= \log p({\bf t} | \boldsymbol w) + \log p(\boldsymbol w) + \text{const.} \\
&= - \frac{\beta}{2} \| {\bf t} - \Phi \boldsymbol w \|^2 - \frac{\alpha}{2} \| \boldsymbol w \|^2 + \text{const.} \tag{3}
\end{align}
となります。ここで $ \boldsymbol w $ に関係ない定数項は $ \text{const.} $ にまとめました。
式 (3) を最大化して $\boldsymbol w$ のMAP推定値を求めるわけですが、式 (3) の最大化は L2 正則化付き最小二乗法に対応することがわかります。
つまり、L2正則化付きの最小二乗法でパラメータを求めるということは目標変数が正規分布に従うと仮定してパラメータをMAP推定していることになるということです(L2正則化項をなくせば最尤推定になる)。
式 (3) を最大化するためには $\boldsymbol w$ で微分して $\boldsymbol 0$ とおけばよいので、
-\beta(\Phi^{\top}\Phi\boldsymbol w - \Phi^{\top}{\bf t}) - \alpha\boldsymbol w = \boldsymbol 0
より
\boldsymbol w_{\text{MAP}} = (\beta\Phi^{\top}\Phi + \alpha I)^{-1}\Phi^{\top}{\bf t}
となります。
よって $ t $ の予測分布は
\begin{align}
p(t | {\bf t})
&= p(t | \boldsymbol w_{\text{MAP}}) \\
&= \mathcal{N}(\boldsymbol w_{\text{MAP}}^{\top}\boldsymbol \phi(\boldsymbol x),\ \beta^{-1})
\end{align}
となります。 $ t $ の予測値が欲しい場合はこの平均 $ w_{\text{MAP}}^{\top}\boldsymbol \phi(\boldsymbol x) $ を用います。
カーネル法を導入するために、次の等式を証明しておきます。
(\beta\Phi^{\top}\Phi + \alpha I)^{-1}\Phi^{\top} = \Phi^{\top}(\beta\Phi\Phi^{\top} + \alpha I)^{-1} \tag{4}
これを示すためには、等式の両辺に左から $ \beta\Phi^{\top}\Phi + \alpha I\ , $ 右から $ \beta\Phi\Phi^{\top} + \alpha I $ をかけた以下の等式
\Phi^{\top}(\beta\Phi\Phi^{\top} + \alpha I) = (\beta\Phi^{\top}\Phi + \alpha I)\Phi^{\top}
が成り立つことを示せばよい（なぜなら左右からかける行列はいずれも正定値、よって正則なので）ですが、この等式は成り立つことがすぐわかるので、それと等価な式 (4) も成り立つことがわかります。
カーネル法
上で述べたMAP推定の結果を、カーネルを用いて書き直してみます。
先ほどMAP推定で求めた $ t $ の予測値（ = 予測分布の平均）は次のように書けます。
\begin{align}
\boldsymbol w_{\text{MAP}}^{\top} \boldsymbol \phi(\boldsymbol x)
&= ((\beta\Phi^{\top}\Phi + \alpha I)^{-1}\Phi^{\top}{\bf t})^{\top}\boldsymbol \phi(\boldsymbol x) \\
&= (\Phi^{\top}(\beta\Phi\Phi^{\top} + \alpha I)^{-1}{\bf t})^{\top}\boldsymbol \phi(\boldsymbol x) \\
&= {\bf t}^{\top}(\beta\Phi\Phi^{\top} + \alpha I)^{-1}\Phi\boldsymbol \phi(\boldsymbol x) \tag{5}
\end{align}
ここで、1行目から2行目の変形に式 (4) を使いました。
最終行に現れた $ \Phi\Phi^{\top} $ と $ \Phi\boldsymbol \phi(\boldsymbol x) $ の成分を書いてみると
\begin{align}
\Phi\Phi^{\top}
&=
\begin{pmatrix}
\boldsymbol \phi(\boldsymbol x_1)^{\top} \\
\vdots \\
\boldsymbol \phi(\boldsymbol x_N)^{\top}
\end{pmatrix}
\begin{pmatrix}
\boldsymbol \phi(\boldsymbol x_1) & \cdots & \boldsymbol \phi(\boldsymbol x_N)
\end{pmatrix} \\
&=
\begin{pmatrix}
& \vdots & \\
\cdots & \boldsymbol \phi(\boldsymbol x_i)^{\top}\boldsymbol \phi(\boldsymbol x_j) & \cdots \\
& \vdots &
\end{pmatrix}\ , \\
\Phi\boldsymbol \phi(\boldsymbol x)
&=
\begin{pmatrix}
\boldsymbol \phi(\boldsymbol x_1)^{\top} \\
\vdots \\
\boldsymbol \phi(\boldsymbol x_N)^{\top}
\end{pmatrix}
\boldsymbol \phi(\boldsymbol x) \\
&=
\begin{pmatrix}
\boldsymbol \phi(\boldsymbol x_1)^{\top}\boldsymbol \phi(\boldsymbol x) \\
\vdots \\
\boldsymbol \phi(\boldsymbol x_N)^{\top}\boldsymbol \phi(\boldsymbol x)
\end{pmatrix}
\end{align}
となり、すべての成分が $ \phi(\boldsymbol x_i) $ と $ \phi(\boldsymbol x_j) $ の内積、あるいは $ \phi(\boldsymbol x_i) $ と $ \phi(\boldsymbol x) $ の内積で書けることがわかります。
よってこの内積をカーネルで置き換えれば、予測値をカーネルのみで記述することができます。
カーネル関数 $ k $ を $ k(\boldsymbol x, \boldsymbol x') = \boldsymbol \phi(\boldsymbol x)^{\top}\boldsymbol \phi(\boldsymbol x') $ で定義し、行列 $ K $ とベクトル $ \boldsymbol k(\boldsymbol x) $ を
\begin{align}
K
&=
\begin{pmatrix}
& \vdots & \\
\cdots & k(\boldsymbol x_i, \boldsymbol x_j) & \cdots \\
& \vdots &
\end{pmatrix}\ = \Phi\Phi^{\top}, \\
\boldsymbol k(\boldsymbol x)
&=
\begin{pmatrix}
k(\boldsymbol x_1, \boldsymbol x) \\
\vdots \\
k(\boldsymbol x_N, \boldsymbol x)
\end{pmatrix} = \Phi\boldsymbol \phi(\boldsymbol x)
\end{align}
で定めると、式 (5) で表される $ t $ の予測値は
{\bf t}^{\top} (\beta K + \alpha I )^{-1} \boldsymbol k(\boldsymbol x)
と書けます。
この式にはもはや $ \boldsymbol \phi $ は現れないので、特徴空間への写像を明示的に書く必要がなくなっています。正定値カーネルを使えば $ k(\boldsymbol x, \boldsymbol x') = \boldsymbol \phi(\boldsymbol x)^{\top}\boldsymbol \phi(\boldsymbol x') $ を満たす写像 $ \boldsymbol \phi $ が存在することが知られています。
カーネル関数として例えばガウシアンカーネルを用いると対応する特徴空間は無限次元となるため、$ \boldsymbol \phi $ を明示的に書く定式化のままでは扱えません。
カーネル化の利点の一つは、カーネル関数を定めることで特徴空間の次元 $ d $ に依存しない計算量で高次元空間上のモデル化ができることですが、$ N \times N $ 行列の逆行列が現れることからわかるようにサンプルサイズ $N$ が大きいときは計算量が大きくなるという問題点もあります。
ガウス過程回帰
次は、ガウス過程を用いて線形回帰を扱ってみます。
書いている途中で気づいたのですが、ガウス過程回帰ではもはやパラメータ $\boldsymbol w $ が存在しないので「線形」回帰ではありませんでした、すみません。ガウス過程回帰という表現が正しいと思います。
ガウス過程回帰では式 (2) の統計モデルを以下のように変更します。
\begin{align}
{\bf t} &= \boldsymbol y + \boldsymbol \varepsilon\ , \\
\boldsymbol \varepsilon &\sim \mathcal{N}(\boldsymbol 0,\ \beta^{-1}I)\ , \\
\boldsymbol y &\sim \mathcal{N}(\boldsymbol 0, K)\ , \tag{6}
\end{align}
ただし
\begin{align}
K
&=
\begin{pmatrix}
& \vdots & \\
\cdots & k(\boldsymbol x_i, \boldsymbol x_j) & \cdots \\
& \vdots &
\end{pmatrix}
\end{align}
です。
周辺分布は
\begin{align}
p({\bf t}) &= \mathcal{N}(\boldsymbol 0,\ C_N)\, \\
C_N &= K + \beta^{-1}I \tag{7}
\end{align}
となります。後のために行列 $ C_N $ を定義しました。
式 (6) で定義した $K$ は（すぐ下で述べるように）カーネル法で定義した $K$ と完全には対応しないので注意してください。
式 (2) では $ \boldsymbol y $ のところを $ \Phi \boldsymbol w $ と書いており、$ \boldsymbol w \sim \mathcal{N}(\boldsymbol 0,\ \alpha^{-1}I) $ としていました。
したがって $ \Phi \boldsymbol w \sim \mathcal{N}(\boldsymbol 0,\ \alpha^{-1}\Phi\Phi^{\top}) $ です。ガウス過程回帰ではパラメータ $ \boldsymbol w $ を明示的に書く代わりに $\alpha^{-1}\Phi\Phi^{\top}$ を正定値行列 $K$ で置き換えていることになります。
$ \boldsymbol y $ がガウス過程に従うと仮定したとき、新たな入力変数 $ \boldsymbol x $ に対する目標変数 $t$ と今までに得られた $ N $ 個の目標変数 $ {\bf t} $ をつなげた $ N+1 $ 次元ベクトル $ {\bf t}_{N+1} = ({\bf t}^{\top},\ t)^{\top} $ は $ N+1 $ 次元正規分布 $ \mathcal{N}(\boldsymbol 0, C_{N+1}) $ に従います。
ここで
\begin{align}
C_{N+1}
&=
\begin{pmatrix}
C_N & \boldsymbol k(\boldsymbol x) \\
\boldsymbol k(\boldsymbol x)^{\top} & k(\boldsymbol x, \boldsymbol x) + \beta^{-1}
\end{pmatrix}\ , \\
\boldsymbol k(\boldsymbol x)
&=
\begin{pmatrix}
k(\boldsymbol x_1, \boldsymbol x) \\
\vdots \\
k(\boldsymbol x_N, \boldsymbol x)
\end{pmatrix}
\end{align}
です。
多次元正規分布の条件付き分布の公式を用いると、予測分布 $ p(t | {\bf t}) $ は
p(t | {\bf t}) = \mathcal{N}\left(\boldsymbol k(\boldsymbol x)^{\top}C_N^{-1}{\bf t},\ k(\boldsymbol x, \boldsymbol x) + \beta^{-1} - \boldsymbol k(\boldsymbol x)^{\top}C_N^{-1}\boldsymbol k(\boldsymbol x)\right)
と書けます。
先ほどのカーネル法で求めた予測分布と比べると、平均は一致していますが、分散はカーネル法では定数 $ \beta^{-1} $ としているのに対して、ガウス過程回帰では $ \boldsymbol y $ がガウス過程に従うと仮定したことで適応的に推定できています。
ハイパーパラメータの(点)推定
ガウス過程回帰におけるハイパーパラメータ $\beta$ を周辺尤度 $p({\bf t}; \beta)$ を最大化することで推定してみます。
ここで、周辺尤度がハイパーパラメータ $\beta$ に依存することを明示的に表すために $p({\bf t}; \beta)$ と表記しています。
対数尤度は
\log p({\bf t}; \beta) = -\frac{1}{2} \log |C_N| - \frac{1}{2} {\bf t}^{\top}C_N^{-1}{\bf t} + \text{const.}
です（$\beta$に依存しない項は$\text{const.}$にまとめています）。
これを$\beta^{-1}$で微分すると
\begin{align}
\frac{\partial}{\partial \beta^{-1}} \log p({\bf t}; \beta)
&= -\frac{1}{2}\text{tr}\left(C_N^{-1}\frac{\partial C_N}{\partial \beta^{-1}}\right) + \frac{1}{2} {\bf t}^{\top}C_N^{-1}\frac{\partial C_N}{\partial \beta^{-1}}C_N^{-1}{\bf t} \\
&= -\frac{1}{2} \text{tr}(C_N^{-1}) + \frac{1}{2} {\bf t}^{\top}C_N^{-2}{\bf t}
\end{align}
となります。ここで行列の微分に関する公式（The Matrix Cookbook 式 (40), (43)）と、$C_N$の式 (7) を使いました。
これが陽に解ければ簡単なのですが、自分はこれを解く方法がわからなかったので、実用上は勾配法や準ニュートン法などを用いて数値的に対数尤度を最大化することで $\beta$ を最適化すればよいと思います。
RVM
次に、線形回帰 (2) をRVM（関連ベクトルマシン）を用いて解いてみます。
RVMは、SVMのように疎な解が得られることが特徴で、カーネル法でサンプルサイズ $ N $ が大きいと計算量が大きくなるという問題を、予測に寄与しないサンプルに対する係数が自動的に $0$ になる定式化を行うことで回避しています。
式 (2) はカーネル化する前の式なのですが、この式を用いて得た結果をあとで $K = \Phi\Phi^{\top} $ と置き換えればよいので、式 (2) を用いて話を進めます。
RVMでは、式 (2) で $\boldsymbol w$ の事前分布を等方的な正規分布 $\mathcal{N}(\boldsymbol 0,\ \alpha^{-1}I)$ から $\mathcal{N}(\boldsymbol 0,\ \text{diag}(\alpha_i^{-1}))$ に変更します。
$\text{diag}(\alpha_i^{-1})$ は対角成分に $\alpha_i^{-1}$ を持つ $d\times d$ 対角行列です。
RVMはSVMと名前が似ていますが、SVMは対応する統計モデルが存在せず、ロス関数最小化により分類基準を学習するのに対して、RVMは目標変数とパラメータに正規分布を仮定したベイズ推定の一種であることに注意してください。
回帰のためのSVM、分類のためのRVMも構成することができますがここでは触れません（詳細はPRML下巻を参照してください）。いずれにしてもSVMには統計モデルが対応せず、RVMはベイズ推定であることは同じです（と自分は理解しています）。
RVMの定式化を書いておくと
\begin{align}
{\bf t} &= \Phi\boldsymbol w + \boldsymbol \varepsilon\ , \\
\boldsymbol \varepsilon &\sim \mathcal{N}(\boldsymbol 0,\ \beta^{-1}I)\ , \\
\boldsymbol w &\sim \mathcal{N}(\boldsymbol 0,\ \text{diag}(\alpha_i^{-1})) \tag{8}
\end{align}
です。上で述べたように、式 (2) から $\boldsymbol w$ の事前分布のみ変更しています。
このとき、$\boldsymbol w$ で条件づけた ${\bf t}$ の条件付き分布と ${\bf t}$ の周辺分布は
\begin{align}
p({\bf t}| \boldsymbol w) &= \mathcal{N}(\Phi\boldsymbol w, \beta^{-1}I)\ , \\
p({\bf t}) &= \mathcal{N}(\boldsymbol 0, \Phi A^{-1}\Phi^{\top} + \beta^{-1}I)
\end{align} \tag{9}
です。ここで、行列 $A = \text{diag}(\alpha_i) $ を定義しました。つまり $A^{-1} = \text{diag}(\alpha_i^{-1}) $ です。
式 (2) を最初に扱った時は $\boldsymbol w$ をMAP推定しましたが、ベイズ的に取り扱うためには
$\boldsymbol w$ の事後分布 $p(\boldsymbol w | {\bf t}) $ を求める必要があります。
というのも予測分布 $p(t|{\bf t}) $ を求めるためには
p(t|\boldsymbol t) = \int p(t|\boldsymbol w)p(\boldsymbol w|{\bf t}) \text{d}\boldsymbol w
という計算が必要で、ここに $\boldsymbol w$ の事後分布 $p(\boldsymbol w | {\bf t}) $ が現れるためです。
事後分布 $p(\boldsymbol w | {\bf t}) $ を求めます。
対数尤度を計算すると
\begin{align}
\log p(\boldsymbol w| {\bf t})
&= \log p({\bf t} | \boldsymbol w) + \log p(\boldsymbol w) + \text{const.} \\
&= - \frac{\beta}{2} \| {\bf t} - \Phi \boldsymbol w \|^2 - \frac{1}{2} \boldsymbol w^{\top}A\boldsymbol w + \text{const.} \\
&= - \frac{1}{2} \boldsymbol w^{\top} (A + \beta \Phi^{\top}\Phi) \boldsymbol w + \boldsymbol w^{\top}\beta\Phi^{\top}{\bf t} + \text{const.} \tag{10}
\end{align}
となります。これは $\boldsymbol w$ の2次式なので事後分布はふたたび正規分布になり、その共分散行列を $ S $ とおくと $ S = (A + \beta \Phi^{\top}\Phi)^{-1} $ であることがわかります。
事後分布の平均を $\boldsymbol m$ とおくと、式 (10) と
-\frac{1}{2} (\boldsymbol w - \boldsymbol m)^{\top}(A + \beta \Phi^{\top}\Phi)(\boldsymbol w - \boldsymbol m)
の $\boldsymbol w$ に関する1次の項が一致することから、
\boldsymbol m = \beta S \Phi^{\top} {\bf t}
とわかります。
事後分布 $p(\boldsymbol w | {\bf t})$ がわかったので予測分布 $p(t|{\bf t})$ を求めることができます。積分を直接計算してもよいですが、予測分布が（実は）正規分布になることからその平均と共分散行列だけ求めればよいです。
\begin{align}
p(t|\boldsymbol w) &= \mathcal{N}(\boldsymbol w^{\top}\boldsymbol \phi(\boldsymbol x),\ \beta^{-1})\ , \\
p(\boldsymbol w | {\bf t}) &= \mathcal{N}(\boldsymbol m,\ S)
\end{align}
より、予測分布は
p(t|{\bf t}) = \mathcal{N}(\boldsymbol m^{\top}\boldsymbol \phi(\boldsymbol x), \ \beta^{-1} + \boldsymbol \phi(\boldsymbol x)^{\top}S\boldsymbol \phi(\boldsymbol x))
となります。
ハイパーパラメータの(点)推定
ここまでは単に尤度と事前分布を定めて事後分布と予測分布を求めるという普通のベイズ推定を行っただけで、疎な解は全く得られていません。
疎な解が得られる理由は、ハイパーパラメータ $(\alpha_i),\ \beta$ を周辺尤度最大化によって推定するといくつかの $ \alpha_i $ は非常に大きな値になり、対応する $\boldsymbol w$ は予測に寄与しなくなるためです。
なぜそのようになるかの詳細はPRML下巻を参照してください。
ここからは、ハイパーパラメータ $(\alpha_i),\ \beta$ を周辺尤度 $p({\bf t})$ を最大化することで点推定します。
（この部分は理解が怪しいので誤っているかもしれません）
このように、ハイパーパラメータを周辺尤度最大化によって点推定する手法を経験ベイズ、または第二種の最尤推定と呼びます。
ハイパーパラメータにも事前分布を設定し、完全にベイズ的に取り扱う手法は階層ベイズと呼ばれます。
後で扱うEMアルゴリズムは経験ベイズに分類され、変分ベイズは階層ベイズに分類されます。
$ C = \Phi A^{-1}\Phi^{\top} + \beta^{-1}I $ とおくと、
式 (9) より対数周辺尤度は
\begin{align}
\log p({\bf t}; (\alpha_i), \beta)
&= - \frac{1}{2} \log |C| - \frac{1}{2}{\bf t}^{\top}C^{-1}{\bf t} + \text{const.}
\end{align}
となります。
これを $\alpha_i^{-1},\ \beta^{-1}$ で微分して……とやりたいところですが、以下で述べる $\alpha_i^{-1},\ \beta^{-1}$ の反復更新式を求めるには式変形の工夫が必要です。
PRML演習問題の解答の7.12の結果を用いると、
対数周辺尤度の微分を $0$ とおいた式から
\begin{align}
\alpha_i^{-1} &= \frac{m_i^2}{1 - \alpha_i S_{ii}}\ , \\
\beta^{-1} &= \frac{\| {\bf t} - \Phi\boldsymbol w \|^2}{N - \sum_i(1 - \alpha_i S_{ii})}
\end{align}
が得られます。ただし $m_i$ は $ \boldsymbol w $ の事後平均 $\boldsymbol m$ の第 $i$ 成分です。
上の2つの式は、右辺も $\alpha_i,\ \beta$ に依存しているため、ハイパーパラメータを陽に求める式にはなっていません。
そこで、$\alpha_i,\ \beta$ の値を初期化して $\boldsymbol m,\ S$ を求め、上式の右辺を用いて $\alpha_i,\ \beta$ の値を更新し……という反復計算を行うことでハイパーパラメータを求めます。
後で述べるEMアルゴリズムによって導かれるハイパーパラメータの更新式を用いれば周辺尤度が単調に増加することが保証される一方で、いま述べた反復計算は収束する保証がない（と思う）のですが、実用上はこちらの式を用いるほうが収束が速いようです（PRML下巻本文より）。
ここまでのまとめ
ここまで、線形回帰をカーネル法・ガウス過程・RVMで解く方法をまとめました。
長くなってしまったのでEMアルゴリズムと変分ベイズは別の記事にします。
まずパラメータ $\boldsymbol w$ のMAP推定から始めて、カーネル法では特徴空間への写像 $\boldsymbol \phi$ を陽に書かずカーネル関数 $k(\boldsymbol x,\ \boldsymbol x')$ のみで予測値を表現しました。
ガウス過程ではパラメータ $\boldsymbol w$ を用いるのではなく $\boldsymbol y$ がガウス過程に従うという仮定を置くことで、条件付き分布の計算によって予測分布を求めました。周辺尤度最大化によるハイパーパラメータの点推定も行いました。
RVMでは、$\boldsymbol w$ の事後分布を求めることで予測分布を計算したあと、周辺尤度最大化によるハイパーパラメータの点推定を行いました。これにより多くの $\alpha_i$ が無限大に発散することで疎な解が得られるようです。
次の記事では、$\boldsymbol w$ を潜在変数とみなしてEMアルゴリズムを適用することでハイパーパラメータの点推定を行います。
変分ベイズでは、ハイパーパラメータにも事前分布（hyperprior）を仮定して、階層ベイズの手法によって予測分布を求めます。
できるだけ正確に書いたつもりですが、誤っているところや怪しいところをみつけたら、コメント等で教えていただけると非常に助かります。