StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - 10.3 再パラメータ化

10822 ワード
PyStan Stan Python Python テキストリンク

実行環境

インポート

import numpy as np
from scipy import stats
import pandas as pd
import pystan
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
%matplotlib inline

データ読み込み

salary2 = pd.read_csv('./data/data-salary-2.txt')

10.3 再パラメータ化

10.3.1 Nealの漏斗

fit_a = pystan.stan('./stan/model10-8a.stan', chains=1, iter=4000, seed=1234)
fit_b = pystan.stan('./stan/model10-8b.stan', chains=1, iter=4000, seed=1234)

ms_a = fit_a.extract()
ms_b = fit_b.extract()

xx, yy = np.mgrid[-5:5:30j, -5:5:30j]
points = np.c_[xx.ravel(), yy.ravel()]
x, y = points[:, 0], points[:, 1]
lp = np.log(stats.norm.pdf(y, loc=0, scale=3)) + np.log(stats.norm.pdf(x, loc=0, scale=np.exp(y/2)))
lp[lp < -15] = -15

_, axes = plt.subplots(1, 2, figsize=figaspect(1/2), sharex=True, sharey=True)
for ms, ax in zip([ms_a, ms_b], axes):
    cs = ax.contourf(xx, yy, lp.reshape(xx.shape), vmin=-15, vmax=0)
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    ax.scatter(ms['r'][:, 0], ms['a'], s=1, c='k')
    plt.setp(ax, xlabel='r[1]', ylabel='a', xlim=xlim, ylim=ylim)
plt.show()

10.3.2 階層モデル

data = {col: salary2[col] for col in salary2.columns}
data['N'] = salary2.index.size
data['K'] = salary2['KID'].max()
fit = pystan.stan('./stan/model8-4c.stan', data=data, seed=1234)