Bayesian線形回帰のGibbs Sampler - おめさんの備忘録

これは何

Bayesian線形回帰のGibbs sampler.
よくわすれちゃうから

Probabilistic Model

$\boldsymbol{y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}, ~~~ \boldsymbol{\epsilon} \sim \mathcal{N}(\boldsymbol{0}_N, \boldsymbol{\Sigma}_{\epsilon}); ~~~ i=1,\dots,N.$

なお、
$\boldsymbol{y}=(y_1,\dots,y_N)^{\prime},$
$\boldsymbol{X}:\text{design matrix},$
$\boldsymbol{\beta}=(\beta_0, \beta_1, \dots, \beta_k)^{\prime},$
$\boldsymbol{\Sigma}_{\epsilon}=\sigma^2_{\epsilon}\boldsymbol{I}_N$ だけど必ずしもそうでなくてよい。

Priorは
$\begin{align*} p(\boldsymbol{\beta})&=(2\pi)^{-N/2}\abs{\boldsymbol{V}_{\beta}}^{-1/2}\exp{\set{-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{M}_{\beta})^{\prime}\boldsymbol{V}_{\beta}^{-1}(\boldsymbol{\beta}-\boldsymbol{M}_{\beta})}}=_d\mathcal{N}(\boldsymbol{M}_{\beta},\boldsymbol{V}_{\beta}), \\ p(\sigma^2_{\epsilon})&=\frac{(\beta_{\sigma})^{\alpha_{\sigma}}}{\Gamma(\alpha_{\sigma})}\bracket{\sigma^2_{\epsilon}}^{-(\alpha_{\sigma}+1)}\exp{\set{-\frac{\beta_{\sigma}}{\sigma^2_{\epsilon}}}}=_d\text{InvGamma}(\alpha_{\sigma},\beta_{\sigma}). \end{align*}$

Likelihoodは、 $\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{X\beta},\sigma^2_{\epsilon}\boldsymbol{I}_N)$ より、

$\begin{align*} p\Bracket{\boldsymbol{y}}{\boldsymbol{\beta},\sigma^2_{\epsilon}} &= \abs{2\pi\boldsymbol{\Sigma}_{\epsilon}}^{-1/2}\exp{\set{-\frac{1}{2}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}\boldsymbol{\Sigma}_{\epsilon}^{-1}(\boldsymbol{y}-\boldsymbol{X\beta})}} \\ &= (2\pi\sigma^2_{\epsilon})^{-N/2}\exp{\set{-\frac{1}{2\sigma^2_{\epsilon}}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{y}-\boldsymbol{X\beta})}}. \end{align*}$

Full Conditionals

$\begin{align*} p\Bracket{\sigma^2_{\epsilon}}{\boldsymbol{y},\boldsymbol{\beta}} & \propto p(\boldsymbol{\beta})P(\sigma^2_{\epsilon})p\Bracket{\boldsymbol{y}}{\boldsymbol{\beta},\sigma^2_{\epsilon}} \\ & \propto (\sigma^2_{\epsilon})^{-(N/2+\alpha_{\sigma}+1)}\exp{\set{-\frac{1}{\sigma^2_{\epsilon}}\bracket{\frac{1}{2}(\boldsymbol{y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{y}-\boldsymbol{X\beta})+\beta_{\sigma}}}} \\ &= \text{InvGamma}\bracket{\frac{N}{2}+\alpha_{\sigma},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X\beta})^{\prime}(\boldsymbol{Y}-\boldsymbol{X\beta})+\beta_{\sigma}}. \end{align*}$

$\begin{align*} p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},\sigma^2_{\epsilon}} &\propto \exp{\set{ -\frac{1}{2}\bracket{ \boldsymbol{\beta}^{\prime}\boldsymbol{\Sigma}_{\beta}^{-1}\boldsymbol{\beta} -2\boldsymbol{\beta}^{\prime}\boldsymbol{\Sigma}_{\beta}^{-1}\boldsymbol{\mu}_{\beta} }}}. \end{align*}$

ここで、
$\boldsymbol{\Sigma}_{\beta}=\bracket{\boldsymbol{V}_{\beta}^{-1}+\boldsymbol{X}^{\prime}\boldsymbol{X}/\sigma^2_{\epsilon}}^{-1}, ~~~ \boldsymbol{\mu}_{\beta}=\boldsymbol{\Sigma}_{\beta}\bracket{\boldsymbol{V}_{\beta}^{-1}\boldsymbol{M}_{\beta} + \boldsymbol{X}^{\prime }\boldsymbol{y}/\sigma^2_{\epsilon}}.$

よって、 $p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},\sigma^2_{\epsilon}}\sim\mathcal{N}(\boldsymbol{\mu}_{\beta},\boldsymbol{\Sigma}_{\beta}).$

上のFC $p\Bracket{\boldsymbol{\beta}}{\text{他}}$ は下の特別な場合。
$\boldsymbol{\Sigma}_{\beta}=\bracket{\boldsymbol{V}_{\beta}^{-1}+\boldsymbol{X}^{\prime}\boldsymbol{\Sigma}_{\epsilon}^{-1}\boldsymbol{X}}^{-1}, ~~~ \boldsymbol{\mu}_{\beta}=\boldsymbol{\Sigma}_{\beta}\bracket{\boldsymbol{V}_{\beta}^{-1}\boldsymbol{M}_{\beta} + \boldsymbol{X}^{\prime }\boldsymbol{\Sigma}_{\epsilon}^{-1}\boldsymbol{y}}.$

