信号検出理論（Signal Detection Theory: SDT）
信号検出理論に基づくデータ分析プログラムをStan(PyStan)を使ってPythonで作成してみた。信号検出理論については＜岡本安晴「心理学データ分析と測定」勁草書房、第１２章＞で説明している。等分散を仮定するモデルと仮定しないモデルのそれぞれについて作成した。Stan(PyStan)については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
モデルの比較をWAICで行う例も用意した。
等分散モデルでは、シグナル刺激とノイズ刺激の感覚の分散は等しく１であると仮定される。この仮定の下で、ノイズとシグナルの感覚の平均値の差が次式
で算出される（Kingdom & Prins, “Psychophysics: A practical introduction, 2nd ed.” 2016, p. 176; Macmillan & Creelman, “Detection Thoery: A user’s guide.” 2005, p. 8）。
しかし、式（１）の右辺は、岡本「心理学データ分析と測定」（2014、p.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”では、式（１）が等分散ガウス型モデルにおけるものであることが明示されており、この訳本＜岡本安晴（訳）「信号検出理論の基礎」＞は、現在（2022.3）、電子書籍版がYahooショッピングのebookjapanなどから入手可能である。このWickens (2002；岡本（訳）)では、Hautus, Macmillan, & Creelman (2022)”Detection theory: A user’s guide”においてbasic bias measureとして提示されている次式
に、sensitivity 
が含まれていることが明示されている（岡本（訳）、p. 41、式（2.6））。Wickensでは、上式の
は
で表されているが、次式で与えられている。
Hautusらのp. 31に挙げられている例では、Trainee 1は訓練後も同じcriterionを維持している（
）と説明されているが、これは上のバイアス基準を用いた場合である。sensitivity の値も変化しているので、刺激
（いわゆるノイズのみの刺激）を基準にした座標（ノイズ刺激の平均値＝０）では、ほぼ
から
に変化している。訓練前はノイズ刺激の左側に設定されていた基準値が訓練後はノイズ刺激の右側に移動している。
上の式（Hautus at al. (2022), eq.(2.1)）の右辺を不等分散モデルの場合で表してみる。ここで、判断の基準値をHautusらに従って
で表し、シグナル刺激の平均値と標準偏差を
と
で表すと、以下のように式が展開される。
すなわち、不等分散モデルにおいては、式（Hautus at al. (2022),
eq.(2.1)）の右辺には、判断の基準値の他に、シグナル刺激の標準偏差と平均値が含まれている。ノイズ刺激の平均値と標準偏差を基準としたとき、判断の基準値の他にシグナル刺激の平均値と標準偏差の要因も含まれた値である。ここで注意しなければならないことは、biasとcriterionの区別であろう。ノイズ刺激とシグナル刺激の中点
に置かれたcriterionがbiasの無いcriterionであると考えるならば、
との差がbiasとなる。しかし、ノイズ刺激に対する正答率とシグナル刺激に対する正答率が等しくなる点
をbiasのない点と考えると、不等分散モデルの下では以下のようになる。
これより、
となる。Biasが０の点は、シグナルの標準偏差が大きくなるとノイズ刺激に近付き、が小さくなるとシグナル刺激に近付く。
なお、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を含まないのは、上で説明したようにシグナル刺激の分散がノイズ刺激の分散に等しいときである。
なお、P(F)とP(H)の組に対して、signal刺激の感覚の分布がの値に応じて変わることを示す図を、https://osf.io/5pzby/のファイルMuSigma.pdfに示した。
等分散を仮定するモデルでは、ノイズに対する感覚量の平均値を０、分散を１として、判断の基準値
（上記、Hautusらのではないことに注意。ノイズ刺激の平均値を０とする座標における値である）に対して、次式で確率が与えられる。
ここで、信号検出理論では、「シグナル＋ノイズ」刺激に対する感覚量の平均値
を
で表すことが多い。「ノイズのみ」刺激に対する感覚量の平均が０であるので、は、「ノイズのみ」刺激と「シグナルを含む」刺激の感覚量の平均の差を表す。
いま、「No|Noise」（Correct Rejection）の反応数をN(CR)で、というように表したとき、データが表１のように得られたとする。
|
表１ |
|
|
|
|
No （「シグナルなし」判断） |
Yes （「シグナルあり」判断） |
|
ノイズのみ |
N(CR) |
N(FA) |
|
シグナル＋ノイズ |
N(Miss) |
N(Hit) |
このとき、尤度関数は次式で与えることができる。
上式に基づいて、２項分布を用いたStanのコードをリストA.1のように用意し、リストA.1のStanスクリプトを利用するPythonスクリプトをリストA.2のように用意した。リストA.1とリストA.2のファイルは圧縮ファイルSDT_EqSgm.zipにまとめた。ダウンロード展開して自由に利用されたい。
いま、次のような表のデータが与えられたとする。
|
|
シグナルではない |
シグナルです |
|
ノイズ刺激 |
69 |
31 |
|
シグナル刺激 |
31 |
69 |
上の圧縮ファイルを展開して作成されるフォルダSDT_EqSgmにおいてPythonスクリプトmain.pyを次のコマンドで実行する。
実行は、Anaconda/Windows 10のPython 3.7仮想環境においてPyStan 2.17.1.0を利用して行った。
(py37) PS
D:\SDT\SDT_EqSgm> python main.py
まず、ノイズ刺激に対する各回答の頻度の入力が求められる。
(py37) 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のフォームを閉じると、プログラムの実行終了である。
実行終了後、フォルダの中を調べると以下のようになっている。
(py37) 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
(py37) PS D:\SDT\SDT_EqSgm>
図A.1のグラフの保存ファイルFigureMu.pngと図A.2のグラフの保存ファイルFigureC.pngが作成保存されおり、出力ファイルOutput.txtも作成されている。
出力ファイルOutput.txtの内容は以下のようである。パラメータの事後分布の推定値や95％CIが出力されている。
(py37) 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.1のStanスクリプトを利用する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直線上でデータ点が散らばるようにするが、これは判断カテゴリの境界値の散らばりが大きくなるようにすることである。判断カテゴリ数を多く取ることが１つの方法である。３カテゴリの場合であっても、中央の「わからない」カテゴリの頻度が少なくならないように、「判りにくいときは『わかりません』と答えて下さい」などと教示しておくとよいと思われる。
ノイズのみの刺激に対する感覚量の分散と「シグナル＋ノイズ」の刺激に対する感覚量の分散が等しいことを仮定しない場合は、ノイズのみの刺激に対する感覚量の平均を０、分散を１とおき、「シグナル＋ノイズ」刺激に対する感覚量の平均を、分散を
とおいて、モデルが設定される。このとき、ＲＯＣ曲線上のデータ点として２点以上必要であるが、評定法を用いることによって１回の実験から複数のデータ点を得ることができる。
いま、評定カテゴリ数がＫ個であるとする。「ノイズのみである」の確信が最も強いカテゴリをカテゴリ１とし、「シグナルを含む」の確信が最も強いカテゴリをカテゴリＫとして、カテゴリ１からカテゴリＫまで順番に「シグナルを含む」とする確信の強い方向に評定カテゴリが用意がされているとする。カテゴリkとカテゴリk+1の感覚量の境界値を判断基準値
で表す。このとき、以下のモデルを設定する。
「ノイズのみ」刺激に対する感覚量の平均０と「シグナルを含む」刺激に対する感覚量の平均との差は
で表される。すなわち、
である（Swets, “Signal detection theory and ROC analysis in psychology and diagnostics: Collected papers” 1996, p. 71; 岡本「心理学データ分析と測定」、2014、p.180）。
いま、「ノイズのみ」刺激に対するカテゴリｋの評定数をN(nk)で、「シグナルを含む」刺激に対するカテゴリｋの評定数をN(sk)で表すと、表３のようにまとめることができる。
|
表３ |
|
|
|
|
|
カテゴリ１ |
・・・ |
カテゴリＫ |
|
ノイズのみ |
N(n1) |
・・・ |
N(nK) |
|
シグナルを含む |
N(s1) |
・・・ |
N(sK) |
表３のデータに対して尤度関数を次のように設定する。
上式に基づいて多項分布を用いたStanのスクリプトをリストB.1のように用意した。リストB.1のStanスクリプトを用いるPythonスクリプトをリストB.2に用意した。
リストB.1とリストB.2のファイルは圧縮ファイルSDT_NESgm.zipにまとめた。圧縮ファイルをダウンロード展開して自由に利用されたい。
ダウンロード展開したフォルダ内にあるPythonスクリプトファイルmain.pyを次のコマンドで実行する。
実行は、Anaconda/Windows 10におけるPython 3.7仮想環境においてPyStan 2.17.1.0を利用して行った。
(py37) PS
D:\SDT_NESgm> python main.py
実行すると、まず評定カテゴリ数Kの入力が次のように求められる。
(py37) 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 |
上の表の場合、カテゴリ数が７であるので、７を入力してEnterキーを押す。
ノイズ刺激に対する回答カテゴリの頻度の入力が求められる。
(py37) 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のフォームを閉じると、プログラムの実行終了である。
プログラムの実行終了後、フォルダを開くと以下のようになっている。
(py37) 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
(py37) PS D:\SDT_NESgm>
図B.1から図B.4までのグラフがFigure_mu.pngからFigure_Cs.pngまでのファイルに保存されていることが分かる。
計算結果の出力はOutput.txtに保存されている。開くと以下のようになっている。
(py37) 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.323
med. = 1.318
mean = 1.324
95% CI = [0.955, 1.727]
sigma:
MAP est. = 1.156
med. = 1.177
mean = 1.185
95% CI = [0.914, 1.511]
C1:
MAP est. = -1.548
median = -1.565
mean = -1.571
95% CI = [-1.966, -1.210]
C2:
MAP est. = -0.711
median = -0.723
mean = -0.724
95% CI = [-0.988, -0.464]
C3:
MAP est. = 0.177
median = 0.183
mean = 0.184
95% CI = [-0.047, 0.419]
C4:
MAP est. = 1.181
median = 1.187
mean = 1.190
95% CI = [0.900, 1.489]
C5:
MAP est. = 2.274
median = 2.295
mean = 2.303
95% CI = [1.842, 2.814]
C6:
MAP est. = 3.160
median = 3.158
mean = 3.177
95% CI = [2.487, 3.988]
パラメータの事後分布のMAP推定値（モード）、中央値、平均値、95％CIが出力されている。
リスト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, 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);
}
リストB.2 リストB.1のStanスクリプトを用いる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
= 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')
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()
不等分散モデルと等分散モデルのデータに対する適合度の比較をWAICによって行ってみた。
不等分散モデルと等分散モデルによるベイズ分析のためのStanスクリプトをリストW1およびリストW2に示すように作成した。PyStan 3対応である。これらのStanスクリプトを利用してWAICによる比較を行うPythonスクリプトをリストW3のように用意した。これらのスクリプトファイルはまとめて圧縮ファイルCompModels.zipとした。ダウンロード・展開して自由に利用されたい。スクリプトファイルは、WindowsにインストールされたUbuntu上で実行した。実行環境は、Ubuntu/WSL/WindowsにインストールされたAnacondaのPython 3.10仮想環境である。WSLなどについては、このウェブサイトを参照されたい。
モデルの比較を行うデータを次のように用意した。
K = 7
Data_Noise = np.array([6, 18, 34, 31,
8, 3, 0])
Data_Signal = np.array([1, 3, 12, 28,
38, 12, 6])
リストW3のスクリプトを実行すると、まず、不等分散モデルによるMCMCサンプリングが実行され、サンプリング後、
の事後分布が表示される（図W1）。
図W1
「
」の等分散モデルの可能性が高いように思われる分布である。
図W1のフォームを閉じると、図W2のグラフが表示される。
図W2
不等分散モデル（UnEq.と表記）および等分散モデル（Eq.と表記）のWAICとそのブートストラップ分布である。
WAICは次のコードで算出している。
def calcWAIC(log_ps):
p_WAIC = np.var(log_ps, ddof=1)
# Gelman
et al. (2014), eq.(7.12)
ps =
np.exp(log_ps)
lppd
= np.log(np.mean(ps)) # Gelman et al. (2014), eq.(7.5)
WAIC
= -2 * (lppd - p_WAIC) # Gelman et al. (2014), p.174
return WAIC
ブートストラップ分布は以下のようなコードで求めている。ブートストラップについては、岡本（2022）、pp. 167-170で説明している。
log_ps = d_frm['LL']
WAIC = calcWAIC(log_ps)
WAICs = [WAIC]
for i in range(999):
t_log_ps = np.random.choice(log_ps, size = len(log_ps), replace = True)
WAICs.append(calcWAIC(t_log_ps))
不等分散モデルのWAIC（＝58.74）より等分散モデルのWAIC（＝56.54）の方が小さく、またブートストラップ分布も明瞭に分離されている。等分散モデルを選ぶことが示唆されているようである。
次に、比率はそのままにして、総数を10倍にしてみる。すなわち、
Data_Noise *= 10
Data_Signal *= 10
である。
図W2のフォームを閉じると、上の10倍にした条件で計算が始まる。
不等分散モデルによるの事後分布は図W3のようになる。
図W3
「」の位置からかなり右に移動していることが認められる。
図W3のフォームを閉じると、図W4のグラフが表示される。
図W4
不等分散モデル（UnEq.で表記）のWAIC（＝120.16）の方が等分散モデル（Eq.で表記）のWAIC（＝122.75）より小さく、ブートストラップの分布も明瞭に不等分散モデルの方が小さい方に（左側に）に分布している。等分散モデルを採用することは難しいと考えられる。
データ数が十分に多くないときは、データとモデルとの関係は弱い、すなわちモデルの情報が十分にデータに含まれていないので、データとモデルとの関係は大まかなものとなり、小さい相違は反映されないと考えられる。データが十分に多くなるとデータとモデルの関係は強くなり、データにモデルの情報が十分に含まれ、データとモデルとの小さい違いも反映されるようになると考えられる。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian Data Analysis, third edition. CRC Press.
岡本安晴（2022） データを使いこなすための統計入門．22世紀アート．
リストW1 不等分散モデルのPyStan3対応Stanスクリプト（YN_rating.stan）
//
// The unequal-variance Gaussian
model
//
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 mu;
real<lower = 0.0>
p_sigma;
simplex[K] pre_c;
}
transformed parameters {
real sigma;
vector[K-1] c;
vector[K] theta_n;
vector[K] theta_s;
real v;
sigma = p_sigma + 0.001;
v = pre_c[1];
for (k in 1:K-1) {
c[k]
= logit(v);
v +=
pre_c[k+1];
}
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 {
pre_c ~ dirichlet(a);
mu ~ normal(0.0, 10.0);
p_sigma ~
exponential(0.1);
Freq_Noise ~
multinomial(theta_n);
Freq_Signal ~
multinomial(theta_s);
}
generated quantities {
real LL;
LL =
multinomial_lpmf(Freq_Noise | theta_n) +
multinomial_lpmf(Freq_Signal | theta_s);
}
リストW2 等分散モデルのPyStan3対応Stanスクリプト（YN_rating_eq.stan）
//
// The
equal-variance Gaussian model
//
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 mu;
simplex[K] pre_c;
}
transformed parameters {
real sigma;
vector[K-1] c;
vector[K] theta_n;
vector[K] theta_s;
real v;
sigma = 1.0;
v = pre_c[1];
for (k in 1:K-1) {
c[k]
= logit(v);
v +=
pre_c[k+1];
}
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 {
pre_c ~ dirichlet(a);
mu ~ normal(0.0, 10.0);
Freq_Noise ~
multinomial(theta_n);
Freq_Signal ~ multinomial(theta_s);
}
generated quantities {
real LL;
LL =
multinomial_lpmf(Freq_Noise | theta_n) +
multinomial_lpmf(Freq_Signal | theta_s);
}
リストW3 モデル比較をWAICで行うPythonスクリプト（CompModels.py）
#
# Yasuharu Okamoto,
2023.07
#
import matplotlib.pyplot as plt
import numpy as np
import stan
import seaborn as sb
import scipy.stats as ss
import arviz as az
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]
K = 7
Data_Noise = np.array([6, 18, 34, 31,
8, 3, 0])
Data_Signal = np.array([1, 3, 12, 28,
38, 12, 6])
for t in range(2):
if t == 1:
Data_Noise *= 10
Data_Signal *= 10
Data = {'K': K,
'Freq_Noise': Data_Noise, 'Freq_Signal': Data_Signal}
def f_init():
return {'mu': 0.0, 'p_sigma': 1.0, 'pre_c':[1.0/K]*K}
with open('YN_rating.stan',
'r') as fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_chains =
4, init = [f_init()]*4)
print(fit)
i_data =
az.from_pystan(posterior = fit, log_likelihood = {'LogL':'LL'},
posterior_model
= sm)
smmry = az.summary(i_data,
var_names = ['mu', 'sigma'])
print(smmry)
d_frm = fit.to_frame()
sb.kdeplot(d_frm['sigma'])
plt.xlabel('$\sigma$',
fontsize = 14)
plt.title(f'N = {np.sum(Data_Noise)
+ np.sum(Data_Signal)}')
plt.show()
def calcWAIC(log_ps):
p_WAIC = np.var(log_ps, ddof=1)
# Gelman et al. (2014),
eq.(7.12)
ps =
np.exp(log_ps)
lppd
= np.log(np.mean(ps)) # Gelman et al. (2014), eq.(7.5)
WAIC
= -2 * (lppd - p_WAIC) # Gelman et al. (2014), p.174
return WAIC
log_ps = d_frm['LL']
WAIC = calcWAIC(log_ps)
WAICs = [WAIC]
for i in range(999):
t_log_ps = np.random.choice(log_ps, size = len(log_ps), replace = True)
WAICs.append(calcWAIC(t_log_ps))
def f_init_eq():
p_c
= [1.0/K]*K
return {'mu': 0.0, 'pre_c': p_c}
with
open('YN_rating_eq.stan', 'r') as fstan:
sm_eq = stan.build(fstan.read(), data = Data)
fit_eq =
sm_eq.sample(num_chains = 4, init = [f_init_eq()]*4)
print(fit_eq)
i_data_eq =
az.from_pystan(posterior = fit_eq, log_likelihood = {'LogL':'LL'},
posterior_model = sm)
smmry_eq = az.summary(i_data_eq,
var_names = ['mu'])
print(smmry_eq)
d_frm_eq = fit_eq.to_frame()
log_ps_eq = d_frm_eq['LL']
WAIC_eq =
calcWAIC(log_ps_eq)
WAICs_eq = [WAIC_eq]
for i in range(999):
t_log_ps_eq = np.random.choice(log_ps_eq,
size = len(log_ps_eq), replace = True)
WAICs_eq.append(calcWAIC(t_log_ps_eq))
plt.title('Bootstrap
distributions of WAIC\nN = {0}'.
format(np.sum(Data_Noise)+np.sum(Data_Signal)))
sb.kdeplot(WAICs, label = f'UnEq.\nWAIC =
{WAIC:.2f}')
sb.kdeplot(WAICs_eq, label =
f'Eq.\nWAIC = {WAIC_eq:.2f}')
plt.legend()
plt.show()