信号検出理論
階層モデル
2025.11：CmdStanPy対応に改訂
複数の条件における信号検出実験のデータを、階層モデルの下にまとめてベイズ分析を行う。
信号検出理論では、ノイズ刺激とシグナルを含むシグナル刺激の２つの刺激の弁別が扱われる（岡本、2025）。ノイズ刺激およびシグナル刺激によって惹き起こされる感覚を
および
で表し、それぞれ平均が
および
、分散が
および
の正規分布に従うものとする。すなわち、
である。ノイズ刺激とシグナル刺激の感覚の分散が等しい
等分散モデルについて考える。
感覚の次元の単位と原点を決めるために、
とおく。
等分散モデルのもとでは、平均値の差は、
で表される。すなわち、
である。
シグナル刺激を正しく「シグナル刺激である」と答える反応はHitと呼ばれ、ノイズ刺激を誤って「シグナル刺激である」と答える反応はFalse Alarmと呼ばれている。感覚がある基準値
を越えたとき、「シグナル刺激である」と判断されるとすれば、HitとFalse Alarmの確率は、次式で与えられる。
いま、信号検出理論の実験を４条件、Ａ，Ｂ，Ｃ，Ｄ，で行い、表１の結果を得たとする。
|
表１ |
||||||||
|
|
条件A |
条件B |
条件C |
条件D |
||||
|
|
No |
Yes |
No |
Yes |
No |
yes |
No |
Yes |
|
ノイズ刺激 |
７７ |
２３ |
５０ |
５０ |
７７ |
２３ |
６０ |
４０ |
|
シグナル刺激 |
２３ |
７７ |
２３ |
７７ |
１０ |
９０ |
１１ |
８９ |
各条件においてノイズ刺激が100回、シグナル刺激が100回、ランダムに提示されたときのNo反応（「ノイズ刺激である」という判断）とYes反応（「シグナル刺激である」という判断）の頻度がまとめられている。条件Aを基準にして考えると、条件Bではシグナル刺激に対する反応の比率は条件Aと同じであるが、ノイズ刺激に対するYes反応が２３回から５０回に増えている。条件Cでは、ノイズ刺激に対する判断の比率は条件Aと同じであるが、シグナル刺激に対するYes反応が７７回から９０回に増えている。条件Dでは、条件Aと比べて、ノイズ刺激とシグナル刺激いずれに対してもYes反応が増えている。
表１のデータを等分散ガウス型モデルで分析するが、４つの条件をまとめて一括して分析する。
条件A、B、C、Dにおけるをda、db、dc、ddで表すとき、これらは共通の分布
に従うものとする。
条件A、B、C、Dにおける判断の基準値も、共通の分布
に従うものとする。
ハイパーパラメータ
、
、
、
は、次のhyperprior分布に従うものとする。
上の階層モデルをリスト１のようにStanスクリプトで表した。Stanについての説明は、岡本（2019）で行っている。
リスト１ 信号検出理論の階層モデルのStanスクリプト（sdt_h.stan）
data {
int NFAa;
int NCRa;
int NHa;
int NMa;
int NFAb;
int NCRb;
int NHb;
int NMb;
int NFAc;
int NCRc;
int NHc;
int NMc;
int NFAd;
int NCRd;
int NHd;
int NMd;
}
parameters {
real dmu_h;
real Cmu_h;
real<lower = 0> dsgm_h;
real<lower = 0> Csgm_h;
real<lower=0,
upper=10.0> da;
real Ca;
real<lower=0,
upper=10.0> db;
real Cb;
real<lower=0,
upper=10.0> dc;
real Cc;
real<lower=0,
upper=10.0> dd;
real Cd;
}
model {
dmu_h
~ normal(0, 5);
Cmu_h
~ normal(0, 5);
dsgm_h
~ exponential(0.5);
Csgm_h
~ exponential(0.5);
da ~ normal(dmu_h, dsgm_h);
db
~ normal(dmu_h, dsgm_h);
dc ~ normal(dmu_h, dsgm_h);
dd ~ normal(dmu_h, dsgm_h);
Ca ~ normal(Cmu_h, Csgm_h);
Cb
~ normal(Cmu_h, Csgm_h);
Cc ~ normal(Cmu_h, Csgm_h);
Cd ~ normal(Cmu_h, Csgm_h);
NFAa
~ binomial(NFAa+NCRa, 1.0 - normal_cdf(Ca
| 0.0, 1.0));
NHa
~ binomial(NHa+NMa, 1.0 - normal_cdf(Ca
| da, 1.0));
NFAb
~ binomial(NFAb+NCRb, 1.0 - normal_cdf(Cb | 0.0, 1.0));
NHb
~ binomial(NHb+NMb, 1.0 - normal_cdf(Cb | db, 1.0));
NFAc
~ binomial(NFAc+NCRc, 1.0 - normal_cdf(Cc
| 0.0, 1.0));
NHc
~ binomial(NHc+NMc, 1.0 - normal_cdf(Cc
| dc, 1.0));
NFAd
~ binomial(NFAd+NCRd, 1.0 - normal_cdf(Cd
| 0.0, 1.0));
NHd
~ binomial(NHd+NMd, 1.0 - normal_cdf(Cd
| dd, 1.0));
}
リスト１のStanスクリプトを用いて表１のデータを分析するPythonスクリプトを、リスト２のように用意した。これらのスクリプトファイルは圧縮ファイルfilesSDT_h.zipにまとめたので、ダウンロード・展開すれば自由に利用できる。
リスト２ 表１のデータを分析するスクリプト（sdt_eqvar_h.py）
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
from cmdstanpy
import CmdStanModel
model = CmdStanModel(stan_file='sdt_h.stan')
NFAa = 23; NCRa
= 77;
NHa = 77; NMa = 23;
NFAb = 50; NCRb
= 50;
NHb = 77; NMb = 23;
NFAc = 23; NCRc
= 77;
NHc = 90; NMc = 10;
NFAd = 40; NCRd
= 60;
NHd = 89; NMd = 11;
sdt_data = {'NFAa' :
NFAa, 'NCRa': NCRa,
'NHa': NHa, 'NMa': NMa,
'NFAb': NFAb, 'NCRb': NCRb,
'NHb': NHb, 'NMb': NMb,
'NFAc': NFAc, 'NCRc':NCRc,
'NHc': NHc, 'NMc': NMc,
'NFAd': NFAd, 'NCRd': NCRd,
'NHd': NHd, 'NMd':NMd}
fit = model.sample(data=sdt_data, adapt_delta=0.9999,max_treedepth=20,
iter_warmup=10_000)
print(fit.diagnose())
print(fit.summary())
inf_data = az.from_cmdstanpy(fit) # Arviz DataFrame型に変換
az.plot_trace(inf_data, figsize=(12,8),
var_names=['da','Ca','db','Cb','dc','Cc','dd','Cd'])
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型に変換
da_med = np.median(fit['da'])
db_med = np.median(fit['db'])
dc_med = np.median(fit['dc'])
dd_med = np.median(fit['dd'])
sb.kdeplot(fit['da'], label = 'da')
sb.kdeplot(fit['db'],
label = 'db')
sb.kdeplot(fit['dc'], label = 'dc')
sb.kdeplot(fit['dd'], label = 'dd')
plt.legend()
plt.title("Distributions of d'" +
'\nda_med = {0:.2f}, db_med
= {1:.2f}, dc_med = {2:.2f}, dd_med
= {3:.2f}'.format(
da_med, db_med, dc_med, dd_med))
plt.savefig('Figdprime.png')
plt.show()
Ca_med = np.median(fit['Ca'])
Cb_med = np.median(fit['Cb'])
Cc_med = np.median(fit['Cc'])
Cd_med = np.median(fit['Cd'])
sb.kdeplot(fit['Ca'], label = 'Ca')
sb.kdeplot(fit['Cb'],
label = 'Cb')
sb.kdeplot(fit['Cc'], label = 'Cc')
sb.kdeplot(fit['Cd'], label = 'Cd')
plt.legend()
plt.title('Distributions of C' +
'\nCa_med = {0:.2f}, Cb_med
= {1:.2f}, Cc_med = {2:.2f}, Cd_med
= {3:.2f}'.format(
Ca_med, Cb_med, Cc_med, Cd_med))
plt.savefig('FigC.png')
plt.show()
dmu_h_med = np.median(fit['dmu_h'])
Cmu_h_med = np.median(fit['Cmu_h'])
sb.kdeplot(fit['dmu_h'],
label = 'dmu_h')
sb.kdeplot(fit['Cmu_h'],
label = 'Cmu_h')
plt.legend()
plt.title('dmu_h_med
= {0:.2f}, Cmu_h_med = {1:.2f}'.format(dmu_h_med, Cmu_h_med))
plt.savefig('Figdmu_h.png')
plt.show()
dsgm_h_med = np.median(fit['dsgm_h'])
Csgm_h_med = np.median(fit['Csgm_h'])
sb.kdeplot(fit['dsgm_h'],
label = 'dsgm_h')
sb.kdeplot(fit['Csgm_h'],
label = 'Csgm_h')
plt.legend()
plt.title('dsgm_h_med
= {0:.2f}, Csgm_h_med = {1:.2f}'.format(dsgm_h_med, Csgm_h_med))
plt.savefig('Figdsgm_h.png')
plt.show()
plt.hist(fit['dsgm_h'],
density=True, label = 'dsgm_h', alpha = 0.3)
plt.hist(fit['Csgm_h'],
density=True, label = 'Csgm_h', alpha = 0.3)
plt.ylabel('Density', fontsize=14)
plt.legend()
plt.title('dsgm_h_med
= {0:.2f}, Csgm_h_med = {1:.2f}'.format(dsgm_h_med, Csgm_h_med))
plt.savefig('Figdsgm_h_hist.png')
plt.show()
リスト１のStanスクリプトファイル（sdt_h.stan）とリスト２のPythonスクリプトファイル（sdt_eqvar_h.py）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、次のコマンドをCmdStanPyのインストールされた環境において実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesSDT_h$ python sdt_eqvar_h.py
上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まり、終了後、図１のトレース図が表示される。
図１
左に事後分布のグラフが表示され、右にMCMCサンプリングのトレース図が表示されている。
図１のWindowを閉じると、図２のグラフが表示される。
図２
、すなわちda、db、dc、ddの事後分布のグラフと、事後分布の中央値が表示されている。Gelmanら（2021）は、事後分布の代表値として、平均値より安定した中央値を勧めている。
条件AとＤの事後分布が重なっており、条件Bより右側、条件Ｃより左側に位置している。条件Cでの弁別力d’が最も高いと言える。
図２のWindowを閉じると、図３のグラフが表示される。
図３
判断の基準値の事後分布のグラフと中央値が表示されている。条件Bの事後分布が最も左側にあり、その右に条件Dの事後分布が位置している。条件AとCの基準値の事後分布が重なって一番右側に位置している。条件Bでの判断基準が最も厳しいと言える。
図３のWindowを閉じると、図４のグラフが表示される。
図４
ハイパーパラメータdmu_hとCmu_hの事後分布のグラフと中央値が表示されている。基準値のハイパー分布の平均値は、弁別力のハイパー分布の平均値より負の側に寄っている。
図４のWindowを閉じると、図５のグラフが表示される。
図５
ハイパーパラメータdsgmとCsgmの事後分布のグラフと中央値が表示されている。これらのパラメータは負の値を取らないはずであるが、グラフが０より左側、負の領域にはみ出しているのは、滑らかなグラフを描くためにKDEを用いているからである。。
図５のフォームを閉じると、ヒストグラムのグラフが表示される（図６）。
図６
事後分布が正の領域に収まっていることが判る。
図６の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.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
岡本安晴（2025） 感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房