Up

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

 

信号検出理論に基づくデータ分析プログラムをStan(PyStan)を使ってPythonで作成してみた。信号検出理論については<岡本安晴「心理学データ分析と測定」勁草書、第12章>で説明している。等分散を仮定するモデルと仮定しないモデルのそれぞれについて作成した。

 

等分散モデル

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

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

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

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

 

表1

 

 

 

No

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

Yes

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

ノイズのみ

N(CR)

N(FA)

シグナル+ノイズ

N(Miss)

N(Hit)

 

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

上式に基づいて、2項分布を用いたStanのコードを以下のように用意した。

 

============ YN_EqSgm.stan ===================

data {

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

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

}

 

parameters {

    real mu;

    real<lower = 0> c;

}

 

transformed parameters {

    real pyes_n;

    real pyes_s;

    pyes_n = 1.0 - Phi(c);          //      Probability of FA

    pyes_s = 1.0 - Phi(c - mu);     //      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);

}

====================================

 

いま、表1の形式のデータとして具体的に表2に示すものを考える。

 

表2

 

 

 

No

Yes

ノイズのみ

69

31

シグナル+ノイズ

31

69

 

表2のデータに対して上のStanスクリプトを次のPythonスクリプトで実行した。「ノイズのみ」に対するデータはリストData_Noiseに、「シグナルを含む」に対するデータはリストData_Signalに設定している。

 

=====================================

import pystan

import matplotlib.pyplot as plt

import numpy as np

 

 

Data_Noise = [69, 31]       #   [Correct Rejection, False Alarm]

Data_Signal = [31, 69]      #   [Miss, Hit]

 

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

 

fit = pystan.stan(file = 'YN_EqSgm.stan', data = Data, iter = 11000, warmup = 1000, seed = 999, n_jobs = 1)

print(fit)

 

mu = fit['mu']

c = fit['c']

mu_s = np.sort(mu)

c_s = np.sort(c)

mu_med = mu_s[int(len(mu_s)/2)]

c_med = c_s[int(len(c_s)/2)]

 

plt.hist(mu)

plt.title("mu (d') = {0:<.3f} (med)".format(mu_med))

plt.show()

 

plt.hist(c)

plt.title("Criterion = {0:<.3f} (med)".format(c_med))

plt.show()

=====================================

 

上のプログラムを実行すると、MCMCによるサンプリング終了後、図1のグラフが表示される。

図1

 

図1は、パラメータの事後分布とその中央値を表示するものである。このフォームの右上のX印をクリックして閉じると、図2のフォームが表示される。

 

図2

 

図2は、判断の基準値の事後分布と中央値を示すものである。

図2を閉じるとプログラムの実行終了となる。

プログラムのファイルは圧縮ファイルSDT_EqSgm.zipとしてまとめた。ファイル名SDT_EqSgm.zipをクリックしてダウンロードすることができる。

 

 

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

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

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

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

である。

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

 

表3

 

 

 

 

カテゴリ1

・・・

カテゴリK

ノイズのみ

N(n1)

・・・

N(nK)

シグナルを含む

N(s1)

・・・

N(sK)

 

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

上式に基づいて多項分布を用いた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> sigma;

    ordered[K-1] c;

}

 

transformed parameters {

    vector[K] theta_n;

    vector[K] theta_s;

 

    theta_n[1] = Phi(c[1]);

    theta_n[K] = 1.0 - Phi(c[K-1]);

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

        theta_n[k] = Phi(c[k]) - Phi(c[k-1]);

   

    theta_s[1] = Phi((c[1] - mu) / sigma);

    theta_s[K] = 1.0 - Phi((c[K-1] - mu) / sigma);

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

        theta_s[k] = Phi((c[k] - mu) / sigma) - Phi((c[k-1] - mu) / sigma);

}

 

model {

    Freq_Noise ~ multinomial(theta_n);

    Freq_Signal ~ multinomial(theta_s);

}

=======================================

 

表3の形式のデータとして表4のものを用意した。

 

表4

 

 

 

 

 

カテゴリ1

「ノイズのみ」

カテゴリ2

「多分、

ノイズのみ」

カテゴリ3

「多分、

シグナル含む」

カテゴリ4

「シグナル含む」

ノイズのみ

10

9

7

1

シグナルあり

2

1

19

5

 

