対応のあるデータの平均値のベイズ分析
対応のあるデータの平均値の差のベイズ分析を行う。データは
と対の形で得られるので、その差
が正規分布に従うとして分析する。すなわち、
とおき、事後分布を求める。
いま、健康法の血圧に対する効果を見るために、ランダムに２０人を選び、健康法を始める前の血圧と健康法を１か月間行った後の血圧を測ると図１に示す結果を得たとする。
図１ Excelデータ（DiffData.xlsx）
図１のデータでは、１行目に変数名が書かれ、データは２行目以降に書かれている。1列名は、データの識別用で、図１では通し番号であるが、文字列でもよい。２列目と３列目に各対のデータ値
と
が書かれている。
図１の形式のデータを式（１）の確率モデルによりベイズ分析を行うため、モデル（１）に対して本ウェブサイト後半に示すリスト１のStanスクリプトを用意した。本ウェブサイトのスクリプトファイルなどは圧縮ファイルfilesDiffMeanRel.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトを使ってベイズ分析するPythonスクリプトを、本ウェブサイトの後半に示すリスト２のように用意した。
リスト２のプログラムを実行すると、まず、図２の入力データファイル名の設定を求めるダイアログウィンドウが表示される。
図２
図１のウィンドウを閉じると、次に画像保存用のファイル名の設定が求められる（図３）。
図３
ベイズ分析による事後分布のグラフ画像を保存するファイル名であるが、ファイル拡張子が「.png」であることを除いて任意である。図３では「MyFigure.png」と設定されている。画像は２枚保存されるが、１枚は設定ファイル名で保存され、他の１枚はファイル名の最後が「_sgm.png」に変更されて保存される。
画像保存用ファイル名を設定して「保存」ボタンをクリックすると、データの読み込みが始まり、Stanスクリプトのコンパイル・ビルドが行われる。このコンパイル・ビルドに多少時間が掛かるが、その後、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、まず、差の平均
の事後分布が表示される（図４）。
図４
図４の場合は、事後分布が0の右側にあり、95％確信区間は0を含んでいない。。差の平均は正の値であり、健康法を１カ月行うと血圧が下がったと言える結果である。
図４のウィンドウを閉じると、差の標準偏差
の事後分布が表示される（図５）。
図５
図５のウィンドウを閉じると、プログラムの実行終了である。
リスト１ Stanスクリプト（DifRelModel.stan）
data {
int
N;
real X[N, 2];
}
transformed data {
real Dif[N];
for (i
in 1:N)
Dif[i] = X[i][1]
- X[i][2];
}
parameters {
real mu;
real<lower = 0.0> sgm;
}
model {
mu ~ normal(0.0, 1000.0);
sgm
~ exponential(0.0001);
for (i
in 1:N)
Dif[i] ~ normal(mu, sgm);
}
リスト２ リスト１のStanスクリプトを用いたPythonスクリプト（MainDiffRel.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
VarName1 = data[0, 1]
VarName2 = data[0, 2]
print(VarName1, VarName2)
X = data[1:, 1:3]
X = [[v[0], v[1]] for v in X]
print(X)
DataStan = {'N': len(X),
'X': X}
sm = pystan.StanModel(file
= 'DifRelModel.stan')
fit = sm.sampling(data
= DataStan, n_jobs = 1)
print(fit)
# 日本語表記の設定
plt.rcParams['font.family']
= ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
mu_025, mu_50, mu_975 = np.percentile(fit['mu'], [2.5, 50, 97.5])
plt.title(('差({0} - {1})の平均の事後分布'.format(VarName1,
VarName2)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
mu_50, mu_025, mu_975), fontsize = 20)
sb.kdeplot(fit['mu'])
plt.xlabel('{0:} - {1:}の平均値'.format(VarName1, VarName2),
fontsize = 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm)
plt.show()
sgm_025, sgm_50, sgm_975 = np.percentile(fit['sgm'], [2.5,
50, 97.5])
plt.title(('差({0} - {1})の標準偏差の事後分布'.format(VarName1,
VarName2)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
sgm_50, sgm_025, sgm_975), fontsize = 20)
sb.kdeplot(fit['sgm'])
plt.xlabel('{0:} - {1:}の標準偏差'.format(VarName1,
VarName2), fontsize = 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_sgm.png')
plt.show()