拙著「Pythonでデータ分析」の第10章「PyStanによる２項分分析――Stan入門――」更新と補足
PyStan 3入門
PyStan 3をWindowsで利用するときは、WindowsにWSL (Windows Subsystem for Linux)をインストールして、その上にUbuntuをインストールして用いる。WSLとWSLでのAnacondaのインストールについては、このウェブサイトを用意して説明しているので、参照されたい。
WindowsにWSLをインストールすると、WSLにUbuntuがインストールされる。Ubuntu/WSL/WindowsにAnacondaをインストールして、Python3.11の仮想環境を構築する。Python 3.11の仮想環境においてpipコマンドを
pip install
pystan
と実行してPyStan 3をインストールする。
Pythonスクリプト内でPyStan 3をimportするときは
import stan
である。
Anaconda/Ubuntu/WSL/Windowsにおいて、PyStan 3のインストールでエラーが出る、あるいはスクリプトを実行したときに「CompileError」が表示されたときは、Ubuntuに適当なC++コンパイラが見つからなかったということなので、例えば、コマンド「g++」を実行してみるとよい。
「g++が見つからない」というようなメッセージとともに、g++をインストールするコマンドが示されるので、その通りに実行してみる。
「g++」のインストールが済めば、PyStan 3のインストールおよび実行が出来るはずである。
PyStan 3では、関数buildによってStanスクリプトで書かれたモデルとデータの組み合わせに対してビルドが行われ、その後、関数sampleによってMCMCサンプリングが行われる。Stanスクリプトは文字列として関数buildの第１パラメータに設定する。Stanスクリプトがファイルとして用意されているときは、ファイルの読み込み関数readによりファイルの中身を文字列として渡せばよい。本サイトにおけるサンプルスクリプトでは、Stanスクリプトのファイルを用意して、ファイルからStanスクリプトを読み込んでいる。
以下の２つの場合について説明する。
（１）MCMCサンプリングとサンプルの分析を１つのスクリプトファイルで行う場合。
（２）MCMCサンプリングとサンプルの分析を分けてスクリプトを用意する場合。MCMCサンプリングを１回行いファイルに保存すれば、サンプルの分析はファイルに保存されたものを使ってできる。
サンプルスクリプトは、まとめて圧縮ファイルc10files.zipとした。ダウンロード・展開して自由に利用されたい。
MCMCサンプリングと分析を１つのスクリプトファイルで行う場合。
２項分布モデルのStanスクリプトをリスト１のように用意した。
リスト１．2項分布モデルのStanスクリプト（stan_binomial.stan）
data {
int N;
int k;
}
parameters {
real<lower=0, upper=1>
theta;
}
model {
k ~ binomial(N, theta);
}
リスト１のStanスクリプトを用いるPythonスクリプトをリスト２のように用意した。
リスト２．リスト１のStanスクリプト（sta_binomial.stan）を用いるPythonスクリプト（sampling_analysis.py）
import stan
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import pandas as pd
N = 10
k = 7
with open('stan_binomial.stan') as f:
sm = stan.build(f.read(),
data = {'N':N, 'k':k})
fit = sm.sample(num_samples=10000)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
d_fm = fit.to_frame()
print('N =', N)
print('k =', k)
with open('Results.txt',
"w") as f:
f.write(az.summary(i_data).__str__())
f.write('\n\n')
f.write(d_fm['theta'].describe(percentiles =
[0.025, 0.25, 0.5, 0.75, 0.975]).__str__())
print('Results.txt was saved.')
az.plot_trace(i_data)
plt.savefig('FigTrace.png')
plt.show()
print('FigTrace.png was saved.')
sb.kdeplot(d_fm['theta'])
plt.xlabel(r'$\theta$', fontsize = 16)
plt.savefig('FigTheta.png')
plt.show()
print('FigTheta.png was saved.')
リスト２のスクリプトをコマンド
(py311) ・・・/c10files$ python
sampling_analysis.py
で実行すると、サンプリングが始まる。サンプリングが終了すると、トレース図が表示される（図１）。
図１
図１のフォームを閉じると、図２のθの事後分布のグラフが表示される。
図２
図２のフォームを閉じると、実行終了となり、下記のメッセージが表示される。
Results.txt was saved.
FigTrace.png was saved.
FigTheta.png was saved.
出力ファイルResults.txtを開くと、以下のようである。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.667 0.131 0.426 0.901 0.001 0.001 14393.0 16592.0 1.0
count 40000.000000
mean
0.666952
std
0.131293
min
0.146446
2.5%
0.391165
25%
0.577323
50%
0.676214
75%
0.766141
97.5%
0.890249
max
0.977678
Name: theta, dtype: float64
MCMCサンプリングとサンプルの分析を別のスクリプトとして用意する場合。
MCMCサンプリングを行うPythonスクリプトをリスト３に示す。Stanスクリプトは、リスト１のものを用いている。
コマンド
(py311) ・・・/c10files$ python
mcmcsampling.py
でスクリプトファイルを実行するとMCMCサンプリングが始まり、MCMCサンプリングが終了してサンプリングの結果がファイルに保存されると、次のメッセージが表示されて終了となる。
temp.pkl was
saved.
リスト３のスクリプトの実行後、リスト４のスクリプトをコマンド
(py311) ・・・/c10files$ python analysis.py
で実行する。ファイルtemp.pklが読み込まれて、分析が始まる。
まず、トレース図が表示される（図３）。
図３
図３のフォームを閉じると、θの事後分布のグラフが表示される（図４）。
図４
図４のフォームを閉じると、実行終了となり、端末には下記のメッセージが表示される。
Results.txt was saved.
FigTrace.png was saved.
FigTheta.png was saved.
出力ファイルResults.txtを開くと、以下のようになっている。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.667 0.13 0.426 0.902 0.001 0.001 15543.0 17000.0 1.0
count 40000.000000
mean
0.667451
std
0.130159
min
0.148430
2.5%
0.391304
25%
0.580336
50%
0.677678
75%
0.763537
97.5%
0.891808
max
0.988770
Name: theta, dtype: float64
リスト３．MCMCサンプリングPythonスクリプト（mcmcsampling.py）
import stan
import arviz as az
import pickle
N = 10
k = 7
with open('stan_binomial.stan') as f:
sm = stan.build(f.read(),
data = {'N':N, 'k':k})
fit = sm.sample(num_samples=10000)
print(fit)
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
d_fm = fit.to_frame()
"""
サンプリングデータの保存
"""
with open('temp.pkl', 'wb') as f:
pickle.dump((N, k,
inference_data, d_fm), f)
print('temp.pkl was saved.')
リスト４．サンプルの分析Pythonスクリプト（analysis.py）
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pickle
import pandas as pd
"""
サンプリングデータの読み出し
"""
with open('temp.pkl', 'rb') as f:
N, k, i_data, d_fm =
pickle.load(f)
print('N =', N)
print('k =', k)
print('Summary...\n',
az.summary(i_data))
print('\n',
d_fm['theta'].describe(percentiles = [0.025, 0.25, 0.5, 0.75, 0.975]))
with open('Results.txt',
"w") as f:
f.write(az.summary(i_data).__str__())
f.write('\n\n')
f.write(d_fm['theta'].describe(percentiles =
[0.025,
0.25, 0.5, 0.75, 0.975]).__str__())
print('Results.txt was saved.')
az.plot_trace(i_data)
plt.savefig('FigTrace.png')
plt.show()
print('FigTrace.png was saved.')
sb.kdeplot(d_fm['theta'])
plt.xlabel(r'$\theta$', fontsize = 16)
plt.savefig('FigTheta.png')
plt.show()
print('FigTheta.png was saved.')