やってみよう
学びたてのころにPyで書いたコードをおいてみる。
期初において、初期値 を設定する。で、
とサンプリングを行う。
ライブラリなど
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))
sigma2_r0 = ((Y - X.dot(Beta_r0)) ** 2).mean()
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)
fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(
xlabel='X',
ylabel='Y',
)
ax.scatter(X[:, 1], Y, label='Obs')
_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]
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.