２項分布のStanによるベイズ分析
StanのPython用モジュールPyStan 3は、Windows上では直接利用できない。WindowsにWSL（Windows Subsystem for Linux）をインストールして用いることができる。WindowsでのPyStan 3の準備については、このウェブサイトなどで説明している。
WindowsでPyStan 3を利用するときは、PyStan 3によるMCMCサンプリングをUbuntu/WSL/Windowsで行うが、WSLでは描画ができなかったので、描画はWindows上で行った。
2項分布のStanスクリプト例をリスト１に示す。
リスト１ MCMCサンプリング
import stan
import arviz as az
import pickle
StanCode = """
data {
int
n;
int
k;
}
parameters {
real<lower=0.0, upper=1.0> p;
}
model {
k ~
binomial(n, p);
}
"""
sm = stan.build(StanCode,
data={'n':10, 'k':7})
fit = sm.sample()
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
print(az.summary(inf_data))
with open('inf_data.pkl', 'wb') as f:
pickle.dump(inf_data, f)
df_sample = fit.to_frame()
with open('df_sample.pkl', 'wb') as f:
pickle.dump(df_sample, f)
Stanスクリプトは、文字列StanCodeとして記述されている。Stan構文はPyStan 2の場合と同様である（拙著ではPyStan2について説明している）。
次のコードでStanスクリプトのコンパイル・ビルドが行われる。
sm = stan.build(StanCode,
data={'n':10, 'k':7})
MCMCサンプリングは、次のコードで行われる。
fit = sm.sample()
次のコードで、MCMCサンプリングがArvizのInferenceData型のオブジェクトinf_dataとして構成され、az.summary(inf_data)によって統計が計算される。
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
print(az.summary(inf_data))
with open('inf_data.pkl', 'wb') as f:
pickle.dump(inf_data, f)
オブジェクトinf_dataは、pickleによってファイルに保存される。
次のコードは、MCMCサンプルデータfitをDataFrame型のオブジェクトdf_sampleとしてファイルに保存するものである。
df_sample = fit.to_frame()
with open('df_sample.pkl', 'wb') as f:
pickle.dump(df_sample, f)
リスト１のスクリプトは、Miniconda/Ubuntu24/WSL/Windows10に用意したPython 3.12仮想環境において実行した。
リスト１のスクリプトの実行によってファイルに保存されたデータは、リスト２のスクリプトで読み出されてグラフ描画が行われる。リスト１のスクリプトファイルとリスト２のスクリプトファイルは同じフォルダに置いておく。
リスト２は、Windows上でPython 3.13.2により実行した。
リスト２ 描画スクリプト
import arviz as az
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
with open('inf_data.pkl', 'rb') as f:
inf_data = pickle.load(f)
az.plot_trace(inf_data)
plt.show()
with open('df_sample.pkl', 'rb') as f:
df_sample = pickle.load(f)
sb.kdeplot(df_sample['p'])
plt.show()
次のコードで、inf_data.pklが読み込まれ、関数az.plot_trace(inf_data)によってトレース図が描画される（図１）。
with open('inf_data.pkl', 'rb') as f:
inf_data = pickle.load(f)
az.plot_trace(inf_data)
plt.show()
図１ az.plot_trace(inf_data)による描画
次のコードで、事後分布が描かれる（図２）。
with open('df_sample.pkl', 'rb') as f:
df_sample = pickle.load(f)
sb.kdeplot(df_sample['p'])
plt.show()
図２ sb.kdeplot(df_sample['p'])による描画
モジュールpickleについては、拙著 pp.197-200を参照されたい。