拙著「Pythonでデータ分析」の第13章「PyMCによる２項分布分析」の補足説明
スクリプトの実行は、Python
3.10/Anaconda/Ubuntu/WSL/Windowsで行った。すなわち、WindowsにWSL (Windows Subsystem for Linux)をインストールするとUbuntuもインストールされるので、そのUbuntuにAnacondaをインストールしてPython 3.10仮想環境を用意した。Python 3.10仮想環境においてコマンド
conda install
pymc3
によりPyMC3をインストールした。詳しくは、このウェブサイトを参照されたい。
MCMCサンプリングを用いる場合のスクリプトの実行を
（２）MCMCサンプリングとサンプルの分析に分けて行う場合。MCMCサンプリングを１回行えば、そのサンプリング結果を用いて分析スクリプトはいろいろ工夫できる。
の３つの場合に分けて説明する。
スクリプトファイルなどは、まとめて圧縮ファイルc13files.zipとした。ダウンロード・展開すれば自由に利用できる。
スクリプトをリスト１のように用意した。Anaconda/Ubuntu/WSL/WindowsのPython 3.10仮想環境において、リスト１のスクリプトファイルBinOneScrpt.pyを実行するためにコマンド
(py310) ・・・/c13files$ python
BinOneScrpt.py
を実行すると、MCMCサンプリングが始まり、サンプリングの終了後、トレース図（図１）が表示される。
図１
図１のフォームを閉じると、θの事後分布のグラフが表示される（図２）
図２
図２のフォームを閉じると、端末に次のメッセージが表示されて実行終了となる。
FigTheta.png was saved.
Summary.txt was saved.
出力ファイルSummary.txtを開くと、以下のようになっている。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.666 0.131 0.42 0.889 0.003 0.002 1576.0 2395.0 1.0
MCMCサンプリングとサンプルの分析とに、スクリプトを分ける場合
MCMCサンプリングのスクリプトをリスト２に、サンプルの分析スクリプトをリスト３に示す。これらのスクリプトをAnaconda/Ubuntu/WSL/WindowsのPython 3.10仮想環境で実行した。
まず、リスト２のスクリプトファイルBin2_1.pyを次のコマンドで実行する。
(py310) ・・・/c13files$ python Bin2_1.py
スクリプトの実行によりMCMCサンプリングが始まる。サンプリングが終了すると、次のメッセージが端末に表示されて実行終了となる。
temp.pkl was
saved.
MCMCサンプリングデータがtemp.pklに保存されている。
次に、リスト３のスクリプトをコマンド
(py310) ・・・/c13files$ python Bin2_2.py
で実行して、上記のファイルtemp.pklを読み込み、サンプルの分析を行う。
まず、トレース図が表示される（図３）。
図３
図３のフォームを閉じると、図４の事後分布のグラフが表示される（図４）。
図４
図４のフォームを閉じると、実行終了となり、以下のメッセージが端末に表示されている。
Summary.txt was saved.
FigTrace.png was saved.
FigTheta.png was saved.
出力ファイルSummary.txtを開くと、以下のような内容である。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.668 0.134 0.422 0.907 0.003 0.002 1833.0 2457.0 1.0
Anaconda/Ubuntu/WSL/WindowsのPython 3.10仮想環境において、Jupyter NotebookでPyMC3を利用することができる。
Jupyter Notebookを用いるために、Ubuntu/WSL/Windowsに適当なブラウザーをインストールする必要がある。コマンド
sudo apt –y install
firefox
を実行して、firefoxをインストールした。
また、jupyterも使用するPython 3.10仮想環境において、コマンド
conda install
jupyter
によりインストールした。
Jupyter Notebookでの実行例を、pdfファイルBinJN.pdfに示す。
リスト１．２項分布モデルの分析（１つのファイル：BinOneScrpt.py）
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az
import seaborn as sb
import numpy as np
N = 10
k = 7
with pm.Model() as BinModel:
theta = pm.Uniform('theta',
lower = 0.0, upper = 1.0)
k_obs = pm.Binomial('k_obs',
n = N, p = theta, observed = k)
trace = pm.sample()
pm.plot_trace(trace)
plt.savefig('FigTrace.png')
plt.show()
smry = pm.summary(trace)
print(smry)
with open('Summary.txt',
'w') as f:
f.write(smry.__str__())
sb.kdeplot(trace['theta'])
plt.xlabel(r'$\theta$', fontsize = 16)
plt.savefig('FigTheta.png')
plt.show()
print('FigTheta.png was saved.')
print('Summary.txt was saved.')
リスト２．2項分布モデルの分析：その１（MCMCサンプリング：Bin2_1.py）
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az
import pickle
N = 10
k = 7
with pm.Model() as BinModel:
theta = pm.Uniform('theta',
lower = 0.0, upper = 1.0)
k_obs = pm.Binomial('k_obs',
n = N, p = theta, observed = k)
trace = pm.sample()
smry = pm.summary(trace)
print(smry)
with open('Summary.txt',
'w') as f:
f.write(smry.__str__())
"""
MCMCサンプルの保存
"""
inference_data =
az.from_pymc3(trace = trace, model = BinModel)
with open('temp.pkl', 'wb')
as f:
pickle.dump((N, k, inference_data), f)
print('temp.pkl was saved.')
リスト３．2項分布モデルの分析：その２（サンプルの分析：Bin2_2.py）
import matplotlib.pyplot as plt
import arviz as az
import pickle
import seaborn as sb
import pickle
import numpy as np
#
# 保存されたデータの読み込み
#
with open('temp.pkl', 'rb') as f:
N, k, i_data =
pickle.load(f)
#
# MCMCサンプルの分析
#
smry = az.summary(i_data)
print('Summary =\n',smry)
with open('Summary.txt', 'w') as f:
f.write(smry.__str__())
print('Summary.txt was saved.')
az.plot_trace(i_data)
plt.savefig('FigTrace.png')
plt.show()
print('FigTrace.png was saved.')
v = np.array(i_data.posterior.theta)
v = np.reshape(v, len(v)*len(v.T))
sb.kdeplot(v)
plt.xlabel(r'$\theta$', fontsize = 16)
plt.savefig('FigTheta.png')
plt.show()
print('FigTheta.png was saved.')