Home

信号検出理論

Signal Detection Theory (SDT)

１AFC評定課題において評定カテゴリ数は多い方が望ましいことを報告した（岡本安晴「検出力指標d′の推定における実験上の問題」基礎心理学研究）。評定カテゴリ数が７個の場合の１AFC評定課題実験のプログラム例を作成した。そのプログラムによる実験結果の分析例を示し、さらにシミュレーションによる検討を行った。

なお、信号検出理論の基礎的説明を（岡本安晴「心理学データ分析と測定」勁草書房）で行っている。

実験

実験プログラムをVisual C++/CLRで作成した。プログラムは、圧縮ファイルSDT_Exp7Cat.zipをダウンロード展開して作成されるフォルダ内のSDT_Exp.slnをVisual Studioで開いて自由に利用できる。Visual C++/CLRについては、＜岡本安晴「大学生のための心理学VC＋＋」、勁草書房＞で説明している。

提示刺激は、ランダムに一様分布するドット・パターンをノイズとし、その上に正規分布する少数のドット・パターンをシグナルとして用いた。刺激パターンの例は、実験用プログラムを起動すると提示される（図１）。

図１

「Noise ->」および「Signal ->」の右のテキストボックス内の数字はノイズ刺激におけるドット（一様分布）の数およびシグナルとして追加されるドット（正規分布）の数を表している。いま、分かりやすくするために、シグナルとして追加するドット数を300に設定して（図２）「Show samples」ボタンをクリックすると、図３のフォームが提示される。

図２

図３

図３のフォームの左側にノイズのみのパターン、右側にシグナルを追加したパターンが示されている。左と右では、ノイズパターンのドットの配置が異なることに注意。これは、ノイズパターンを作成するときの乱数が、左と右とで異なる系列が用いられているからである。図２のフォームの「Show samples」ボタンを再度クリックすると、別のサンプルが表示される。

ノイズとシグナルのドット数を適当に設定した後、「GO」ボタン（図１）をクリックすると実験が始まる。図１の「GO」ボタンのクリックで、図４のフォームが表示される。試行数は全部で420試行あり、その第１試行目であることが表示されている。

図４

図４の「GO」ボタンのクリックで、図５のようにドットパターンが表示される。

図５

その後、すぐに図６のフォームになる。「確かにシグナルあり」、「シグナルあり」、「たぶん、シグナルあり」、「？」、「たぶん、ノイズのみ」、「ノイズのみ」、「確かにノイズのみ」のボタンがイネーブルになっている。

図６

提示されたパターンがシグナルを含むものであると思うときは「シグナルあり」ボタンを、よくわからないが多分ノイズのみであろうと思うときは「たぶん、ノイズのみ」というように、提示パターンに対してもっともよくあてはまる判断を表すボタンをクリックする。判断のボタンをクリックすると、提示されたパターンが「シグナルを含むもの」であったか、「ノイズのみ」のものであったか、正解が図７のように示される（図７）。

図７

図７の画面の後、図４のフォーム（試行数は増える）が提示される。図４から図７までの試行を420試行繰り返して実験が終了すると、図８のフォームが表示される。

図８

実験についてのコメントなどをリッチ・テキスト・ボックス内に記入した後「Close」ボタンをクリックすると、図９の出力ファイル名設定用ダイアログボックスが表示される。

図９

適当なファイル名（テキストファイル名）を設定して「保存」ボタンをクリックすると、実験データがファイルに書き出されてプログラムの終了となる。データの書き出されたファイルは、次に説明するデータ分析用プログラムで分析することができる。

データは以下のような形式で書き出されている。

n_noise = 500 n_signal = 30

Comments:

Experimenter

trial Stimulus Response

1 N 7

2 N 1

3 S 1

4 S 5

5 N 4

6 S 5

/　　第１ブロック（練習試行）終了

7 S 3

8 S 5

9 N 4

10 N 2

11 N 3

・

・

・

415 S 5

416 N 2

417 N 3

418 N 4

419 S 5

420 S 5

/

Responseは「確かにノイズのみ」が１，「確かにシグナルあり」が７で表されている。

