対応のあるデータの平均値のベイズ分析
2025.11：CmdStanPy対応に改訂
ベイズ分析のソフトとしてStanは良く知られている。CmdStanPyは、StanのPythonインターフェースである。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
対応のあるデータの平均値の差のベイズ分析を行う。データは
と対の形で得られるので、その差
が正規分布に従うとして分析する。すなわち、
とおき、事後分布を求める。
なお、岡本（2025）は、対応のあるデータの分析を回帰分析で行うことを勧めている。
いま、健康法の血圧に対する効果を見るために、ランダムに２０人を選び、健康法を始める前の血圧と健康法を１か月間行った後の血圧を測ると図１に示す結果を得たとする。
図１ Excelデータ（DiffData.xlsx）
図１のデータでは、１行目に変数名が書かれ、データは２行目以降に書かれている。1列名は、データの識別用で、図１では通し番号であるが、文字列でもよい。２列目と３列目に各対のデータ値
と
が書かれている。
図１の形式のデータを式（１）の確率モデルによりベイズ分析を行うため、モデル（１）に対して本ウェブサイト後半に示すリスト１のStanスクリプトを用意した。本ウェブサイトのスクリプトファイルなどは圧縮ファイルfilesDiffMeanRel.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトを使ってベイズ分析を行うPythonスクリプトを、本ウェブサイトの後半に示すリスト２のように用意した。
リスト１のStanスクリプトファイル（DifRelModel..stan）、リスト２のPythonスクリプトファイル（MainDiffRel.py）、およびデータファイル（例えば、DiffData.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境において、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesDiffMeanRel$ python MainDiffRel.py
上のコマンドを実行すると、図２の入力データファイル名の設定を求めるダイアログウィンドウが表示される。
図２
図２では、図１のデータDiffData.xlsxが設定されている。「Open」ボタンをクリックすると、次に画像保存用のファイル名の設定が求められる（図３）。
図３
ベイズ分析による事後分布のグラフ画像を保存するファイル名であるが、ファイル拡張子が「.png」であることを除いて任意である。図３では「MyFig.png」と設定されている。画像は複数枚保存されるが、図３で設定された名前から最後の４文字「.png」を画像に合わせて変更された名前で保存される。
画像保存用ファイル名を設定して「Save」ボタンをクリックすると、出力テキストファイル名の設定を求めるダイアログボックスが表示される（図４）。
図４
適当なテキストファイル名（*.txt）を設定する。
「Save」ボタンをクリックすると、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが行われ、終了後、トレース図が表示される（図５）。
図５
左側にパラメータの事後分布のグラフ、右側にMCMCサンプリングのトレース図が表示されている。
図５のWindowを、右上角のXアイコンをクリックして閉じると、図６のグラフが表示される。
図６
左のグラフは差の平均の事後分布である。事後分布が0の右側にあり、95％確信区間は0を含んでいない。差の平均は正の値であり、健康法を１カ月行うと血圧が下がったと言える結果である。
右のグラフは、差の標準偏差の事後分布のグラフである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
差の事後予測分布のグラフである。差のパラメータ、平均muと標準偏差sgm、の事後分布のサンプルに基づいて、差のサンプル値を次のスクリプトで求めている。
prdctn = []
for diff_mu,
diff_sgm in zip(fit['mu'], fit['sgm']):
v = ss.norm.rvs(diff_mu, diff_sgm)
prdctn.append(v)
prdctn = np.array(prdctn)
pp = np.mean(prdctn>0)
prdctn_med = np.median(prdctn)
予測事後サンプル値が正である事後予測確率は0.98であり、事後予測値の中央値は9.1である。
図７のWindowを閉じると、スクリプトの実行終了である。
リスト１ Stanスクリプト（DifRelModel.stan）
data {
int N;
array[N, 2] real X;
}
transformed data {
array[N] real Dif;
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 scipy.stats
as ss
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
import arviz
as az
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
Data File', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg
= 'Output Image File',
types = [('Image(*.png)', '*.png')],
deftype = '.png')
txtoutnm = askoutfilename(msg='Output
Text File',
types=[('Text(*.txt)',
'*.txt')],
deftype='.txt')
txtout = open(txtoutnm,
'w')
outfile_names = [txtoutnm]
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}
model = CmdStanModel(stan_file='DifRelModel.stan')
fit = model.sample(data=DataStan)
print(fit.summary())
txtout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_trace.png')
outfile_names.append(foutnm[:-4]
+ '_trace.png')
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
plt.figure(figsize=(12,5))
plt.subplot(121)
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.subplot(122)
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]
+ '_mu_sgm.png')
outfile_names.append(foutnm[:-4]
+ '_mu_sgm.png')
plt.show()
prdctn = []
for diff_mu,
diff_sgm in zip(fit['mu'], fit['sgm']):
v = ss.norm.rvs(diff_mu, diff_sgm)
prdctn.append(v)
prdctn = np.array(prdctn)
pp = np.mean(prdctn>0)
prdctn_med = np.median(prdctn)
sb.kdeplot(prdctn)
plt.title(f'予測事後確率（データ＞０）= {pp:.5f}' +
f'\n中央値.={prdctn_med:.3f}',
fontsize= 14)
plt.xlabel('Post Pred. P(Data>0)', fontsize=14)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_predp.png')
outfile_names.append(foutnm[:-4]
+ '_predp.png')
plt.show()
txtout.close()
for nm in outfile_names:
print(nm, 'was saved.')