東京大学・株式会社Nopareの菅澤です．今回から縮小事前分布を使ったベイズ的変数選択の方法について，背景の原理やRでの実装について数回に分けて紹介していこうと思います．
今回は正則化法として有名なLassoのベイズ版であるBayesian Lassoについて紹介していきます．
線形回帰モデル
以下のような線形回帰モデルを考えます．
$$
y_i=x_{i}^T\beta+\varepsilon_i, \ \ \ \ i=1,\ldots,n.
$$
ここで$y_i$は被説明変数，$x_i=(x_{i1},\ldots,x_{ip})^T$は説明変数のベクトル，$\beta=(\beta_1,\ldots,\beta_p)^T$は回帰係数のベクトル，$\varepsilon_i$は誤差項を表します．以下では簡単のために，$\varepsilon_i$は独立かつ$\varepsilon_i\sim N(0,\sigma^2)$を仮定しておきます．
線形回帰モデルは統計学における基本的なモデルの一つですが，見た目以上に扱えるモデリングの範囲は広く，多項式回帰や基底関数を使ったノンパラメトリック回帰など，変数間の非線形的な関連まで表現することも可能です．
現実のデータ分析で線形回帰モデルを適用する際のよくある問題点として，説明変数の数$p$が大きいことが挙げられます．その背景としては
- 多様なデータが収集できる昨今では$y_i$に関連しそうな変数が多く得られる．
- 説明変数の非線形な構造を柔軟に表現するために多くの基底関数が必要である．
などが典型的なものとして考えられます．
その場合，通常の最小二乗推定で$\beta$を推定するとあまり精度の高い推定値が得られない（場合によっては存在しない）ことが知られています．主な原因としては，モデルの中に（ほとんど）影響を与えていない項が存在することで無駄なパラメータを推定することになり効率性が下がってしまうことが挙げられます．
そこで，精度の高い$\beta$の推定値を得るために無駄な説明変数を削除する（変数選択をする）ことが重要になってきます．
Lasso（Least Absolute Shrinkage and Selection Operator）
以下のようなベクトル表現を用いて議論を進めていきます．
$$
y= X\beta + \varepsilon, \ \ \ \varepsilon\sim N(0, \sigma^2 I_n)
$$
ここで$1_n$は1が並んだ$n$次元ベクトル，$X$は$n\times p$の計画行列です．また$y$, $\varepsilon$, $\beta$は$y_i$, $\varepsilon_i$, $\beta_k$をそれぞれ並べたベクトルです．
変数選択と$\beta$の推定を同時に実行する方法として以下のようなLassoと呼ばれる方法が非常に有名です．
\hat{\beta}={\rm argmin}_{\beta} \Big\{\|y- X\beta\|^2 + \lambda \sum_{k=1}^p |\beta_k| \Big\}
第2項に現れる係数の絶対値の罰則の効果で最小解のいくつかの要素がexactに0になるため，変数選択とパラメータ推定が同時に実行できます．また$\lambda$はチューニングパラメータで，大きい値に設定するほど0と推定される要素の数が増加（選ばれる変数の数が減少）します．通常はcross validationや情報量規準などを用いて適応的に選択されます．
Lassoのベイズ的解釈
Lassoで登場した目的関数に対してベイズ的な解釈を与えることができます．
線形回帰モデル$y\sim N(X\beta, \sigma^2I_n)$には$\beta$および$\sigma^2$のパラメータがありますが，ひとまず$\sigma^2$は固定し，$\beta$に対する事前分布として以下のようなラプラス分布を導入します．
\pi(\beta\ | \ \sigma^2)=\prod_{k=1}^p\frac{\lambda}{2\sigma}\exp\left(-\frac{\lambda|\beta_k|}{\sigma}\right)
ラプラス分布は原点上で尖った密度関数をもつ分布です．（密度関数の形については例えばWikipediaを参照してください．）ラプラス事前分布は原点付近に値が出やすい分布になっているため，「回帰係数$\beta$の多くは原点付近にあるだろう」という考えを反映した分布になっています．このような事前分布を一般に 縮小事前分布 と言います．
$\beta$の事後分布は以下のようになります．
\pi(\beta \ | \ y,\sigma^2)\propto \phi_n(y;X\beta, \sigma^2 I_n) \ \pi(\beta\ | \ \sigma^2) \\
\propto \exp\left\{-\frac{1}{2\sigma^2}(y-X\beta)^T(y-X\beta)-\frac{\lambda\sum_{k=1}^p|\beta_k|}{\sigma}\right\} \tag{1}
この表現から$\beta$の最大事後確率推定（MAP推定）は以下の関数の最小値と等価であることがわかります．
(y-X\beta)^T(y-X\beta)+\lambda^{\ast}\sum_{k=1}^p|\beta_k|, \ \ \ \ \lambda^{\ast}=2\sigma\lambda
$\lambda^{\ast}$を新しいチューニングパラメータとみなすと，結論として 「Lasso推定量は$\beta$にラプラス事前分布を導入した際のMAP推定とみなせる」 ことがわかります．
Bayesian Lasso
$\beta$にラプラス事前分布を導入し，MAPだけではなく事後分布全体を用いて$\beta$の統計的推測を行う方法を Bayesian Lasso と言います．特徴としては以下の3点が挙げられます．
- 信用区間を用いることで不確実性評価が容易に行える．
- $\lambda$に事前分布を導入することでチューニングパラメータの推定をMCMCに組み込める．
- 事後平均や事後中央値を点推定値として用いた場合はLasso推定量のようなexactな0の推定値は得られない．
事後分布に基づく統計的推測を行うために事後分布から乱数を生成することを考えます．しかし，(1)式で与えられた$\beta$の事後分布を見ると，あまり馴染みのない分布になっています．その原因はラプラス事前分布に由来する第2項の部分です．（この部分がなければ$\beta$の事後分布は多変量正規分布になります．）
もちろんMetropolis-Hastings法などの一般的なMCMCアルゴリズムを使って(1)式から直接$\beta$の乱数を生成することは可能ですが，（$\beta$の次元が高いときには特に）計算効率の観点から好ましいアプローチではありません．
そこで重要になってくるのがPark and Casella (2008)で提案された 階層表現によるパラメータ拡大法（parameter augmentation） です．
ラプラス分布の階層表現
一般にラプラス分布は以下のような正規分布のスケール混合で表現することができます．
\frac{a}{2}e^{-a|z|}=\int_0^{\infty}\phi(z;0,s)\cdot \frac{a^2}{2}e^{-a^2s/2} \ ds, \ \ \ a>0.
ここで$\phi(z;0,s)$は平均0, 分散$s$の正規分布の密度関数を表します．
この積分表現により，ラプラス分布に従う確率変数$Z$は以下のような階層表現を持つことがわかります．
Z\ |\ S\sim N(0, S), \ \ \ S\sim {\rm Exp}(a^2/2)
ただし，${\rm Exp}(\psi)$は平均$1/\psi$の指数分布を表します．
Bayesian Lassoの階層表現
ラプラス分布の階層表現を用いることで，Bayesian Lassoのモデルは以下のように表せます．
y \ | \ \beta, \sigma^2\sim N(X\beta, \sigma^2 I_n), \\
\beta_k \ | \ \sigma^2, u_k \sim N(0, \sigma^2 u_k), \ \ \ \
u_k \ | \ \lambda \sim {\rm Exp}(\lambda^2/2), \ \ \ k=1,\ldots,p \tag{2}
残るは$\sigma^2$および$\lambda$の事前分布ですが，今回は条件付き共役なものとして以下を用います．
\pi(\sigma^2)\propto \frac{1}{\sigma}, \ \ \ \ \ \lambda^2\sim {\rm Ga}(r, \delta)
すなわち，$\sigma^2$に対してはimproperな事前分布，$\lambda^2$については形状パラメータ$r$，レートパラメータ$\delta$のガンマ分布を設定します．
ギブスサンプリングによるMCMC
(2)式の階層表現の利点としてはギブスサンプリングによって簡単にMCMCが実行できる点です．まず全てのパラメータの同時事後分布は以下のような形になります．
\pi(\beta,\sigma^2,\lambda^2,u \ | \ y)\propto
\exp\left\{-\frac{1}{2\sigma^2}(y-X\beta)^T(y-X\beta)-\frac{1}{2\sigma^2}\sum_{k=1}^p\frac{\beta_k^2}{u_k}-\frac{\lambda^2}{2}\sum_{k=1}^pu_k\right\}\\
\times (\sigma^2)^{-(n+p-1)/2-1}\left\{\prod_{k=1}^p u_k^{-1/2}\right\}(\lambda^2)^{r+p}\exp(-\delta\lambda^2)
ギブスサンプリングを実行するためにそれぞれのパラメータの完全条件付き分布がどのような分布になるか調べる必要があります．例えば$\beta$の完全条件付き分布を求めたい場合，上の表現で$\beta$に依存する部分だけを抜き出してどんな分布になっているかを調べます．
最終的に各パラメータの完全条件付き分布が以下のように求まります．
-
$\beta$の完全条件付き分布: $N(A^{-1}X^Ty, \sigma^2 A^{-1})$．ただし$A=X^TX+D$, $D={\rm diag}(u_1^{-1},\ldots,u_p^{-1})$
-
$\sigma^2$の完全条件付き分布: ${\rm IG}(\ (n+p-1)/2, \ (y-X\beta)^T(y-X\beta)/2+\beta^TD\beta/2 \ )$．
-
$u_k^{-1}$の完全条件付き分布: ${\rm InvG}(\sqrt{\lambda^2\sigma^2/\beta_k^2}, \lambda^2)$．ただし${\rm InvG}(\mu, \eta)$は以下の密度関数を持つ逆正規分布です．
f(x; \mu, \eta)=\sqrt{\frac{\eta}{2\pi}}x^{-3/2}
\exp\left\{-\frac{\eta(x-\mu)^2}{2\mu^2x}\right\}, \ \ \ \ x>0.
- $\lambda^2$の完全条件付き分布: ${\rm Ga}(r+p,\delta+\sum_{k=1}^p u_k/2 )$．
このように各パラメータの完全条件付き分布は全て標準的な分布になり，容易にMCMCが実行できることがわかります．
ポイントとしては 「(1)式で与えられる$\beta$の事後分布からの直接的な乱数生成は容易ではないが，潜在変数$u_k$を用いた階層表現を導入することで簡単なギブスサンプリングを導出できて簡単に$\beta$の事後サンプルを生成できる」 ということです．一見難しく見える問題に対して潜在変数を導入することで簡単な細かい問題に分解しているイメージです．
ギブスサンプリングによって生成された$\beta$の事後サンプルを用いることで，事後平均や事後中央値などの点推定値だけでなく信用区間などの不確実性を考慮した量を計算することも可能です．
潜在変数$u_k$はギブスサンプリングを導出するために機械的に導入したもので，それ自体に関心があるパラメータではないですが，(2)式を見ると「$u_k$の値が大きいと$\beta_k$の事前分布の影響が小さく，逆に$u_k$の値が小さいと$\beta_k$の事前分布の分散が小さくなり0方向に縮小された事後分布が導出されやすい」といった解釈を与えることができます．
最後に
今回はベイズ的変数選択の手法としてBayesian Lassoについて紹介しました．具体的な数値例や実装については次回以降で紹介していきます．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．データ利活用支援や統計アドバイザリー・ビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
