PyStan 3によるMCMCサンプルのtickleによる保存と読み出し
Stan（PyStan 3）によるMCMCサンプリングの結果をファイルに保存して、そのファイルを別のプログラムで読み出して分析することを試みた。MCMC
サンプリングに時間が掛かるときは、サンプルを保存しておいて別のプログラムで描画や事後分析を行うとサンプルの分析を効率よく行うことができる。
以下の例では、PyStan 3の実行をWindowsにWSLによりインストールされたUbuntuで行い、MCMCサンプルの分析プログラムはWindows上のPythonで実行した。WindowsへのUbuntuとPythonのインストールについては、このウェブサイトで説明している。
PyStan 3のスクリプト例としてリスト１のものを用意した。
リスト１ MCMCサンプリング（PyStan 3）
import stan
import matplotlib.pyplot
as plt
import arviz
as az
import pickle
stan_code = """
data {
int N;
int k;
}
parameters {
real<lower=0, upper=1> p;
}
model {
k ~
binomial(N, p);
}
"""
sm = stan.build(stan_code, data={'N':10, 'k':7})
fit = sm.sample()
#
# MCMCサンプリングの結果fitをPandasのDataFrame型に変換して保存
#
dframe = fit.to_frame()
with open('dframe.pkl',
'wb') as f:
pickle.dump(dframe, f)
plt.hist(dframe['p'])
plt.show()
#
# MCMCサンプリングの結果をArvizのInferenceData型として保存
#
infdata = az.from_pystan(posterior
= fit, posterior_model = sm)
with open('infdata.pkl',
'wb') as f:
pickle.dump(infdata, f)
summary = az.summary(infdata)
print('\nSummary...\n',
summary)
az.plot_trace(infdata)
plt.show()
次のコード
sm = stan.build(stan_code, data={'N':10, 'k':7})
fit = sm.sample()
で作成したサンプルfitを2種類のファイルに保存している。
１つはPandasのDataFrameの形式で、次のコード
dframe = fit.to_frame()
with open('dframe.pkl',
'wb') as f:
pickle.dump(dframe, f)
で保存されている。この形式のデータはPandasのDataFrameとして分析できる。
次のコードは、パラメータpの事後分布を描くものである。
plt.hist(dframe['p'])
plt.show()
dframeから'p'をキーとしてパラメータpの事後分布のサンプルを取り出している。
図１
サンプルfitは、ArvizのInferenceData型として次のコード
infdata = az.from_pystan(posterior
= fit, posterior_model = sm)
with open('infdata.pkl',
'wb') as f:
pickle.dump(infdata, f)
で保存することができる。InferenceData型のデータはArvizの関数で分析することができる。
ＭＣＭＣサンプルの統計量を次のコードで求めている。
summary = az.summary(infdata)
print('\nSummary...\n',
summary)
上のコードが実行されると、以下の出力が端末に表示される。
Summary...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.662 0.128 0.422 0.89 0.003 0.002 1498.0 1683.0 1.0
次のコード
az.plot_trace(infdata)
plt.show()
を実行すると、MCMCサンプリングのトレース図が描画される。
図２
リスト１を実行すると、次のコード
with open('dframe.pkl',
'wb') as f:
pickle.dump(dframe, f)
により、dframeがファイルdframe.pklに保存され、
次のコード
with open('infdata.pkl',
'wb') as f:
pickle.dump(infdata, f)
により、infdataがファイルinfdata.pklに保存される。
リスト１の実行によって保存されたファイルdframe.pklとinfdata.pklを読み込んで処理するスクリプト例をリスト２に示す。
リスト２ ファイルの読み込みと処理
import matplotlib.pyplot
as plt
import arviz
as az
import pickle
#
# dframe.pklからの読み出し
#
with open('dframe.pkl',
'rb') as f:
dframe
= pickle.load(f)
plt.hist(dframe['p'])
plt.show()
#
# infdata.pklからの読み出し
#
with open('infdata.pkl',
'rb') as f:
infdata
= pickle.load(f)
summary = az.summary(infdata)
print('\nSummary...\n',
summary)
az.plot_trace(infdata)
plt.show()
リスト２は、WindowsにインストールされたPythonで実行することができる。
以下は、WindowsにインストールされたPython 3.12.2の実行により表示されるグラフと端末出力である。
図３
図４
図５
Jupyter Notebookによるリスト２の実行例を以下に示す。PythonのスクリプトファイルをNotebookにコピー・貼り付けを行う例については、このウェブサイトで説明している。
なお、PyStan2の場合のMCMCサンプルの保存と読み出しについては、岡本（2019）のpp. 197-200で説明している。