線形回帰の解がスパースである(0が多い)ことが事前にわかっているとき、 ベイズ線形回帰(Bayesian Linear Regression) において重みの事前分布に対して、 馬蹄分布(Horseshoe Distribution) を用いるのが一般的です。
正規分布などを事前分布に用いたベイズ線形回帰と比べて、馬蹄分布を事前分布に用いたベイズ線形回帰では、事後分布からのサンプリングには比較的複雑なコードが必要になります。
この記事では、 マルコフ連鎖モンテカルロ法(Markov chain Monte Calro method) のアルゴリズムの一つである ギブスサンプリング(Gibbs Sampling) を用いた、最も基本的な馬蹄事後分布からのサンプリング方法を紹介します。
ベイズ線形回帰とスパース性
次のような線形回帰モデルにおいて、
\mathbf{y}=\mathbf{X}\boldsymbol \beta + \boldsymbol \epsilon
\boldsymbol \epsilon\sim\mathbb{N}(\mathbf{0},\sigma^2\mathbf{I})
最も一般的に設定される事前分布は正規分布です。
\boldsymbol \beta \sim \mathbb{N}(\boldsymbol \mu,\Sigma)
このとき、 $\mathbf{y}$ の分布は
\mathbf{y}|\mathbf{X},\boldsymbol \beta \sim \mathbb{N}(\mathbf{X}\boldsymbol \beta,\sigma^2\mathbf{I})
と表すことができるので、ベイズの定理を用いることで比較的簡単に $\boldsymbol \beta$ の事後分布を求めることができます。
p(\boldsymbol \beta|\mathbf{y},\mathbf{X})=\frac{p(\mathbf{y}|\mathbf{X},\boldsymbol \beta)p(\boldsymbol \beta)}{\int p(\mathbf{y}|\mathbf{X},\boldsymbol \beta)p(\boldsymbol \beta) d\boldsymbol \beta}
ここで、線形回帰の各係数 $\beta_1,\beta_2,\cdots,\beta_n$ がスパースであるという情報をモデルの設計者が持っている場合を考えます。
y_i = \beta_1 x_{i1}+\beta_2 x_{i2} + \cdots +\beta_n x_{in} + \epsilon_i
例えば高次元の問題の解は一般的にスパースとなる傾向にあることが知られています。
このとき、 $\beta_1,\beta_2,\cdots,\beta_n$ にスパース性を考慮した事前分布を設定すると、事前知識をうまくモデルに組み込むことができます。
スパース性を誘導する事前分布として、 馬蹄事前分布(Horseshoe Prior) が知られています。
馬蹄分布(Horseshoe Distribution)
馬蹄分布は次のような数式で記述されます。
\beta_i|\lambda_i,\tau\sim\mathbb{N}(0,\lambda_i^2\tau^2)
\lambda_i\sim C^+(0,1)
$C^+(0,1)$ は半コーシー分布です。半コーシー分布は次のような数式で記述されます。馬蹄分布は、$\lambda_i$ が半コーシー分布から生成されているとき、平均 $0$ で分散 $\lambda_i^2\tau^2$ のガウス分布から生成されます。
\begin{split} \begin{align}
f(x;a, b) = \left\{\begin{array}{cll}
\frac{2}{\pi b}\,\frac{1}{1 + (x-a)^2/b^2} & & x \ge a \\[1em]
0 & & \text{otherwise}.
\end{array}\right.
\end{align}\end{split}
馬蹄分布から10000個サンプリングした結果を次図に示します。馬蹄分布は分布がスパースであるときに用いられる分布なので、サンプリングの結果もほとんど0となっています。
いくつかのサンプルは0以外の値となっていますが、0のサンプル数と比べて極端に少ないためヒストグラムでは確認できません。
別な記事で馬蹄分布の基本的な事柄についても解説したので、よろしければ参照してみてください。
ギブスサンプリング(Gibbs Sampling)
確率分布$P(\mathbf{Z})$から直接$\mathbf{Z}$全体をサンプリングすることが難しい場合、$\mathbf{Z}=(\mathbf{Z}_1,\mathbf{Z}_2, \dots ,\mathbf{Z}_M)$のようにM個の部分集合に分けることによって逐次的にサンプリングを行います。
この方法を ギブスサンプリング(Gibbs Sampling) と呼びます。
\mathbf{Z}_1 \sim p(\mathbf{Z}_1|\mathbf{Z}_2,\mathbf{Z}_3, \dots ,\mathbf{Z}_{M-1}, \mathbf{Z}_M) \\
\mathbf{Z}_2 \sim p(\mathbf{Z}_2|\mathbf{Z}_1,\mathbf{Z}_3, \dots ,\mathbf{Z}_{M-1}, \mathbf{Z}_M) \\
\vdots \\
\mathbf{Z}_M \sim p(\mathbf{Z}_M|\mathbf{Z}_1,\mathbf{Z}_2, \dots ,\mathbf{Z}_{M-2}, \mathbf{Z}_{M-1})
ギブスサンプリングはサンプルを得たい変数の数が膨大である場合や、複数の確率分布が組み合わさった巨大な確率モデル化からサンプルを得たい場合に用いられます。
上式の各条件付き分布が解析的に計算できるように共役事前分布をうまく用いてモデルを構築できれば非常に効率的に動作します。
別な記事でMCMCの基本的なアルゴリズムについて、解説や実験を行ったので、よろしければ参照してみてください。
回帰係数の推定
線形回帰モデルにおいて、重み $\boldsymbol \beta$ の事前分布に馬蹄分布を設定します。
\mathbf{y}=\mathbf{X}\boldsymbol \beta + \boldsymbol \epsilon
\boldsymbol \epsilon\sim\mathbb{N}(\mathbf{0},\sigma^2\mathbf{I})
\beta_i|\lambda_i^2,\tau^2,\sigma^2\sim\mathbb{N}(0,\lambda_i^2\tau^2\sigma^2)
\boldsymbol \beta\in\mathbb{R}^p
このとき $\boldsymbol \beta$ の事後分布は次のような条件付き分布として計算することができます。詳細な導出は省略します。気になる人はこの論文を参照してください。
\boldsymbol \beta|\cdot\sim\mathbb{N}(\mathbf{A}^{-1}\mathbf{X}^{\top}\mathbf{y},\sigma^2\mathbf{A}^{-1}),\quad \mathbf{A}=(\mathbf{X}^{\top}\mathbf{X}^{-1}+\mathbf{\Lambda_{*}}^{-1}),\quad \mathbf{\Lambda_{*}}=\tau^2\mathbf{\Lambda}
\mathbf{\Lambda}=\textrm{diag}(\lambda_1^2,\dots,\lambda_p^2)
同様に $\sigma^2,\lambda_i^2,\tau^2$ の条件付き事後分布も次のように計算することができます。
\sigma^2|\cdot\sim\textrm{IG}(\frac{(n+p)}{2},\frac{(\mathbf{y}-\mathbf{X}\boldsymbol\beta)^{\top}(\mathbf{y}-\mathbf{X}\boldsymbol\beta)}{2} + \frac{\boldsymbol\beta^{\top}\mathbf{\Lambda_{*}^{-1}\boldsymbol\beta}}{2})
\lambda_i^2|\cdot\sim\textrm{IG}(1,\frac{1}{\nu_i}+\frac{\beta_i^2}{2\tau^2\sigma^2})
\tau_i^2|\cdot\sim\textrm{IG}(\frac{p+1}{2},\frac{1}{\xi}+\frac{1}{2\sigma^2}\sum_{i=1}^p \frac{\beta_i^2}{\lambda_i^2})
上記の条件付き事後分布には、元々はなかった補助変数 $\nu_i,\xi$ が使われています。これらの事後分布は次のように求められます。
\nu_i|\cdot\sim\textrm{IG}(1,1+\frac{1}{\lambda_i^2})
\xi|\cdot\sim\textrm{IG}(1,1+\frac{1}{\tau^2})
興味深いことに回帰係数以外の全てのパラメータは逆ガンマ分布の形で記述することができました。
これらの数式をプログラムに落とし込むと次のようになります。基本的に数式をそのままプログラムに変換しただけで、難しい処理等は必要ありません。
def bhs(X, y, n_samples=1000, burnin=200):
n, p = X.shape
XtX = X.T @ X
beta = np.zeros((p, n_samples))
sigma2 = 1
lambda2 = np.random.uniform(size=p)
tau2 = 1
nu = np.ones(p)
xi = 1
# Run Gibbs Sampler
for i in range(n_samples + burnin):
Lambda_star = tau2 * np.diag(lambda2)
A = XtX + np.linalg.inv(Lambda_star)
A_inv = np.linalg.inv(A)
b = np.random.multivariate_normal(A_inv @ X.T @ y, sigma2 * A_inv)
# Sample sigma^2
e = y - np.dot(X, b)
shape = (n + p) / 2.
scale = np.dot(e.T, e) / 2. + np.sum(b**2 / lambda2) / tau2 / 2.
sigma2 = 1. / np.random.gamma(shape, 1. / scale)
# Sample lambda^2
scale = 1. / nu + b**2. / 2. / tau2 / sigma2
lambda2 = 1. / np.random.exponential(1. / scale)
# Sample tau^2
shape = (p + 1.) / 2.
scale = 1. / xi + np.sum(b**2. / lambda2) / 2. / sigma2
tau2 = 1. / np.random.gamma(shape, 1. / scale)
# Sample nu
scale = 1. + 1. / lambda2
nu = 1. / np.random.exponential(1. / scale)
# Sample xi
scale = 1. + 1. / tau2
xi = 1. / np.random.exponential(1. / scale)
if i >= burnin:
beta[:, i - burnin] = b
return beta
人口データを用いてスパース性を導入したベイズ線形回帰を実行してみます。
回帰係数は20個でサンプルデータ数は20個で回帰を実行してみます。また、入力は平均ゼロで分散が単位行列の正規分布から生成されるものとして、ギブスサンプリングに回帰係数を求めます。
今回の実験では、係数がスパースであるという状況に合わせるために真の回帰係数のうち $\beta_1,\beta_2,\beta_3$ を0以外の値をとるものとし、他の17個の係数は全て0であるという状況を考えます。
if __name__ == '__main__':
n = 20
p = 20
# true coefficients
coefs = np.zeros((p,))
coefs[0] = 1.0
coefs[1] = 1.5
coefs[2] = 0.5
X = np.random.multivariate_normal(np.zeros((p,)), np.eye(p), size=n)
y = X @ coefs
beta = bhs(X, y)
# traceplot
fig, axes = plt.subplots(3, 3, figsize=(24, 24))
for i in range(3):
for j in range(3):
index = 3 * i + j
if i == 0:
axes[i][j].set_ylim(-5, 5)
axes[i][j].set_title("beta_{}".format(index + 1))
axes[i][j].plot(beta[index, :])
fig.tight_layout()
fig.savefig('bhs-tp.png')
plt.close()
# mean, variance
fig = plt.figure(figsize=(16, 8))
plt.xlabel('coefficients')
plt.ylabel('value')
for i in range(p):
mean = np.mean(beta[i, :])
var = np.var(beta[i, :])
std = np.sqrt(np.abs(var))
plt.errorbar(i, mean, yerr=std * 4, capsize=5, fmt='o', markersize=10,)
fig.savefig('bhs-ci.png')
plt.close()
このプログラムを実行すると次のような結果が得られます。
下の図は $\beta_1\sim\beta_9$ のトレースプロットです。今回の実験では、 $\beta_1,\beta_2,\beta_3$ に対して0以外の値を設定しました。トレースプロットによると $\beta_1,\beta_2,\beta_3$ の値の推定には成功していることが見て取れます。
また、回帰係数の事前分布に馬蹄分布を用いたため、 $\beta_4\sim\beta_9$ の
サンプリングの結果は0に近い値を取ることが多くなっています。
下の図は各回帰係数の平均とその信頼区間を描画したものです。エラーバーも描画したのですが、各回帰係数の分散は非常に小さい値となっていたため、エラーバーは視認できませんでした。
まとめ
ベイズ線形回帰の回帰係数の事前分布に馬蹄分布を設定し、ギブスサンプリングを用いて回帰係数の推定を行いました。
実験では人口データを用意し、回帰係数に対して馬蹄事前分布を設定したベイズ線形回帰を行いました。。問題が簡単すぎたためか、非常に良い結果が得られました。
今回の実験コードはGitHub上にあります。
参考
[1] CARLOS M. CARVALHO, O, NICHOLAS G. POLSON. The horseshoe estimator for sparse signal. 2010.
[2] Carlos M. Carvalho, Nicholas G. Polson, James G. Scott. Handling Sparsity via the Horseshoe. 2009.
[3] Enes Makalic, Daniel F. Schmidt. A simple sampler for the horseshoe estimator. 2015.
[4] Ricardo Baptista, Matthias Poloczek. Bayesian Optimization of Combinatorial Structures. 2018.