１変量正規分布によるベイズ分析
Stanの使用
１変量のデータ、例えば20人の学生の国語の点数が図１のように与えられているとする。
図１ Excelデータファイル（OneVar.xlsx）
図１のデータは、１行目に変数名が書かれている。１列目はデータの識別用で、図１では通し番号が設定されているが文字列でもよい。データ値は２列目に書かれている。
このデータを正規分布モデル
によってベイズ分析するためのStanスクリプトを本ウェブサイト後半に示すリスト１のように用意した。スクリプトファイルなどは圧縮ファイルfilesbysonevar.zipにまとめたので、ダウンロード解凍すれば自由に用いることができる。
リスト１のStanスクリプトを利用して図１のようなデータをベイズ分析するPythonスクリプトとして、本ウェブサイト後半に示すリスト２を用意した。
リスト２のスクリプトを実行すると、まず入力データ名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
図２では、図１のExcelデータファイル名が設定されている。入力データファイルは図１の形式のExcelファイル（互換ソフトで作成されたものでもよい）を指定する。
入力データファイル名を設定して「開く」ボタンをクリックすると、次に画像ファイルの保存先ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
保存先ファイル名はファイル拡張子が「.png」であることを除いて任意である。図３では「MyFigure.png」が設定されている。画像ファイルは２枚保存され、一枚は設定通りの名前、もう一枚はファイル名の最後が「_sgm.png」に変更されたものである。
保存先ファイル名を設定して図３の「保存」ボタンをクリックすると計算が始まるが、Stanスクリプトのコンパイルとビルドにかなりの時間が掛かる。ビルド後にMCMCサンプリングが始まり、サンプリングが終了すると平均
の事後分布を描いたグラフが表示される（図４）。
図４
図４のウィンドウを閉じると、標準偏差
の事後分布のグラフが表示される（図５）。
図５
図５のウィンドウを閉じるとプログラムの実行終了である。表示されたグラフは、ファイルに保存されている。
リスト１ 1変量正規分布用Stanスクリプト（OneVarNormal.stan）
data {
int
N;
real X[N];
}
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
import pystan
import matplotlib.pyplot
as plt
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('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.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}
sm = pystan.StanModel(file
= 'OneVarNormal.stan')
fit = sm.sampling(data
= DataStan, n_jobs = 1)
print(fit)
mu_025, mu_50, mu_975 = np.percentile(fit['mu'], [2.5, 50, 97.5])
# 日本語表記の設定
plt.rcParams['font.family']
= ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
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()