試行数は６×70試行である。ノイズのみ試行３試行とシグナルあり試行３試行の計６試行をランダマイズしたものを１ブロックとし、70ブロックよりなる。最初の１ブロックは練習試行として扱う。

なお、出力ファイルはs-jisコードで書き出されている。Pythonのエディタで開くときは、文字コードが確認されるので、s-jisを指定すればよい。

分析

上の実験で得られたデータをベイズ分析するスクリプトをリストA.1およびリストA.2のように用意した。リストA.1は、１AFC評定課題のモデルに基づくStanスクリプトであるが、このウェブサイトを参照されたい。カテゴリ数が４個の場合のモデルは、拙著pp.176-177で説明している。リストA.1のStanスクリプトを用いてベイズ分析をリストA.2のPythonスクリプトで行っている。PythonとStanについては、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。スクリプトファイルは圧縮ファイルSDT_Anal7Cat.zipとして用意しているが、ダウンロード展開すれば自由に利用できる。

リストA.2のスクリプトは、Anaconda/Windowsに用意したPython 3.7仮想環境にpystan2.17.1.0をインストールして実行した。

スクリプトファイルsdt_7cat_anal.pyを次のように実行すると、分析用データのファイル名の設定が求められるので、データファイル名を設定する（次の例ではResults.txt）。

(py37) PS D:\... \SDT_Anal7Cat> python sdt_7cat_anal.py

INFO:numexpr.utils:NumExpr defaulting to 8 threads.

Input data = Results.txt

入力データファイル名を設定すると、MCMCサンプリングが始まる。サンプリング後、まず(Δm)の事後分布のグラフが表示される（図A.1）。

図A.1

パラメータの点推定値としてmode (MAP推定値)を示しているのは、MAP推定値が点推定値として最も良いという報告に基づいている（岡本安晴、2022 https://doi.org/10.14947/psychono.41.1）。

図A.1のフォームを閉じると、シグナル刺激のσの事後分布が表示される（図A.2）。

図A.2

図A.2のフォームを閉じると、しばらくの計算時間の後、ROC曲線が表示される（図A.3）。

図A.3

図A.3のフォームを閉じると、判断基準のパラメータの事後分布が表示される。

図A.4

図A.4のフォームを閉じると、スクリプトの実行終了である。

スクリプトの実行終了後、出力ファイルOutput.txtを見ると、以下のようになっている。

Data for Noise:

[13 47 62 34 43 8 0]

Data for Signal:

[ 1 4 12 10 100 74 6]

mu:

MAP est. = 1.726

med. = 1.719

mean = 1.724

95% CI = [1.417, 2.062]

sigma:

MAP est. = 1.101

med. = 1.099

mean = 1.105

95% CI = [0.897, 1.356]

C1:

MAP est. = -1.489

median = -1.514

mean = -1.518

95% CI = [-1.791, -1.262]

C2:

MAP est. = -0.554

median = -0.552

mean = -0.553

95% CI = [-0.731, -0.378]

C3:

MAP est. = 0.198

median = 0.203

mean = 0.204

95% CI = [0.039, 0.369]

C4:

MAP est. = 0.612

median = 0.613

mean = 0.613

95% CI = [0.439, 0.787]

C5:

MAP est. = 1.983

median = 1.995

mean = 2.000

95% CI = [1.693, 2.335]

C6:

MAP est. = 3.751

median = 3.810

mean = 3.833

95% CI = [3.146, 4.647]

なお、StanによるMCMCサンプリングは乱数を用いるものであるので、分析を繰り返すごとに推定値は少し異なるものが得られる、

リストA.1　１AFC評定課題データのベイズ分析用Stanスクリプト（YN_rating.stan）

data {

int<lower = 3> K;

int<lower = 0> Freq_Noise[K];

int<lower = 0> Freq_Signal[K];

}

parameters {

real mu;

real<lower = 0.0, upper = 5.0> sigma;

ordered[K-1] c;

}

