StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - Chapter 4 練習問題

実行環境

インポート

import pystan
import matplotlib.pyplot as plt
%matplotlib inline

データの準備

%load_ext rpy2.ipython

%%R
set.seed(123)
N1 <- 30
N2 <- 20
Y1 <- rnorm(n=N1, mean=0, sd=5)
Y2 <- rnorm(n=N2, mean=1, sd=4)

from rpy2.robjects import r
N1 = 30
N2 = 20
Y1 = r.get('Y1')
Y2 = r.get('Y2')

(1)

plt.hist([Y1, Y2], density=True, label=['Y1', 'Y2'])
plt.legend()
plt.show()

(4)

data = dict(N1=N1, N2=N2, Y1=Y1, Y2=Y2)
fit = pystan.stan('./stan/model4ex1.stan', data=data, seed=1234)
ms = fit.extract()
prob = (ms['mu1'] < ms['mu2']).mean()
print(prob)

0.9305

(5)

data = dict(N1=N1, N2=N2, Y1=Y1, Y2=Y2)
fit = pystan.stan('./stan/model4ex2.stan', data=data, seed=1234)
ms = fit.extract()
prob = (ms['mu1'] < ms['mu2']).mean()
print(prob)

0.935