ポアッソン回帰モデルのベイズ分析
PyStan 3対応
岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版の第12章において、ポアッソン回帰モデルのPyStanスクリプトについて説明している。これはPyStan 2対応であるので、PyStan 3を利用するときは、少し変更が必要である。以下に変更例を示す。サンプルスクリプトのファイルは圧縮ファイルc12_files.zipにまとめているので、ダウンロード展開すれば利用できる。
なお、WindowsでPyStan 3を利用するときは、WindowsにUbuntuをインストールしてUbuntu上で利用する。WindowsにUbuntuをインストールすることはWSL（Windows Subsystem for Linux）により簡単にできるが、このウェブサイトで説明している。
拙著第12章のリスト12.2は次のリスト１のように変更するとPyStan 3が利用できる。
PyStan 2では、StanModelによる翻訳、samplingによるMCMCサンプリングの実行と進むが、PyStan 3ではbuildによるビルド、sampleによる実行と進む。リスト１では、
sm
= stan.build(f_stan.read(),
data = Data)
fit = sm.sample()
となっている。関数buildの第１引数は、Stanスクリプトの文字列であり、パラメータdataにMCMCサンプリングで使用されるデータが設定されている。第１引数の文字列は、上の例ではファイルからの読み込み値f_stan.read()が設定されている。
MCMCサンプリングは関数sampleによって実行される。
MCMCサンプリングの統計は、サンプリング結果fitを次のコードでArvizのInferenceData型に変換したi_dataから求めている。
i_data = az.from_pystan(posterior
= fit, posterior_model = sm)
summary = az.summary(i_data, var_names = ['mu', 'a2',
'a3', 'b2', 'b3',
'g22', 'g23', 'g32', 'g33'])
print(summary)
統計を求めるパラメータは、関数summaryのパラメータvar_namesで指定している。指定を行わないときは、すべてのパラメータについて統計が計算される。
リスト１では、fitをPandasのDataFrame型に変換している。コードは次のようである。
d_frame = fit.to_frame()
事後分布などは、PandasのDataFrame型オブジェクトのデータとして取り出して処理が行われている。
リスト１ 拙著第12章リスト12.2をPyStan 3用に変更したもの（L12_2.py）
import stan
# pystan
--> stan
import matplotlib.pyplot
as plt
import numpy
as np
import arviz
as az
RawData = [[1, 1, 0, 0, 1, 0, 0],
[1, 1, 0, 0, 0, 1, 0],
[10, 1, 0, 0, 0, 0, 1],
[3, 0, 1, 0, 1, 0, 0],
[2, 0, 1, 0, 0, 1, 0],
[10, 0, 1, 0, 0, 0, 1],
[6, 0, 0, 1, 1, 0, 0],
[2, 0, 0, 1, 0, 1, 0],
[4, 0, 0, 1, 0, 0, 1]]
F = []
Xr1 = []
Xr2 = []
Xr3 = []
Xc1 = []
Xc2 = []
Xc3 = []
for v in RawData:
F.append(v[0])
Xr1.append(v[1])
Xr2.append(v[2])
Xr3.append(v[3])
Xc1.append(v[4])
Xc2.append(v[5])
Xc3.append(v[6])
for i in
range(9):
print(' ', F[i], ' ', Xr1[i], ' ', Xr2[i], ' ', Xr3[i],
' ', Xc1[i],
' ', Xc2[i],
' ', Xc3[i])
N = 9
Data = {'N': N, 'F': F, 'Xr1': Xr1,
'Xr2': Xr2, 'Xr3': Xr3,
'Xc1':
Xc1, 'Xc2': Xc2, 'Xc3': Xc3}
#sm = pystan.StanModel(file = 'Poisson.stan')
#fit = sm.sampling(data
= Data,
#
pars = ['mu', 'a2', 'a3', 'b2', 'b3',
#
'g22', 'g23', 'g32', 'g33'], n_jobs = 1)
with open('Poisson.stan',
'r') as f_stan:
sm
= stan.build(f_stan.read(),
data = Data)
fit = sm.sample()
#print(fit)
i_data = az.from_pystan(posterior
= fit, posterior_model = sm)
summary = az.summary(i_data, var_names = ['mu', 'a2',
'a3', 'b2', 'b3',
'g22',
'g23', 'g32', 'g33'])
print(summary)
d_frame = fit.to_frame()
a3 = d_frame['a3']
ck_a3 = a3 > 0
print('P(a3 > 0.0) =', np.sum(ck_a3)/len(ck_a3))
b3 = d_frame['b3']
ck_b3 = b3 > 0
print('P(b3 > 0.0) =', np.sum(ck_b3)/len(ck_b3))
g33 = d_frame['g33']
ck_g33 = g33 < 0
print('P(g33 < 0.0) =',np.sum(ck_g33)/len(ck_g33))
リスト１を実行すると、以下の出力が端末に表示される。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu -0.569 1.234 -2.877 1.414 0.050 0.036 648.0 647.0 1.00
a2 1.502 1.377 -0.873 4.158 0.056 0.039 640.0 748.0 1.00
a3 2.269 1.296 0.041 4.764 0.051 0.038 688.0 735.0 1.00
b2 0.000 1.737 -3.341 3.171 0.066 0.047 684.0 912.0 1.01
b3 2.815 1.278 0.669 5.188 0.052 0.037 644.0 680.0 1.00
g22 -0.504 1.969 -4.029 3.347 0.073 0.052 726.0 1063.0 1.01
g23 -1.493 1.449 -4.325 0.995 0.057 0.040 662.0 787.0 1.00
g32 -1.281 1.942 -4.798 2.481 0.070 0.049 775.0 1237.0 1.01
g33 -3.258 1.443 -5.999 -0.719 0.055 0.040 719.0 738.0 1.00
P(a3 > 0.0) = 0.98725
P(b3 > 0.0) = 0.999
P(g33 < 0.0) = 0.99825