PyStan 3による2項分布の分析
拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章で扱ったPyStan 2による2項分布の分析StanスクリプトをPyStan 3で実行した。Stanスクリプトは変更なしで用いることが出来たが、PythonスクリプトはPyStan 3用に書き換えた。書き換えたスクリプトを、リスト１とリスト２に示す。スクリプトファイルは、圧縮ファイルfiles_bin.zipにまとめた。
リスト１は、Stanスクリプトが文字列で与えられている場合である。Stanスクリプトを別ファイルとして用意する例は、本ウェブサイトの後半に示す。
リスト１ PyStan 3によるMCMCサンプリングを行うスクリプト（BinStr.py）
import stan
import arviz as az
import pickle
Stan_Model = """
data{
int N;
int k;
}
parameters{
simplex[2] theta;
}
model{
k ~ binomial(N, theta[1]);
}
"""
Data = {'N': 10, 'k': 7}
sm = stan.build(Stan_Model, data =
Data)
fit = sm.sample()
with open('v_theta.pkl', 'wb') as f:
pickle.dump(fit['theta'], f)
dframe = fit.to_frame()
print('Keys...\n', dframe.keys())
with open('dframe.pkl', 'wb') as f:
pickle.dump(dframe, f)
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
with open('inf_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print(az.summary(inference_data))
リスト２ リスト１に続いて実行するスクリプト（DrawBinStr.py）
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
with open('dframe.pkl', 'rb') as f:
dframe = pickle.load(f)
print('keys...\n', dframe.keys())
with open('inf_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
summary = az.summary(inference_data)
print('Summary:\n', summary)
az.plot_trace(inference_data)
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
print(dframe)
sb.kdeplot(dframe['theta.1'])
plt.title("KDEPlot of
dframe['theta.1]", fontsize = 16)
plt.savefig('Fig1.png')
plt.show()
with open('v_theta.pkl', 'rb') as f:
v_theta = pickle.load(f)
print(v_theta[:5])
sb.kdeplot(v_theta[0])
plt.title('KDEPlot of v_theta')
plt.savefig('Fig2.png')
plt.show()
リスト１を実行してからリスト２を実行する。ただし、Windowsの場合は、リスト１は直接Windows上では実行できない。WindowsにWSL（Windows Subsystem for Linux）を用意して、その上にインストールされたUbuntu上でPyStan3を用いる。詳しくはこのウェブサイトを参照されたい。PyStan 3は、直接Windows上で実行できないので、リスト１とリスト２に分けたが、Windows以外では分ける必要はない。
リスト１では、PyStan 3は次のコードで実行されている。
sm = stan.build(Stan_Model, data =
Data)
fit = sm.sample()
ここで、PyStan 3は次のimport文
import stan
でimportされていることに注意。
PyStan 3では、まず、関数buildによってStanスクリプトのコンパイルとデータを設定したビルドを行い、Modelオブジェクトを作成する。
次に、Modelオブジェクトの関数sample()の呼び出してMCMCサンプリングを行い、結果をFitオブジェクトとして返す。
パラメータthetaの事後分布のサンプル値fit['theta']を次のコードでファイルに保存している。
ith open('v_theta.pkl', 'wb') as f:
pickle.dump(fit['theta'], f)
Fitオブジェクトの関数to_frame()は、Fitオブジェクトの内容をPandas DataFrameにまとめるものである。
dframe = fit.to_frame()
DataFrame dframeのキーは、パラメータ名に対応していて、リスト１の場合、次のコード
print('keys...\n', dframe.keys())
を実行すると、以下のように出力を得る。
Keys...
Index(['lp__', 'accept_stat__',
'stepsize__', 'treedepth__', 'n_leapfrog__',
'divergent__', 'energy__', 'theta.1', 'theta.2'],
dtype='object',
name='parameters')
次のパラメータ
simplex[2] theta;
は、キーとして
'theta.1', 'theta.2'
で表されていることが分かる。
このDataFrameオブジェクトdframeをpickleによりファイルに保存して、次のリスト２で用いる。
次のコード
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
は、FitオブジェクトをInferenceDataオブジェクトにまとめるものであり、これもpickleでファイルに保存して、次のリスト２で用いる。
リスト２のスクリプトは、直接Windows上で実行できる。リスト１とリスト２のスクリプトファイルは、同じフォルダ内に置いておくと、pickleで扱うファイルも同じフォルダ内で扱うことが出来て便利である。
リスト１の実行後、リスト２を実行すると図１のグラフが表示される。
図１
これは、次のコード
with open('inf_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
で読み込んだものを、次のコード
az.plot_trace(inference_data)
で描画したものである。
図１のフォームを閉じると、図２のグラフが表示される。
図２
これは、次のコード
with open('dframe.pkl', 'rb') as f:
dframe = pickle.load(f)
で読み込んだものを、次のコード
sb.kdeplot(dframe['theta.1'])
で描画したものである。
図２のフォームを閉じると、図３のグラフ表示される。
図３
これは、次のコード
with open('v_theta.pkl', 'rb') as f:
v_theta = pickle.load(f)
で読み込んだものを、次のコード
sb.kdeplot(v_theta[0])
で描画したものである。
図３のフォームを閉じると、プログラムの実行終了である。
StanスクリプトをPythonスクリプトファイルとは別のファイルとして用意する場合
PyStan 3のビルドでは、Stanコードは文字列で与えられるが、Stanコードをファイルとして用意している場合は、そのスクリプトファイルを読み込んで、文字列としてbuildの引数とすればよい。リスト３では、リスト４のStanコードのスクリプトファイルを読み込んで、ビルドを行っている。スクリプトファイルは、圧縮ファイルfiles_bin.zipにまとめているので自由にダウンロード展開して利用されたい。
Stanスクリプトがファイルとして用意されているときは、ファイルを読み込んで文字列としてbuildの第１パラメータに指定すればよい。すなわち、
with open('Bin.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
である。
後は、リスト１の場合と同様であるが、リスト１ではパラメータthetaが
simplex[2] theta;
と、simplex型であるが、リスト４では
real<lower = 0.0, upper =
1.0> theta;
と、real型である。
これに伴い、リスト３の実行によるMCMCサンプリングの結果の分析スクリプト（リスト５）に変更がある。
また、Windowsでは、直接リスト３を実行することはできないので、WSL上のUbuntuで実行することはリスト１の場合と同じである。
リスト５は、直接Windows上で実行できる。
リスト３ PyStan 3のビルドにおいて、Stanコードをファイルから読み込む場合（BinFile.py）
import stan
import arviz as az
import pickle
Data = {'N': 10, 'k': 7}
with open('Bin.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
fit = sm.sample()
with open('v_theta.pkl', 'wb') as f:
pickle.dump(fit['theta'], f)
print('fit...\n', fit['theta'])
dframe = fit.to_frame()
print('Keys...\n', dframe.keys())
with open('dframe.pkl', 'wb') as f:
pickle.dump(dframe, f)
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
with open('inf_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print(az.summary(inference_data))
リスト４のStanスクリプトは、リスト１の場合とパラメータthetaの宣言を少し変えてみた。
data{
int N;
int k;
}
parameters{
real<lower = 0.0, upper =
1.0> theta;
}
model{
k ~ binomial(N, theta);
}
リスト３のMCMCサンプリングのスクリプトの実行後、サンプリングの分析スクリプト（リスト５）を実行する。
Windowsでは、リスト３のスクリプトは、Anaconda/Ubuntu/WSL/Windowsに用意したPython 3.10仮想環境において実行したが、リスト５は直接Windows上で実行できる。
リスト５ リスト３によるMCMCサンプリングの結果の処理（DrawBinFile.py）
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
with open('dframe.pkl', 'rb') as f:
dframe = pickle.load(f)
print('Keys...\n', dframe.keys())
with open('inf_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
summary = az.summary(inference_data)
print('Summary:\n', summary)
az.plot_trace(inference_data)
plt.tight_layout()
plt.show()
print(dframe)
sb.kdeplot(dframe['theta'])
plt.title("KDEPlot
dframe['theta']", fontsize = 16)
plt.savefig('Fig1.png')
plt.show()
with open('v_theta.pkl', 'rb') as f:
v_theta = pickle.load(f)
print(v_theta)
sb.kdeplot(v_theta[0])
plt.title("KDEPlot
v_theta[0]", fontsize = 16)
plt.savefig('Fig2.png')
plt.show()
リスト３のコード
with open('v_theta.pkl', 'wb') as f:
pickle.dump(fit['theta'], f)
で保存したデータは、次のコード
print('fit...\n', fit['theta'])
の出力
fit...
[[0.73891671 0.44881966 0.68348178 ...
0.75439311 0.77459035 0.65838827]]
に示されるように、（１，サンプルサイズ）型の配列になっているので、リスト５におけるthetaのサンプルの描画コードも、それに合わせて
sb.kdeplot(v_theta[0])
となっている。
リスト５を実行すると、まず図４のグラフが表示される。
図４
図１のグラフと異なるのは、図１ではsimplex型のパラメータtheta[1]とtheta[2]の2つのパラメータの分布が描かれているが、図４ではreal型のパラメータthetaの分布が描かれているからである。