複数の判断基準を用いるYes-No（１AFC）課題実験のベイズ分析
信号検出理論（Signal Detection Theory: Green & Swets, 1988; Wickens, 2002）の実験データからシグナル刺激に対するsensitivityを算出するとき、分散が不明であれば、複数の判断基準におけるデータが必要である（岡本、2014）。図１に示すデータは、4つの判断基準条件、C1、C2、C3、C4、においてそれぞれノイズ刺激を１００試行、シグナル刺激を１００試行提示したものである（仮想データ）。「fa」はFlse Alarm、「cr」はCorrect Rejection、「hit」はHit、「miss」はMissの回数である。
図１ sample.xlsx
図１のデータをベイズ分析（Gelman et al., 2014；岡本、2019）するStanスクリプトを以下のように用意した。
data {
int n_c;
// Number of
conditions
array[n_c]
int n_f; // Number of false alarms
array[n_c]
int n_cr; // Number of correct rejections
array[n_c]
int n_h; // Number of hits
array[n_c]
int n_miss; // Number of misses
}
parameters {
real mu_s;
real<lower=0.0001> sgm_s;
array[n_c]
real c;
}
transformed parameters {
array[n_c]
real p_f;
array[n_c]
real p_h;
for (i
in 1:n_c) {
p_f[i] = normal_cdf(-c[i] | 0, 1);
p_h[i] = normal_cdf((-c[i] + mu_s) / sgm_s
| 0, 1);
}
}
model {
mu_s
~ uniform(-5, 10);
sgm_s
~ uniform(0.0001, 10);
for (i
in 1:n_c) {
c[i] ~ uniform(-10, 10);
}
for (i
in 1:n_c) {
n_f[i] ~ binomial(n_f[i]+n_cr[i], p_f[i]);
n_h[i] ~ binomial(n_h[i]+n_miss[i], p_h[i]);
}
}
上のStanスクリプトにより図１の形式のデータを分析するPythonスクリプトをリスト２のように用意した。ファイルは、multiCfiles.zipにまとめた。
上のStanスクリプト（リスト１にも示した）のファイル、リスト２のPythonスクリプトファイル、および入力データファイル（例えば、図１のファイル）を同じフォルダに置き、リスト２のスクリプトをCmdStanPyのインストールされた環境で実行する。CmdStanPyの簡単な説明のウェブサイトも用意している。
リスト２のスクリプトを実行すると、入力データファイル名を聞いてくる。
(stan) ****/multiCfiles$
python est_param.py
Data file (*.xlsx) = sample.xlsx
上の例では、図１のファイル名sample.xlsxが設定されている。
入力データファイルは、図１の形式のExcelファイルであればよい。１行目に、実験条件のラベルを書く。実験条件は２個以上であればよいが、３個以上が望ましい。１列目は反応のカテゴリを書く。名前は任意であるが、False Alarm、Correct Rejection、Hit、Missの順とする。
入力データファイル名を設定すると、ファイルが読み込まれ、Stanスクリプトのビルド、ビルドに続いてMCMCサンプリングが始まる。MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、ROC曲線の事後予測が計算されて、図３のグラフが表示される。
図３
事後分布の中央値をパラメータの推定値（Gelman et al., 2021）としたROC曲線が青の実線で描かれている。MCMCサンプルから50組をランダムに選び、それらの値により描いたROC曲線が橙色の細い曲線で描かれている。データは、赤色の小円で、条件のラベルとともに描かれている。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータmu_sの事後分布である。中央値と95％CIがタイトル欄に表示されている。
図４のWindowを閉じると、図５のグラフが表示される。
図５
パラメータsgm_sの事後分布のグラフである。タイトル欄に、中央値と95％CIが表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
判断の基準値Cの事後分布である。各判断の基準値の事後分布の中央値がmed.(中央値)として表示されている。
図６のWindowを閉じると、スクリプトの実行終了である。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press. (free download)
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.
Green, D. M. & Swets, J. A. (1988). Signal detection theory and psychophysics. Peninsula Publishing.
岡本安晴（2014） 心理学データ分析と測定．勁草書房（電子書籍版あり）
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房
Wickens, T. D. (2002) Elementary signal detection theory. Oxford University Press. (訳：岡本安晴「信号検出理論の基礎」2005、共同出版：電子書籍版)
リスト１ Stanスクリプト（multiC1afc.stan）
data {
int n_c;
// Number of
conditions
array[n_c]
int n_f; // Number of false alarms
array[n_c]
int n_cr; // Number of correct rejections
array[n_c]
int n_h; // Number of hits
array[n_c]
int n_miss; // Number of misses
}
parameters {
real mu_s;
real<lower=0.0001> sgm_s;
array[n_c]
real c;
}
transformed parameters {
array[n_c]
real p_f;
array[n_c]
real p_h;
for (i
in 1:n_c) {
p_f[i] = normal_cdf(-c[i] | 0, 1);
p_h[i] = normal_cdf((-c[i] + mu_s) / sgm_s
| 0, 1);
}
}
model {
mu_s
~ uniform(-5, 10);
sgm_s
~ uniform(0.0001, 10);
for (i
in 1:n_c) {
c[i] ~ uniform(-10, 10);
}
for (i
in 1:n_c) {
n_f[i] ~ binomial(n_f[i]+n_cr[i], p_f[i]);
n_h[i] ~ binomial(n_h[i]+n_miss[i], p_h[i]);
}
}
リスト２ Pythonスクリプト（est_param.py）
import pandas as pd
import numpy
as np
import scipy.stats
as ss
from cmdstanpy
import CmdStanModel
import arviz
as az
import matplotlib.pyplot
as plt
import seaborn as sb
f = input('Data file (*.xlsx) =
')
dframe = pd.read_excel(f)
print(dframe)
conds = dframe.keys()[1:]
values = dframe.values[:,1:]
print(conds)
print(values)
sm = CmdStanModel(stan_file = 'multiC1afc.stan')
n_f = values[0]
n_cr = values[1]
n_h = values[2]
n_miss = values[3]
Data = {'n_c':len(conds), 'n_f':values[0], 'n_cr':values[1],
'n_h':values[2], 'n_miss':values[3]}
fit = sm.sample(data=Data)
fit.diagnose()
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# ArvizのInferenceData型に変換
az.plot_trace(inf_data, figsize=(10,8))
plt.tight_layout()
plt.show()
fit_dframe = fit.draws_pd()
# PandasのDataFrame型に変換
mu_s_hat = np.median(fit_dframe['mu_s'])
sgm_s_hat = np.median(fit_dframe['sgm_s'])
def my_cdf(z):
if z < -9:
return 0.0
elif
z > 9:
return 1.0
else:
return ss.norm.cdf(z)
class Calc_p:
def __init__(self,
p):
self.p = p
def f(self, z):
diff
= my_cdf(z) - self.p
return diff
def calc_z(p):
calc_p
= Calc_p(p)
z = so.brentq(calc_p.f, -10, 10)
return z
d_f = n_f / (n_f + n_cr)
d_h = n_h / (n_h + n_miss)
plt.figure(figsize=(7,
7))
rng = np.random.default_rng()
idx_rnd = rng.choice(range(len(fit_dframe['mu_s'])), size=50, replace=False)
for i in idx_rnd:
mu_v
= fit_dframe['mu_s'][i]
sgm_v
= fit_dframe['sgm_s'][i]
xcoord
= np.linspace(-3, 3, 1000)
t_pf
= []
t_ph
= []
for x in xcoord:
t_pf.append(my_cdf(-x))
t_ph.append(my_cdf((-x + mu_v) / sgm_v))
plt.plot(t_pf, t_ph, lw=1,
c='orange', alpha=0.5)
plt.plot([], lw=1,
c='orange', label='post. distri')
xcoord = np.linspace(-3,
3, 1000)
pf = []
ph = []
for x in xcoord:
pf.append(my_cdf(-x))
ph.append(my_cdf((-x + mu_s_hat) / sgm_s_hat))
plt.plot(pf, ph,
c='g', lw=3, label='med. param.')
plt.plot(d_f, d_h, 'o', c='r', label='data')
for i, cs in
enumerate(conds):
plt.text(d_f[i], d_h[i], conds[i],
ha='left', va='top', fontsize=14)
plt.xlabel('P(FA)', fontsize=14)
plt.ylabel('P(Hit)', fontsize=14)
plt.title(r'$\widehat{\mu}_{s}$'+f'={mu_s_hat:.3f}'
+ ', ' +
r'$\widehat{\sigma}_{s}$' +
f'={sgm_s_hat:.3f}', fontsize=18)
plt.legend(fontsize=16)
plt.show()
p025, p50, p975 = np.percentile(fit_dframe['mu_s'], [2.5, 50,
97.5])
sb.kdeplot(fit_dframe['mu_s'])
plt.xlabel(r'$\widehat{\mu}_s$',
fontsize=16)
plt.title(f'Med.={p50:.3f}, 95%CI=[{p025:.3f}, {p975:.3f}]', fontsize=16)
plt.show()
p025, p50, p975 = np.percentile(fit_dframe['sgm_s'], [2.5, 50,
97.5])
sb.kdeplot(fit_dframe['sgm_s'])
plt.xlabel(r'$\widehat{\sigma}_s$',
fontsize=16)
plt.title(f'Med.={p50:.3f}, 95%CI=[{p025:.3f}, {p975:.3f}]', fontsize=16)
plt.show()
for i, cs in
enumerate(conds):
v_med
= np.median(fit_dframe[f'c[{i+1}]'])
sb.kdeplot(fit_dframe[f'c[{i+1}]'], label=conds[i] +
f': med.({v_med:.3f})')
plt.xlabel('c', fontsize=14)
plt.legend(fontsize=12)
plt.show()