１変量正規分布によるベイズ分析
CmdStanPy対応
StanのPythonインターフェースCmdStanPyの説明のウェブサイトを用意した。
１変量のデータ、例えば20人の学生の国語の点数が図１のように与えられているとする。
図１ Excelデータファイル（OneVar.xlsx）
図１のデータは、１行目に変数名が書かれている。１列目はデータの識別用で、図１では通し番号が設定されているが文字列でもよい。データ値は２列目に書かれている。
このデータを正規分布モデル
によってベイズ分析するためのStanスクリプトを本ウェブサイト後半に示すリスト１のように用意した。スクリプトファイルなどは圧縮ファイルfilesbysonevar.zipにまとめたので、ダウンロード解凍すれば自由に用いることができる。
リスト１のStanスクリプトを利用して図１のようなデータをベイズ分析するPythonスクリプトとして、本ウェブサイト後半に示すリスト２を用意した。
リスト１のStanスクリプトファイル（OneVarNormal.stan）、リスト２のPythonスクリプトファイル（MainOneVar.py）、入力データファイル（OneVar.xlsx）を同じフォルダに置き、そのフォルダ（ディレクトリ）にカレントディレクトリを移す。CmdStanPyのインストールされた環境にいて、リスト２の宿里オプトファイルMainOneVar.pyを実行する。例えば、以下のようである。
(stan) *****/filesbysonevar$
ls
MainOneVar.py
OneVar.xlsx OneVarNormal.stan
(stan) *****/filesbysonevar$
python MainOneVar.py
コマンド
python MainOneVar.py
を実行すると、まず入力データ名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
図２では、図１のExcelデータファイル名が設定されている。入力データファイルは図１の形式のExcelファイル（互換ソフトで作成されたものでもよい）を指定する。
入力データファイル名を設定して「Open」ボタンをクリックすると、次に画像ファイルの保存先ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
保存先ファイル名はファイル拡張子が「.png」であることを除いて任意である。図３では「MyFigure.png」が設定されている。画像ファイルは２枚保存され、一枚は設定通りの名前、もう一枚はファイル名の最後が「_sgm.png」に変更されたものである。
保存先ファイル名を設定して図３の「Save」ボタンをクリックすると計算が始まるが、Stanスクリプトのコンパイルにかなりの時間が掛かる。ビルド後にMCMCサンプリングが始まり、サンプリングが終了すると平均
の事後分布を描いたグラフが表示される（図４）。
図４
図４のWindowの右上のXアイコンをクリックして閉じると、標準偏差
の事後分布のグラフが表示される（図５）。
図５
図５のWindowを閉じるとプログラムの実行終了である。表示されたグラフは、ファイルに保存されている。
リスト１ 1変量正規分布用Stanスクリプト（OneVarNormal.stan）
data {
int N;
array[N] real X;
}
parameters {
real mu;
real<lower = 0.0> sgm;
}
model {
mu ~ normal(0.0, 10000.0);
sgm
~ exponential(0.0001);
for (i
in 1:N) {
X[i] ~ normal(mu, sgm);
}
}
リスト２ ベイズ分析Pythonスクリプト（MainOneVar.py）
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
def askinfilename(msg
= None, types = [('', '*.*')]):
"""
入力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.askopenfilename(title = msg, filetypes = types)
rt.destroy()
return flnm
def askoutfilename(msg
= None, types = [('', '*.*')], deftype = '.txt'):
"""
出力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.asksaveasfilename(title = msg, filetypes =
types,
defaultextension = deftype)
rt.destroy()
return flnm
finnm = askinfilename('Input
Excel File', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg
= 'Output Figure File',
types = [('Image(*.png)', '*.png')],
deftype = '.png')
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
VarName = data[0,1]
X = data[1:, 1]
X = list(X)
print(X)
DataStan = {'N': len(X),
'X': X}
model = CmdStanModel(stan_file='OneVarNormal.stan')
fit = model.sample(data=DataStan)
print(fit.summary())
fit = fit.draws_pd()
mu_025, mu_50, mu_975 = np.percentile(fit['mu'], [2.5, 50, 97.5])
plt.title('平均値の事後分布' + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
mu_50, mu_025, mu_975), fontsize = 20)
sb.kdeplot(fit['mu'])
plt.xlabel(VarName + 'の平均値', fontsize
= 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '.png')
plt.show()
sgm_025, sgm_50, sgm_975 = np.percentile(fit['sgm'], [2.5,
50, 97.5])
plt.title('標準偏差の事後分布' + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f},
{2:.2f}]'.format(
sgm_50, sgm_025, sgm_975), fontsize = 20)
sb.kdeplot(fit['sgm'])
plt.xlabel(VarName + 'の標準偏差', fontsize
= 14)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_sgm.png')
plt.show()