２項分布のPyStan 3サンプルスクリプト
岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版の第10章において、2項分布のPyStanスクリプトについて説明している。PyStan 3を利用するときは、少し変更が必要である。以下に変更例を示す。サンプルスクリプトのファイルは圧縮ファイルc10_files.zipにまとめているので、ダウンロード展開すれば利用できる。
なお、WindowsでPyStan 3を利用するときは、WindowsにUbuntuをインストールしてUbuntu上で利用する。WindowsにUbuntuをインストールすることはWSL（Windows Subsystem for Linux）により簡単にできるが、このウェブサイトで説明している。
拙著第10章のリスト10.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型に変換してArvizの関数summaryにより得ている。コードは以下のようになっている。
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
print(az.summary(i_data))
MCMCサンプリングのトレースはArvizの関数plot_traceによって以下のコードで行っている。
az.plot_trace(i_data)
リスト１ 拙著第10章リスト10.2をPyStan 3用に変更したもの（L10_2.py）
import stan
#
oysta --> stan
import matplotlib.pyplot as plt
import arviz as az
Data = {'N': 10, 'k': 7}
with open('Bin.stan', 'r') as f_stan:
sm =
stan.build(f_stan.read(), data = Data)
# pystan.StanModel(file = 'Bin.stan')
fit = sm.sample()
# sm.sampling(data = Data, n_jobs = 1)
#print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
print(az.summary(i_data))
f = open('fit.txt', 'w')
f.write('fit =
\n{}'.format(az.summary(i_data)))
f.close()
#fit.plot()
az.plot_trace(i_data)
plt.show()
リスト１を実行すると、統計として以下の出力が画面に表示される。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.668 0.133 0.426 0.911 0.004 0.003 1366.0 1557.0 1.0
また、図１のトレース図も表示される。
図１
拙著のリスト10.3をPyStan 3用に変更したものをリスト２に示す。関数buildの第１引数がStanスクリプトの文字列Stan_Modelで与えられている。
リスト２ リスト10.3をPyStan 3用に編集したもの（L10_3.py）
import stan
# pystan
import matplotlib.pyplot as plt
import arviz as az
Stan_Model = """
data{
int N;
int k;
}
parameters{
simplex[2] theta;
}
model{
k ~ binomial(N, theta[1]);
}
"""
Data = {'N': 10, 'k': 7}
with open('Bin.stan', 'r') as f_stan:
sm = stan.build(Stan_Model,
data = Data) # sm = pystan.StanModel(model_code =
Stan_Model)
fit = sm.sample()
# sm.sampling(data = Data, n_jobs = 1)
#print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
print(az.summary(i_data))
f = open('fit.txt', 'w')
f.write('fit =
\n{}'.format(az.summary(i_data)))
f.close()
#fit.plot()
az.plot_trace(i_data)
plt.show()
リスト２を実行すると、MCMCサンプリングの統計が画面に以下のように出力される。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta[0] 0.671 0.129 0.430 0.895 0.004 0.003 1247.0 1948.0 1.0
theta[1] 0.329 0.129 0.105 0.570 0.004 0.003 1247.0 1948.0 1.0
トレース図は、図２のようである。
図２
拙著リスト10.4をPyStan 3用にリスト３のように変更した。
PyStan 3では、MCMCサンプリングの初期値は各マルコフ鎖に対して行う。リスト３では、マルコフ鎖の本数を指定して、各マルコフ鎖の初期値をリスト[f_init()]*4で与えている。コードは以下の通りである。
fit = sm.sample(num_chains = 4, init = [f_init()]*4)
リスト３ リスト10.4をPyStan 3用に変更したもの（L10_4.py）
import stan
#
oysta --> stan
import matplotlib.pyplot as plt
import arviz as az
Data = {'N': 10, 'k': 7}
with open('Bin.stan', 'r') as f_stan:
sm =
stan.build(f_stan.read(), data = Data)
# pystan.StanModel(file = 'Bin.stan')
def f_init():
return dict({'theta': 0.5})
fit = sm.sample(num_chains =
4, init = [f_init()]*4)
# sm.sampling(data = Data, n_jobs = 1)
#print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
print(az.summary(i_data))
f = open('fit.txt', 'w')
f.write('fit = \n{}'.format(az.summary(i_data)))
f.close()
#fit.plot()
az.plot_trace(i_data)
plt.show()
リスト３を実行すると、以下の出力が表示される。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.67 0.13 0.421 0.894 0.003 0.002 1502.0 1568.0 1.0
トレース図は図３のようである。
図３