Home

等分散の仮定を置かないモデルによる１AFC評定課題データのベイズ分析

信号検出理論において、ノイズ刺激とシグナル刺激の感覚の分散が同じであるという等分散の仮定が、一般には成り立たないことが指摘されている（Bradyら、2023；Hautusら、2022；Swets、1996；Wixted、2020）。したがって、ノイズ刺激とシグナル刺激の感覚の分散が同じであるとする根拠が無い場合には、この等分散の仮定を置かないモデルに基づいてデータの分析を行わなければならない。等分散の仮定が成り立たないときのパラメータ推定値の問題については、岡本（2023）の報告が参考になると思われる。信号検出理論については、岡本（2014）で説明している。２AFC課題における同様の問題（等分散モデルを前提としない分析）を、Okamoto (2023) で論じた。

ノイズ刺激の感覚の分散とシグナル刺激の感覚の分散が等しいことを仮定しない場合は、ノイズ刺激に対する感覚量の平均を０、分散を１とおき、シグナル刺激の感覚の平均を、分散をとおいて、モデルが設定される（Wickens, 2002）。このとき、ＲＯＣ曲線上のデータ点として２点以上必要であるが、評定法を用いることによって１回の実験から複数のデータ点を得ることができる。

いま、評定カテゴリ数がＫ個であるとする。「ノイズのみである」の確信が最も強いカテゴリをカテゴリ１とし、「シグナルを含む」の確信が最も強いカテゴリをカテゴリＫとして、カテゴリ１からカテゴリＫまで順番に「シグナルを含む」とする確信の強い方向に評定カテゴリが設定されているとする。カテゴリkとカテゴリk+1の感覚の境界値を判断基準値で表す。このとき、以下のモデルを設定する。

ノイズ刺激の感覚の平均０とシグナル刺激の感覚の平均との差はで表される。すなわち、

である（Swets, 1996; 岡本、2014）。

いま、ノイズ刺激に対するカテゴリｋの評定数をN(nk)で、シグナル刺激に対するカテゴリｋの評定数をN(sk)で表すと、表１のようにまとめることができる。

表１
	カテゴリ１	・・・	カテゴリＫ
ノイズのみ	N(n1)	・・・	N(nK)
シグナルを含む	N(s1)	・・・	N(sK)

表１のデータに対して尤度関数を次のように設定する。

上式に基づいて多項分布によるベイズ分析を行うStanスクリプトをリスト１のように用意した。リスト１のStanスクリプトを用いるPythonスクリプトをリスト２に用意した。

リスト１とリスト２のファイルは圧縮ファイルanal1afcr.zipにまとめた。圧縮ファイルをダウンロード展開して自由に利用されたい。

Pythonスクリプトは、PyStan 3対応であるので、Windowsの場合は、UbuntuをWindowsにインストールして、UbuntuにインストールされたAnacondaのPython 3.11仮想環境においてPyStan 3をインストールして実行する。Ubuntuは、WSL (Windows Subsystem for Linux)としてWindowsにインストールできる。WSLなどは、このウェブサイトを参照されたい。

リスト２のスクリプトファイルmain.pyをAnaconda/Ubuntu/WSL/WindowsのPython3.11仮想環境において次のコマンドにより実行した。

(py311) ・・・/anal1afcr$ python main.py

実行すると、まず評定カテゴリ数Kの入力が次のように求められる。

(py311) ・・・/anal1afcr$ python main.py

K =

データとして次の表のものを入力する場合は、カテゴリ数Kは７である。

	ノイズ刺激	シグナル刺激
間違いなくノイズである	6	1
ノイズである	18	3
多分、ノイズである	34	12
わからない	31	28
多分、シグナルである	8	38
シグナルである	3	12
間違いなくシグナルである	0	6

カテゴリー数を入力すると、次にノイズ刺激に対する判断カテゴリーの頻度の入力が求められる。

(py311) ・・・/anal1afcr$ python main.py

K = 7

Data for Noise =

ノイズ刺激に対する判断カテゴリーの頻度を入力すると、次にシグナル刺激に対するカテゴリー判断の頻度の入力が求められる。

(py311) ・・・/anal1afcr$ python main.py

K = 7

Data for Noise = 6 18 34 31 8 3 0

Data for Signal =

シグナル刺激に対するカテゴリー判断の頻度を入力する。

(py311) ・・・/anal1afcr$ python main.py

K = 7

Data for Noise = 6 18 34 31 8 3 0

Data for Signal = 1 3 12 28 38 12 6

シグナル刺激に対する各カテゴリの判断頻度数を入力してEnterキーを押すと計算が始まる。

Stanによる計算が終了すると、まずトレース図が表示される（図１）。

図１

図１のフォームを閉じると、すなわちの事後分布のグラフが表示される（図２）。

図２

図２のフォームを閉じると、次にの事後分布のグラフが表示される（図３）。

図３

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

図４

