Bayesian Ridge Regression (BRR)
BRRの定式化
ここではギブスサンプリングを用いたBRR（Bayesian Ridge Regression）の実装を試みる。BRRはその名の通り、通常のリッジ回帰をベイズ統計の枠組みに拡張したモデルであり、マーカー効果を確率変数として扱う点に特徴がある。
通常のリッジ回帰では、表現型とマーカー情報の関係は以下の線形モデルによって表される。
\begin {align*}
\mathbf{y} = \mathbf{X}\mathbf{w} + \mathbf{e}
\end {align*}
ここで，$\mathbf{y}$は$n \times 1$の表現型値ベクトル、$\mathbf{X}$は$n \times L$のマーカー行列、$\mathbf{w}$は$L \times 1$のマーカー効果ベクトル、$\mathbf{e}$は$n \times 1$の誤差ベクトルである。今回は説明簡略化のため、表現型値の数と遺伝子型の数は等しいものとする（つまり各遺伝子型に対して反復が存在しない）。
通常のリッジ回帰では、マーカー効果 $\mathbf{w}$は以下の式により推定される。
\begin {align*}
\hat{\mathbf{w}} = (\mathbf{X}^T\mathbf{X} + \lambda \mathbf{I} )^{-1}\mathbf{X}^T\mathbf{y}
\end {align*}
ここで、$\lambda$は正則化パラメータである。上記の式ではベイズ的な要素、つまり事前知識は全く入っておらず、解は一意に決まる。BRRでは、マーカー効果は正規分布に従うはずであるという事前知識を入れ込んだ形でマーカー効果の推定を行う。つまり、
\begin {align*}
p(\mathbf{w}|\sigma_w^2) = N(\mathbf{w} | 0, \mathbf{I}\sigma_w^2)
\end {align*}
となる。ここで、$\sigma_w^2$はマーカー効果の分散である。そのため、$\sigma_w^2$がパラメータとなるが、さらに、$\sigma_w^2$は逆カイ二乗分布に従うはずだと仮定する。
\begin {align*}
p(\sigma_w^2) = \chi^{-2}(\sigma_w^2 | \nu_w, \tau_w)
\end {align*}
ここで、$\nu_w$、$\tau_w$はパラメータ$\sigma_w^2$のパラメータであるということから超パラメータとよばれる。マーカー効果と同じように同じように残差$\mathbf{e}$が平均0、誤差分散$\sigma_e^2$の正規分布に従うと仮定すると、
\begin {align*}
p(\mathbf{e}|\sigma_e^2) = N(\mathbf{e} | 0, \mathbf{I}\sigma_e^2) \\
p(\sigma_e^2) = \chi^{-2}(\sigma_e^2 | \nu_e, \tau_e)
\end {align*}
となる。ここで、$\nu_e$、$\tau_e$はパラメータ$\sigma_e^2$の超パラメータである。まず、データ（$\mathbf{y}$）をもとにこれらの3パラメータ、$\mathbf{w}$、$\sigma_w^2$、$\sigma_e^2$を直接推定するためには、これら3変数の同時事後分布$p(\mathbf{w},\sigma_w^2,\sigma_e^2|\mathbf{y})$をベイズの定理をもとに考える。
\begin {align*}
p(\mathbf{w},\sigma_w^2,\sigma_e^2|\mathbf{y})
&\propto p(\mathbf{y}|\mathbf{w},\sigma_e^2)p(\mathbf{w} | \sigma_{w}^2) p(\sigma_{w}^2 | \nu_{w}, \tau_{w}) p(\sigma_{e}^2 | \nu_{e}, \tau_{e})\\
&=N(\mathbf{y}|\mathbf{X}\mathbf{w}, \mathbf{I}\sigma_e^2)N(\mathbf{w} | 0, \mathbf{I}\sigma_w^2)\chi^{-2}(\sigma_w^2 | \nu_w, \tau_w)\chi^{-2}(\sigma_e^2 | \nu_e, \tau_e)
\end {align*}
上記の事後分布$p(\mathbf{w},\sigma_w^2,\sigma_e^2|\mathbf{y}) = p(\theta | \mathbf{y})$を計算することができれば、周辺化（$p(\mathbf{w}|\mathbf{y}) = \iint p(\mathbf{w},\sigma_w^2,\sigma_e^2|\mathbf{y})d\sigma_w^2 d\sigma_e^2$）を行うことによって事後分布を求めたり（実際は周辺化も逆カイ二乗分布の存在により難しい気がする）、新たな入力データに対して予測分布を求めたりすることができる。
しかし、上記の事後分布を考えるにあたっては、$p(\mathbf{y})=\iiint p(\mathbf{y} | \mathbf{w},\sigma_w^2,\sigma_e^2)p(\mathbf{w},\sigma_w^2,\sigma_e^2)d\mathbf{w} d\sigma_w^2 d\sigma_e^2$の正規化定数を考える必要があるが、この時、逆カイ二乗分布の存在により解析的には求められない（はず）。そのため、ギブスサンプリングを活用し、各パラメータの事後分布からサンプリングを行うことで、パラメータの推定を行っていく。
ギブスサンプリング
ギブスサンプリングを行うにあたっては、各パラメータの事後分布を考える必要がある。つまり、$p(\mathbf{w}|\sigma_w^2,\sigma_e^2,\mathbf{y})$、$p(\sigma_w^2|\mathbf{w}, \sigma_e^2,\mathbf{y})$、$p(\sigma_e^2|\mathbf{w},\sigma_w^2,\mathbf{y})$である。ここについてもベイズの定理を使って考えていく。
マーカー効果の事後分布
マーカー効果の事後分布は、
\begin {align*}
p(\mathbf{w}|\sigma_w^2,\sigma_e^2,\mathbf{y})
&\propto p(\mathbf{y}|\mathbf{w}, \sigma_e^2)p(\mathbf{w} | \sigma_{w}^2) \\
&=N(\mathbf{y}|\mathbf{X}\mathbf{w}, \mathbf{I}\sigma_e^2)N(\mathbf{w} | 0, \mathbf{I}\sigma_w^2)
\end {align*}
となる。ここで、マーカー効果の事前分布として設定した正規分布は共役事前分布であるため、事後分布$p(\mathbf{w}|\sigma_w^2,\sigma_e^2,\mathbf{y})$も正規分布となる。共役事前分布である特徴を生かして式変形を行うと、
\begin {align*}
p(\mathbf{w}|\sigma_w^2,\sigma_e^2,\mathbf{y}) = N(\mathbf{w} | \mathbf{m}, \mathbf{V})
\end {align*}
と表すことができる。ここで、$\mathbf{m}$はマーカー効果の平均値ベクトル、$\mathbf{V}$はマーカー効果の分散共分散行列であり、
\begin {align*}
\mathbf{m} = \left( \mathbf{X}^{T} \mathbf{X} + \frac{\sigma_e^2}{\sigma_w^2} \mathbf{I} \right)^{-1} \mathbf{X}^{T} \mathbf{y}\\
\mathbf{V} = \left( \mathbf{X}^{T} \mathbf{X} + \frac{\sigma_e^2}{\sigma_w^2} \mathbf{I} \right)^{-1} \sigma_e^2
\end {align*}
となる。
マーカー効果の分散の事後分布
マーカー効果の分散の事後分布は、
\begin {align*}
p(\sigma_w^2|\mathbf{w}, \sigma_e^2,\mathbf{y})
&=p(\sigma_w^2|\mathbf{w})\\
&\propto p(\mathbf{w} | \sigma_{w}^2) p(\sigma_{w}^2 | \nu_{w}, \tau_{w})\\
&=N(\mathbf{w} | 0, \mathbf{I}\sigma_w^2)\chi^{-2}(\sigma_w^2 | \nu_w, \tau_w)
\end {align*}
となる。ここで、マーカー効果の分散の事前分布として設定した逆カイ二乗分布は共役事前分布であるため、事後分布$p(\sigma_w^2|\mathbf{w})$も逆カイ二乗分布となる。逆カイ二乗分布についてはもう少し詳しく更新後の分布について見ていく。ここで、尤度関数$p(\mathbf{w} | \sigma_{w}^2)$は以下のように表される。
\begin {align*}
p(\mathbf{w} \mid \sigma_w^2)
&= \prod_{\ell=1}^{L} \mathcal{N}(w_\ell \mid 0, \sigma_w^2)\\
&= \prod_{\ell=1}^{L} \left( \frac{1}{\sqrt{2\pi \sigma_w^2}} \exp\left\{ -\frac{w_\ell^2}{2\sigma_w^2} \right\} \right)\\
&\propto (\sigma_w^2)^{-L/2} \exp\left( -\frac{1}{2\sigma_w^2} \sum_{\ell=1}^{L} w_\ell^2 \right)
\end {align*}
ここでは、$L$をマーカー数としている。事前分布$p(\sigma_{w}^2 | \nu_{w}, \tau_{w})$については、
\begin {align*}
p(\sigma_{w}^2 | \nu_{w}, \tau_{w}) = (\sigma_w^2)^{-(\nu_w/2 + 1)} \exp\left( -\frac{\nu_w s_w^2}{2 \sigma_w^2} \right)
\end {align*}
であるため、
\begin {align*}
p(\sigma_w^2|\mathbf{w})
&\propto p(\mathbf{w} | \sigma_{w}^2) p(\sigma_{w}^2 | \nu_{w}, \tau_{w}) \\
&=N(\mathbf{w} | 0, \mathbf{I}\sigma_w^2)\chi^{-2}(\sigma_w^2 | \nu_w, \tau_w)\\
&\propto (\sigma_w^2)^{-(L + \nu_w)/2 - 1} \exp\left( -\frac{1}{2\sigma_w^2} \left( \sum w_\ell^2 + \nu_w s_w^2 \right) \right)
\end {align*}
となる。そのため、マーカー効果の分散の事後分布は、
\begin {align*}
p(\sigma_w^2 \mid \mathbf{w}) &= \chi^{-2}(\sigma_w^2 \mid \nu_L, s_L^2)\\
\nu_L &= L + \nu_w\\
s_L^2 &= \frac{ \sum w_\ell^2 + \nu_w s_w^2 }{L + \nu_w}
\end {align*}
となる。
誤差分散の事後分布
誤差分散の事後分布は、
\begin {align*}
p(\sigma_e^2|\mathbf{w}, \sigma_w^2,\mathbf{y})
&=p(\sigma_e^2|\mathbf{w},\mathbf{y})\\
&\propto p(\mathbf{y} \mid \mathbf{w} , \sigma_e^2) \, p(\sigma_e^2)\\
&=N(\mathbf{y} | \mathbf{X}\mathbf{w}, \mathbf{I}\sigma_e^2)\chi^{-2}(\sigma_e^2 | \nu_e, \tau_e)
\end {align*}
となる。誤差分散の事前分布が逆カイ二乗分布であることから、マーカー効果の分散の更新方法と全く同じように更新できるため、
\begin {align*}
p(\sigma_e^2 \mid \mathbf{w}, \mathbf{y}) &= \chi^{-2}(\sigma_e^2 \mid \nu_N, s_N^2)\\
\nu_N &= N + \nu_e\\
s_N^2 &= \frac{ \sum (y_i - \mathbf{x}_i^T \mathbf{w})^2 + \nu_e s_e^2 }{N + \nu_e}
\end {align*}
となる。
完全事後分布のまとめ
上記をまとめると各パラメータの事後分布は以下のようになる。
\begin {align*}
p(\mathbf{w}|\sigma_w^2,\sigma_e^2,\mathbf{y}) = N(\mathbf{w} | \mathbf{m}, \mathbf{V})\\
p(\sigma_w^2 \mid \mathbf{w}) = \chi^{-2}(\sigma_w^2 \mid \nu_L, s_L^2)\\
p(\sigma_e^2 \mid \mathbf{w}, \mathbf{y}) = \chi^{-2}(\sigma_e^2 \mid \nu_N, s_N^2)
\end {align*}
これらの完全事後分布をもとにして、ギブスサンプリングを行っていく。ひとまず今回はここまで。