1AFC rating task for signal detection theory (SDT)
The scripts have been revised for CmdStanPy (2026.01). A simple installation of CmdStanPy is shown at this website.
In One-alternative forced choice (1AFC)
task, one of the two stimulus, e.g., noise only
stimulus and noise + signal stimulus, is presented, and the observer is asked
to judge which one was presented. The sensation of the noise only stimulus is
represented by 
, and that of the signal stimulus by 
. They are assumed to have the following distribution
Origin and unit of the sensation are
assumed so that 
. Hence we have
The observer is required to respond with rating about his/her judgment. For example, “Certainly noise only stimulus”, “Probably signal stimulus”, and so on.
In K category rating task, let’s the judgment be
represented by integer from 1 to K. “1” represents the most certain judgment of the noise only, and K the
most certain judgment of the stimulus. Corresponding to the category judgment,
the sensation dimension is partitioned by 
category boundaries.
When the sensation is between 
and 
, judgment “k” is
done.
So, we have the following equations:
where 
represents the
cumulative distribution function of the standard normal distribution.
The difference between the means of
sensations and is denoted by 
, that is,
Data of 1AFC rating task can be summarized like Table 1.
|
Tabnle 1 Data of 1AFC rating task |
|||
|
|
Category-1 |
… |
Category-K |
|
Noise only |
N(n1) |
… |
N(nK) |
|
Noise + Signal |
N(s1) |
… |
N(sK) |
N(nk) and N(sk) denote the number of judgment category-k for the noise stimulus and the noise+signal stimulus, respectively.
The likelihood function for the data of Table 1 is given by
(1)
Stan script for eq,(1) is shown in Listing 1.
Python script using the Stan script in Listing 1 is shown in Listing 2.
Files of these scripts are archived in the file sdt_1afc_rating.zip, which can be freely downloaded and used.
An example of data is shown in Figure 1. The data is stored in an Excel file data.xlsx (Figure 1).
Figure 1
Frequences of judgment of categories are set from the second row. In the second row, frequences of most sure judgment of noise stimulus are set. In the last row, frequences of most sure judgment of signal are set. In the second column, judgments for noise stimulus are set, and in the third column, judgment for signal stimulus are set.
Put the script files of Listings 1 and 2, and the input data file in the same folder, to which the current directory is changed.
Run the script main.py in an environment with cmdstanpy installed.
(stan) ****/sdt_1afc_rating$ ls
YN_rating.stan data.xlsx main.py
(stan) ****/sdt_1afc_rating$ python
main.py
Input data file (*.xlsx) = data.xlsx
Input data file name is asked. In the above example, the data file of Figure 1 is set.
When the input data file name is set, building (compilation) of the Stan script starts. After building, MCMC sampling starts. Trace plots are shown like Figure 2
Figure 2
Close the window of Figure 2, graph of the
posterior distribution of 
is displayed (Figure 3).
Figure 3
Close the window of the posterior
distribution of , graph of the posterior distribution of 
is displayed (Figure 4).
Figure 4
After closing the window of the posterior
distribution of , it takes some time for calculation.
After the calculation, graph of ROC curve is displayed (Figure 5).
Figure 5
Close the window of ROC curve, graph of posterior distributions of criterion boundaries is displayed (Figure 6).
Figure 6
Close the window, then the script ends.
Execution of the script produces a text file Output .txt, in which some information of the analysis is stored as follows:
Data for Noise:
[6 18 34 31 8 3 0]
Data for Signal:
[1 3 12 28 38 12 6]
mu:
med. = 1.318
mean = 1.323
95% CI = [0.955, 1.718]
sigma:
med. = 1.176
mean = 1.185
95% CI = [0.915, 1.524]
C1:
median = -1.564
mean = -1.569
95% CI = [-1.960, -1.205]
C2:
median = -0.725
mean = -0.723
95% CI = [-0.981, -0.470]
C3:
median = 0.182
mean = 0.183
95% CI = [-0.052, 0.416]
C4:
median = 1.187
mean = 1.188
95% CI = [0.902, 1.479]
C5:
median = 2.291
mean = 2.300
95% CI = [1.850, 2.811]
C6:
median = 3.157
mean = 3.178
95% CI = [2.481, 4.005]
Listing
1. Stan script for 1AFC rating task data (file name: 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);
}
Listing
2. Python script using YN_rating.stan
in Listing 1 (file name: main.py)
import pandas as pd
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot as plt
import numpy
as np
import seaborn as sb
import scipy.stats as ss
import arviz
as az
fout = open('Output.txt',
'w')
inflnm = input('Input
data file (*.xlsx) =
')
data_pf = pd.read_excel(inflnm)
print('data_pf
=\n', data_pf)
data = data_pf.values
print('data =\n', data)
Data_Noise = data[:,1]
Data_Signal = data[:,2]
K = len(Data_Noise)
print('Data_Noise
=\n', Data_Noise)
print('Data_Signal
=\n', Data_Signal)
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.diagnose())
print(fit.summary())
infdata = az.from_cmdstanpy(fit) # Transforming to Arviz InfereceData
az.plot_trace(infdata, var_names=['mu', 'sigma'])
plt.tight_layout()
plt.savefig('Fig_trace.png')
plt.show()
fit = fit.draws_pd()
#
Transforming to Pandas DataFrame
print(fit.keys())
mu = fit['mu']
sigma = fit['sigma']
mu_mean = mu.mean()
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_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(' 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.savefig('Figure_mu.png')
plt.show()
sb.kdeplot(sigma)
plt.title(r"$\sigma$ = {0:<.3f}(med) ".format(sgm_med), fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize
= 14)
plt.xlabel(r'$\sigma$', fontsize
= 18)
plt.tight_layout()
plt.savefig('Figure_sgm.png')
plt.show()
c_smpls = []
for k in range(K-1):
c_smpls.append(fit[f'c[{k+1}]'])
c_smpls = np.array(c_smpls).T
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 =
(5,5))
plt.plot(x_values, y_values, 'b-')
plt.title(fr"ROC-curve($\Delta$m={mu_med:.3f}, $\sigma_s$={sgm_med:.3f})",
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_q025s = np.zeros(K-1)
c_q975s = np.zeros(K-1)
for k in range(K-1):
c_means[k]
= c_smpls.T[k].mean()
c_q025s[k], c_meds[k], c_q975s[k] = \
np.percentile(c_smpls.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.legend(fontsize =
12)
plt.tight_layout()
plt.savefig('Figure_ROC.png')
plt.show()
for k in range(K-1):
sb.kdeplot(c_smpls.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.savefig('Figure_Cs.png')
plt.show()
fout.close()
print('Output.txt was saved.')