Bayesian スパースモデリング

• Idea
• 機械学習

Bayes 統計に基づいた線形回帰における変数選択について、その理論と Python による簡単な実装をまとめます。

## LASSO の Bayes 的解釈

こちらの記事で示したように、LASSO はラプラス事前分布を課した回帰係数の MAP 推定と見なすことができます。

まず、これを示します。ここでは、データ $\mathcal{D} = (X, \bm{y})$：

に対する線形回帰モデル

を考えます。線形回帰モデルを構成する特徴量 $\bm{x}_j$ について、一部の係数はゼロ ($\beta_j = 0$) であるとみなす仮定をスパース性の仮定といい、スパース性を仮定した統計モデリングを総称してスパースモデリングといいます。

### Laplace 分布

Laplace 分布 $\mathrm{Laplace}(\mu, b)$ は、以下の確率密度関数：

で与えられる確率分布で、位置パラメータ $\mu$ と尺度パラメータ $b$ を持ちます。
$\mathrm{Laplace}(0, 1)$ の可視化は以下のようになります。正規分布と比べて裾が重く、ゼロの位置に不連続なピークを持つことがわかります。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import laplace

x = np.linspace(-6, 6, 1024)

# Laplace(0, 1)
laplace_pdf = laplace.pdf(x, loc=0, scale=1)

plt.figure(figsize=(8, 6))
plt.plot(x, laplace_pdf, label="Laplace")
plt.xlim(-6, 6)
plt.ylim(0, None)
plt.xlabel(r"$x$")
plt.ylabel("Propability density")
plt.legend()
plt.tight_layout()
plt.show()

### MAP 推定量としての LASSO

回帰係数 $\bm{\beta} = [\beta_1, \ldots, \beta_p]^\top$ の事前分布として、i.i.d な Laplace 分布：

を課します。確率密度関数とその対数は、

です。ここで、線形回帰モデルの対数尤度関数は

ですから、LASSO 推定量 $\hat{\bm{\beta}}_\mathrm{LASSO}$ は適当な $\lambda > 0$ を用いた MAP 推定量として、

と表すことができます。ここで

を $L_1$-ノルムといいます。

### Bayesian LASSO

LASSO が推定する回帰係数は MAP 推定量、すなわち事後確率が最大となるときのパラメータの値であるため、事後分布 $P(\bm{\beta}|X)$ の点推定 (最頻値) となります。
これに対し、$P(\bm{\beta}|X)$ そのものを推定する方法は Bayesian LASSO と呼ばれます。この場合の $P(\bm{\beta}|X)$ は通常解析的でないため、Markov 連鎖モンテカルロ (MCMC) をはじめとするサンプリング法によって推定されます。

### Bayesian LASSO の階層表現

Laplace 分布は $\mathrm{Laplace}(0, b)$ は指数分布と正規分布の混合によって表すことができます。

これを (やんわり) 示します。Laplace 分布の密度関数 $p(\beta_j)$ を指数分布のパラメータ $\tau$ で拡大すると、

となります。正規分布と指数分布の密度関数を明示すると、

であるので、Gauss 積分公式を使用して

となります。

指数分布と正規分布からのサンプルで以上の階層表現を構成してみると、確かに $\mathrm{Laplace}(0, 1)$ に近づくことがわかります。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import laplace

b = 1.0
N = 200000

# τ ~ Exp(1/(2b²))
tau = np.random.exponential(scale=2 * b**2, size=N)

# beta | τ ~ N(0,τ)
beta= np.random.normal(loc=0, scale=np.sqrt(tau))

x = np.linspace(-8, 8, 1000)

plt.figure(figsize=(8, 6))
plt.hist(beta, bins=200, density=True, alpha=0.6, label="mixture samples")

plt.plot(x, laplace.pdf(x, scale=b), linewidth=3, label="Laplace(0,1)", color="dodgerblue")

plt.legend()
plt.tight_layout()
plt.show()

この事実を利用すれば、Bayesian LASSO を以下のように階層表現することができます。

実際に用いる場合は、$\sigma^2, \tau$ にも適当な事前分布をおきます (共役事前分布を用いる場合が多いのでしょうか？)。

## 縮小事前分布

次の問題設定を考えてみます。

そのような $\beta_j ~ (j = 1,\ldots, p)$ の事前分布として、以下の形式の階層表現を導入します。

$P(\cdot)$ は適当な確率分布で、たとえば $\mathrm{Exp}$ とすれば Bayesian LASSO に相当します。この形式で表現できる事前分布は線形回帰モデルに $\beta_j \simeq 0$ となるような事前知識を与えます。そのような理由から、縮小事前分布と呼ばれるようです。

なお、$\sigma^2, \tau, \lambda_j$ が所与である場合、事後分布は正規分布になります (こちらの記事↗を参照)。

### 縮小係数

$\sigma^2 = 1$ とし、$\tau, \lambda_j$ が所与のもとで $\beta_j$ の事後平均は以下のようになります：

ここで、

を縮小係数などと呼び、$\tau^2\lambda_j^2 > 0$ であることから $0 < \kappa_j < 1$ の範囲を取ります。事後平均を縮小係数で書き直すと、

となることから、縮小係数 $\kappa_j$ はシグナル $y_j$ をどれくらい圧縮するかを決めるパラメータであることがわかります。また、

• $\tau$ は $\kappa_1, \ldots, \kappa_p$ 全てに寄与するパラメータ
• $\lambda_j$ は $\kappa_j$ のみに寄与するパラメータ

