信号検出理論のシミュレーションとベイズ分析
サンプルサイズが小さいときのベイズ分析の有用性
信号検出理論による分析において、サンプルサイズが小さい場合のsensitivity d’の推定への影響を調べてみた。モデルとして、もっとも簡単な１AFC（Yes-No）課題に対するガウス等分散モデルを取り上げた。ファイルは、scriptfiles.zipにまとめた。
１選択肢強制選択（One Alternative Forced Choice: 1AFC）課題
ノイズ刺激の感覚を
、シグナル刺激の感覚を
で表し、それぞれ次の正規分布に従うとする（Wickens, 2002; 岡本, 2025）。
感覚の次元の原点と単位を、次の条件を満たすものとして設定する（Wickens, 2002）。
平均値の差は
で表される（Green & Swets, 1988, p.96）。
等分散モデル（
）の下では、は
で表される（Green & Swets, 1988, p.96）。
このとき、次式が成り立つ（Hautus et al., 2022, p.7）。
シミュレーションのためのモデル
提示された刺激の感覚が、基準値
より大きいときシグナル刺激が提示されたと判断され、そうでなければノイズ刺激が提示されたと判断されるとする。基準値に等しい場合は確率が０であるので計算上では無視する。
シグナル刺激が提示されたときに、シグナル刺激であると判断する（Hit: H）確率は
で与えられ、ノイズ刺激が提示されたときにシグナル刺激であると判断する（False Alarm: F）は
で与えられる（Wickens, 2002、岡本,2025）。
シグナル刺激を提示する試行をn_trials試行、ノイズ刺激を提示する試行もn_trials試行行う実験のシミュレーションをn_simu回行うプログラムをリスト１（付録に掲載）のように作成した。シミュレーション総数はn_simu=1000であり、n_trialsは25試行と200試行である。25試行は、Hautus et al.(2022）の説明（p.4）における試行数である。これに対して、心理物理学実験では100試行以上の試行数は珍しくないのでn_trials=200の条件も設定した。
モデルのパラメータは、Hautus et al. (2022)のTable 1.5に示されているものを用い、次のように設定した。
以上の条件の下でリスト１のスクリプトを実行して、図１の結果を得た。
図１ 推定値d’（式（１）により算出）の分布。
１回の実験ごとに式（１）によりを求め、n_trials条件別にヒストグラムを描いたものである。図中nとあるのがn_trialsを表している。中央の黒の縦線は、パラメータ
の位置を示している。n_trials=25およびn_trials200のいずれの条件においてもは真値
を中心として分布しているが、n_trials=200の条件に対して、n_trials=25の条件では分布の広がりが大きい。
実験からを求めた場合、その真の値
は不明である。の推定値との関係は、n_trials=25の場合のように推定値の分散が大きい場合は不確かなものとなる。ベイズ分析を行うと、推定値と真値の関係を調べることができる。次にベイズ分析の場合のシミュレーションを行う。
ベイズ分析のシミュレーション
ベイズ分析（Gelman et al. 2014；岡本、2019）を行うためのStanスクリプト（CmdStan対応）を次のように用意した。
data {
int n;
int n_f;
int n_h;
}
parameters {
real<lower=-5,
upper=5> lmd;
real<lower=-5,
upper=5> d_prime;
}
model {
lmd ~ uniform(-5, 5);
d_prime ~ uniform(-5, 5);
n_f ~ binomial(n,
normal_cdf(-lmd | 0.0, 1.0));
n_h ~ binomial(n,
normal_cdf(-lmd+d_prime | 0.0, 1.0));
}
上のStanスクリプトを用いてベイズ分析を行うPythonスクリプト（CmdStanPy対応）をリスト３（付録に掲載）のように用意した。シミュレーションの条件は、上のリスト１によるを式（１）により算出するものと同じである。なお、CmdStanPyの簡単な説明のウェブサイトも用意している。
リスト３を実行して、図２と図３の結果を得た。
図２ 推定値d’（事後分布の中央値）の分布。
図３は、ベイズ分析で得られたの分布である。推定値は事後分布の中央値である（Gelman et al., 2021）。式（１）による推定（図１）とほぼ同じ結果である。
図２の表示Windowを閉じると、図３のグラフが表示される。
図３ 95％CIの幅の分布。med.は分布の中央値、
prop.は真値が95％CIに含まれた比率を表す。
真値
と推定値（事後分布の中央値）の関係を見るために、事後分布の95%CIの幅と真値が95%CIに含まれた比率を調べたものである。
グラフの横軸は、95%CIの幅を表している。n_trialが25のとき（n25で表記）は95%CIの幅は1.5ほどであり（中央値をmed.に示す）、真値
と比べると極めて広いと言える。しかし、これらのn_trials=25における95%CIに真値が含まれた比率は0.936（prop.として表記）であり、95%CIのcredibility 95%にほぼ対応している。
n_trialsが200のとき（n200で表記）のときは、95％CIの幅はn_trials=25のときの約１／３である。95%CI区間に真値が含まれた比率は、ほぼcredibility 95%に対応している。
サンプルサイズの大きさによりCIの幅は変動するが、それらが真値を含む頻度はcredibilityの値にほぼ対応していると考えられる。CIと真値との関係は、Okamoto (2022)で調べられているが参考になると思う。
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, reprint revision edition. Peninsula Publishing.
Hautus, M. J., Macmillan, N. A., & Creelman, C. D. (2022) Detection theory: A user’s guide, third edition. Routledge, Taylor & Francis group.
岡本安晴（2019） いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析（PyStan、PyMC）．丸善出版
Okamoto, Y. (2022) How Do Posterior Distributions Correspond to Facts? https://doi.org/10.17605/OSF.IO/AN72U
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房
Wickens, T. D. (2002) Elementary signal detection theory. Oxford University Press. (訳：岡本安晴「信号検出理論の基礎」2005、共同出版：電子書籍版)
リスト１ 式（１）による計算のシミュレーション（simu_dprime.py）
import scipy.stats as ss
import scipy.optimize as so
import numpy as np
import matplotlib.pyplot as plt
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
mu_n = 0.0
sgm_n = 1.0
mu_s = 1.0950
sgm_s = 1.0
lmd = 0.2533
rng = np.random.default_rng()
n_simu = 1000
d_prime_lst = [[],[]]
n_trials_lst = [25, 200]
for i, n_trials in
enumerate(n_trials_lst):
for i_simu in range(n_simu):
if
i_simu % 100 == 0:
print(f'\n{i_simu}/{i}\n')
x_n
= rng.normal(mu_n, sgm_n, size=n_trials)
x_s
= rng.normal(mu_s, sgm_s, size=n_trials)
n_h
= np.sum(x_s > lmd)
n_f
= np.sum(x_n > lmd)
if
n_h == n_trials:
n_h -= 1 / (2*n_trials)
elif
n_h == 0:
n_h += 1 / (2*n_trials)
if
n_f == n_trials:
n_f -= 1 / (2*n_trials)
elif
n_f == 0:
n_f += 1 / (2*n_trials)
d_prime = calc_z(n_h / n_trials) - calc_z(n_f / n_trials)
d_prime_lst[i].append(d_prime)
n25_005, n25_995 =
np.percentile(d_prime_lst[0], [0.5, 99.5])
n200_005, n200_995 = np.percentile(d_prime_lst[1],
[0.5, 99.5])
plt.hist(d_prime_lst[0], density=True,
color='g', alpha=0.7,
label=f'n={n_trials_lst[0]}' +
f'\n99%I=[{n25_005:.2f}, {n25_995:.2f}]')
plt.hist(d_prime_lst[1], density=True,
color='orange', alpha=0.5,
label=f'n={n_trials_lst[1]}' +
f'\n99%I=[{n200_005:.2f}, {n200_995:.2f}]')
plt.plot([mu_s, mu_s], [0, 3], c='k',
label=f"d'={mu_s}")
plt.xlabel("estimate d'",
fontsize= 16)
plt.title("Distributions of
estimate d'", fontsize=18)
plt.legend(fontsize=14)
plt.show()
リスト２ １AFC課題のStanスクリプト（d_prime_1afc.stan）
data {
int n;
int n_f;
int n_h;
}
parameters {
real<lower=-5,
upper=5> lmd;
real<lower=-5,
upper=5> d_prime;
}
model {
lmd ~ uniform(-5, 5);
d_prime ~ uniform(-5, 5);
n_f ~ binomial(n,
normal_cdf(-lmd | 0.0, 1.0));
n_h ~ binomial(n,
normal_cdf(-lmd+d_prime | 0.0, 1.0));
}
リスト３ ベイズ分析のシミュレーション（simu_dprime_ba.py）
import scipy.stats as ss
import scipy.optimize as so
import numpy as np
import matplotlib.pyplot as plt
from cmdstanpy import CmdStanModel
import time
mu_n = 0.0
sgm_n = 1.0
mu_s = 1.0950
sgm_s = 1.0
lmd = 0.2533
rng = np.random.default_rng()
n_trials_lst = [25, 200]
n_simu = 1000
d_prime_lst = [[],[]]
CI_p025_lst = [[],[]]
CI_p975_lst = [[],[]]
ck_CI_lst = [[],[]]
for i, n_trials in
enumerate(n_trials_lst):
for i_simu in range(n_simu):
if
i_simu % 100 == 0:
print(f'\n{i}
{i_simu}/{n_simu}\n')
ck_t = time.time();
while time.time() < ck_t+3:
pass
x_n
= rng.normal(mu_n, sgm_n, size=n_trials)
x_s
= rng.normal(mu_s, sgm_s, size=n_trials)
n_h
= np.sum(x_s > lmd)
n_f
= np.sum(x_n > lmd)
model = CmdStanModel(stan_file = 'd_prime_1afc.stan')
Data
= {'n':n_trials, 'n_f':n_f, 'n_h':n_h}
fit
= model.sample(data=Data, adapt_delta=0.99, max_treedepth=20,
show_progress=False)
fit_dframe = fit.draws_pd() # PandasのDataFrame型に変換
d_prime_lst[i].append(np.median(fit_dframe['d_prime']))
p025, p975 = np.percentile(fit_dframe['d_prime'], [2.5, 97.5])
CI_p025_lst[i].append(p025)
CI_p975_lst[i].append(p975)
if
(p025 <= mu_s) and (mu_s <= p975):
ck_ci = 1
else:
ck_ci = 0
ck_CI_lst[i].append(ck_ci)
plt.hist(d_prime_lst[0], density=True,
label='n=25', alpha=0.50)
plt.hist(d_prime_lst[1], density=True,
label='n200', alpha=0.50)
plt.plot([mu_s, mu_s], [0, 3],
color='k', label=f"d'={mu_s}")
plt.xlabel("estimate of d'",
fontsize=14)
plt.title("Distributions of
estimate d'", fontsize=16)
plt.legend()
plt.show()
w_CI_lst = [[],[]]
for i in range(2):
w_CI_lst[i].append(np.array(CI_p975_lst[i]) - np.array(CI_p025_lst[i]))
plt.hist(w_CI_lst[0], density=True,
label=f'n25: med.({np.median(w_CI_lst[0]):.3f})' +
f', prop.({np.mean(ck_CI_lst[0]):.3f})')
plt.hist(w_CI_lst[1], density=True,
label=f'n200: med.({np.median(w_CI_lst[1]):.3f})' +
f', prop.({np.mean(ck_CI_lst[1]):.3f})')
plt.title(f'Widths of 95%Intervals',
fontsize=14)
plt.xlabel('Width of 95% CI',
fontsize=14)
plt.legend(fontsize=14)
plt.show()