（Exponential familyのcanonical formにおける自然パラメータを彷彿とさせる形。）

IGもIWと同様の関係にある。

やってみよう

学びたてのころにPyで書いたコードをおいてみる。

期初 $r=0$ において、初期値 $\boldsymbol{\beta}_{(r)}\in\mathbb{R}^{k+1},(\sigma^2_{\epsilon})_{(r)}\in\mathbb{R}_{++}$ を設定する。 $r=1,\dots, R$ で、

$\boldsymbol{\beta}_{(r)} \sim p\Bracket{\boldsymbol{\beta}}{\boldsymbol{y},(\sigma^2_{\epsilon})_{(r-1)}}$
$(\sigma^2_{\epsilon})_{(r)} \sim p\Bracket{\sigma^2_{\epsilon}}{\boldsymbol{y},\boldsymbol{\beta}_{(r)}}$

とサンプリングを行う。

ライブラリなど

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from tqdm import tqdm

R = 3_000
BURN_IN = 1_000
N = 500

trueBeta = np.array([[2, 4]]).transpose()
K = len(trueBeta)

trueSigma2 = 2

"""Artificial Data"""
X = [np.ones((N, 1)), np.sort(np.random.uniform(size=(N, 1)), axis=0)]
X += [np.random.uniform(size=(N, 1)) for _i in range(K - 2)]
X = np.hstack(X)

U = np.random.normal(loc=0, scale=trueSigma2, size=(N, 1))

Y = X.dot(trueBeta) + U

"""r = 0"""
Beta_r0 = np.linalg.solve(X.transpose().dot(X), X.transpose().dot(Y)) # ols
sigma2_r0 = ((Y - X.dot(Beta_r0)) ** 2).mean() # mle

Theta = dict()
Theta['Beta'] = np.zeros((BURN_IN + R + 1, K))
Theta['Beta'][0] = Beta_r0.flatten()
Theta['sigma2'] = np.zeros((BURN_IN + R + 1, 1))
Theta['sigma2'][0] = sigma2_r0

実行

for r in tqdm(range(1, R + BURN_IN + 1)):
    Sigma_beta = np.linalg.solve(np.linalg.inv(V_Beta_prior()) + X.transpose().dot(X) / Theta['sigma2'][r-1], np.eye(K))
    Mu_beta = Sigma_beta.dot(np.linalg.inv(V_Beta_prior()).dot(M_Beta_prior()) + X.transpose().dot(Y) / Theta['sigma2'][r-1])
    
    C = np.linalg.cholesky(Sigma_beta)
    Beta_r = Mu_beta + C.dot(np.random.randn(K, 1))
    Theta['Beta'][r] = Beta_r.flatten()
    
    e = Y - X.dot(Beta_r)
    sigma2_r = 1 / np.random.gamma(N/2 + alpha_sigma2_prior, 1 / (e.transpose().dot(e) / 2 + beta_sigma2_prior))
    Theta['sigma2'][r] = sigma2_r

col = 'Beta'
ax = pd.DataFrame(Theta[col]).add_prefix(f'{col}_').plot(figsize=(14,6), alpha=0.6)
ax.vlines(BURN_IN, Theta[col].min(), Theta[col].max())
ax.hlines(trueBeta, 0, BURN_IN + R, linestyles='--', label='True')
ax.legend()

col = 'Beta'
ax = pd.DataFrame(Theta[col]).add_prefix(f'{col}_').plot.hist(bins=R // 10, figsize=(14,6), alpha=0.6)
ax.vlines(Theta[col][BURN_IN:].mean(axis=0), 0, 60)
ax.vlines(trueBeta, 0, 60, linestyles='--', label='True')
ax.legend()

Theta['Beta'][BURN_IN:][::100].mean(axis=0)
# array([2.03793947, 4.02184397])

fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(
    xlabel='X',
    ylabel='Y',
    # title='Reg Lines over Draws',
)

ax.scatter(X[:, 1], Y, label='Obs')#, c=colors, cmap=plt.get_cmap("viridis"), lw=0)

_X = [np.ones((N, 1)), np.sort(np.random.uniform(size=(N, 1)), axis=0)]
_X += [np.linspace(0, 1, N).reshape(-1,1) for _i in range(K - 2)]
_X = np.hstack(_X)

Betas = Theta['Beta'][BURN_IN:][::30]#.transpose()
for _r in range(len(Betas)):
    ax.plot(_X[:,1], _X.dot(Betas[[_r]].transpose()), linestyle='solid', alpha=0.4, color='r', label=_r*'_' + 'Reg')

ax.plot(_X[:,1], _X.dot(trueBeta), color='black', label="'True'")
ax.legend()

Chan, J., Koop, G., Poirier, D., & Tobias, J. (2019). Bayesian Econometric Methods (2nd ed., Econometric Exercises). Cambridge: Cambridge University Press.