であることから、$\tau$ を大域縮小パラメータ、$\lambda_j$ を局所縮小パラメータなどと呼びます。

### Laplace 分布と縮小係数

Laplace 分布の問題点は、任意の回帰係数 $\beta_j$ で縮小が起きやすくなってしまう点にあります。この問題が起きる理由は、Laplace 分布に対応する $\kappa_j$ の事前分布を見ることで理解できます。

Laplace 分布の階層表現における局所縮小パラメータ $\lambda_j$ に対して、

とすると縮小係数を定義できます。$\lambda_j \sim \mathrm{Exp}(\tau^2/2)$ であるので、その確率密度関数は

と表すことができます。変数変換 $\lambda_j \to \kappa_j$ に伴う Jacobian は

であるので、

を得ます。$\tau = 1$ としてこれをプロットしたのが以下の図です。

import numpy as np
import matplotlib.pyplot as plt

# kappa
kappa = np.linspace(0, 1, 1024)

# p(kappa)
prob = np.exp((kappa - 1)/(2 * kappa)) / (2 * kappa**2)

plt.figure(figsize=(8, 6))
plt.plot(kappa, prob, linewidth=3, label=r"$p(\kappa)$", color="dodgerblue")
plt.legend()
plt.tight_layout()
plt.show()

図から明らかなように、Laplace 分布は $\kappa_j \simeq 0$ の近辺に密度をほとんど持たず、それゆえ多くの係数がゼロに縮退してしまいます。

### 馬蹄分布

強いスパース性を担保しつつも過剰縮小を抑えるためには、$\kappa_j \simeq 0$ および $\kappa_j \simeq 1$ に密度を集中させれば良いことがわかります。$\kappa_j$ の事前分布としてベータ分布 $\mathrm{Beta}(1/2, 1/2)$ をとると、そのような性質を満足することができます。確率密度関数は以下で与えられます：

これをプロットしたものは以下のようになります。

import numpy as np
import matplotlib.pyplot as plt

# kappa
kappa = np.linspace(0, 1, 1024)

# p(kappa)
prob = (1 / np.pi) * (1 / np.sqrt(kappa * (1 - kappa)))

plt.figure(figsize=(8, 6))
plt.plot(kappa, prob, linewidth=3, label=r"$p(\kappa)$", color="dodgerblue")
plt.legend()
plt.tight_layout()
plt.show()

関数形が馬の蹄 (ひづめ) のように見えることから、これに対応する $\beta_j$ の事前分布は馬蹄分布と呼ばれます。局所縮小パラメータ $\lambda_j$ と大域縮小パラメータ $\tau$ をもつ馬蹄分布は以下の階層表現で表すことができます。

$\mathrm{HalfCauchy}(0, 1)$ は位置パラメータ $0$、尺度パラメータ $1$ の (標準) 半 Cauchy 分布で、確率密度関数は以下で与えられます：

これを利用して、馬蹄事前分布をおいた Bayes 線形回帰を以下のような階層表現で定式化できます。

### 馬蹄分布の逆ガンマ表現

馬蹄分布のデメリットとして、半 Cauchy 事前分布の性質上、解析的な事後分布を得にくい点が挙げられます。この問題に対処するため、解析的な事後分布を得やすい逆ガンマ分布 $\mathrm{InvGamma}(a, b)$ で半 Cauchy 分布を書き直すことができます。

形状パラメータ $a$、尺度パラメータ $b$ を持つ逆ガンマ関数の確率密度関数は以下です：

なお、$\Gamma(t)$ はガンマ関数です。逆ガンマ分布を用いた半 Cauchy 分布の表現は以下のようになります。

これを用いると、$\tau, \lambda_j$ にそれぞれ対応する補助変数 $\xi, \nu_j$ を用意して、

という階層表現を得ることができます。$\sigma^2$ の事前分布も逆ガンマ分布を使えば、馬蹄分布をすべて逆ガンマ分布で表現できたことになります。正規分布の分散パラメータに対応する共役事前分布は逆ガンマ分布ですから、このように構成した馬蹄分布は事後分布の一部もまた逆ガンマ分布で表現できます。

## References

### 2次情報

• 縮小事前分布、特に馬蹄分布に関する導入を参考にしました。

• 続き物の記事です。大変参考になりました。

縮小事前分布によるベイズ的変数選択1: Bayesian Lasso - Qiita

東京大学・株式会社Nopareの菅澤です．今回から縮小事前分布を使ったベイズ的変数選択の方法について，背景の原理やRでの実装について数回に分けて紹介していこうと思います． 今回は正則化法として有名なLassoのベイズ版であるBayesian Lassoについて紹介していきま...

https://qiita.com

また、サンプルコードや数学的補足は一部 ChatGPT から生成しました。

### 論文等

1. T. Park and G. Casella, The Bayesian Lasso, J. Am. Stat. Assoc., 103, 681–686 (2008).
2. C. M. Carvalho et al., Handling Sparsity via the Horseshoe, Proc. Mach. Learn. Res. (PMLR), 5, 73–80 (2009).
3. C. M. Carvalho, J. G. Scott. The horseshoe estimator for sparse signals, Biometrika, 97, 465–480 (2010).
4. E. Makalic and D. F. Schmidt, A simple sampler for the horseshoe estimator, IEEE Signal Process. Lett., 23, 179–182 (2016).

Discussions

Related Articles