Yasuharu Okamoto, 2017/07
Stan/Python Scripts for Analysis in Signal Detection
Theory (SDT)
Two programs were developed, one is for a Yes/No type with equal variance assumption and the other is for a Yes/No method using a rating method without equal variance assumption.
Equal Variance Model
Under the condition of equal variances for
noise and signal stimuli, sensation by the noise is assumed to have mean 0 and
variance 1, and sensation by the signal is assumed to have mean 
and variance 1. With a
criterion 
of judgment, we have the
following equations.
The difference of means for signal and
noise is designated by 
, so we have
Suppose that we have data shown in Table 1, where N(CR) denotes the number of responses “No” to “Noise” stimulus (Correct Rejection), and so on.
|
Table 1 |
|
|
|
|
No (No “Signal”) |
Yes (Yes, “Signal” displayed) |
|
Noise |
N(CR) |
N(FA) |
|
Signal |
N(Miss) |
N(Hit) |
We have the following likelihood function.
Stan script for the above model is given as follows
================== 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);
}
====================================================
An example data of Table 1 is shown in Table 2.
|
Table 2 |
|
|
|
|
No |
Yes |
|
Noise |
６９ |
３１ |
|
Signal |
３１ |
６９ |
For the data of Table 2, the following Python script which uses the above Stan script was written. Numbers of responses to Noise are set in Data_Noise list, and those to Signal in Data_Signal list.
======================================================
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()
=======================================================
Run this script, a histogram of samples
from the posterior distribution of (Figure 1) is shown
after MCMC sampling by Stan.
Figure 1
Above the histogram, the median is shown.
When the form is closed by clicking on the icon X at the upper right corner, the next histogram is shown (Figure 2).
Figure 2
On the form of Figure 2, a histogram of
samples from the posterior distribution of criterion and the median are shown.
When the form of Figure 2 is closed, the program ends.
The program files are archived in SDT_EqSgm.zip . Click on the file name SDT_EqSgm.zip to down load.
Model, which allows unequal variances
To set an origin and a unit, mean and
variance for noise are assumed to be 0 and 1, respectively. Mean and variance
for signal are denoted by and 
, respectively. This type of model with variances, which may be unequal
each other, requires data which gives more than one points with respect to ROC
curve. We can get more than one points by employing a rating method.
Suppose that an observer responds to the
stimulus choosing one from K rating categories, Category-1 (the strongest
confidence in judgment of “no signal”) to Category-K (the strongest confidence in judgment of “signal”). Denote the criterion between
Category-k and Category-(k+1) by 
. Then, we have the following model.
The difference of means for noise and
signal is denoted by 
, i.e.,
When we denote the number of responses of Category-k to noise by N(nk) and to signal by N(sk), we have Table 3.
|
Table 3 |
|
|
|
|
|
Category-1 |
・・・ |
Category-K |
|
Noise |
N(n1) |
・・・ |
N(nK) |
|
Signal |
N(s1) |
・・・ |
N(sK) |
A likelihood function for Table 3 can be set as follows.
A Stan script for the above model is given as follows.
================== 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);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
An example data in format of Table 3 is shown in Table 4, where “Category-1” to “Category-4” represent judgment of “Noise only”, “Probably Noise only”, “Probably Signal present” to “Signal present”, respectively.
|
Table 4 |
|
|
|
|
|
|
Category-1 |
Category-2 |
Category-3 |
Category-4 |
|
Noise |
10 |
9 |
7 |
1 |
|
Signal |
2 |
1 |
19 |
5 |
To analyze the data of Table 4, the following Python script using the above Stan script is given as follows. The frequencies of responses of each category to “Noise” and “Signal” are set in lists Data_Noise and Data_Signal, respectively.
==========================================
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()
===========================================================
Run this program, then a form of Figure 3 is shown after MCMC sampling by Stan.
Figure 3
In Figure 3, a histogram and the median of
the samples from the posterior distribution of are shown. “Delta-m” denotes .
Click the icon X at the upper right corner of the form, then the form is closed and the next form is presented (Figure 4).
Figure 4
Figure 4 shows a histogram and the median
of the samples from the posterior distribution of 
.
Close the form of Figure 4, then the next form is presented (Figure 5).
Figure 5
Figure 5 shows the ROC curve. Small red circles represent data points, and small green squares represent the corresponding points estimated by the model using the median values of the samples.
Close the form of Figure 5, then the next form is shown (Figure 6).
Figure 6
In Figure 6, histograms drawn by lines and
medians of samples from the posterior distributions of criterions s are shown.
When the form is closed, the program ends.
The program files are archived into file SDT_NESgm.zip ,which can be down loaded by clicking the file name SDT_NESgm.zip.