Up

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

注意:ブラウザーの表示が乱れる、ダウンロードしたファイルの中身が以前の古いものであるなどの場合は、ブラウザーの履歴をクリアしてから、再読み込みを行って下さい。Microsoft Edgeの場合は、「履歴|閲覧データをクリア」の操作です。

 

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

 

 

等分散モデル

 

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

で算出される(Kingdom & Prins, Psychophysics: A practical introduction, 2nd ed. 2016, p. 176; Macmillan & Creelman, Detection Thoery: A users guide. 2005, p. 8)。

しかし、式(1)の右辺は、岡本「心理学データ分析と測定」2014p.181)に示されているように、等分散の仮定が成り立たないときは、感覚以外の要因、すなわち判断基準を含むので感覚の差の指標として用いることができない。等分散の仮定が成り立たない場合については、例えば、Swets, Signal Detection Theory and ROC Analysis in Psychology and Diagnostics: Collected papers (1996, pp. 52-54)で論じられている。T.M.Wickens (2002) Elementary Signal Detection Theory”では、式(1)が等分散ガウス型モデルにおけるものであることが明示されており、この訳本<岡本安晴(訳)「信号検出理論の基礎」>は、現在(2022.3)、電子書籍版がYahooショッピングのebookjapanなどから入手可能である。このWickens (2002;岡本(訳))では、Hautus, Macmillan, & Creelman (2022)Detection theory: A users guide”においてbasic bias measureとして提示されている次式

に、sensitivity が含まれていることが明示されている(岡本(訳)、p. 41、式(2.6))。Wickensでは、上式ので表されているが、次式で与えられている。

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 user’s 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を含まないのは、上で説明したようにシグナル刺激の分散がノイズ刺激の分散に等しいときである。

 

等分散を仮定するモデルでは、ノイズに対する感覚量の平均値を0、分散を1として、判断の基準値(上記、Hautusらのではないことに注意。ノイズ刺激の平均値を0とする座標における値である)に対して、次式で確率が与えられる。

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

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

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

 

表1

 

 

 

No

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

Yes

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

ノイズのみ

N(CR)

N(FA)

シグナル+ノイズ

N(Miss)

N(Hit)

 

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

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

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

 

シグナルではない

シグナルです

ノイズ刺激

69

31

シグナル刺激

31

69

 

上の圧縮ファイルを展開して作成されるフォルダSDT_EqSgmにおいてPythonスクリプトmain.pyを次のコマンドで実行する。

 

(py39) PS D:\SDT\SDT_EqSgm> python main.py

 

まず、ノイズ刺激に対する各回答の頻度の入力が求められる。

 

(py39) PS D:\SDT\SDT_EqSgm> python main.py

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

Data for noise =

 

上の表の数値を任意個の空白(スペース)で区切って入力する。

 

Data for noise = 69 31

 

ノイズ刺激に対する回答頻度を入力してEnterキーを押すと、次にシグナル刺激に対する回答数の入力が求められる。

 

Data for noise = 69 31

Data for signal =

 

シグナル刺激に対する回答数を入力する。

 

Data for noise = 69 31

Data for signal = 31 69

 

シグナル刺激に対する回答数を入力して、Enterキーを押すと、計算が始まる。

計算が進むと、まずシグナル刺激とノイズ刺激の感覚差d’の事後分布のグラフが表示される(図A.1)。

A.1

 

A.1のフォームを右上角のX印をクリックして閉じると、次に判断基準Cの事後分布のグラフが表示される(図A.2)。

A.2

 

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

実行終了後、フォルダの中を調べると以下のようになっている。

 

(py39) PS D:\SDT\SDT_EqSgm> ls

 

 

    ディレクトリ: D:\SDT\SDT_EqSgm

 

 

Mode                 LastWriteTime         Length Name

----                 -------------         ------ ----

-a----        2021/08/21     20:48          23251 FigureC.png

-a----        2021/08/21     20:45          22338 FigureMu.png

-a----        2021/08/21     15:27           2469 main.py

-a----        2021/08/21     20:50            233 Output.txt

-a----        2021/08/21     11:06            578 YN_EqSgm.stan

 

 

(py39) PS D:\SDT\SDT_EqSgm>

 

A.1のグラフの保存ファイルFigureMu.pngと図A.2のグラフの保存ファイルFigureC.pngが作成保存されおり、出力ファイルOutput.txtも作成されている。

出力ファイルOutput.txtの内容は以下のようである。パラメータの事後分布の推定値や95CIが出力されている。

 

(py39) PS D:\SDT_EqSgm> cat Output.txt

Data for Noise:  [69, 31]

Data for Signal: [31, 69]

 

 

mu:

MAP est. = 1.005

median   = 0.998

mean     = 0.996

95% CI   = [0.630, 1.359]

 

 

C:

MAP est. = 0.512

median   = 0.500

mean     = 0.497

95% CI   = [0.236, 0.749]

(py39) PS D:\SDT_EqSgm>

 

 

リストA.1 等分散SDTモデルの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 - 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);

}

 

 

