StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - 12.6 1次元の空間構造

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

実行環境

インポート

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

データ読み込み

kubo11a = pd.read_csv('./data/data-kubo11a.txt')

12.6 1次元の空間構造

data = dict(
    I=kubo11a.index.size,
    Y=kubo11a['Y']
)
stanmodel = pystan.StanModel('./stan/model12-11.stan')
fit = stanmodel.sampling(data=data, seed=1234)
stanmodel_b = pystan.StanModel('./stan/model12-12.stan')
fit_b = stanmodel_b.sampling(data=data, iter=4000, seed=1234)

probs = (10, 25, 50, 75, 90)
columns = ['p{}'.format(p) for p in probs]
_, axes = plt.subplots(1, 2, figsize=figaspect(3/8), sharex=True, sharey=True)

for i, f in enumerate([fit, fit_b]):
    ax = axes[i]
    ms = f.extract()

    d_est = pd.DataFrame(np.percentile(ms['Y_mean'], probs, axis=0).T, columns=columns)
    d_est['x'] = d_est.index + 1

    ax.plot('x', columns[2], data=d_est, color='k')
    for j in range(2):
        ax.fill_between('x', columns[j], columns[-j-1], data=d_est, color='k', alpha=0.2*(j+1))
    ax.scatter(kubo11a.index+1, kubo11a['Y'], facecolor='w', edgecolor='k')
    plt.setp(ax, xlabel='i', ylabel='Y[i]')

plt.show()