正則化と縮小事前分布
回帰分析で回帰係数の推定値を求める場合、最小二乗推定値はよく使われています。
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},
\qquad
\boldsymbol{\varepsilon} \sim \mathcal{N}\!\left(\mathbf{0},\,\sigma^2\mathbf{I}\right)
\widehat{\boldsymbol{\beta}}_{\text{OLS}}
= \left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{y}
\tag{1}
この最小二乗推定値は、$\left(\mathbf{X}^{\top}\mathbf{X}\right)$は正則行列（可逆行列）であることを前提にして計算されています。
しかし、データの次元数pがデータ数nより多い場合、または変数同士に強い多重共線性がある場合、この行列は正則とは限らない。
もっと具体的に言うと、1個の従属変数に対し、説明変数の候補が90個もある一方、サンプルが60人しかない場合、説明変数の数はデータ数よりも多いため、$\left(\mathbf{X}^{\top}\mathbf{X}\right)$はフルランク行列ではなくなり、逆行列を求めることはできない。
このような状況で推定を安定化し、過学習を抑えて予測性能を 向上させるために用いられるのが 正則化（regularization） です。
頻度論の文脈では、誤差関数に正則化項を加えることで、正則化を実現する。代表的正則化手法として、Lasso回帰（L1正則化）とリッジ回帰（L2正則化）が知られている。
Lasso回帰の誤差関数
\min_{\boldsymbol{\beta}} \;\bigl\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\|^{2} \;+\;\lambda\,\,\|\boldsymbol{\beta}\|_{2}^{2}, \quad \lambda \ge 0.
リッジ回帰の誤差関数
\min_{\boldsymbol{\beta}} \;\bigl\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\bigr\|^{2} \;+\;\lambda\,\|\boldsymbol{\beta}\|_{1}, \quad \lambda \ge 0.
ベイズ統計学の文脈では、正則化項の代わりに、縮小事前分布（Shrinkage Prior）を使います。代表的な縮小事前分布は、Bayesian Lasso、馬蹄事前分布（Horseshoe Prior）、Spike and Slab Priorなどがあります。本記事は、馬蹄事前分布の設定、事後分布の導出およびStanでの実装について解説します。
馬蹄事前分布（Horseshoe Prior）
馬蹄事前分布は、Carvalho et al(2010)で提案された縮小事前分布です。回帰分析の文脈において、回帰係数の馬蹄事前分布の設定方法はこちら：
\beta_j \mid \lambda_j,\tau,\sigma \;\sim\; \mathcal{N}\!\left(0,\; \tau^{2}\lambda_j^{2}\right),\quad j=1,\dots,p
行列表現 \quad \boldsymbol{\beta}\mid\boldsymbol{\lambda},\tau,\sigma \sim \mathcal{N}\!\left(\mathbf{0},\; \tau^{2}\,\mathrm{diag}(\lambda_1^{2},\dots,\lambda_p^{2})\right)
\lambda_j \;\sim\; \mathcal{C}^{+}(0,1),
ここでの$\mathcal{C}^{+}(0,1)$は半コーシー分布を指します。
注意すべきなのは、ここに切片$\beta_0$が含まれていない。一般的に、切片を0まで縮小させる必要はないからです。切片$\beta_0$に対して、flatな事前分布（弱情報事前分布）を別途に設定します。
事後分布の導出
Piironen& Vehtari(2017)では、馬蹄事前分布を使う場合の回帰係数ベクトル$\boldsymbol{\beta}$の事後分布を以下のように示した。
p(\boldsymbol{\beta}\mid\boldsymbol{\lambda},\tau,\sigma^{2},\mathbf{X},\mathbf{y}) = N(\boldsymbol{\beta}\mid\bar{\boldsymbol{\beta}},\boldsymbol{\Sigma})
\bar{\boldsymbol{\beta}} = \tau^{2}\boldsymbol{\Lambda}\ (\tau^{2}\boldsymbol{\Lambda} + \sigma^{2}\ \left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1})^{-1}\widehat{\boldsymbol{\beta}_{\text{OLS}}}
\boldsymbol{\Sigma} = (\tau^{-2}\boldsymbol{\Lambda}^{-1} + \frac{1}{\sigma^{2}}\left(\mathbf{X}^{\top}\mathbf{X}\right))^{-1}
事後分布の導出自体は一般的な回帰モデルの事後分布の導出とほぼ同じです。具体的な計算はピーター・D・ホフの「標準ベイズ統計学」のp171-p173を参照されたい。しかし、$\boldsymbol{\beta}$の事後平均の導出にはいくつかの線形代数の恒等式を使う必要があるため、ここで記します。
(C^{-1} + B)^{-1}B = (I + CB)^{-1}CB
(I + CB)^{-1}C = C(I + BC)^{-1}
C + \sigma^{2}A^{-1} = \sigma^{2}A^{-1}(I + \sigma^{2}AC)
逆行列をとれば
(I + \sigma^{2}AC)^{-1}\sigma^{-2}A = (C + \sigma^{2}A^{-1})^{-1}
さらに、$\tau$を調整することでパラメータへの縮小効果も変化する。$\tau$を大きくすると 事前分散が大きく なり、縮小効果は弱くなります。反対に$\tau$が0に近づくと、回帰係数の事前分布の分散も0に近づき、最終的にすべての回帰係数が0まで縮小される。
Stanによる実装
#必要なライブラリを導入
library(tidyverse)
library(dplyr)
library(tidyverse)
library(cmdstanr)
library(bayesplot)
library(bridgesampling)
#重回帰分析のシミュレーションデータを生成
set.seed(123)
n<- 500
beta0 <- 2.3 #切片
beta1 <- 3.1
beta2 <- 1.7
sigma <- 2 #残差の標準偏差
#説明変数の生成
x1 <- runif(n,-1,1)
x2 <- runif(n,-1.5,1.5)
#残差の生成
e <- rnorm(n,0,sigma)
#従属変数の生成
y <- beta0 + beta1 * x1 + beta2 * x2 + e
#関係のない説明変数を生成する
x3 <- runif(n,-1,1)
x4 <- runif(n,-1,1)
#dataframeとしてまとめる
df1<- data.frame(y = y, x1 = x1, x2 = x2,x3 = x3, x4 = x4)
#Stanによる実装
data {
int<lower=0> N;
array[N] real y;
array[N] real x1;
array[N] real x2;
array[N] real x3;
array[N] real x4;
real<lower=0> tau; #tauはデータリストから事前に渡す
}
parameters {
real b0;
real b1;
real b2;
real b3;
real b4;
real<lower=0> sigma;
real<lower=0> lambda1;
real<lower=0> lambda2;
real<lower=0> lambda3;
real<lower=0> lambda4;
}
transformed parameters{
array[N] real mu;
for(n in 1:N){
mu[n] = b0 + b1*x1[n] + b2*x2[n] + b3*x3[n] + b4*x4[n];
}
}
model {
for(n in 1:N){
y[n] ~ normal(mu[n], sigma);
}
sigma ~ exponential(2);
b0 ~ normal(0,10); #b0の事前分布を別途に設置
b1 ~ normal(0, tau * lambda1);
b2 ~ normal(0, tau * lambda2);
b3 ~ normal(0, tau * lambda3);
b4 ~ normal(0, tau * lambda4);
lambda1 ~ cauchy(0,1);
lambda2 ~ cauchy(0,1);
lambda3 ~ cauchy(0,1);
lambda4 ~ cauchy(0,1);
}
#データリストを作成する
N<- nrow(df1)
tau = 10 #tauを10に設定
data_horseshoe <- list(
N = N,
y = df1$y,
x1 = df1$x1,
x2 = df1$x2,
x3 = df1$x3,
x4 = df1$x4,
tau = tau
)
#Stanコードをcmdstan_modelに回す
model_horseshoe<-cmdstan_model("stanmodel/horseshoe_prior.stan",force_recompile = TRUE)
#MCMCを実行
chains = 4
fit_horseshoe<- model_horseshoe$sample(data = data_horseshoe, chains = chains,
parallel_chains = chains,
iter_sampling = 5000,
iter_warmup = 5000,
refresh = 200,
output_dir = "cmdstanr/horseshoe_prior")
#結果の確認
fit_horseshoe$summary(c("b0", "b1", "b2", "b3", "b4"))
fit_horseshoe$summary(c("sigma"))
結果は以下のような感じです。
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b0 2.30 2.30 0.0901 0.0900 2.15 2.45 1.00 15956. 12322.
2 b1 2.95 2.96 0.160 0.159 2.69 3.22 1.00 12141. 11764.
3 b2 1.64 1.64 0.107 0.105 1.47 1.82 1.00 4889. 2465.
4 b3 0.227 0.225 0.154 0.161 -0.0132 0.482 1.00 7053. 12735.
5 b4 0.0216 0.0138 0.142 0.130 -0.209 0.264 1.00 15859. 11455.
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 sigma 2.02 2.02 0.0633 0.0633 1.92 2.13 1.00 14543. 11566.
多少違いはあるものの、b0, b1, b2の結果は事前に設定された数値と近い結果になりましたね。
一方で b3・b4 の事後分布は 0 付近に強く縮小され、95% 信用区間はいずれも 0 を含むことが確認できました。
考察
Stan 実装を通して、馬蹄事前分布（Horseshoe prior）の持つ縮小効果を確認できました。
また、縮小事前分布は回帰に限りません。たとえば因子分析における交差負荷（cross-loading)の推定を安定化する目的でも有効で、Horseshoe や Regularized Horseshoe、Bayesian Lasso などが用いられます。今後はこの応用についても取り上げたいと思います。
引用文献
Carvalho, C. M., Polson, N. G., & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 465-480.
Piironen, J., & Vehtari, A. (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors.
日本統計学会.統計学実践ワークブック
ピーター・D・ホフ[著], 入江 薫,菅澤 翔之助, & 橋本 真太郎[訳] (2022). 標準ベイズ統計学
小杉 考司, 紀ノ定 保礼, 清水裕士 (2023). 数値シミュレーションで読み解く統計の仕組み