transformed parameters {

vector[K] theta_n;

vector[K] theta_s;

theta_n[1] = normal_cdf(c[1], 0.0, 1.0);

theta_n[K] = 1.0 - normal_cdf(c[K-1], 0.0, 1.0);

for (k in 2:(K-1))

theta_n[k] = normal_cdf(c[k], 0.0, 1.0) - normal_cdf(c[k-1], 0.0, 1.0);

theta_s[1] = normal_cdf((c[1] - mu) / sigma, 0.0, 1.0);

theta_s[K] = 1.0 - normal_cdf((c[K-1] - mu) / sigma, 0.0, 1.0);

for (k in 2:(K-1))

theta_s[k] =normal_cdf((c[k] - mu) / sigma, 0.0, 1.0) -

normal_cdf((c[k-1] - mu) / sigma, 0.0, 1.0);}

model {

mu ~ normal(0.0, 10.0);

sigma ~ uniform(0, 5);

Freq_Noise ~ multinomial(theta_n);

Freq_Signal ~ multinomial(theta_s);

}

リストA.2　１AFC評定課題データの分析用Pythonスクリプト（sdt_7cat_anal.py）

import pystan

import matplotlib.pyplot as plt

import numpy as np

from pystan import StanModel

import seaborn as sb

import scipy.stats as ss

with open(input('Input data = '), 'r') as f:

raw_data = f.readlines()

s_data = [ s.split() for s in raw_data]

i = 0

while True:

s = s_data[i]

if len(s) > 0:

if s[0] == '/':

i += 1

break

print(s)

i += 1

K = 7

Data_Noise = np.zeros(K, dtype = int)

Data_Signal = np.zeros(K, dtype = int)

while True:

s = s_data[i]

if s[0] == '/':

break

else:

if s[1] == 'N':

Data_Noise[int(s[2])-1] += 1

else:

Data_Signal[int(s[2])-1] += 1

i += 1

print('Data_Noise =', Data_Noise)

print('Data_Signal =', Data_Signal)

def CalcMAPEst(samples, a = 0.05, n_points = 2000):

"""

Calculatte a MAP estimate from a KDE graph on [Lp, Up]

Lp and Up are 100*a/2 and 100(1-a/2) percentile points of samples

"""

Lp, Up = np.percentile(samples, [100 * a/2, 100 * (1 - a/2)]) # import numpy as np

coord = np.linspace(Lp, Up, n_points)

est_pdf = ss.gaussian_kde(samples).pdf(coord) # import scipy.stats as ss

map_idx = np.argmax(est_pdf)

MAP_Est = coord[map_idx]

return MAP_Est, est_pdf[map_idx]

fout = open('Output.txt', 'w')

fout.write('\nData for Noise:\n')

fout.write(' {}\n'.format(Data_Noise))

fout.write('\nData for Signal:\n')

fout.write(' {}'.format(Data_Signal))

Data = {'K': K, 'Freq_Noise': Data_Noise, 'Freq_Signal': Data_Signal}

sm = pystan.StanModel(file = 'YN_rating.stan')

fit = sm.sampling(data = Data, iter = 10000, n_jobs = 1)

print(fit)

mu = fit['mu']

sigma = fit['sigma']

mu_mean = mu.mean()

mu_MAP = CalcMAPEst(mu)[0]

mu_q025, mu_med, mu_q975 = np.percentile(mu, [2.5, 50, 97.5]) # mu_s[int(len(mu_s) / 2)]

sgm_mean = sigma.mean()

sgm_MAP = CalcMAPEst(sigma)[0]

sgm_q025, sgm_med, sgm_q975 = np.percentile(sigma, [2.5, 50, 97.5]) # sgm_s[int(len(sgm_s) / 2)]

fout.write('\nmu:\n')

fout.write(' MAP est. = {0:.3f}\n'.format(mu_MAP))

fout.write(' med. = {0:.3f}\n'.format(mu_med))

fout.write(' mean = {0:.3f}\n'.format(mu_mean))

fout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(mu_q025, mu_q975))

fout.write('\nsigma:\n')

fout.write(' MAP est. = {0:.3f}\n'.format(sgm_MAP))

