Extending 2AFC rating task
Estimation of both the mean and variance of a signal
stimulus
Yasuharu Okamoto, 2023.05
2AFC (two alternative forced choice) task is one of the popular methods of SDT (signal detection theory) (Wixted, 2019). As to measuring memory, Brady et al. (2023) recommend 2AFC, but they point out that, in 2AFC task, equal variance is assumed, and unequal variance models provide a good account of ROC curves. The unequal variance SDT model has been widely used (DeCarlo, 2010). Hence, we need to estimate mean and variance for a signal stimulus in 2AFC task. To estimate mean and variance, extension of 2AFC task was proposed by Okamato (2023; https://osf.io/ts2eq/)
In the proposed extension of 2AFC task, two pairs of stimuli, (noise,noise) and (signal,signal), were added to the pairs, (noise,signal) and (signal,noise) , of the standard 2AFC task, and rating of judgment was introduced. From rating data of the pair (noise,noise) , boundaries of categories of rating can be estimated. From rating data of the pair (signal,signal) and estimates of category boundaries, variance of sensation of the signal stimulus can be estimated. From these estimates and data of pairs of (noise,signal) and (signal,noise), mean of sensation of the signal stimulus can be estimated. These are intuitive explanation, which justifies the procedure. For a given data, statistical analysis is done by Bayesian method using a stochastic model as follows.
Stochastic Model
Let sensations of the noise stimulus 
and the signal one 
be represented by 
and 
, respectively. They are assumed to have normal distributions as
follows:
Set the origin and the unit so that
In two alternative forced choice (2AFC)
task, two stimuli are presented in one trial, at positon 1 and position 2.
These positons are temporal or spacial ones, e.g., the first one and the second
one, or the left side and the right side. Denote these positons using angle
brackets. For example, <S1,S2> means that S1 is presented at the position1
and S2 is at the positon 2. <
> means that sensation of the stimulus at positon 1 is 
and sensation of the
stimulus at positon 2 is 
.
When the bias of position 2 relative to
positon 1 is denoted by 
, we have:
If 
, then 
If 
, then 
, then 
If 
, then 
The bias parameter may include sensory ones, e.g., constant time or space error, and response inclination.
Now, set category boundaries of rating as follows:
where number of categories is 
.
When the following condition holds
rating 
for the pair <> is 
., i.e.
For example, ratings in the case of 
may be as follows:
denotes “It is sure that the signal stimulus was presented at positon 1”
denotes “I think that the signal stimulus was presented at positon 1”
denotes “Maybe, the signal stimulus was presented at positon 1”
denotes “I don’t know/Both stimuli seemed to be the
same”
denotes “Maybe, the signal stimulus was presented at positon 2”
denotes “I think that the signal stimulus was presented at positon 2”
denotes “It is sure that the signal stimulus was presented at positon 2”
Denote a cumulative normal distribution
with mean 
and variance 
by
As a convention, set
Then we have the following equations
When the stimuli 
are presented, we have
Hence,
When the stimuli 
are presented, we have
Hence,
When the stimuli 
are presented, we have
Hence,
When the stimuli 
are presented, we have
Hence,
Category boundaries are set to be symmetric about the origin 0.
In case of K even, we set as follows
In case of K odd, we set as follows
According to this difference, Stan script are prepared as in Listing 1 (K even) and in Listing 2 (K odd). The script files used in this website are archived into the file ex2afcrfiles.zip, which can be downloaded and used freely.
The script file sampling.py in Listing 3 uses the script files in Listngs 1 and 2. This file was run in Python 3.11 virtual environment of Anaconda on Ubuntu. So, when you use the script in Windows, you should use it on WSL (Windows Subsystem for Linux). In a Python 3.11 virtual environment, PyStan 3 (Python interface to Stan) can be installed by the command
pip install
pystan
Before installing PyStan 3 on Ubuntu, g++ need to be installed. In Ubuntu, run the code
sudo apt
install g++
Then, install PyStan 3 in the Python 3.11 virtual environment.
After MCMC sampling by the script sampling.py in Listing 3, analysis of the sample can be done by the script analysis.py in Listing 4. This script can be run in Python 3.11 virtual environment of Anaconda in Windows, not in /WSL/Windows, in which graphs would not be presented.
For WSL and PyStan, check websites
http://y-okamoto-psy1949.la.coocan.jp/Python/en2/wsl_ubuntu/.
and
http://y-okamoto-psy1949.la.coocan.jp/Python/en3/PyStan3WSL/
Example of Analysis of the extended 2AFC rating data
I show how to use the script files separately in cases where K is even and where K is odd.
Number of categories K is even
Analyze the data in Table 1.
|
Table 1. Data of 6 categories |
||||
|
Categories |
n-n |
n-s |
s-n |
s-s |
|
1: surely position 1 |
2 |
1 |
6 |
6 |
|
2: positon 1 |
5 |
1 |
10 |
3 |
|
3: probably position 1 |
15 |
13 |
24 |
12 |
|
4: probably position 2 |
21 |
21 |
7 |
21 |
|
5: position 2 |
6 |
7 |
3 |
6 |
|
6: surely position 2 |
1 |
7 |
0 |
2 |
Run the script sampling.py in Listing 3 in Python 3.11 virtual environment of Anaconda in Ubuntu/WSL/Windows
(py311) …/ex2afcrfiles$ python
sampling.py
Promt, which asks you to input data for <n, n> condition, is presented.
n_n =
Set the data from the column ‘n-n’ of Table 1 as follows
n_n = 2 5 15
21 6 1
After putting the data, which are separated by whitespaces,, push down the Enter key.
Then the data for <n, s> is required.
n_n = 2 5 15
21 6 1
n_s =
Repeat the above process.
After setting the data for <s, s> condition, the terminal is as follows.
n_n = 2 5 15 21 6 1
n_s = 1 1 13 21 7 7
s_n = 6 10 24 7 3 0
s_s = 6 3 12 21 6 2
After setting the data for <s, s>, MCMC sampling starts.
The program ends as shown below.
i_data.pkl was saved.
d_frame.pkl was saved.
K.pkl was saved.
End.
(py311) …/ex2afcrfiles$
MCMC sampling data were saved in the files i_data.pkl, d_frame.pkl, and K.pkl.
These files are used by analysis.py in the next step.
Run the script analysis.py in Python 3.11 virtual environment of Anaconda/Windows, not Anaconda/Ubuntu/WSL/Windows, as follows. In WSL, graphs would not be displayed.
(py311) …\ex2afcrfiles> python analysis.py
Posterior distribution of 
is displayed (Figure 1)
Figure 1
A point estimate of is shown as MAP
estimate.
Close the form window of Figure 1, then the
posterior distribution of 
is displayed (Figure 2).
Figure 2
Close the form window of Figure 2, then posterior distributions of category boundaries is shown (Figure 3).
Figure 3
In case of K even, the middle criterion is fixed at the origin 0, and shown as a red semicircle.
Close the form window of Figure 3, then the posterior distribution
of is shown (Figure 4).
Figure 4
Close the form window of Figure 4, then a graph which shows relation between the cumulative proportion of category responses and the predicted values is shown (Figure 5).
Figure 5
Close the form window of Figure 5, the program ends.
On the terminal window, the following message is shown.
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \
mu_s 0.947 0.192 0.583 1.306 0.001 0.001 25679.0 25362.0
sgm_s 1.269 0.184 0.931 1.618 0.001 0.001 19923.0 23961.0
r_hat
mu_s 1.0
sgm_s 1.0
I am calculating.
End.
(py311) …\Ex2AFCR\ex2afcrfiles>
Graphs 1 to 5 are saved in the files.
Number of categories K is odd.
Analyze the data in Table 2.
|
Table 2. Data of 7 categories |
||||
|
Categories |
n-n |
n-s |
s-n |
s-s |
|
1: surely position 1 |
2 |
0 |
5 |
0 |
|
2: positon 1 |
7 |
3 |
9 |
7 |
|
3: probably position 1 |
9 |
10 |
14 |
13 |
|
4: don’t know/same |
17 |
14 |
11 |
7 |
|
5: probably position 2 |
11 |
9 |
10 |
13 |
|
6: position 2 |
1 |
10 |
1 |
9 |
|
7: surely position 2 |
3 |
4 |
0 |
1 |
Run the script sampling.py in Listing 3 in Python 3.11 virtual environment of Anaconda in Ubuntu/WSL/Windows
(py311) …/ex2afcrfiles$ python
sampling.py
Promt, which asks you to input data for <n, n> condition, is presented.
n_n =
Set the data from the column ‘n-n’ of Table 2 as follows
n_n = 2 7 9
17 11 1 3
After putting the data, which are separated by whitespaces,, push down the Enter key.
Then the data for <n, s> is required.
n_n = 2 7 9
17 11 1 3
n_s =
Repeat the above process.
After setting the data for <s, s> condition, the terminal is as follows.
n_n = 2 7 9 17 11 1 3
n_s = 0 3 10 14 9 10 4
s_n = 5 9 14 11 10 1 0
s_s = 0 7 13 7 13 9 1
After setting the data for <s, s>, MCMC sampling starts.
The program ends as shown below.
i_data.pkl was saved.
d_frame.pkl was saved.
K.pkl was saved.
End.
(py311) …/ex2afcrfiles$
MCMC sampling data are saved in the files i_data.pkl, d_frame.pkl, and K.pkl.
These files are used by analysis.py in the next step.
Run the script analysis.py in Python 3.11 virtual environment of Anaconda/Windows, not Anaconda/Ubuntu/WSL/Windows, as follows. In WSL, graphs would not be displayed.
(py311) …\ex2afcrfiles> python
analysis.py
Posterior distribution of is displayed (Figure 6)
Figure 6
A point estimate of is shown as MAP
estimate.
Close the form window of Figure 6, then the
posterior distribution of is displayed (Figure 7).
Figure 7
Close the form window of Figure 7, then posterior distributions of category boundaries is shown (Figure 8).
Figure 8
Close the form window of Figure 8, then the
posterior distribution of is shown (Figure 9).
Figure 9
Close the form window of Figure 9, then a graph which shows relation between the cumulative proportion of category responses and the predicted values is shown (Figure 10).
Figure 10
Close the form window of Figure 10, the program ends.
On the terminal window, the following message is shown.
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \
mu_s 0.651 0.159 0.357 0.949 0.001 0.001 36605.0 28818.0
sgm_s 1.026 0.153 0.742 1.311 0.001 0.001 27605.0 28474.0
r_hat
mu_s 1.0
sgm_s 1.0
I am calculating.
End.
(py311) …\ex2afcrfiles>
Graphs 6 to 10 are saved in the files.
Brady, T. F., Robinson, M. M., Williams, J. R., & Wixted, J. T. (2023). Measuring memory is harder than you think: How to avoid problematic measurement practices in memory research. Psychonomic Bulletin & Review, 30, 421-449.
DeCarlo, L. T. (2010). On the statistical and theoretical basis of signal detection theory and extensions: Unequal variance, random coefficients, and mixture models. Journal of Mathematical Psychology, 54, 304-313.
Wickens, T. D. (2002). Elementary Signal Detection Theory. Oxford University Press.
Wixted, J. T. (2019). The forgotten history of signal detection theory. Journal of Experimental Psychology: Learning, Memory, & Cognition, 46, 201-233.
Listings
The script files below are archived into the file ex2afcrfiles.zip ,
which can be downloaded and used freely.
Listing 1. Stan script for the extended 2AFC rating task in case of K even (SDT_2AFC_KCat_even.stan)
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_cond;
array[K] int ss_cond;
}
transformed data {
int hK;
vector[K%/%2] alpha;
hK = K %/% 2;
for (i in 1:hK) {
alpha[i] = 1.0;
}
}
parameters {
real mu_s;
real<lower = 0.0>
sgm_s;
real b;
simplex[hK] preC;
}
transformed parameters {
vector[K] nn_theta;
vector[K] ns_theta;
vector[K] sn_theta;
vector[K] ss_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K+1] ss_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
real sgm_ss;
C[hK] = 0.0;
vsum = 0.0;
for (i in 1:hK-1) {
vsum
+= preC[i];
C[hK+i] = vsum / (1.0 - vsum);
C[hK-i] = -C[hK+i];
}
sgm_nn = sqrt(1.0 + 1.0);
sgm_ns = sqrt(1.0 +
square(sgm_s));
sgm_ss =
sqrt(square(sgm_s)*2);
nn_cum_p[1] = 0.0;
ns_cum_p[1] = 0.0;
sn_cum_p[1] = 0.0;
ss_cum_p[1] = 0.0;
nn_cum_p[K+1] = 1.0;
ns_cum_p[K+1] = 1.0;
sn_cum_p[K+1] = 1.0;
ss_cum_p[K+1] = 1.0;
for (i in 2:K) {
nn_cum_p[i] = normal_cdf(C[i-1] | b, sgm_nn);
ns_cum_p[i] = normal_cdf(C[i-1] | mu_s + b, sgm_ns);
sn_cum_p[i] = normal_cdf(C[i-1] | -mu_s + b, sgm_ns);
ss_cum_p[i] = normal_cdf(C[i-1] | b, sgm_ss);
}
for (i in 1:K) {
nn_theta[i] = nn_cum_p[i+1] - nn_cum_p[i];
ns_theta[i] = ns_cum_p[i+1] - ns_cum_p[i];
sn_theta[i] = sn_cum_p[i+1] - sn_cum_p[i];
ss_theta[i] = ss_cum_p[i+1] - ss_cum_p[i];
}
}
model {
mu_s ~ normal(0.0, 10.0);
sgm_s ~ exponential(0.1);
preC ~ dirichlet(alpha);
b ~ normal(0.0, 10.0);
nn_cond ~
multinomial(nn_theta);
ns_cond ~
multinomial(ns_theta);
sn_cond ~
multinomial(sn_theta);
ss_cond ~
multinomial(ss_theta);
}
Listing 2 . Stan script for the extended 2AFC rating task in case of K odd (SDT_2AFC_KCat_odd.stan)
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_cond;
array[K] int ss_cond;
}
transformed data {
int hK;
vector[(K%/%2)+1] alpha;
hK = K %/% 2;
for (i in 1:hK+1) {
alpha[i] = 1.0;
}
}
parameters {
real mu_s;
real<lower = 0.0>
sgm_s;
real b;
simplex[hK+1] preC;
}
transformed parameters {
vector[K] nn_theta;
vector[K] ns_theta;
vector[K] sn_theta;
vector[K] ss_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K+1] ss_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
real sgm_ss;
vsum = 0.0;
for (i in 1:hK) {
vsum
+= preC[i];
C[hK+i] = vsum / (1.0 - vsum);
C[hK+1-i] = -C[hK+i];
}
sgm_nn = sqrt(1.0 + 1.0);
sgm_ns = sqrt(1.0 +
square(sgm_s));
sgm_ss =
sqrt(square(sgm_s)*2);
nn_cum_p[1] = 0.0;
ns_cum_p[1] = 0.0;
sn_cum_p[1] = 0.0;
ss_cum_p[1] = 0.0;
nn_cum_p[K+1] = 1.0;
ns_cum_p[K+1] = 1.0;
sn_cum_p[K+1] = 1.0;
ss_cum_p[K+1] = 1.0;
for (i in 2:K) {
nn_cum_p[i] = normal_cdf(C[i-1] | b, sgm_nn);
ns_cum_p[i] = normal_cdf(C[i-1] | mu_s + b, sgm_ns);
sn_cum_p[i] = normal_cdf(C[i-1] | -mu_s + b, sgm_ns);
ss_cum_p[i] = normal_cdf(C[i-1] | b, sgm_ss);
}
for (i in 1:K) {
nn_theta[i] = nn_cum_p[i+1] - nn_cum_p[i];
ns_theta[i] = ns_cum_p[i+1] - ns_cum_p[i];
sn_theta[i] = sn_cum_p[i+1] - sn_cum_p[i];
ss_theta[i] = ss_cum_p[i+1] - ss_cum_p[i];
}
}
model {
mu_s ~ normal(0.0, 10.0);
sgm_s ~ exponential(0.1);
preC ~ dirichlet(alpha);
b ~ normal(0.0, 10.0);
nn_cond ~
multinomial(nn_theta);
ns_cond ~
multinomial(ns_theta);
sn_cond ~
multinomial(sn_theta);
ss_cond ~
multinomial(ss_theta);
}
Listing 3. Python script for MCMC sampling using the Stanscripts in Listings 1 and 2 (sampling.py)
import numpy as np
import stan
import arviz as az
import pandas as pd
import pickle
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
n_n = [int(v) for v in input('n_n =
').split()]
K = len(n_n)
if K < 3:
print('Number of categories
should be larger than 2.')
sys.exit()
n_s = [int(v) for v in input('n_s =
').split()]
if len(n_s) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
s_n = [int(v) for v in input('s_n =
').split()]
if len(s_n) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
s_s = [int(v) for v in input('s_s =
').split()]
if len(s_s) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
print('n_n =', n_n)
print('n_s =', n_s)
print('s_n =', s_n)
print('s_s =', s_s)
Data = {'K':K, 'nn_cond':n_n,
'ns_cond':n_s, 'sn_cond':s_n, 'ss_cond':s_s}
with open('SDT_2AFC_KCat_odd.stan' if
K % 2 == 1 else
'SDT_2AFC_KCat_even.stan',
'r') as f:
sm = stan.build(f.read(),
data = Data)
fit = sm.sample(num_samples=10000)
print('fit\n', fit)
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
with open('i_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print('i_data.pkl was saved.')
df = fit.to_frame()
with open('d_frame.pkl', 'wb') as f:
pickle.dump(df, f)
print('d_frame.pkl was saved.')
with open('K.pkl', 'wb') as f:
pickle.dump((K, n_n, n_s,
s_n, s_s), f)
print('K.pkl was saved.')
print('\nEnd.\n')
Listing 4. Python script for analysis of MCMC sample (analysis.py)
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pandas as pd
import pickle
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
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]
with open('i_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
with open('d_frame.pkl', 'rb') as f:
fit = pickle.load(f)
with open('K.pkl', 'rb') as f:
(K, n_n, n_s, s_n, s_s) =
pickle.load(f)
smry = az.summary(inference_data,
var_names = ['mu_s', 'sgm_s'])
print('Summary =\n',smry)
mu_s_map = CalcMAPEst(fit['mu_s'])[0]
sb.kdeplot(fit['mu_s'])
plt.xlabel('$\mu_s$')
plt.title('$\mu_s$(MAP) =
{0:.3f}'.format(mu_s_map))
plt.savefig('Fig_mu.png')
plt.show()
sgm_s_map =
CalcMAPEst(fit['sgm_s'])[0]
sb.kdeplot(fit['sgm_s'])
plt.title('$\sigma_s$(MAP) =
{0:.3f}'.format(sgm_s_map))
plt.xlabel('$\sigma_s$')
plt.savefig('Fig_sgm.png')
plt.show()
print('\n I am calculating.')
s_title = ''
if K % 2 == 1:
for k in range(K-1):
Csmpl = fit[f'C.{k+1}']
sb.kdeplot(Csmpl)
c_map = CalcMAPEst(Csmpl)[0]
s_title += f'C{k+1}={c_map:.2f}'
if k
< K-2:
s_title += ', '
else:
for k in range(K-1):
if k
!= (K//2) - 1:
Csmpl = fit[f'C.{k+1}']
sb.kdeplot(Csmpl)
c_map = CalcMAPEst(Csmpl)[0]
s_title += f'C{k+1}={c_map:.2f}'
if k < K-2:
s_title += ', '
else:
plt.plot([0], [0], marker = 'o', markersize = 15, c = 'r')
s_title += f'C{k+1}=0, '
plt.title('MAP estimates\n' + s_title)
plt.savefig('FigC.png')
plt.show()
b_map = CalcMAPEst(fit['b'])[0]
sb.kdeplot(fit['b'])
plt.title('b(MAP) =
{0:.3f}'.format(b_map))
plt.savefig('Fig_b.png')
plt.show()
cum_n_n = np.cumsum(n_n)
pcum_n_n = cum_n_n/cum_n_n[-1]
cum_n_s = np.cumsum(n_s)
pcum_n_s = cum_n_s / cum_n_s[-1]
cum_s_n = np.cumsum(s_n)
pcum_s_n = cum_s_n / cum_s_n[-1]
cum_s_s = np.cumsum(s_s)
pcum_s_s = cum_s_s / cum_s_s[-1]
est_pcum_n_n = np.empty(K-1)
for k in range(K-1):
est_pcum_n_n[k] =
np.median(fit[f'nn_cum_p.{k+2}'])
est_pcum_n_s = np.empty(K-1)
for k in range(K-1):
est_pcum_n_s[k] =
np.median(fit[f'ns_cum_p.{k+2}'])
est_pcum_s_n = np.empty(K-1)
for k in range(K-1):
est_pcum_s_n[k] =
np.median(fit[f'sn_cum_p.{k+2}'])
est_pcum_s_s = np.empty(K-1)
for k in range(K-1):
est_pcum_s_s[k] =
np.median(fit[f'ss_cum_p.{k+2}'])
plt.plot(est_pcum_n_n, pcum_n_n[:-1],
label = 'n_n')
plt.plot(est_pcum_n_s, pcum_n_s[:-1],
label = 'n_s')
plt.plot(est_pcum_s_n, pcum_s_n[:-1],
label = 's_n')
plt.plot(est_pcum_s_s, pcum_s_s[:-1],
label = 's_s')
plt.plot([0,1], [0,1], c = 'k', ls =
'--', label = 'Obs.=Est.')
plt.xlabel('Est.Cum.P')
plt.ylabel('Obs.Cum.P')
plt.legend()
plt.title('Cumulative Proportions')
plt.savefig('FigCumEstObs.png')
plt.show()
print('End.')