4のデータに対して以下のPythonスクリプトを用意した。「ノイズのみ」に対するデータをリストData_Noiseに、「シグナルを含む」に対するデータをリストData_Signalに設定している。

 

======================================

import pystan

import matplotlib.pyplot as plt

import numpy as np

import pickle

from pystan import StanModel

from CumNormalDistri import *

 

 

K = 4

Data_Noise = [10, 9, 7, 1]

Data_Signal = [2, 1, 19, 5]

 

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

 

fit = pystan.stan(file = 'YN_rating.stan', data = Data, iter = 10000, warmup = 1000, seed = 12345, n_jobs = 1)

 

print(fit)

 

mu = fit['mu']

# print('mu = \n', mu)

sigma = fit['sigma']

# print('sigma = \n', sigma)

mu_s = np.sort(mu)

sgm_s = np.sort(sigma)

mu_med = mu_s[int(len(mu_s) / 2)]

sgm_med = sgm_s[int(len(sgm_s) / 2)]

print('\nmu_med = ', mu_med)

print('\nsgm_sgm = ', sgm_med)

 

plt.hist(mu)

plt.title("mu (Delta-m) = {0:<.3f} (med)".format(mu_med))

plt.show()

 

plt.hist(sigma)

plt.title("Sigma = {0:<.3f} (med)".format(sgm_med))

plt.show()

 

 

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 = []

y_values = []

for z in z_values:

    x = 1.0 - cum_normal(z)

    x_values.append(x)

    y = 1.0 - cum_normal((z - mu_med) / sgm_med)

    y_values.append(y)

 

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

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

plt.title("ROC-curve")

plt.xlabel("P(False Alarm)")

plt.ylabel("P(Hit)")

 

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

 

c_s = []

c_med = []

for k in range(K-1):

    c_s.append([])

    c_s[k] = np.sort(c_comp[k])

    c_med.append(c_s[k][int(len(c_s[k]) / 2)])

 

print('\nCriterions = ')

for k in range(K-1):

    print('c-{0} = {1:<.3f} (med)'.format(k + 1, c_med[k]))

 

   

x_fa = []

y_ht = []

 

for k in range(K-1):

    x_fa.append(1.0 - cum_normal(c_med[k]))

    y_ht.append(1.0 - cum_normal((c_med[k] - mu_med) / sgm_med))

plt.plot(x_fa, y_ht, 'gs')

 

plt.show()

 

n_smpl = len(c_s[0])

 

min_c = c_s[0][0]

max_c = c_s[0][0]

for k in range(K-1):

    for i in range(n_smpl):

        if min_c > c_s[k][i]:

            min_c = c_s[k][i]

        if max_c < c_s[k][i]:

            max_c = c_s[k][i]

 

c_cnt = []

for k in range(K-1):

    c_cnt.append([])

    for i in range(20):

        c_cnt[k].append(0)

    for i in range(n_smpl):

        pos = 20.0 * (c_s[k][i] - min_c) / (max_c - min_c)

        pos = int(pos) if pos < 20.0 else 19

        c_cnt[k][pos] += 1

 

x_v = [0]*20

for i in range(20):

    x_v[i] = min_c + (max_c - min_c) * (i + 0.5) / 20.0

 

for k in range(K-1):

    plt.plot(x_v, c_cnt[k], 'b-')

title_str = "Criterions (medians)\nc{0}({1:<.2f})".format(1, c_med[0])

for k in range(1, K-1):

    title_str += "  c{0}({1:<.2f})".format(k+1, c_med[k])

plt.title(title_str)

 

plt.show()

======================================

 

上のプログラムを実行すると、図3のフォームが表示される。

図3

 

パラメータの事後分布とその中央値が表示されている。Delta-mのことである。

フォーム右上のX印をクリックしてフォームを閉じると、次に図4のフォームが表示される。

図4

 

図4には、パラメータの事後分布と中央値が表示されている。

図4のフォームを閉じると、次に図5のフォームが表示される。

図5

 

図5は、ROC曲線を描いたものである。赤丸はデータ点であり、緑の小正方形は理論値としての点である。

図5を閉じると、次は図6のフォームが表示される。

図6

 

判断の境界基準値の事後分布と中央値が表示されている。

図6を閉じるとプログラムの実行終了となる。

プログラムのファイルは圧縮ファイルSDT_NESgm.zipにまとめた。ファイル名SDT_NESgm.zipをクリックしてダウンロードすることができる。

 

 

Up