１変量正規分布によるベイズ分析
Stanの使用
１変量のデータ、例えば20人の学生の国語の点数が図１のように与えられているとする。
図１ Excelデータファイル（OneVar.xlsx）
図１のデータは、１行目に変数名が書かれている。１列目はデータの識別用で、図１では通し番号が設定されているが文字列でもよい。データ値は２列目に書かれている。
このデータを正規分布モデル
によってベイズ分析するためのStanスクリプトを本ウェブサイト後半に示すリスト１のように用意した。スクリプトファイルなどは圧縮ファイルfilesbysonevar_nfd.zipにまとめたので、ダウンロード解凍すれば自由に用いることができる。
リスト１のStanスクリプトを利用して図１のようなデータをベイズ分析するPythonスクリプトとして、本ウェブサイト後半に示すリスト２を用意した。
リスト２のスクリプトを実行すると、まず入力データ名の設定が求められる（図２）。
図２
図２では、図１のExcelデータファイル名が設定されている。入力データファイルは図１の形式のExcelファイル（互換ソフトで作成されたものでもよい）を指定する。
入力データファイル名を設定すると、次に画像ファイルの保存先ファイル名の設定が求められる（図２）。保存先ファイル名はファイル拡張子が「.png」であることを除いて任意である。図２では「MyFigure.png」が設定されている。画像ファイルは２枚保存され、一枚は設定通りの名前、もう一枚はファイル名の最後が「_sgm.png」に変更されたもので保存される。
保存先ファイル名を設定して「Enter」キーを押すと、データが読み込まれて計算が始まる。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 pandas as pd
import numpy
as np
import pystan
import matplotlib.pyplot
as plt
import seaborn
as sb
finnm = input('入力データファイル（*.xlsx） = ')
foutnm = input('画像保存ファイル（*.png）= ')
savedfiles = []
# 保存ファイル名のリスト
#
# Excelファイルの読み込み
#
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')
savedfiles.append(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')
savedfiles.append(foutnm[:-4]
+ '_sgm.png')
plt.show()
print('保存ファイル：')
for v in savedfiles:
print(v)