fout.write(' med. = {0:.3f}\n'.format(sgm_med))

fout.write(' mean = {0:.3f}\n'.format(sgm_mean))

fout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(sgm_q025, sgm_q975))

sb.kdeplot(mu)

plt.title((r'$\mu$') + (r'($\Delta$') + 'm) = ' +

'{0:<.3f}(mode)'.format(mu_MAP),

fontsize = 18)

plt.yticks([])

plt.ylabel('Density', fontsize = 14)

plt.xlabel('$\mu$', fontsize = 16)

plt.tight_layout()

plt.savefig('Figure_mu.png')

plt.show()

sb.kdeplot(sigma)

plt.title("$\sigma$ = {0:<.3f}(mode) ".format(sgm_MAP), fontsize = 18)

plt.yticks([])

plt.ylabel('Density', fontsize = 14)

plt.xlabel('$\sigma$', fontsize = 18)

plt.tight_layout()

plt.savefig('Figure_sgm.png')

plt.show()

print('\nI am calculating...')

c = fit['c']

c_comp = []

for i in range(3):

c_comp.append([])

for v in c:

for j in range(3):

c_comp[j].append(v[j])

z_values = np.arange(-4.0, 4.001, 0.1)

x_values = 1.0 - ss.norm(0.0, 1.0).cdf(z_values)

y_values = 1.0 - ss.norm(mu_MAP, sgm_MAP).cdf(z_values)

plt.figure(figsize = (5,5))

plt.plot(x_values, y_values, 'b-')

plt.title("ROC-curve", fontsize = 18)

plt.xlabel("P(False Alarm)", fontsize = 14)

plt.ylabel("P(Hit)", fontsize = 14)

cum_freq_n = []

cum_freq_s = []

for k in range(K):

cum_freq_n.append(Data_Noise[K - 1 - k])

cum_freq_s.append(Data_Signal[K - 1 - k])

for k in range(1, K):

cum_freq_n[k] += cum_freq_n[k-1]

cum_freq_s[k] += cum_freq_s[k-1]

x_fa = []

y_ht = []

for k in range(K-1):

x_fa.append(cum_freq_n[k] / cum_freq_n[K-1])

y_ht.append(cum_freq_s[k] / cum_freq_s[K-1])

plt.plot(x_fa, y_ht, 'ro', label = 'Data')

c_means = np.zeros(K-1)

c_meds = np.zeros(K-1)

c_MAPs = np.zeros(K-1)

c_q025s = np.zeros(K-1)

c_q975s = np.zeros(K-1)

for k in range(K-1):

c_means[k] = c.T[k].mean()

c_MAPs[k] = CalcMAPEst(c.T[k])[0]

c_q025s[k], c_meds[k], c_q975s[k] = np.percentile(c.T[k], [2.5, 50, 97.5])

fout.write('\n')

for k in range(K-1):

fout.write('\nC{}:\n'.format(k+1))

fout.write(' MAP est. = {0:.3f}\n'.format(c_MAPs[k]))

fout.write(' median = {0:.3f}\n'.format(c_meds[k]))

fout.write(' mean = {0:.3f}\n'.format(c_means[k]))

fout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(

c_q025s[k], c_q975s[k]))

x_fa = []

y_ht = []

for k in range(K-1):

x_fa.append(1.0 - ss.norm().cdf(c_MAPs[k]))

y_ht.append(1.0 - ss.norm(mu_MAP, sgm_MAP).cdf(c_MAPs[k]))

plt.plot(x_fa, y_ht, 'gs', label = 'Prediction')

plt.legend(fontsize = 12)

plt.tight_layout()

plt.savefig('Figure_ROC.png')

plt.show()

for k in range(K-1):

sb.kdeplot(c.T[k], label = '$C_{0}$: mode = {1:.3f}'.format(k+1, c_MAPs[k]))

title_str = "Criterions"

plt.title(title_str, fontsize = 18)

plt.yticks([])

plt.ylabel('Density', fontsize = 14)

