前置き
投稿論文の査読を待っている間、ベイズがなんだかよくわからないままになっていたことに気づきちゃんと勉強することにした。
リッジ回帰をベイズ的に解こうとするなら、変数が正規分布(事前分布)に従うことを仮定すれば良い。この場合正規分布のパラメータ(分散)は実はなんでもよいので、決め打ちでやるしかない。実はリッジのペナルティを決めるのと、正規分布のパラメータを決めることは同義になる(https://qiita.com/yoyoyo11131113/items/e03831863f746d8ea0b7)
正規分布パラメータを決め打ちするのが嫌なら、事前分布の前の階層分布Xを考えて、このXのパラメータ(つまりモデルのハイパーパラメータ)を考えることとする。これがBRR (Bayean Ridge Regression)。
BRR設計では事前分布パラメータ(正規分布のパラメータ)は自動で更新されるようにgibbsサンプリングするので、パラメータ選択がある程度柔軟になるのだと信じる。
線形回帰モデルの設計
線形回帰モデルは，以下のように表せます：
$$
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},
$$
ここで，
$$
\boldsymbol{\varepsilon} \sim \mathrm{N}\left(\mathbf{0}, \sigma^2\mathbf{I}_n\right)
$$
さらに事前分布として
$$
\boldsymbol{\beta} \mid \tau^2
\sim \mathrm{N}\left(\mathbf{0}, \tau^2 \mathbf{I}_p\right).
$$
です。
MAP 推定では $\sigma^2$ と $\tau^2$ を定数として扱いますが，実際にはこれらのハイパーパラメータに分布を仮定することもできます。
Ridge 回帰では
$$
\lambda=\sigma^2/\tau^2
$$
を外部から与える必要がありますが，
$\sigma^2$ や $\tau^2$ にも事前分布を置くことで，
両者をデータから同時に推定することができます。
この完全ベイズ化が Bayesian Ridge Regression (BRR) です。
BRR (Bayesian Ridge Regression)
ここでは $\sigma^2$ と $\tau^2$ にも事前分布を置き，
完全ベイズ的に $\boldsymbol{\beta}$ を推定します。
これにより，Ridge（MAP）では外部から与える必要があった
$\lambda = \sigma^2/\tau^2$ を，データから同時に学習できるようになります。
階層ベイズモデル
BRR は次の階層ベイズモデルとして書けます。
$$
\mathbf{y}\mid\boldsymbol{\beta}, \sigma^2
\sim \mathrm{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I}_n)
$$
$$
\boldsymbol{\beta}\mid\tau^2
\sim \mathrm{N}(\mathbf{0}, \tau^2\mathbf{I}_p)
$$
$$
\sigma^2
\sim \mathrm{IG}(a_\sigma, b_\sigma)
$$
$$
\tau^2
\sim \mathrm{IG}(a_\tau, b_\tau)
$$
ここで $\mathrm{IG}(a,b)$ は形パラメータ $a$，尺度パラメータ $b$ をもつ逆ガンマ分布です。
逆ガンマ分布 $\mathrm{IG}(a,b)$ は
$$
\mathrm{IG}(s \mid a,b) =
\frac{b^{a}}{\Gamma(a)} s^{-(a+1)} \exp\left(-\frac{b}{s}\right),\qquad s>0.
$$
という形を持ちます。$\Gamma(a)$ はガンマ関数です。
$$
\Gamma(x) =
\int_{0}^{\infty}
t^{,x-1},e^{-t},dt.
$$
実際にガンマ関数含めて、定数は無視して
$$
p(s) \propto s^{-(a+1)}\exp\left(-\frac{b}{s}\right),\qquad s>0
$$
と考えると計算が容易になります。これを使います。
正規分布と逆ガンマ分布は共役な組み合わせになっているため，
条件付き事後分布が解析的に求まり，Gibbs サンプリングで容易に推定できます。
1. beta導出
$\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2$ の導出
条件付き事後分布は
$$
p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)
\propto
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2),
p(\boldsymbol{\beta}\mid\tau^2)
$$
から得られます。
1.1 尤度と事前の形
尤度：
$$
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2)
\propto
(\sigma^2)^{-n/2}
\exp\left(
-\frac{1}{2\sigma^2}\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right)
$$
事前：
$$
p(\boldsymbol{\beta}\mid\tau^2)
\propto
(\tau^2)^{-p/2}
\exp\left(
-\frac{1}{2\tau^2}\bigl\lVert\boldsymbol{\beta}\bigr\rVert^2
\right)
$$
$\sigma^2,\tau^2$ は条件付き事後で固定とみなせるので，$\boldsymbol{\beta}$ に依らない係数はまとめて無視すると，
$$
p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)
\propto
\exp\left(
-\frac{1}{2\sigma^2}\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
-\frac{1}{2\tau^2}\bigl\lVert\boldsymbol{\beta}\bigr\rVert^2
\right)
$$
1.2 2次形式の展開と平方完成
ノルムを展開すると，
$$
|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}|^2
= (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^{\mathsf{T}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})
= \mathbf{y}^{\mathsf{T}}\mathbf{y} - 2\boldsymbol{\beta}^{\mathsf{T}}\mathbf{X}^{\mathsf{T}}\mathbf{y} + \boldsymbol{\beta}^{\mathsf{T}}\mathbf{X}^{\mathsf{T}}\mathbf{X}\boldsymbol{\beta}
$$
また，
$$
\bigl\lVert\boldsymbol{\beta}\bigr\rVert^2 = \boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta}
$$
したがって，指数の中身は
$$
-\frac{1}{2\sigma^2}
\left(\mathbf{y}^{\mathsf{T}}\mathbf{y} - 2\boldsymbol{\beta}^{\mathsf{T}}\mathbf{X}^{\mathsf{T}}\mathbf{y} + \boldsymbol{\beta}^{\mathsf{T}}\mathbf{X}^{\mathsf{T}}\mathbf{X}\boldsymbol{\beta}\right) -\frac{1}{2\tau^2}\boldsymbol{\beta}^{\mathsf{T}}\boldsymbol{\beta}.
$$
$\boldsymbol{\beta}$ に依存しない $\mathbf{y}^{\mathsf{T}}\mathbf{y}$ は定数なので捨ててよく，
$$
-\frac{1}{2}
\left[
\boldsymbol{\beta}^{\mathsf{T}}\left(
\frac{1}{\sigma^2}\mathbf{X}^{\mathsf{T}}\mathbf{X}
+\frac{1}{\tau^2}\mathbf{I}_p
\right)\boldsymbol{\beta} - 2\frac{1}{\sigma^2}\boldsymbol{\beta}^{\mathsf{T}}\mathbf{X}^{\mathsf{T}}\mathbf{y}
\right].
$$
ここで
$$
\mathbf{A} = \frac{1}{\sigma^2}\mathbf{X}^{\mathsf{T}}\mathbf{X} +\frac{1}{\tau^2}\mathbf{I}_p,\qquad \mathbf{c} = \frac{1}{\sigma^2}\mathbf{X}^{\mathsf{T}}\mathbf{y}
$$
とおくと，
$$
-\frac{1}{2}\left(
\boldsymbol{\beta}^{\mathsf{T}}\mathbf{A}\boldsymbol{\beta} - 2\boldsymbol{\beta}^{\mathsf{T}}\mathbf{c}
\right)
$$
となります。
ここで平方完成すると，
$$
\boldsymbol{\beta}^{\mathsf{T}}\mathbf{A}\boldsymbol{\beta}
{}- 2\boldsymbol{\beta}^{\mathsf{T}}\mathbf{c}
$$
ただし，
$$
\boldsymbol{\mu}_\beta = \mathbf{A}^{-1}\mathbf{c}
$$
$\mathbf{A}\boldsymbol{\mu}_\beta=\mathbf{c}$ より，
$$
\boldsymbol{\mu}_\beta = \mathbf{A}^{-1}\mathbf{c}
$$
$$
\Sigma_\beta = \mathbf{A}^{-1}
$$
すなわち，
$$
\Sigma_\beta
= \left(
\frac{1}{\sigma^2}\mathbf{X}^{\mathsf{T}}\mathbf{X}
{} + \frac{1}{\tau^2}\mathbf{I}_p
\right)^{-1}
$$
\boldsymbol{\mu}_\beta
= \frac{1}{\sigma^2}\Sigma_\beta\,\mathbf{X}^{\mathsf{T}}\mathbf{y}.
よって
p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2)
\propto
\exp\left( -\frac{1}{2}
(\boldsymbol{\beta}-\boldsymbol{\mu}_\beta)^{\mathsf{T}}
\Sigma_\beta^{-1}
(\boldsymbol{\beta}-\boldsymbol{\mu}_\beta)
\right)
であり，
\boldsymbol{\beta}\mid\mathbf{y},\sigma^2,\tau^2
\sim \mathrm{N}\!\left(\boldsymbol{\mu}_\beta,\ \Sigma_\beta\right)
となります。
2. sigmaの導出
$\sigma^2\mid\mathbf{y},\boldsymbol{\beta}$ の導出
$$
p(\sigma^2\mid\mathbf{y},\boldsymbol{\beta})
\propto
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2),
p(\sigma^2)
$$
2.1 尤度の sigma 依存部分
尤度の $\sigma^2$ 依存部分
$$
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2) =
(2\pi\sigma^2)^{-n/2}
\exp\left(
-\frac{1}{2\sigma^2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right)
$$
個人的には，
正規化定数の中の $|\sigma^2 \mathbf{I}_n|$ が $\sigma^2$ に依存しているため，今回の $\sigma^2$ の事後分布の導出では無視できない という点が注意。
逆に $\beta$ の式では定数として消えていました。
$\pi$ はいらないということで
$$
p(\mathbf{y}\mid\boldsymbol{\beta},\sigma^2)
\propto
(\sigma^2)^{-n/2}
\exp\left(
-\frac{1}{2\sigma^2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right).
$$
ここで $(\sigma^2)^{-n/2}$ が自然に出てきます。
なお行列式の基本性質：
$$
|c\mathbf{A}| = c^{n}|\mathbf{A}|
$$
（$c$: スカラー，$\mathbf{A}$: $n\times n$ 行列）
ここでは $n$ 次元正規分布において
$$
\mathbf{A}=\mathbf{I}_n,\quad c = \sigma^2,
$$
なので，
$$
|\sigma^2\mathbf{I}_n |=(\sigma^2)^n|\mathbf{I}_n|=(\sigma^2)^n.
$$
2.2 事前分布
$$
p(\sigma^2)
\propto
(\sigma^2)^{-(a_\sigma+1)}
\exp\left(-\frac{b_\sigma}{\sigma^2}\right).
$$
2.3 かけ合わせて逆ガンマの形を確認
\begin{aligned}
p(\sigma^2\mid\mathbf{y},\boldsymbol{\beta})
&\propto
(\sigma^2)^{-n/2}
\exp\!\left(
-\frac{1}{2\sigma^2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right)
(\sigma^2)^{-(a_\sigma+1)}
\exp\!\left(-\frac{b_\sigma}{\sigma^2}\right) \\
&=
(\sigma^2)^{-(a_\sigma+1+n/2)}
\exp\!\left(
-\frac{1}{\sigma^2}
\left[
b_\sigma + \frac{1}{2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right]
\right).
\end{aligned}
これは，
$$
\mathrm{IG}!\left(
a_\sigma + \frac{n}{2},
b_\sigma + \frac{1}{2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right)
$$
の核と一致するので，
$$
\sigma^2\mid\mathbf{y},\boldsymbol{\beta}
\sim \mathrm{IG}!\left(
a_\sigma + \frac{n}{2},\
b_\sigma + \frac{1}{2}
\bigl\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\rVert^2
\right).
$$
3. tau2の導出
$\tau^2\mid\boldsymbol{\beta}$ の導出
$$
p(\tau^2\mid\boldsymbol{\beta})
\propto
p(\boldsymbol{\beta}\mid\tau^2),
p(\tau^2)
$$
3.1 beta依存部分
$\boldsymbol{\beta}\mid\tau^2$ の $\tau^2$ 依存部分
$$
\boldsymbol{\beta}\mid\tau^2\sim\mathrm{N}(\mathbf{0},\tau^2\mathbf{I}_p)
$$
その密度は
$$
p(\boldsymbol{\beta}\mid\tau^2)
\propto
(\tau^2)^{-p/2}
\exp!\left(
-\frac{1}{2\tau^2}
\lVert\boldsymbol{\beta}\rVert^2
\right).
$$
ここで $(\tau^2)^{-p/2}$ があとで形状パラメータの $p/2$ に対応します。
3.2 事前分布
$$
p(\tau^2)
\propto
(\tau^2)^{-(a_\tau+1)}
\exp!\left(-\frac{b_\tau}{\tau^2}\right).
$$
3.3 かけ合わせて逆ガンマの形を確認
\begin{aligned}
p(\tau^2\mid\boldsymbol{\beta})
&\propto
(\tau^2)^{-p/2}
\exp\!\left(
-\frac{1}{2\tau^2}
\lVert\boldsymbol{\beta}\rVert^2
\right)
(\tau^2)^{-(a_\tau+1)}
\exp\!\left(-\frac{b_\tau}{\tau^2}\right) \\
&=
(\tau^2)^{-(a_\tau+1+p/2)}
\exp\!\left(
-\frac{1}{\tau^2}
\left[
b_\tau + \frac{1}{2}\lVert\boldsymbol{\beta}\rVert^2
\right]
\right).
\end{aligned}
これは
$$
\mathrm{IG}\left(
a_\tau + \frac{p}{2},
b_\tau + \frac{1}{2}\lVert\boldsymbol{\beta}\rVert^2
\right)
$$
の核と一致するので，
$$
\tau^2\mid\boldsymbol{\beta}
\sim \mathrm{IG}\left(
a_\tau + \frac{p}{2},\
b_\tau + \frac{1}{2}
\bigl\lVert\boldsymbol{\beta}\bigr\rVert^2
\right).
$$
このように，すべての条件付き事後分布が既知の分布で書けるため，
Gibbs サンプリングで繰り返しサンプルすることで
$\boldsymbol{\beta},\sigma^2,\tau^2$ の事後分布を近似できます。
Gibbs サンプラーの流れ
Gibbs サンプリングは，複数の未知パラメータ
$$
\theta = (\theta_1, \theta_2, \dots, \theta_K)
$$
の同時事後分布
$$
p(\theta \mid \text{data})
$$
からサンプリングするための MCMC (Markov Chain Monte Carlo) アルゴリズムの一種です。
Gibbs サンプラーの更新ステップは次のようになります：
- 初期値 $\boldsymbol{\beta}^{(0)}, \sigma^{2(0)}, \tau^{2(0)}$ を用意する。
- $m=1,\dots,M$ について以下を繰り返す：
- $\boldsymbol{\beta}^{(m)} \sim p(\boldsymbol{\beta}\mid\mathbf{y},\sigma^{2(m-1)},\tau^{2(m-1)})$
- $\sigma^{2(m)} \sim p(\sigma^2\mid\mathbf{y},\boldsymbol{\beta}^{(m)})$
- $\tau^{2(m)} \sim p(\tau^2\mid\boldsymbol{\beta}^{(m)})$
バーンイン後のサンプル ${\boldsymbol{\beta}^{(m)}}$ に基づき，
事後平均や信用区間を推定できます。MAP（Ridge）が一点の推定値しか与えないのに対して，
BRR では 事後分布全体を扱える点が特徴です。
R による簡単な実装例
以下に，最もシンプルな BRR の Gibbs サンプラーの例を示します。
set.seed(123)
brr_gibbs <- function(y, X,
n_iter = 5000,
burn_in = 1000,
a_sigma = 1,
b_sigma = 1,
a_tau = 1,
b_tau = 1) {
n <- length(y)
p <- ncol(X)
# 初期値
beta <- rep(0, p)
sigma2 <- 1
tau2 <- 1
# 保存用
n_save <- n_iter - burn_in
beta_save <- matrix(NA_real_, nrow = n_save, ncol = p)
sigma2_save <- numeric(n_save)
tau2_save <- numeric(n_save)
XtX <- crossprod(X)
Xty <- crossprod(X, y)
for (iter in 1:n_iter) {
# 1. beta | y, sigma2, tau2 ~ N(mu_beta, Sigma_beta)
Sigma_beta <- solve(XtX / sigma2 + diag(1 / tau2, p))
mu_beta <- Sigma_beta %*% (Xty / sigma2)
beta <- as.vector(MASS::mvrnorm(1, mu = mu_beta, Sigma = Sigma_beta))
# 2. sigma2 | y, beta ~ IG(a_sigma + n/2, b_sigma + 0.5 * ||y - X beta||^2)
resid <- y - X %*% beta
shape_sigma <- a_sigma + n / 2
rate_sigma <- b_sigma + 0.5 * sum(resid^2)
sigma2 <- 1 / rgamma(1, shape = shape_sigma, rate = rate_sigma)
# 3. tau2 | beta ~ IG(a_tau + p/2, b_tau + 0.5 * ||beta||^2)
shape_tau <- a_tau + p / 2
rate_tau <- b_tau + 0.5 * sum(beta^2)
tau2 <- 1 / rgamma(1, shape = shape_tau, rate = rate_tau)
# 保存
if (iter > burn_in) {
idx <- iter - burn_in
beta_save[idx, ] <- beta
sigma2_save[idx] <- sigma2
tau2_save[idx] <- tau2
}
}
list(
beta = beta_save,
sigma2 = sigma2_save,
tau2 = tau2_save
)
}
n <- 100
p <- 5
X <- matrix(rnorm(n * p), n, p)
beta_true <- c(1, -1, 0.5, 0, 0)
sigma_true <- 0.5
y <- X %*% beta_true + rnorm(n, sd = sigma_true)
res_brr <- brr_gibbs(y = as.numeric(y), X = X,burn_in = 1)
#事後平均による beta 推定
beta_hat_brr <- colMeans(res_brr$beta)
beta_hat_brr
beta_true
$\beta$推定結果の例
[1] 0.98610167 -0.92287089 0.49895582 -0.04017238 0.10070413
[1] 1.0 -1.0 0.5 0.0 0.0
plot(res_brr$beta[,1], type = "l")
plot(res_brr$sigma2, type = "l")
sigmaに関しては0.2-0.3くらいの予測になった。しかし本当は0.5. ややoverfit ?
Hyper parameter は適当に全部1にした。小さい値を設定すると、裾が広くなり、無情報に近くなるはず。
結論
線形回帰モデルにおいて、係数$\beta$に正規分布の事前分布を置き、さらにその分散に逆ガンマ分布を事前分布として与えると、正規—逆ガンマの共役性により各パラメータの事後分布が解析的に求めやすくなる。このとき、
- $\beta \mid \sigma_{\beta}^2, y$ は正規分布、
- $\sigma_{\beta}^2 \mid \beta, y$ は逆ガンマ分布
となり、Gibbsサンプリングによって $\beta$ と分散パラメータを交互にサンプリングしながら事後分布に収束させることができる。高校数学で出てくる逆像法、ファクシミリ論法とかいうのに近い気がする。ただしここで得られるのは点推定ではなく、事後分布からのサンプル列であり、更新過程には確率的な要素が含まれる。
今回の設定では、すべての係数 $\beta_j$ が同じ事前分散を共有しているが、実際には係数ごとに異なる逆ガンマ分布を与えることもできる。その場合、階層モデルによって SNP ごとの効果分散を推定することになり、これがいわゆる BayesA に対応する。BayesB, Cでもやはり変数ごとに別々の事前分布を考える。
ただしこの全てで、実はハイパーパラメータ(事前分布のパラメータのパラメータ)は一位に定まるのでやはりある程度の決めうちは避けられない。(多分かなりマイナーな影響にはなると思う)、この辺は階層ベイズの設計でどこまで深くすれば良いのかという一般的な問題だが、どこまで行っても決め打ちするんだし二層くらいで大体いいだろうという感じがする(というか、階層の数より形のほうが重要な感じもする：BRRでinverseガンマと正規分布の掛け算をしてsigmaの周辺分布を考えると結局Bがt分布に従うのと同じ話になる。。)
3月卒業までにbayes A,B,Cの実装までやってみたい。
参考文献
PRML: Pattern Recognition and Machine Learning
https://www.microsoft.com/en-us/research/wp-content/uploads/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf
BGLR: Genome-Wide Regression and Prediction with the BGLR Statistical Package
https://academic.oup.com/genetics/article/198/2/483/5935898
OLS は 最尤推定、Ridge 回帰は MAP 推定、 ― ベイズの備忘録
https://qiita.com/yoyoyo11131113/items/e03831863f746d8ea0b7
同僚の書いた良質な記事として、こちらも参考に
BRRのギブスサンプリングによるパラメータ推定
https://qiita.com/Kazutoshi08/items/34807483be952c7279dd
線形混合モデルのギブスサンプリングによるパラメータ推定
https://qiita.com/KimuraSoh/items/f33d57859bd101cb487c