リストA.2 リストA.1Stanスクリプトを利用するPythonスクリプト(ファイル名main.py

 

import pystan

import matplotlib.pyplot as plt

import numpy as np

import scipy.stats as ss

import seaborn as sb

 

 

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

    """

        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]

 

 

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 = open('Output.txt', 'w')

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}

 

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

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

print(fit)

 

mu = fit['mu']

c = fit['c']

 

mu_mean = mu.mean()

mu_MAP = CalcMAPEst(mu)[0]

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

 

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

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

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("$\mu$(d') = {0:<.3f}(med),  {1:<.3f}(mode)".format(mu_med, mu_MAP))

plt.yticks([])

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

plt.tight_layout()

plt.savefig('FigureMu.png')

plt.show()

 

c_mean = c.mean()

c_MAP = CalcMAPEst(c)[0]

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

 

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

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

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.yticks([])

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

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

plt.savefig('FigureC.png')

plt.show()

 

fout.close()

 

 

 

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

 

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

 

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

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

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

である(Swets, Signal detection theory and ROC analysis in psychology and diagnostics: Collected papers 1996, p. 71; 岡本「心理学データ分析と測定」、2014p.180)。

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

 

表3

 

 

 

 

カテゴリ1

・・・

カテゴリK

ノイズのみ

N(n1)

・・・

N(nK)

シグナルを含む

N(s1)

・・・

N(sK)

 

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

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

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

ダウンロード展開したフォルダ内にあるPythonスクリプトファイルmain.pyを次のコマンドで実行する。

 

(py39) PS D:\SDT_NESgm> python main.py

 

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

 

(py39) PS D:\SDT_NESgm> python main.py

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

K =

 

データとして、次の表のものを入力してみる。

 

 

ノイズ刺激

シグナル刺激

間違いなくノイズである

6

1

ノイズである

18

3

多分、ノイズである

34

12

わからない

31

28

多分、シグナルである

8

38

シグナルである

3

12

間違いなくシグナルである

0

6

 

上の表の場合、カテゴリ数が7であるので、7を入力してEnterキーを押す。

ノイズ刺激に対する回答カテゴリの頻度の入力が求められる。

 

(py39) PS D:\SDT_NESgm> python main.py

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

K = 7

Data for Noise =

 

上の表のデータを、空白(スペース)で区切って並べて入力する。区切りの空白文字の個数は任意である。

ノイズ刺激に対する各カテゴリの回答頻度を入力してEnterキーを押すと、次にシグナル刺激に対する回答頻度の入力が求められる。

 

Data for Noise = 6 18 34 31 8 3 0

Data for Signal =

 

上の表のデータを以下のように入力する。

 

Data for Noise = 6 18 34 31 8 3 0

Data for Signal = 1 3 12 28 38 12 6

 

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

Stanによる計算が終了すると、まずすなわちの事後分布のグラフが表示される(図B.1)。

B.1

 

B.1のフォームの右上角のX印をクリックして閉じると、次にの事後分布のグラフが表示される(図B.2)。

B.2

 

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

B.3

 

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

B.3のフォームを閉じると、カテゴリ境界の事後分布のグラフが表示される(図B.4)。

B.4

 

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

プログラムの実行終了後、フォルダを開くと以下のようになっている。

 

(py39) PS D:\SDT_NESgm> ls

 

 

    ディレクトリ: D:\SDT_NESgm

 

 

Mode                 LastWriteTime         Length Name

----                 -------------         ------ ----

-a----        2021/08/22     11:27          51765 Figure_Cs.png

-a----        2021/08/22     11:21          23660 Figure_mu.png

-a----        2021/08/22     11:24          24082 Figure_ROC.png

-a----        2021/08/22     11:22          22583 Figure_sgm.png

-a----        2021/08/22     11:16           5026 main.py

-a----        2021/08/22     11:28            946 Output.txt

-a----        2021/08/21     14:52            871 YN_rating.stan

 

 

(py39) PS D:\SDT_NESgm>

 

B.1から図B.4までのグラフがFigure_mu.pngからFigure_Cs.pngまでのファイルに保存されていることが分かる。

計算結果の出力はOutput.txtに保存されている。開くと以下のようになっている。

 

(py39) PS D:\SDT_NESgm> cat Output.txt

 

Data for Noise:

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

 

Data for Signal:

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

mu:

    MAP est. = 1.305

    med.     = 1.319

    mean     = 1.324

    95% CI   = [0.951, 1.731]

 

sigma:

    MAP est. = 1.155

    med.     = 1.174

    mean     = 1.185

    95% CI   = [0.916, 1.514]

 

 

C1:

    MAP est. = -1.563

    median   = -1.566

    mean     = -1.571

    95% CI   = [-1.968, -1.210]

 

C2:

    MAP est. = -0.704

    median   = -0.720

    mean     = -0.722

    95% CI   = [-0.989, -0.460]

 

C3:

    MAP est. = 0.199

    median   = 0.184

    mean     = 0.184

    95% CI   = [-0.052, 0.416]

 

C4:

    MAP est. = 1.195

    median   = 1.189

    mean     = 1.191

    95% CI   = [0.900, 1.493]

 

C5:

    MAP est. = 2.283

    median   = 2.294

    mean     = 2.302

    95% CI   = [1.844, 2.816]

 

C6:

    MAP est. = 3.095

    median   = 3.148

    mean     = 3.174

    95% CI   = [2.490, 3.997]

(py39) PS D:\SDT_NESgm>

 

パラメータの事後分布のMAP推定値(モード)、中央値、平均値、95CIが出力されている。

 

 

リストB.1 SDT不等分散モデルの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] = normal_cdf(c[1], 0.0, 1.0);  // Phi(c[1]);

    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 {

    Freq_Noise ~ multinomial(theta_n);

    Freq_Signal ~ multinomial(theta_s);

}

 

 

リストB.2 リストB.1Stanスクリプトを用いるPythonスクリプト(ファイル名main.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

 

 

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

    """

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

 

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}

 

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

 

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

 

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

 

 

Up