plt.xlabel('C', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

plt.savefig('Figure_Cs.png')

plt.show()

fout.close()

シミュレーション

実験とデータ分析によって得られた推定値の特性をシミュレーションによって確かめる。

パラメータの推定値を、上の分析で得られたMAP推定値（モード）として、

, ,

とおく。

上の実験デザインに合わせて、ノイズ刺激およびシグナル刺激の提示回数を練習試行を除いた207試行とした。

このデザインでのシミュレーション実験を1000回行い、得られた1000個のパラメータの推定値の分布を調べるシミュレーションプログラムをリストS.2のように作成した。リストS.2におけるベイズ分析のStanスクリプトは、リストS.1に示すものである。スクリプトファイルは、圧縮ファイルSDT_Simu_7Cat.zipとして纏めたので、自由にダウンロード展開して利用されたい。リストS.2のスクリプトは、Anaconda/WindowsにおけるPython 3.7仮想環境においてPyStan 2.17.1.0をインストールして実行した。

リストS.2のスクリプトを実行すると、1000回のシミュレーション実験後、図S.1のグラフが表示される。

図S.1

パラメータの推定値1000個のバイオリン図である。横破線は四分位数の位置を示す。点推定値が事後分布の平均値（mean）、中央値（median）、MAP推定値（map）いずれの場合も、パラメータの真値は第1四分位数と第３四分位数の間にあり、四分位範囲（第３四分位数と第１四分位数の差）は約0.2である。MAP推定値の場合は、他の平均値および中央値と比べて、第２四分位数の位置が真値に近い。

図S.1のフォームを閉じると、図S.2のグラフが表示される。

図S.2

パラメータの推定値1000個のバイオリン図である。横破線は四分位数の位置を示す。点推定値が事後分布の平均値（mean）、中央値（median）、MAP推定値（map）いずれの場合も、パラメータの真値は第1四分位数と第３四分位数の間にあり、四分位範囲（第３四分位数と第１四分位数の差）は約0.25である。MAP推定値の場合は、他の平均値および中央値と比べて、第２四分位数の位置が真値に近い。

図S.2のフォームを閉じると、スクリプトの実行終了となる。

リストS.1　Stanスクリプト（SDT.stan）

data {

int K;

int ResNoise[K];

int ResSignal[K];

}

transformed data {

vector[K] a;

for (k in 1:K) {

a[k] = 1.0;

}

}

parameters {

real mu_s;

real<lower = 0.0> sgm_s;

simplex[K] preCs;

}

transformed parameters {

real C[K-1];

real v;

vector[K] theta_n;

vector[K] theta_s;

v = 0.0;

for (k in 1:(K-1)) {

v = v + preCs[k];

C[k] = -5 + v * 15.0; // -5 < C[k] < 10

}

theta_n[1] = normal_cdf(C[1], 0.0, 1.0);

theta_s[1] = normal_cdf(C[1], mu_s, sgm_s);

for (k in 2:(K-1)) {

theta_n[k] = normal_cdf(C[k], 0.0, 1.0) - normal_cdf(C[k-1], 0.0, 1.0);

theta_s[k] = normal_cdf(C[k], mu_s, sgm_s) - normal_cdf(C[k-1], mu_s, sgm_s);

}

theta_n[K] = 1.0 - normal_cdf(C[K-1], 0.0, 1.0);

theta_s[K] = 1.0 - normal_cdf(C[K-1], mu_s, sgm_s);

}

model {

preCs ~ dirichlet(a);

mu_s ~ normal(0.0, 1000.0);

sgm_s ~ exponential(0.001);

ResNoise ~ multinomial(theta_n);

ResSignal ~ multinomial(theta_s);

}

リストS.2　シミュレーションスクリプト（simu_n207C7.py）

import pystan

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats as ss

import seaborn as sb

import pickle

import tkinter as tk

sm = pystan.StanModel(file = 'SDT.stan')

with open('sm.pkl', 'wb') as f:

pickle.dump(sm, f)

print('OK')

def CalcMAPEst(samples, a = 0.05, n_points = 2000):

"""

Calculatte a MAP estimate from a KDE graph on [Lp, Up]

Lp and Up are 100*a/2 and 100(1-a/2) percentile points of samples

"""

Lp, Up = np.percentile(samples, [100 * a/2, 100 * (1 - a/2)]) # import numpy as np

coord = np.linspace(Lp, Up, n_points)

est_pdf = ss.gaussian_kde(samples).pdf(coord) # import scipy.stats as ss

map_idx = np.argmax(est_pdf)

MAP_Est = coord[map_idx]

return MAP_Est, est_pdf[map_idx]

n_simu = 1000

mu_means = []

mu_meds = []

mu_maps = []

sgm_means = []

sgm_meds = []

sgm_maps = []

s_sgm = 1.101

s_mu = 1.726

c = [-1.489, -0.554, 0.198, 0.612, 1.983, 3.751]

K = len(c) + 1

n_st = 207 # 3 * 69

def judge(v):

if v < c[0]:

return 0

elif v < c[1]:

return 1

elif v < c[2]:

return 2

elif v < c[3]:

return 3

elif v < c[4]:

return 4

elif v < c[5]:

return 5

else:

return 6

root = tk.Tk()

root.title(f'simu_n100C7.py')

cnvs = tk.Canvas(width = 700, height = 300, background = '#00ffff')

cnvs.pack()

cnvs.create_text(350, 150, text = '', font = ('', 100), fill = 'blue',

tags = 'MyText')

for t in range(n_simu):

#print(f't = {1+t}/{n_simu} ', end = '\r' )

cnvs.itemconfig('MyText', text = f'{t+1}/{n_simu}')

cnvs.update()

d_noise = np.zeros(K, dtype = int)

d_signal = np.zeros(K, dtype = int)

for t in range(n_st):

v = ss.norm.rvs(0, 1)

d_noise[judge(v)] += 1

v = ss.norm.rvs(s_mu, s_sgm)

d_signal[judge(v)] += 1

Data = {'K':K, 'ResNoise':d_noise, 'ResSignal':d_signal}

with open('sm.pkl', 'rb') as f:

sm = pickle.load(f)

fit = sm.sampling(data = Data, iter = 10000, n_jobs = 1)

mu_means.append(fit['mu_s'].mean())

mu_meds.append(np.median(fit['mu_s']))

mu_maps.append(CalcMAPEst(fit['mu_s'])[0])

sgm_means.append(fit['sgm_s'].mean())

sgm_meds.append(np.median(fit['sgm_s']))

sgm_maps.append(CalcMAPEst(fit['sgm_s'])[0])

root.destroy()

print('\nSimulation ended.')

Y = sgm_means + sgm_meds + sgm_maps

Grp = ['mean']*len(sgm_means) + ['median']*len(sgm_meds) + ['map']*len(sgm_maps)

BPData = {'X':Grp, 'Y':Y}

Xmin, Xmax, Ymin, Ymax = sb.violinplot(x = 'X', y = 'Y', data = BPData, inner = 'quartile').axis()

plt.hlines(y = s_sgm, xmin = Xmin, xmax = Xmax, label = f's_sgm({s_sgm})')

plt.legend()

plt.title(f'n_noise=n_signal={n_st}, 7 Categories\nDistributions of Estimates of s_sgm')

plt.tight_layout()

plt.savefig(f'FigSgms_n{n_st}C7.png')

plt.show()

Y = mu_means + mu_meds + mu_maps

Grp = ['mean']*len(mu_means) + ['median']*len(mu_meds) + ['map']*len(mu_maps)

BPData = {'X':Grp, 'Y':Y}

Xmin, Xmax, Ymin, Ymax = sb.violinplot(x = 'X', y = 'Y', data = BPData, inner = 'quartile').axis()

plt.hlines(y = s_mu, xmin = Xmin, xmax = Xmax, label = f's_mu({s_mu})')

plt.legend()

plt.title(f'n_noise=n_signal={n_st}, 7 Categories\nDistributions of Estimates of s_mu')

plt.tight_layout()

plt.savefig(f'FigMus_n{n_st}C7.png')

plt.show()

Home