赤の小円でデータ値、緑の小正方形でパラメータの点推定値（mode、すなわちMAP推定値）による予測値が表示されている。

図４のフォームを閉じると、カテゴリ境界の事後分布のグラフが表示される（図５）。

図５

図５のフォームを閉じると、プログラムの実行終了である。

プログラムの実行終了後、出力ファイルResults.txtを開くと以下のようになっている。

Data for Noise:

[6, 18, 34, 31, 8, 3, 0]

Data for Signal:

[1, 3, 12, 28, 38, 12, 6]

Summary...

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat

pre_mu 0.623 0.093 0.450 0.800 0.001 0.001 10673.0 12729.0 1.0

pre_sigma 0.555 0.071 0.426 0.686 0.001 0.000 10083.0 13165.0 1.0

ic[0] 0.179 0.027 0.126 0.229 0.000 0.000 17740.0 14491.0 1.0

ic[1] 0.151 0.028 0.101 0.205 0.000 0.000 17816.0 14953.0 1.0

ic[2] 0.213 0.027 0.165 0.266 0.000 0.000 18635.0 14380.0 1.0

ic[3] 0.213 0.026 0.161 0.259 0.000 0.000 16489.0 15162.0 1.0

ic[4] 0.140 0.020 0.105 0.178 0.000 0.000 23638.0 14828.0 1.0

ic[5] 0.053 0.013 0.031 0.077 0.000 0.000 19102.0 14536.0 1.0

ic[6] 0.051 0.016 0.023 0.083 0.000 0.000 9166.0 11369.0 1.0

mu 1.246 0.187 0.900 1.600 0.002 0.001 10673.0 12729.0 1.0

sigma 1.120 0.141 0.862 1.383 0.001 0.001 10083.0 13165.0 1.0

prec[0] 0.179 0.027 0.126 0.229 0.000 0.000 17740.0 14491.0 1.0

prec[1] 0.330 0.029 0.276 0.385 0.000 0.000 24354.0 15039.0 1.0

prec[2] 0.543 0.029 0.488 0.596 0.000 0.000 21614.0 15826.0 1.0

prec[3] 0.755 0.027 0.706 0.807 0.000 0.000 12802.0 14091.0 1.0

prec[4] 0.896 0.021 0.856 0.935 0.000 0.000 9561.0 13066.0 1.0

prec[5] 0.949 0.016 0.917 0.977 0.000 0.000 9166.0 11369.0 1.0

c[0] -1.537 0.189 -1.894 -1.187 0.001 0.001 17740.0 14491.0 1.0

c[1] -0.711 0.133 -0.962 -0.467 0.001 0.001 24354.0 15039.0 1.0

c[2] 0.171 0.117 -0.051 0.388 0.001 0.001 21614.0 15826.0 1.0

c[3] 1.133 0.147 0.854 1.406 0.001 0.001 12802.0 14091.0 1.0

c[4] 2.174 0.232 1.749 2.613 0.002 0.002 9561.0 13066.0 1.0

c[5] 2.967 0.347 2.332 3.620 0.004 0.003 9166.0 11369.0 1.0

theta_n[0] 0.065 0.024 0.026 0.111 0.000 0.000 17740.0 14491.0 1.0

theta_n[1] 0.175 0.035 0.111 0.240 0.000 0.000 20683.0 15020.0 1.0

theta_n[2] 0.327 0.040 0.252 0.404 0.000 0.000 18371.0 14510.0 1.0

theta_n[3] 0.301 0.038 0.232 0.374 0.000 0.000 19624.0 15816.0 1.0

theta_n[4] 0.114 0.025 0.068 0.163 0.000 0.000 15949.0 14643.0 1.0

theta_n[5] 0.015 0.008 0.003 0.029 0.000 0.000 10356.0 13412.0 1.0

theta_n[6] 0.002 0.003 0.000 0.007 0.000 0.000 9166.0 11369.0 1.0

theta_s[0] 0.008 0.006 0.000 0.019 0.000 0.000 19026.0 17264.0 1.0

theta_s[1] 0.034 0.013 0.011 0.058 0.000 0.000 22655.0 17466.0 1.0

theta_s[2] 0.128 0.024 0.085 0.176 0.000 0.000 22371.0 17769.0 1.0

theta_s[3] 0.291 0.037 0.224 0.362 0.000 0.000 21893.0 16492.0 1.0

theta_s[4] 0.335 0.041 0.257 0.411 0.000 0.000 28970.0 16668.0 1.0

theta_s[5] 0.139 0.033 0.079 0.199 0.000 0.000 21222.0 16057.0 1.0

theta_s[6] 0.065 0.025 0.023 0.112 0.000 0.000 25100.0 17887.0 1.0

mu:

MAP est. = 1.228

med. = 1.241

mean = 1.246

95% CI = [0.895, 1.631]

sigma:

MAP est. = 1.098

med. = 1.112

mean = 1.120

95% CI = [0.871, 1.423]

C1:

