Up

 

信号検出理論(Signal Detection Theory: SDT

CmdStanPy対応

 

 

信号検出理論に基づくベイジアン・データ分析プログラムを作成してみた。信号検出理論については、岡本安晴(2025で説明している。等分散を仮定するモデル仮定しないモデルのそれぞれについて作成した。CmdStanPyの準備についての簡単な説明のウェブサイトも用意した。

 

 

 

等分散モデル

 

 

等分散モデルでは、シグナル刺激とノイズ刺激の感覚の分散は等しく1であると仮定される。この仮定の下で、ノイズとシグナルの感覚の平均値の差が次式

 

 

で算出される(Kingdom & Prins, 2016, p. 176; Macmillan & Creelman, 2005, p. 8)。

しかし、式(A1)の右辺は、岡本「心理学データ分析と測定」2014p.181)に示されているように、等分散の仮定が成り立たないときは、感覚以外の要因、すなわち判断基準を含むので感覚の差の指標として用いることができない。等分散について、本ウェブサイトの後半において若干の議論を行った。

 

等分散を仮定するモデルでは、ノイズに対する感覚量の平均値を0、分散を1とするとき、判断の基準値に対して、次式で確率が与えられる。

 

 

ここで、信号検出理論では、「シグナル+ノイズ」刺激に対する感覚量の平均値

 

 

で表すことが多い。「ノイズのみ」刺激に対する感覚量の平均が0であるので、は、「ノイズのみ」刺激と「シグナルを含む」刺激の感覚量の平均の差を表す。

ここで、等分散の仮定の下で

 

 

である。上式は、等分散モデルの下で成り立つが、一般には成り立たない(、2023p.107)。

 

等分散モデルのときにを用いるが、等分散を仮定しない場合は、が用いられる(岡本、2025

 

いま、「No|Noise」(Correct Rejection)の反応数をN(CR)で、というように表したとき、データが表A1のように得られたとする。

 

A

 

 

 

No

(「シグナルなし」判断)

Yes

(「シグナルあり」判断)

ノイズのみ

N(CR)

N(FA)

シグナル+ノイズ

N(Miss)

N(Hit)

 

このとき、尤度関数は次式で与えることができる。

 

 

上式に基づいて、2項分布を用いたStanのコードをリストA1のように用意し、リストA1Stanスクリプトを利用するPythonスクリプトをリストA2のように用意した。リストA1とリストA2のファイルは圧縮ファイルSDT_EqSgm.zipにまとめた。ダウンロード展開して自由に利用されたい。

 

いま、図A1のようなデータが与えられたとする。

 

A1

 

A1Excelデータファイル(Data.xlsx)をリストA1のスクリプトファイル(YN_EqSgm.stan)、リストA2のスクリプトファイル(main.py)と同じフォルダに置く。そのフォルダにカレントディレクトリ(フォルダ)を移して、次のコマンドをCmdStanPyのインストールされた環境で実行する。CmdStanPyの準備は、このウェブサイトに簡単な説明を用意した。

 

(stan) *****/SDT_EqSgm$ python main.py

 

上のコマンドを実行すると、出力ファイル名と入力データファイル名の設定が求められる。

 

(stan) *****/SDT_EqSgm$ python main.py

Outut file (*.txt) = Results.txt

Input DataFile (*.xlsx) = Data.xlsx

 

出力ファイルは、任意のテキストファイル名(*.txt)である。

入力データファイルは、図A1の形式のExcelファイル名(*.xlsx)である。

1行目に、反応カテゴリ名を設定する。2列目に「ノイズである」など、3列目に「シグナルである」などの順に置く。

1列目に提示刺激条件を選定する。2行目に「ノイズ刺激」、3行目に「シグナル刺激」である。

2行目、2列目以降に、対応する反応数を設定する。

1行目の反応カテゴリ名、1列目の提示刺激条件名は、プログラムの中では用いられないので任意であるが、2行目2列目以降の反応数は、この配列で処理されるので、並べ方の順序は図A1のものに合わせる。

書き込みが終わったら、ファイル拡張子を「.xlsx」として保存する。保存は、プログラムの実行前に行う。

 

ファイル名の設定が終わると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、多少の時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリングが終了すると、トレース図が表示される(図A2)。

 

A2

 

A2Windowを閉じると、図A3のグラフが表示される。

 

A3

 

すなわちの事後分布のグラフである。分布の代表値として中央値0.991が示されている。Gelmanら(2021は、代表値として、平均値より安定した中央値を勧めている。

A3Windowを閉じると、図A4のグラフが表示される。

 

A4

 

判断の基準値Cの事後分布のグラフである。

A4Windowを閉じると、図A5のグラフが表示される。

 

A5

 

事後分布の中央値を点推定値として、ノイズ刺激(橙色)、シグナル刺激(緑色)の感覚の確率分布のグラフ、および判断の基準値(赤色の垂直線)を描いたものである。垂直線の右側の領域が「Yes(シグナル刺激であると判断)」の確率を表し、左側の領域が「No(ノイズ刺激であると判断)」の確率を表している。

 

A5Windowを閉じると、実行終了である。

端末には、以下のようなメッセージが表示されている。

 

 

Results.txt was saved.

 

(stan) *****/SDT_EqSgm$

 

出力ファイルResults.txtを開くと、以下のような内容である。

 

 

Data for Noise:  [69 31]

Data for Signal: [31 69]

 

              Mean      MCSE    StdDev  ...  ESS_tail  ESS_bulk/s    R_hat

lp__   -125.573000  0.033158  1.084000  ...   1407.70     36770.7  1.00100

mu        0.987610  0.005508  0.185902  ...   1381.52     32709.3  1.00028

c         0.491365  0.004023  0.131404  ...   1066.06     30678.8  1.00107

pyes_n    0.313067  0.001424  0.046201  ...   1066.06     30678.8  1.00107

pyes_s    0.688662  0.000877  0.046037  ...   2524.60     81094.2  1.00031

 

[5 rows x 11 columns]

 

mu:

median   = 0.991

mean     = 0.988

95% CI   = [0.632, 1.354]

 

 

C:

median   = 0.492

mean     = 0.491

95% CI   = [0.237, 0.751]

 

 

 

 

リストA1 等分散SDTモデルのStanスクリプト(ファイル名YN_EqSgm.stan

 

data {

    array[2] int<lower = 0> Freq_Noise;       //  [Correct Rejection, False Alarm]

    array[2] int<lower = 0> Freq_Signal;      //  [Miss, Hit]

}

 

parameters {

    real mu;

    real<lower = 0> c;

}

 

transformed parameters {

    real pyes_n;

    real pyes_s;

    pyes_n = 1.0 - normal_cdf(c | 0.0, 1.0);    //      Probability of FA

    pyes_s = 1.0 - normal_cdf(c | mu, 1.0);     //      Probability of Hit

}

 

model {

    Freq_Noise[2] ~ binomial(Freq_Noise[1] + Freq_Noise[2], pyes_n);

    Freq_Signal[2] ~ binomial(Freq_Signal[1] + Freq_Signal[2], pyes_s);

}

 

 

 

 

リストA2 リストA1Stanスクリプトを利用するPythonスクリプト(ファイル名main.py

 

from cmdstanpy import CmdStanModel

import pandas as pd

import matplotlib.pyplot as plt

import numpy as np

import scipy.stats as ss

import seaborn as sb

import arviz as az

 

 

outfnm = input('Outut file (*.txt) = ')

fout = open(outfnm, 'w')

 

inflnm = input('Input DataFile (*.xlsx) = ')

df_data = pd.read_excel(inflnm)

print(df_data)

 

Data_Noise = df_data.values[0,1:]

print(Data_Noise)

 

Data_Signal = df_data.values[1,1:]

print(Data_Signal)

 

fout.write('Data for Noise:  {}\n'.format(Data_Noise))

fout.write('Data for Signal: {}\n'.format(Data_Signal))

 

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

 

model = CmdStanModel(stan_file="YN_EqSgm.stan")    #   コンパイル

fit = model.sample(data=Data)                      #   サンプリング

 

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

 

inf_data = az.from_cmdstanpy(fit)           #   Arviz(az)用のデータに変換

az.plot_trace(inf_data, figsize=(10,7))     #   トレースプロット

plt.tight_layout()

plt.show()

 

fit = fit.draws_pd()                        #   Pandas DataFrameへ変換

 

mu = fit['mu']

c = fit['c']

 

mu_mean = mu.mean()

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

 

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

fout.write('median   = {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))

          

sb.kdeplot(mu)

plt.title(fr"$\mu$(d') = {mu_med:.3f}(med)," +

          f"   95% CI = [{mu_q025:.3f}, {mu_q975:.3f}]", fontsize=16)

plt.yticks([])

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

plt.tight_layout()

plt.show()

 

c_mean = c.mean()

c_q025, c_med, c_q975 = np.percentile(c, [2.5, 50, 97.5])

 

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

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

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

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

        

sb.kdeplot(c)

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

plt.title(f"Criterion = {c_med:.3f}(med)," +

          f"  95% CI = [{c_q025:.3f}, {c_q975:.3f}]", fontsize=16)

plt.show()

 

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

x_min = -3

x_max = 3 + mu_med

xcoord = np.arange(x_min, x_max, 0.01)

y_noise = ss.norm.pdf(xcoord, 0.0, 1.0)

plt.plot(xcoord, y_noise, lw=3, c='orange', label='noise')

y_signal = ss.norm.pdf(xcoord, mu_med, 1.0)

plt.plot(xcoord, y_signal, lw=3, c='g', label='signal')

plt.plot([c_med, c_med], [0.0, 0.4], c='r', label='Criterion')

plt.plot([x_min, x_max], [0, 0], c='k', lw=2)

plt.title(f"C = {c_med:.3f},  mu(d') = {mu_med:.3f}", fontsize=18)

plt.legend(fontsize=14)

plt.tight_layout()

plt.show()

 

fout.close()

print('\n', outfnm, 'was saved.\n')

 

 

 

 

等分散を仮定しないモデル

 

 

 等分散を仮定しないモデルは不等分散モデルと呼ばれているが、これは等分散を仮定しないということである。データが等分散であっても、不等分散モデルによる分析は可能である。分散の推定値が等分散になるだけであり、問題は無い。等分散を仮定しないモデルによる分析を行う場合は、分析において分散の推定、すなわちROC直線の傾きの推定が安定するようにデータの収集を行う必要がある。ROC直線上でデータ点が散らばるようにするが、これは判断カテゴリの境界値の散らばりが大きくなるようにすることである。判断カテゴリ数を多く取ることが1つの方法である。3カテゴリの場合であっても、中央の「わからない」カテゴリの頻度が少なくならないように、「判りにくいときは『わかりません』と答えて下さい」などと教示しておくとよいと思われる。

 

 ノイズのみの刺激に対する感覚量の分散と「シグナル+ノイズ」の刺激に対する感覚量の分散が等しいことを仮定しない場合は、ノイズのみの刺激に対する感覚量の平均を0、分散を1とおき、「シグナル+ノイズ」刺激に対する感覚量の平均を、分散をとおいて、モデルが設定される。このとき、ROC曲線上のデータ点として2点以上必要であるが、評定法を用いることによって1回の実験から複数のデータ点を得ることができる。オーバーフィッティングを避けるためには、ROC曲線に対して3点以上が必要である。

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

 

 

「ノイズのみ」刺激に対する感覚量の平均0と「シグナルを含む」刺激に対する感覚量の平均との差はで表される。すなわち、

 

 

である(Swets, 1996, p. 71)。

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

 

B1

 

 

 

 

カテゴリ1

・・・

カテゴリK

ノイズのみ

N(n1)

・・・

N(nK)

シグナルを含む

N(s1)

・・・

N(sK)

 

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

 

 

上式に基づいて多項分布を用いたStanのスクリプトをリストBのように用意した。リストB1Stanスクリプトを用いるPythonスクリプトをリストBに用意した。

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

リストB1のファイル(YN_rating.stan)、リストB2のファイル(main.py)、および入力データファイルを同じフォルダ内に置く。そのフォルダにカレントディレクトリ(フォルダ)を移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyの準備の簡単な説明を、このウェブサイトに用意した。

 

(stan) *****/SDT_NESgm$ python main.py

 

上のコマンドを実行すると、出力ファイル名と入力データファイル名の設定が求められる。

 

(stan) *****/SDT_NESgm$ python main.py

Output file (*.txt) = Results.txt

Input file (*.xlsx) = Data.xlsx

 

出力ファイル名は、任意のテキストファイル名(*.txt)である。

入力データファイルは、Excelファイルを図B1のように用意する。

 

B1

 

1行目は、2列目に「ノイズ刺激」、3列目に「シグナル刺激」と書く。

1列目は、評定カテゴリを、7件法の場合は、「間違いなくノイズである」から「間違いなくシグナルである」までを順番に書く。

これらの1行目、および1列目に設定された文字列は、分析では用いられない。データを設定するための助けである。

2行目、2列目以降に、対応する提示刺激と評定カテゴリの組み合わせに対する観測度数を設定する。この部分は、列および行が、列が「ノイズ刺激」、「シグナル刺激」の順番、行が「間違いなくノイズである」から「間違いなくシグナルである」のように「ノイズ刺激の確信の強さ」から「シグナル刺激の確信の強さ」への程度の順番としてプログラムでは読み込まれるので、順序はランダムではない。表B1におけるデータの配置と図B1におけるデータの配置は、一方は他方を転置した配置になっているので注意。

データを図B1のように用意したら、プログラムの実行前に保存しておく。

 

ファイル名を設定すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリングが終了すると、トレース図が表示される(図B2)。

 

B2

 

B2Windowを閉じると、図B3のグラフが表示される。

 

B3

 

パラメータ(グラフではで表されている)、すなわちの事後分布のグラフである。代表値(点推定値)として、中央値が表示されている。Gelmanら(2021)は、事後分布の代表値として、平均値より安定した中央値を勧めている。

B3Windowを閉じると、図B4のグラフが表示される。

 

B4

 

シグナル刺激の感覚の標準偏差の事後分布のグラフである。

B4Windowを閉じると、図B5のグラフが表示される。

 

B5

 

ROC曲線である。赤丸がデータ、緑の正方形がモデルの予測である。データにモデルはよく適合していると言えそうである。

B5Windowを閉じると、図B6のグラフが表示される。

 

B6

 

判断のカテゴリ境界の事後分布のグラフである。

図6のWindowを閉じると、図B7のグラフが表示される。

 

B7

 

パラメータの事後分布の中央値を推定値としたときの、ノイズ刺激の感覚(橙色の曲線)、シグナル刺激の感覚(緑色の曲線)の確率分布と判断の境界値(赤色の垂直破線)を表示したものである。

 

B7Windowを閉じると、スクリプトの実行終了である。

 

端末には、以下のようなメッセージが表示されている。

 

 

Results.txt was saved.

 

(stan) *****/SDT_NESgm$

 

出力ファイルResults.txtを開くと、以下のようになっている。

 

 

 

Data for Noise:

    [6 18 34 31 8 3 0]

 

Data for Signal:

    [1 3 12 28 38 12 6]

                  Mean      MCSE    StdDev  ...  ESS_tail  ESS_bulk/s     R_hat

lp__       -313.473000  0.047884  2.005020  ...   2405.26     8793.55  1.000930

mu            1.321310  0.003423  0.193971  ...   3065.50    15749.20  0.999944

sigma         1.181660  0.003011  0.152225  ...   3056.14    12507.10  0.999564

c[1]         -1.564580  0.003565  0.191335  ...   2740.30    14162.60  1.001460

c[2]         -0.719186  0.002033  0.134831  ...   3100.48    21550.90  0.999817

c[3]          0.184657  0.001727  0.120042  ...   3430.33    23508.70  1.000610

c[4]          1.187230  0.002489  0.152542  ...   3296.62    18404.80  1.000400

c[5]          2.299720  0.004671  0.245982  ...   2943.09    13502.90  0.999927

c[6]          3.165970  0.007585  0.387495  ...   3036.90    12766.40  1.000040

theta_n[1]    0.062142  0.000418  0.022780  ...   2740.30    14162.60  1.001460

theta_n[2]    0.175867  0.000514  0.034997  ...   2895.25    22180.40  1.000070

theta_n[3]    0.334728  0.000661  0.042831  ...   2710.05    20241.40  1.000280

theta_n[4]    0.307002  0.000584  0.038670  ...   2765.37    21757.30  0.999950

theta_n[5]    0.107565  0.000395  0.025691  ...   3189.33    20935.60  1.000710

theta_n[6]    0.011231  0.000121  0.006610  ...   2883.48    14083.70  0.999995

theta_n[7]    0.001466  0.000035  0.001818  ...   3036.90    12766.40  0.999776

theta_s[1]    0.008665  0.000111  0.006514  ...   3106.48    19155.80  1.001250

theta_s[2]    0.034950  0.000189  0.013082  ...   3491.87    23418.00  0.999880

theta_s[3]    0.125233  0.000360  0.023924  ...   3499.25    21671.60  0.999812

theta_s[4]    0.286997  0.000566  0.037080  ...   3101.55    20758.90  1.000370

theta_s[5]    0.339537  0.000575  0.040930  ...   3047.62    24698.60  1.000730

theta_s[6]    0.142343  0.000499  0.033452  ...   3384.35    21845.90  1.000920

theta_s[7]    0.062276  0.000342  0.023866  ...   3649.65    24566.10  1.001120

 

[23 rows x 11 columns]

mu:

    med.     = 1.315

    mean     = 1.321

    95% CI   = [0.957, 1.707]

 

sigma:

    med.     = 1.174

    mean     = 1.182

    95% CI   = [0.916, 1.499]

 

 

C1:

    median   = -1.558

    mean     = -1.565

    95% CI   = [-1.950, -1.216]

 

C2:

    median   = -0.717

    mean     = -0.719

    95% CI   = [-0.985, -0.469]

 

C3:

    median   = 0.184

    mean     = 0.185

    95% CI   = [-0.052, 0.419]

 

C4:

    median   = 1.185

    mean     = 1.187

    95% CI   = [0.894, 1.494]

 

C5:

    median   = 2.294

    mean     = 2.300

    95% CI   = [1.853, 2.804]

 

C6:

    median   = 3.137

    mean     = 3.166

    95% CI   = [2.493, 3.983]

 

 

 

 

リストB1 評定法データのベイズ分析のためのStanスクリプト(YN_rating.stan

 

data {

    int<lower = 3> K;

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

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

}

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

}

 

 

 

リストB2 リストB1Stanスクリプトを用いるPythonスクリプト(main.py

 

from cmdstanpy import CmdStanModel

import pandas as pd

import matplotlib.pyplot as plt

import numpy as np

import seaborn as sb

import scipy.stats as ss

import arviz as az

 

 

outfnm = input('Output file (*.txt) = ')

fout = open(outfnm, 'w') 

 

inflnm = input('Input file (*.xlsx) = ')

df_data = pd.read_excel(inflnm)

Data_Noise = df_data.values.T[1]

Data_Signal = df_data.values.T[2]

K = len(Data_Noise)

 

print('Data_Noise =', Data_Noise)

print('Data_Signal =', Data_Signal)

print('K =', K)

 

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}

 

model = CmdStanModel(stan_file="YN_rating.stan")    #   コンパイル

fit = model.sample(data=Data)                       #   サンプリング

 

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

 

inf_data = az.from_cmdstanpy(fit)           #   Arviz(az)用のデータに変換

az.plot_trace(inf_data, var_names=['mu', 'sigma', 'c'], figsize=(10,7))     #   トレースプロット

plt.tight_layout()

plt.show()

 

fit = fit.draws_pd()                        #   Pandas DataFrameへ変換

 

mu = fit['mu']

sigma = fit['sigma']

mu_mean = mu.mean()

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

sgm_mean = sigma.mean()

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

 

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

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('    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}(med)'.format(mu_med),

          fontsize = 18)

plt.yticks([])

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

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

plt.tight_layout()

plt.show()

 

sb.kdeplot(sigma)

plt.title(r"$\sigma_s$ = {0:<.3f}(med) ".format(sgm_med), fontsize = 18)

plt.yticks([])

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

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

plt.tight_layout()

plt.show()

 

c = []

for k in range(K-1):

    c.append(fit[f'c[{k+1}]'])

c = np.array(c).T

 

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_med, sgm_med).cdf(z_values)

 

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

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_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_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('    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_meds[k])) 

    y_ht.append(1.0 - ss.norm(mu_med, sgm_med).cdf(c_meds[k])) 

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

 

plt.plot([], lw=0, label=fr"$\mu_s$={mu_med:.3f}" + "\n" +

                     fr"$\sigma_s$={sgm_med:.3f}")

plt.legend(fontsize = 12)

plt.tight_layout()

plt.show()

 

for k in range(K-1):

    sb.kdeplot(c.T[k], label = '$C_{0}$: med = {1:.3f}'.format(k+1, c_meds[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.show()

 

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

x_min = -3

x_max = 3 + mu_med

xcoord = np.arange(x_min, x_max, 0.01)

y_noise = ss.norm.pdf(xcoord, 0.0, 1.0)

plt.plot(xcoord, y_noise, lw=3, c='orange', label='noise')

y_signal = ss.norm.pdf(xcoord, mu_med, sgm_med)

plt.plot(xcoord, y_signal, lw=3, c='g', label='signal')

for k in range(K-1):

    plt.plot([c_meds[k], c_meds[k]], [0, 0.4], c='r', lw=2, ls='--'),

 

plt.plot([x_min, x_max], [0, 0], c='k', lw=2)

plt.plot([], ls='--', c='r', lw=2, label='Criterion')

plt.legend(fontsize=14)

plt.tight_layout()

plt.show()

 

fout.close()

print(f'\n{outfnm} was saved.\n')

 

 

 

 

等分散についての若干の議論

 

 

等分散の仮定が成り立たない場合については、例えば、Swets (1996, pp. 52-54)で論じられている。Wickens (2002)では、式(A1)が等分散ガウス型モデルにおけるものであることが明示されている。このWickens (2002)では、Hautus(2022)においてbasic bias measureとして提示されている次式

 

 

に、sensitivity が含まれていることが明示されている。Wickens2002)では、上式ので表されているが、次式で与えられている。

 

 

Hautusらのp. 31に挙げられている例では、Trainee 1は訓練後も同じcriterionを維持している()と説明されているが、これは上のバイアス基準を用いた場合である。sensitivity の値も変化しているので、刺激(いわゆるノイズのみの刺激)を基準にした座標(ノイズ刺激の平均値=0)では、ほぼからに変化している。訓練前はノイズ刺激の左側に設定されていた基準値が訓練後はノイズ刺激の右側に移動している。

 上の式(Hautus at al. (2022), eq.(2.1))の右辺を不等分散モデルの場合で表してみる。ここで、判断の基準値をHautusらに従ってで表し、シグナル刺激の平均値と標準偏差をで表すと、以下のように式が展開される。

 

 

すなわち、不等分散モデルにおいては、式(Hautus at al. (2022), eq.(2.1))の右辺には、判断の基準値の他に、シグナル刺激の標準偏差と平均値が含まれている。ノイズ刺激の平均値と標準偏差を基準としたとき、判断の基準値の他にシグナル刺激の平均値と標準偏差の要因も含まれた値である。ここで注意しなければならないことは、biascriterionの区別であろう。ノイズ刺激とシグナル刺激の中点に置かれたcriterionbiasの無いcriterionであると考えるならば、との差がbiasとなる。しかし、ノイズ刺激に対する正答率とシグナル刺激に対する正答率が等しくなる点biasのない点と考えると、不等分散モデルの下では以下のようになる。

 

 

これより、

 

 

となる。Biasが0の点は、シグナルの標準偏差が大きくなるとノイズ刺激に近付き、が小さくなるとシグナル刺激に近付く。

なお、Hautus, Macmillan & Creelman (2022) Detection theory: A users guide, third edition, p.99

A major accomplishment of detection theory is to abstract a measure of sensitivity, such as d, that does not depend on response criterion.

とあるが、d’がresponse criterionを含まないのは、上で説明したようにシグナル刺激の分散がノイズ刺激の分散に等しいときである。

なお、P(F)P(H)の組に対して、signal刺激の感覚の分布がの値に応じて変わることを示す図を、https://osf.io/5pzby/のファイルMuSigma.pdfに示した。

 

 

 

 

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian Data Analysis, third edition. CRC Press.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University

Hautus, Macmillan, & Creelman (2022) Detection theory: A users guide, third edition. Routledge

Kingdom, F.A.A. & Prins, N. (2016)  Psychophysics: A practical introduction, 2nd ed. Elsevier

Macmillan, N.A. & Creelman, C.D. (2005)  Detection Thoery: A users guide. Lawrence Erlbaum Associates

岡本安晴(2014)心理学データ分析と測定.勁草書房

岡本安晴(2023)検出力指標d′の推定における実験上の問題.基礎心理学研究、41107114

岡本安晴(2025)感覚・知覚測定法.和氣典二・重野純・村上郁也(編)感覚・知覚心理学ハンドブック 第三版(第2章)、誠信書房

Swets, J.A. (1996) Signal Detection Theory and ROC Analysis in Psychology and Diagnostics: Collected papers. Lawrence Erlbaum Assoiates.

Wickens, T.D. (2002) Elementary signal detection theory. Oxford University Press〈岡本安晴(訳)信号検出理論の基礎.共同出版(電子書籍、Yahoo)〉

 

 

 

Home