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

実行環境

インポート

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)