MAP est. = -1.520

median = -1.531

mean = -1.537

95% CI = [-1.923, -1.183]

C2:

MAP est. = -0.710

median = -0.711

mean = -0.711

95% CI = [-0.974, -0.454]

C3:

MAP est. = 0.167

median = 0.171

mean = 0.171

95% CI = [-0.059, 0.399]

C4:

MAP est. = 1.124

median = 1.130

mean = 1.133

95% CI = [0.854, 1.428]

C5:

MAP est. = 2.148

median = 2.166

mean = 2.174

95% CI = [1.745, 2.643]

C6:

MAP est. = 2.900

median = 2.948

mean = 2.967

95% CI = [2.349, 3.698]

参考文献

Brady,T.F.,Robinson,M.M.,Williams,J.R.,&Wixted,J.T.(2023). Measuring memory is harder than you think: How to avoid problematic measurement practices in memory research. Psychonomic Bulletin & Review, 30, 421-449.

Hautus, M. J., Macmillan, N. A., & Creelman, C. D (2022). Detection Theory: A user’s guide, third edition. Routledge.

岡本安晴（2023）検出力指標d′の推定における実験上の問題．基礎心理学研究、41、107-114. https://doi.org/10.14947/psychono.41.16

Okamoto, Y. (2023) Extended 2AFC Rating Task. https://osf.io/ts2eq

岡本安晴（2019）いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析（PyStan、PyMC）．丸善出版．

岡本安晴（2014）心理学データ分析と測定．勁草書房．

Swets, J. A. (1996). Signal detection theory and ROC analysis in psychology and diagnostics: Collected papers. Lawrence Erlbaum Associations, Publishers.

Wickens, T. D. （2002）. Elementary signal detection theory. Ox ford: Oxford University Press. （ウィッケンズ，T. D.　岡本安晴（訳）（2005）．信号検出理論の基礎　協同出版）

Wixted, J. T. (2020). The forgotten history of signal detection theory. Journal of Experimental Psychology: Learning, Memory, and Cognition, 2020, 46, 2, 201-233.

リスト１　SDT不等分散モデルのStanスクリプト（ファイル名YN_rating.stan）

data {

int<lower = 3> K;

array[K] int<lower = 0> Freq_Noise;

array[K] int<lower = 0> Freq_Signal;

}

transformed data {

vector[K] a;

for (i in 1:K) {

a[i] = 1.0;

}

}

parameters {

real pre_mu;

real<lower = 0.0> pre_sigma;

simplex[K] ic;

}

transformed parameters {

real mu;

real sigma;

vector[K-1] prec;

vector[K-1] c;

vector[K] theta_n;

vector[K] theta_s;

mu = pre_mu * 2.0;

sigma = pre_sigma * 2.0 + 0.01;

prec[1] = ic[1];

for (i in 2:K-1) {

prec[i] = prec[i-1] + ic[i];

}

for (i in 1:K-1) {

c[i] = logit(prec[i]);

}

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 {

ic ~ dirichlet(a);

pre_mu ~ normal(0.0, 10.0);

pre_sigma ~ exponential(0.1);

Freq_Noise ~ multinomial(theta_n);

Freq_Signal ~ multinomial(theta_s);

}

リスト２　リスト１のStanスクリプトを用いるPythonスクリプト（ファイル名main.py）

import matplotlib.pyplot as plt

import numpy as np

import stan

import seaborn as sb

import scipy.stats as ss

import arviz as az

import pickle

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

K = int(input('K = '))

s = input('Data for Noise = ')

Data_Noise = []

for v in s.split():

Data_Noise.append(int(v))

s = input('Data for Signal = ')

Data_Signal = []

for v in s.split():

Data_Signal.append(int(v))

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}

def f_init():

return {'pre_sigma': 0.5}

with open('YN_rating.stan', 'r') as fstan:

sm = stan.build(fstan.read(), data = Data)

fit = sm.sample(num_samples = 5000, num_chains=4, init=[f_init()]*4)

print(fit)

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

smry = az.summary(i_data)

print('Summary =\n',smry)

fout.write('\n\n')

fout.write('Summary...\n')

fout.write(smry.__str__())

fout.write('\n')

d_fm = fit.to_frame()

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]

az.plot_trace(i_data, var_names = ['mu', 'sigma', 'c'])

plt.show()

mu = d_fm['mu']

sigma = d_fm['sigma']

mu_mean = mu.mean()

mu_MAP = CalcMAPEst(mu)[0]

mu_q025, mu_med, mu_q975 = np.percentile(mu, [2.5, 50, 97.5])

sgm_mean = sigma.mean()

sgm_MAP = CalcMAPEst(sigma)[0]

sgm_q025, sgm_med, sgm_q975 = np.percentile(sigma, [2.5, 50, 97.5])

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()

c = []

for i in range(K-1):

c.append(d_fm[f'c.{i+1}'])

c = np.array(c).T

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')

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

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()

Home