Up

Two-Alternative-Forced-Choice Task of Signal Detection Theory with Confidence Rating

Model and Bayesian Analysis

Yasuharu Okamoto, 2019

Script and data files are archived in a file SDT_2AFC.zip,and the files for Jupyter Notebook in a file SDT_2AFC_jn.zip. These archived files can be freely downloaded and used.

 

Python scripts with Stan scripts for Bayesian analyses of data from the standard 2AFC (alternative-forced-choice) task and ones with confidence rating of signal detection theory (SDT) were presented. Klein (2001, p. 1436) recommends using more than two response categories.

Analysis of the standard 2AFC task needs the assumption on the variance of sensation of the signal stimulus. For extension of 2AFC task, which does not need assumption on the variance of sensation of the signal stimulus, check the website https://osf.io/ts2eq/, and the website http://y-okamoto-psy1949.la.coocan.jp/Python/en2/Ex2AFCR/.

On the standard 2AFC task, an observer chooses between the two response categories, e.g. which one is the signal stimulus?. On a 2AFC task with confidence rating, an observer is asked as follows:

In a three rating category task, e.g., Which one is your judgment?, The left stimulus is the signal., Dont know., or The right stimulus is the signal.’”

In a four rating category task, e.g., Which one is your judgment?, The left stimulus is the signal., Maybe, the left stimulus is the signal., Maybe, the right stimulus is the signal,, or The right stimulus is the signal.’”

In the following, models are presented first, then scripts and examples of their uses are shown.

 

The Basic Model

Sensations of the noise and signal stimuli are denoted by random variables  and , which have the following normal distributions:

To set the origin and unit, it is assumed that

In a 2AFC task, a pair of noise (denoted by N) and signal (denoted by S) stimuli are presented in one of the two conditions, condition <N, S> and condition <S, N>. In conditon <N, S>, N is presented at position 1 and S in position 2. In condition <S, N>, S in position 1 and N in position 2. When positions mean temporal ones, position 1 precedes position 2 in temporal order. When positions mean spatial ones, position 1 is the left side and position 2 the right side, for example. And so on.

There may be biases in psychophysical judgment. Total effect of the biases is represented by , which is measured from position 2.

Hence, in condition <N, S>, sensations have the following distributions

In condition <S, N>,

Psychophysical judgment is assumed to be done based on the difference between the sensations . When the difference is between category boundaries  and , observer chooses category  as his/her judgment. It is assumed that when number of categories is ,

For example, for ,

Category 1 means that signal is perceived to be weaker than noise, so wrongly judged to be noise. Category 2 means that I dont know. Category 3 means that signal is perceived to be stronger than noise, so correctly judged to be signal.

In codition <N, S>,

Denote the probability that rating judgment is category  in codition <N, S> by , then we have

where  is the cumulative standard normal distribution function.

In codition <S, N>,

Denote the probability that rating judgment is category  in codition <S, N> by , then we have

Because the independent variables in function  are given in ratios of parameters, the unit must be given to identify parameter values. The popular method is to set , i.e. equal-variance model . In this case,

This assumption may be replaced more generally by

From either of these assumptions, we have

In codition <N, S>, we have

In codition <S, N>, we have

Parameters  and  appear as difference, so an assumption to set a origin is necessary. To set an origin, the following constraints are set on category boundaries.

For even number of rating categories , set

For odd number of rating categories , set

 

Stan scripts and Python scripts were developed for three cases of number of category boundaries,  i.e. the standard task,  is odd, and  is even(> 2).

 

In case of :

Stan script in this case is shown in Listing A.1.

 

Listing A.1  Stan script for K=2 (SDT2AFC2Cat.stan).

data {

    int n_NS[2];

    int n_SN[2];

}

parameters {

    real mu_s;

    real tau;

}

transformed parameters {

    simplex[2] p_NS;

    simplex[2] p_SN;

    p_NS[1] = Phi(-(mu_s - tau)/sqrt(2.0));

    p_NS[2] = 1 - p_NS[1];

    p_SN[1] = Phi(-(mu_s + tau)/sqrt(2.0));

    p_SN[2] = 1 - p_SN[1];

}

model {

    mu_s ~ normal(0.0, 100.0);

    tau ~ normal(0.0, 100.0);

    n_NS ~ multinomial(p_NS);

    n_SN ~ multinomial(p_SN);

}

 

Python script which uses the Stan script above is shown in Listing A.2.

 

Listing A.2  Python script which uses the Stan script above (SDT2AFC2Cat.py).

import numpy as np

import pystan

import matplotlib.pyplot as plt

import csv

 

flnm = input('Input data file (*.csv) = ')

with open(flnm, 'r') as f:

    data = [v for v in csv.reader(f)]

 

fout_nm = input('Output text file (*.txt) = ')

fout = open(fout_nm, 'w')

fout.write('Input Data File = {}\n'.format(flnm))

   

Freq_NS = []

Freq_NS.append(int(data[1][1]))

Freq_NS.append(int(data[1][2]))

Freq_SN = []

Freq_SN.append(int(data[2][1]))

Freq_SN.append(int(data[2][2]))

print(Freq_NS)

print(Freq_SN)

fout.write('<N, S>: {}'.format(Freq_NS))

fout.write('\n<S, N>: {}'.format(Freq_SN))

TN_NS = sum(Freq_NS)

Rating_NS = np.array(Freq_NS)

 

TN_SN = sum(Freq_SN)

Rating_SN = np.flip(np.array(Freq_SN))

 

Data = {'n_NS': Rating_NS, 'n_SN': Rating_SN}

 

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

 

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

 

print(fit)

fout.write('\n{}\n'.format(fit))

 

mu = fit['mu_s']

tau = fit['tau']

p_NS = fit['p_NS']

p_SN = fit['p_SN']

 

mu_L05 = np.percentile(mu, 5)

mu_med = np.percentile(mu, 50)

mu_U95 = np.percentile(mu, 95)

fout.write("\nd': \nMed. = {0:<.3f},  90%CI = [{1:<.3f}, {2:<.3f}".

           format(mu_med, mu_L05, mu_U95))

 

tau_L05 = np.percentile(tau, 5)

tau_med = np.percentile(tau, 50)

tau_U95 = np.percentile(tau, 95)

fout.write('\ntau: \nMed. = {0:<.3f},  90%CI ~ [{1:<.3f}, {2:<.3f}]'.

           format(tau_med, tau_L05, tau_U95))

 

p1 = np.median(p_NS, axis = 0)

p2 = np.median(p_SN, axis = 0)

p2 = [p2[1], p2[0]]

 

plt.hist(mu)

plt.xlabel("d'($\mu_S$)", fontsize = 14)

plt.title("Posterior Distribution of d'($\mu_S$)" + \

          '\nMed. = {0:<.3f},  90%CI = [{1:M<.3f}, {2:<.3f}]'.

          format(mu_med, mu_L05, mu_U95), fontsize = 16)

plt.show()

 

plt.hist(tau)

plt.xlabel(r'$\tau$', fontsize = 18)

plt.title(r'Posterior Distribution of $\tau$' + \

          '\nMed. = {0:<.3f},  90%CI = [{1:<.3f}, {2:<.3f}]'.

          format(tau_med, tau_L05, tau_U95), fontsize = 16)

plt.show()

 

xcat = [1, 2]

xlabels = ['Cat-1', 'Cat-2']

y1 = np.array(Freq_NS) / TN_NS

y2 = np.array(Freq_SN) / TN_SN

 

plt.plot(xcat, y1, 'b:', linewidth = 10, label = 'Data/<N, S>')

plt.plot(xcat, p1, 'g-', linewidth = 3, label = 'Model/<N, S>')

plt.plot(xcat, y2, 'r:', linewidth = 10, label = 'Data/<S, N>')

plt.plot(xcat, p2, 'm-', linewidth = 3,label = 'Model/<S, N>')

plt.xticks(xcat, xlabels, fontsize = 14)

plt.xlim(0.8, 2.2)

plt.xlabel('Rating Category', fontsize = 14)

plt.ylabel('Probability/Proportion', fontsize = 14)

plt.title('Rating in Conditions <N, S> and <S, N>', fontsize = 18)

plt.legend(loc = 'upper center', fontsize = 10)

plt.show()

     

fout.close()

print('\n' + fout_nm + ' was saved.\n')

 

Files in Listing A1 and A.2 and a sample data file are archived into a zip file Rating2Cat.zip, which can be down loaded.

Run the above script SDT2AFC2Cat.py, names of the input data file and an output text file are required to be set (Figure A.1).

Figure A.1

 

Name of an output text file is arbitrary. Format of the input data file should be CSV file(Figure .A.2).

A.2 データファイルData2Cat.csv

 

In the first row, names for columns are set. The first column shows conditions, <N, S> and <S, N>, the second column shows frequencies of rating judgment of category 1, and the third column shows frequencies of rating judgment of category 2.

In the second row, data from condition <N, S> are set, and in the third row, data from condition <S, N>.

When the file names are set, the program starts calculation. After sampling by Stan ends, posterior distribution of parameter d is presented (Figure A.3).

A.3

 

Closing the form in Figure A.3, posterior distribution of  is shown (Figure A.4).

Figure A.4

 

Closing the form in Figure A.4, graphs of data and model predictions are shown (Figure A.5).

Figure A.5

 

Dotted lines represent proportions of rating categories calculated from the data, and solid lines represent estimated probabilities of rating categories, which are medians of posterior distributions. Cat-1 denotes judgment of category 1, and Cat-2 denotes judgment of category 2.

Closing the form of Figure A.5, the program ends. After the program ends, output text file can be opened. The content is like this:

 

Input Data File = Data2CAt.csv

<N, S>: [30, 70]

<S, N>: [84, 16]

Inference for Stan model: anon_model_88e26578a5e2b024ad50d7db907e822c.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

mu_s      1.08  2.5e-3   0.14   0.79   0.98   1.08   1.17   1.35   3130    1.0

tau       0.33  2.5e-3   0.14   0.06   0.24   0.33   0.42   0.62   3288    1.0

p_NS[0]    0.3  7.0e-4   0.04   0.22   0.27    0.3   0.33   0.39   4000    1.0

p_NS[1]    0.7  7.0e-4   0.04   0.61   0.67    0.7   0.73   0.78   4000    1.0

p_SN[0]   0.16  6.9e-4   0.04    0.1   0.14   0.16   0.19   0.24   2919    1.0

p_SN[1]   0.84  6.9e-4   0.04   0.76   0.81   0.84   0.86    0.9   2919    1.0

lp__    -106.0    0.02   0.98 -108.5 -106.4 -105.7 -105.3 -105.0   1677    1.0

 

Samples were drawn using NUTS at Thu Oct 24 19:29:23 2019.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

d':

Med. = 1.079,  90%CI = [0.842, 1.309

tau:

Med. = 0.327,  90%CI ~ [0.101, 0.568]

 

 

In case of  being odd:

 

Stan script for  being odd is shown in Listing B1.

 

Listing B.1  Stan script for  being odd (SDT2AFCOddCat.stan).

functions {

    real half_normal_lpdf(real y, real sgm){

        if (y > 0.0){

            return normal_lpdf(y | 0.0, sgm);

        } else {

            return log(0.0);

        }

    }

}

data {

    int K;

    int n_NS[K];

    int n_SN[K];

}

 

parameters {

    real mu_s;

    real tau;

    real<lower = 0.0> theta[(K - 1) / 2];

}

transformed parameters {

    real C[K - 1];

    simplex[K] p_NS;

    simplex[K] p_SN;

    C[(K - 1) / 2 + 1] = theta[1];

    C[(K - 1) / 2] = -C[(K - 1) / 2 + 1];

    for (k in 2: ((K - 1) / 2)){

        C[(K - 1) /2 + k] = C[(K - 1) / 2 + k - 1] + theta[k];

        C[(K - 1) / 2 - k + 1] = -C[(K - 1) / 2 + k];

    }

    p_NS[1] = Phi((C[1] - (mu_s - tau))/sqrt(2.0));

    p_NS[K] = 1 - Phi((C[K - 1] - (mu_s - tau))/sqrt(2.0));

 

    for (k in 2:(K-1)) {

        p_NS[k] = Phi((C[k] - (mu_s - tau))/sqrt(2.0)) -

                   Phi((C[k - 1] - (mu_s - tau))/sqrt(2.0));

    }

 

    p_SN[1] = Phi((C[1] - (mu_s + tau))/sqrt(2.0));

    p_SN[K] = 1 - Phi((C[K - 1] - (mu_s + tau))/sqrt(2.0));

   

    for (k in 2:(K-1)) {

        p_SN[k] = Phi((C[k] - (mu_s + tau))/sqrt(2.0)) -

                   Phi((C[k - 1] - (mu_s + tau))/sqrt(2.0));

    }

}

model {

    for (k in 1:((K - 1) / 2)) {

        theta[k] ~ half_normal(1000.0);

    }

    mu_s ~ normal(0.0, 100.0);

    tau ~ normal(0.0, 100.0);

    n_NS ~ multinomial(p_NS);

    n_SN ~ multinomial(p_SN);

}

 

Python script which uses the above Stan script SDT2AFCOddCat.stan is shown in Listing B.2.

 

Listing B.2  Python script which uses Stan script SDT2AFCOddCat.stan (SDT2AFCOddCat.py)

import numpy as np

import pystan

import matplotlib.pyplot as plt

import csv

 

flnm = input('Input data file (*.csv) = ')

with open(flnm, 'r') as f:

    data = [v for v in csv.reader(f)]

fout_nm = input('Output text file (*.txt) = ')

fout = open(fout_nm, 'w')

fout.write('\nInput Data File = {}\n'.format(flnm))

 

K = len(data[0]) - 1

print('K = ', K)

if (K % 2) == 0:

    print('The number of response categories is not odd.')

    import sys

    sys.exit()

   

Freq_NS = []

for k in range(K):

    Freq_NS.append(int(data[1][1 + k]))

Freq_SN = []

for k in range(K):

    Freq_SN.append(int(data[2][1 + k]))

print(Freq_NS)

print(Freq_SN)

fout.write('\n<N, S>\n{}\n'.format(Freq_NS))

fout.write('\n<S, N>\n{}\n'.format(Freq_SN))          

TN_NS = sum(Freq_NS)

Rating_NS = Freq_NS

TN_SN = sum(Freq_SN)

Rating_SN = np.flip(np.array(Freq_SN))

 

Data = {'K': K, 'n_NS': Rating_NS, 'n_SN': Rating_SN}

 

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

 

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

 

print(fit)

fout.write('\n{}\n'.format(fit))

 

mu = fit['mu_s']

tau = fit['tau']

C = fit['C']

p_NS = fit['p_NS']

p_SN = fit['p_SN']

 

mu_L05 = np.percentile(mu, 5)

mu_med = np.percentile(mu, 50)

mu_U95 = np.percentile(mu, 95)

fout.write("\nd':\nMedian = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]\n".

           format(mu_med, mu_L05, mu_U95))

 

 

tau_L05 = np.percentile(tau, 5)

tau_med = np.percentile(tau, 50)

tau_U95 = np.percentile(tau, 95)

fout.write('\ntau:\nMedian = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]\n'.

           format(tau_med, tau_L05, tau_U95))

 

C_L05 = np.zeros(K-1)

C_med = np.zeros(K-1)

C_U95 = np.zeros(K-1)

for k in range(K-1):

    C_L05[k] = np.percentile(C.T[k], 5)

    C_med[k] = np.percentile(C.T[k], 50)

    C_U95[k] = np.percentile(C.T[k], 95)

for k in range(K-1):

    fout.write(('\nC[{0}]:\n' + \

               'Median = {1:<.3f},   90%CI = [{2:<.3f}, {3:<.3f}]\n').

               format(k+1, C_med[k], C_L05[k], C_U95[k]))

 

plt.hist(mu)

plt.xlabel("d'($\mu_S$)", fontsize = 14)

plt.title("Posterior Distribution of d'($\mu_S$)" + \

          '\nMed. = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]'.

          format(mu_med, mu_L05, mu_U95), fontsize = 16)

plt.show()

 

plt.hist(tau)

plt.xlabel(r'$\tau$', fontsize = 18)

plt.title(r'Posterior Distribution of $\tau$' + \

          '\nMed. = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]'.

          format(tau_med, tau_L05, tau_U95), fontsize = 16)

plt.show()

 

for k in range(K-1):

    plt.hist(C.T[k], label = 'C{}'.format(k+1))

plt.title('Posterior Distributions of Category Boundaries', fontsize = 18)

plt.legend()

plt.show()

 

p1 = np.median(p_NS, axis = 0)

p2 = np.median(p_SN, axis = 0)

p2 = np.flip(p2)

 

xcat = np.arange(1, K+0.1, 1)

xlabels = ['{}'.format(int(v)) for v in xcat]

y1 = np.array(Freq_NS) / TN_NS

y2 = np.array(Freq_SN) / TN_SN

   

plt.plot(xcat, y1, 'b--', linewidth = 3, alpha = 0.7, label = 'Data/<N, S>')

plt.plot(xcat, p1, 'g-', linewidth = 3, alpha = 0.7, label = 'Model/<N, S>')

plt.plot(xcat, y2, 'r--', linewidth = 3, alpha = 0.7, label = 'Data/<S, N>')

plt.plot(xcat, p2, 'm-', linewidth = 3, alpha = 0.7, label = 'Model/<S, N>')

 

plt.xticks(xcat, xlabels, fontsize = 14)

plt.xlabel('Rating Category', fontsize = 14)

plt.ylabel('Probability/Proportion', fontsize = 14)

plt.title('Ratings in Conditions <N, S> and <S, N>', fontsize = 18)

plt.legend(loc = 'upper center', fontsize = 10)

plt.show()

 

fout.close()

print('\n{} was saved.\n'.format(fout_nm))

 

Files in Listing B.1 and Listing B.2 and sample data files are archived into a zip file RatingOddCat.zip, which can be downloaded.

Run SDT2AFCOddCat.py in Listing B.2, names of the input data file and an output text file are required to be set (Figure B.1).

Figure B.1

 

Output text file name is arbitrary. The input datafile should be CSV format (Figure B.2).

Figure B.2  Input data file Data3Cat.csv

 

In the first row, names for columns are set. The first column shows conditions, <N, S> and <S, N>, the second column shows frequencies of rating judgment of category 1 (e.g., The signal is on the left), the third column shows frequencies of rating judgment of category 2 (e.g., Dont know), and the fourth column shows frequencies of rating judgment of category 3 (e.g., The signal is on the right).

In the second row, data from condition <N, S> are set, and in the third row, data from condition <S, N>.

When the file names are set, the program starts calculation. After sampling by Stan ends, posterior distribution of parameter d is presented (Figure B.3).

Figure B.3

 

Closing the form in Figure B.3, posterior distribution of  is shown (Figure B.4).

Figure B.4

 

Closing the form in Figure B.4, posterior distribution of category boundaries are shown (Figure B.5).

Figure B.5

 

Closing the form in Figure B,5, graphs of data and model predictions are shown (Figure B.6).

Figure B.6

 

Dotted lines represent proportions of rating categories calculated from the data, and solid lines represent estimated probabilities of rating categories, which are medians of posterior distributions. Rating categories 1 to 3 correspond to the columns 2 to 4 in Figure B.2, respectively.

Closing the form of Figure B.6, the program ends. After the program ends, output text file can be opened. The content is like this:

 

Input Data File = Data3Cat.csv

 

<N, S>

[20, 26, 54]

 

<S, N>

[77, 12, 11]

 

Inference for Stan model: anon_model_4f4c02dce251e874ae19399d5b4bfa08.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

mu_s       1.05  2.2e-3   0.13    0.8   0.96   1.05   1.14    1.3   3528    1.0

tau         0.4  2.0e-3   0.13   0.16   0.31    0.4   0.49   0.64   4000    1.0

theta[0]   0.46  1.1e-3   0.07   0.34   0.41   0.46    0.5    0.6   3652    1.0

C[0]      -0.46  1.1e-3   0.07   -0.6   -0.5  -0.46  -0.41  -0.34   3652    1.0

C[1]       0.46  1.1e-3   0.07   0.34   0.41   0.46    0.5    0.6   3652    1.0

p_NS[0]    0.22  6.2e-4   0.04   0.15   0.19   0.22   0.24    0.3   4000    1.0

p_NS[1]    0.23  5.2e-4   0.03   0.17   0.21   0.23   0.25    0.3   4000    1.0

p_NS[2]    0.55  7.5e-4   0.05   0.46   0.52   0.55   0.58   0.65   4000    1.0

p_SN[0]    0.09  4.5e-4   0.02   0.05   0.07   0.09   0.11   0.15   3004    1.0

p_SN[1]    0.15  4.2e-4   0.03    0.1   0.13   0.15   0.17   0.21   4000    1.0

p_SN[2]    0.76  6.6e-4   0.04   0.67   0.73   0.76   0.79   0.83   3889    1.0

lp__     -181.3    0.03   1.21 -184.5 -181.9 -181.0 -180.4 -179.9   1907    1.0

 

Samples were drawn using NUTS at Fri Oct 25 09:04:47 2019.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

d':

Median = 1.050,   90%CI = [0.838, 1.256]

 

tau:

Median = 0.403,   90%CI = [0.194, 0.607]

 

C[1]:

Median = -0.455,   90%CI = [-0.574, -0.352]

 

C[2]:

Median = 0.455,   90%CI = [0.352, 0.574]

 

 

Next, consider another example, which uses K=5 number of rating categories (Figure B.7).

Figure B.7  Data5Cat.csv

 

Run SDT2AFCOddCat.py in Listing B.2, and set names of the input data file and an output text file as in Figure B.8.

Figure B.8

 

When the file names are set, the program starts calculation. After sampling by Stan ends, posterior distribution of parameter d is presented (Figure B.9).

Figure B.9

 

Closing the form in Figure B.9, posterior distribution of  is shown (Figure B.10).

Figure B.10

 

Closing the form in Figure B.10, posterior distribution of category boundaries are shown (Figure B.11).

Figure B.11

 

Closing the form in Figure B,11, graphs of data and model predictions are shown (Figure B.12).

Figure B.12

 

Dotted lines represent proportions of rating categories calculated from the data, and solid lines represent estimated probabilities of rating categories, which are medians of posterior distributions. Rating categories 1 to 5 correspond to the columns 2 to 6 in Figure B.7, respectively.

Closing the form of Figure B.12, the program ends. After the program ends, output text file can be opened. The content is like this:

 

Input Data File = Data5Cat.csv

 

<N, S>

[21, 10, 5, 18, 46]

 

<S, N>

[66, 10, 7, 12, 5]

 

Inference for Stan model: anon_model_4f4c02dce251e874ae19399d5b4bfa08.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

mu_s       0.93  1.9e-3   0.12   0.69   0.85   0.93   1.01   1.17   4000    1.0

tau        0.38  1.8e-3   0.12   0.16   0.31   0.38   0.46   0.61   4000    1.0

theta[0]   0.14  6.1e-4   0.04   0.07   0.11   0.14   0.16   0.22   4000    1.0

theta[1]   0.58  1.1e-3   0.07   0.44   0.53   0.58   0.63   0.73   4000    1.0

C[0]      -0.72  1.2e-3   0.08  -0.87  -0.77  -0.72  -0.66  -0.57   4000    1.0

C[1]      -0.14  6.1e-4   0.04  -0.22  -0.16  -0.14  -0.11  -0.07   4000    1.0

C[2]       0.14  6.1e-4   0.04   0.07   0.11   0.14   0.16   0.22   4000    1.0

C[3]       0.72  1.2e-3   0.08   0.57   0.66   0.72   0.77   0.87   4000    1.0

p_NS[0]    0.19  5.6e-4   0.04   0.13   0.16   0.19   0.21   0.27   4000    1.0

p_NS[1]    0.13  2.6e-4   0.02    0.1   0.12   0.13   0.14   0.16   4000    1.0

p_NS[2]    0.07  3.2e-4   0.02   0.04   0.06   0.07   0.09   0.12   4000    1.0

p_NS[3]    0.16  3.2e-4   0.02   0.12   0.15   0.16   0.17    0.2   4000    1.0

p_NS[4]    0.45  7.5e-4   0.05   0.36   0.42   0.45   0.48   0.55   4000    1.0

p_SN[0]    0.08  3.6e-4   0.02   0.04   0.06   0.08   0.09   0.12   3523    1.0

p_SN[1]    0.08  2.0e-4   0.01   0.05   0.07   0.08   0.08    0.1   4000    1.0

p_SN[2]    0.05  2.4e-4   0.01   0.03   0.04   0.05   0.06   0.08   4000    1.0

p_SN[3]    0.13  3.0e-4   0.02    0.1   0.12   0.13   0.15   0.17   4000    1.0

p_SN[4]    0.66  7.1e-4   0.05   0.57   0.63   0.66   0.69   0.75   4000    1.0

lp__     -270.3    0.03   1.47 -274.0 -271.1 -270.0 -269.3 -268.5   2012    1.0

 

Samples were drawn using NUTS at Fri Oct 25 09:44:06 2019.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

d':

Median = 0.930,   90%CI = [0.733, 1.131]

 

tau:

Median = 0.383,   90%CI = [0.192, 0.577]

 

C[1]:

Median = -0.716,   90%CI = [-0.852, -0.596]

 

C[2]:

Median = -0.137,   90%CI = [-0.208, -0.084]

 

C[3]:

Median = 0.137,   90%CI = [0.084, 0.208]

 

C[4]:

Median = 0.716,   90%CI = [0.596, 0.852]

 

 

In case of  being even greater than 2:

 

Stan script for K being even greater that 2 is shown in Listing C.1.

 

Listing C.1  Stan script for K being even greater than 2 (SDT2AFCEvenCat.stan)

functions {

    real half_normal_lpdf(real y, real sgm){

        if (y > 0.0){

            return normal_lpdf(y | 0.0, sgm);

        } else {

            return log(0.0);

        }

    }

}

data {

    int K;

    int n_NS[K];

    int n_SN[K];

}

 

parameters {

    real mu_s;

    real tau;

    real<lower = 0.0> theta[K/2 - 1];

}

transformed parameters {

    real C[K - 1];

    simplex[K] p_NS;

    simplex[K] p_SN;

    C[K/2] = 0.0;

    for (k in 1: (K/2 -1)){

        C[K/2 + k] = C[K/2 + k - 1] + theta[k];

        C[K/2 - k] = -C[K/2 + k];

    }

    p_NS[1] = Phi((C[1] - (mu_s - tau))/sqrt(2.0));

    p_NS[K] = 1 - Phi((C[K - 1] - (mu_s - tau))/sqrt(2.0));

 

    for (k in 2:(K-1)) {

        p_NS[k] = Phi((C[k] - (mu_s - tau))/sqrt(2.0)) -

                   Phi((C[k - 1] - (mu_s - tau))/sqrt(2.0));

    }

 

    p_SN[1] = Phi((C[1] - (mu_s + tau))/sqrt(2.0));

    p_SN[K] = 1 - Phi((C[K - 1] - (mu_s + tau))/sqrt(2.0));

   

    for (k in 2:(K-1)) {

        p_SN[k] = Phi((C[k] - (mu_s + tau))/sqrt(2.0)) -

                   Phi((C[k - 1] - (mu_s + tau))/sqrt(2.0));

    }

}

model {

    for (k in 1:(K/2 - 1)) {

        theta[k] ~ half_normal(1000.0);

    }

    mu_s ~ normal(0.0, 100.0);

    tau ~ normal(0.0, 100.0);

    n_NS ~ multinomial(p_NS);

    n_SN ~ multinomial(p_SN);

}

 

Python script which uses the Stan script above is shown in Listing C.2.

 

Listing C.2  Python script which uses the Stan script in Listing C.2 (SDT2AFCEvenCat.py)

#

#       Number of rating categories should be even.

#

#           Yasuharu Okamoto, 2019.10

#

import numpy as np

import pystan

import matplotlib.pyplot as plt

import csv

import pickle

 

flnm = input('Input data file (*.csv) = ')

with open(flnm, 'r') as f:

    data = [v for v in csv.reader(f)]

 

fout_nm = input('Output text file (*.txt) = ')

fout = open(fout_nm, 'w')

fout.write('Input Data File = {}\n'.format(flnm))

 

K = len(data[0]) - 1

print('K = ', K)

if (K % 2) == 1:

    print('The number of rating categories is not even!')

    import sys

    sys.exit()

   

Freq_NS = []

for k in range(K):

    Freq_NS.append(int(data[1][1 + k]))

Freq_SN = []

for k in range(K):

    Freq_SN.append(int(data[2][1 + k]))

print(Freq_NS)

print(Freq_SN)

fout.write('\n<N, S>\n{}\n'.format(Freq_NS))

fout.write('\n<S, N>\n{}\n'.format(Freq_SN))  

TN_NS = sum(Freq_NS)

Rating_NS = Freq_NS

print(Rating_NS)

TN_SN = sum(Freq_SN)

Rating_SN = np.flip(np.array(Freq_SN))

 

Data = {'K': K, 'n_NS': Rating_NS, 'n_SN': Rating_SN}

 

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

 

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

 

print(fit)

fout.write('\n{}\n'.format(fit))

 

mu = fit['mu_s']

tau = fit['tau']

C = fit['C']

p_NS = fit['p_NS']

p_SN = fit['p_SN']

 

mu_L05 = np.percentile(mu, 5)

mu_med = np.percentile(mu, 50)

mu_U95 = np.percentile(mu, 95)

fout.write("\nd':\nMedian = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]\n".

           format(mu_med, mu_L05, mu_U95))

 

tau_L05 = np.percentile(tau, 5)

tau_med = np.percentile(tau, 50)

tau_U95 = np.percentile(tau, 95)

fout.write('\ntau:\nMedian = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]\n'.

           format(tau_med, tau_L05, tau_U95))

 

C_L05 = np.zeros(K-1)

C_med = np.zeros(K-1)

C_U95 = np.zeros(K-1)

for k in range(K-1):

    C_L05[k] = np.percentile(C.T[k], 5)

    C_med[k] = np.percentile(C.T[k], 50)

    C_U95[k] = np.percentile(C.T[k], 95)

for k in range(K-1):

    fout.write(('\nC[{0}]:\n' + \

               'Median = {1:<.3f},   90%CI = [{2:<.3f}, {3:<.3f}]\n').

               format(k+1, C_med[k], C_L05[k], C_U95[k]))

 

plt.hist(mu)

plt.xlabel("d'($\mu_S$)", fontsize = 14)

plt.title("Posterior Distribution of d'($\mu_S$)" + \

          '\nMed. = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]'.

          format(mu_med, mu_L05, mu_U95), fontsize = 16)

plt.show()

 

plt.hist(tau)

plt.xlabel(r'$\tau$', fontsize = 18)

plt.title(r'Posterior Distribution of $\tau$' + \

          '\nMed. = {0:<.3f},   90%CI = [{1:<.3f}, {2:<.3f}]'.

          format(tau_med, tau_L05, tau_U95), fontsize = 16)

plt.show()

 

for k in range(K-1):

    if k + 1 != K/2:

        plt.hist(C.T[k], label = 'C{}'.format(k+1))

plt.title('Posterior Distributions of Category Boundaries' + \

          '\nC{} = 0 fixed'.format(K//2), fontsize = 16)

plt.legend(loc = 'lower center')

plt.show()

 

xcat = np.arange(1, K+0.1, 1)

xlabels = ['{}'.format(int(v)) for v in xcat]

y1 = np.array(Freq_NS) / TN_NS

y2 = np.array(Freq_SN) / TN_SN

 

p1 = np.median(p_NS, axis = 0)

p2 = np.median(p_SN, axis = 0)

p2 = np.flip(p2)

 

plt.plot(xcat, y1, 'b--', linewidth = 3, alpha = 0.7, label = 'Data/<N, S>')

plt.plot(xcat, p1, 'g-', linewidth = 3, alpha = 0.7, label = 'Model/<N, S>')

plt.plot(xcat, y2, 'r--', linewidth = 3, alpha = 0.7, label = 'Data/<S, N>')

plt.plot(xcat, p2, 'm-', linewidth = 3, alpha = 0.7, label = 'Model/<S, N>')

 

plt.xticks(xcat, xlabels, fontsize = 14)

plt.xlabel('Rating Category', fontsize = 14)

plt.ylabel('Probability/Proportion', fontsize = 14)

plt.title('Ratings in Conditions <N, S> and <S, N>', fontsize = 18)

plt.legend(loc = 'upper center', fontsize = 10)

plt.show()

 

fout.close()

print('\n{} was saved.\n'.format(fout_nm))

 

Files in Listing C.1 and Listing C.2 and sample data files are archived into a zip file RatingEvenCat.zip, which can be downloaded.

Run SDT2AFCEvenCat.py in Listing C.2, names of the input data file and an output text file are required to be set (Figure C.1).

Figure C.1

 

Output text file name is arbitrary. The input datafile should be CSV format (Figure C.2).

Figure C.2  Data4Cat.csv

 

In the first row, names for columns are set. The first column shows conditions, <N, S> and <S, N>, the second column shows frequencies of rating judgment of category 1 (e.g., The signal is on the left), the third column shows frequencies of rating judgment of category 2 (e.g., Probably, the signal is on the left), the fourth column shows frequencies of rating judgment of category 3 (e.g., Probably, the signal is on the right), and the fifth column shows frequencies of rating judgment of category 4 (e.g., The signal is on the right).

In the second row, data from condition <N, S> are set, and in the third row, data from condition <S, N>.

When the file names are set, the program starts calculation. After sampling by Stan ends, posterior distribution of parameter d is presented (Figure C.3).

Figure C.3

 

Closing the form in Figure C.3, posterior distribution of τ is shown (Figure C.4).

Figure C.4

 

Closing the form in Figure C.4, posterior distribution of category boundaries are shown (Figure C.5).

Figure C.5

 

Closing the form in Figure C,5, graphs of data and model predictions are shown (Figure C.6).

Figure C.6

 

Dotted lines represent proportions of rating categories calculated from the data, and solid lines represent estimated probabilities of rating categories, which are medians of posterior distributions. Rating categories 1 to 4 correspond to the columns 2 to 5 in Figure C.2, respectively.

Closing the form of Figure C.6, the program ends. After the program ends, output text file can be opened. The content is like this:

 

Input Data File = Data4Cat.csv

 

<N, S>

[24, 9, 12, 55]

 

<S, N>

[74, 13, 7, 6]

 

Inference for Stan model: anon_model_f2c4de04275e32f42265169574310b77.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

mu_s       1.05  2.2e-3   0.13   0.81   0.96   1.05   1.14    1.3   3432    1.0

tau        0.43  2.0e-3   0.13   0.18   0.35   0.43   0.52   0.68   4000    1.0

theta[0]   0.51  1.3e-3   0.07   0.38   0.46   0.51   0.55   0.65   2759    1.0

C[0]      -0.51  1.3e-3   0.07  -0.65  -0.55  -0.51  -0.46  -0.38   2759    1.0

C[1]        0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0   4000    nan

C[2]       0.51  1.3e-3   0.07   0.38   0.46   0.51   0.55   0.65   2759    1.0

p_NS[0]    0.22  6.1e-4   0.04   0.15   0.19   0.21   0.24    0.3   4000    1.0

p_NS[1]    0.12  2.5e-4   0.02   0.09   0.11   0.12   0.13   0.15   4000    1.0

p_NS[2]    0.14  3.6e-4   0.02    0.1   0.12   0.14   0.15   0.18   2950    1.0

p_NS[3]    0.53  7.8e-4   0.05   0.44    0.5   0.53   0.56   0.63   4000    1.0

p_SN[0]    0.08  4.4e-4   0.02   0.04   0.07   0.08    0.1   0.13   2688    1.0

p_SN[1]    0.07  1.9e-4   0.01   0.05   0.06   0.07   0.07   0.09   4000    1.0

p_SN[2]     0.1  2.7e-4   0.02   0.07   0.09    0.1   0.11   0.13   4000    1.0

p_SN[3]    0.75  6.8e-4   0.04   0.67   0.73   0.75   0.78   0.83   3710    1.0

lp__     -210.0    0.03   1.26 -213.3 -210.5 -209.6 -209.0 -208.6   2124    1.0

 

Samples were drawn using NUTS at Fri Oct 25 10:46:34 2019.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

d':

Median = 1.052,   90%CI = [0.844, 1.257]

 

tau:

Median = 0.433,   90%CI = [0.225, 0.643]

 

C[1]:

Median = -0.506,   90%CI = [-0.624, -0.400]

 

C[2]:

Median = 0.000,   90%CI = [0.000, 0.000]

 

C[3]:

Median = 0.506,   90%CI = [0.400, 0.624]

 

 

Next, consider another example, which uses K=6 number of rating categories (Figure C.7).

Figure C.7  Data6Cat.csv

 

Run SDT2AFCEvenCat.py in Listing C.2, and set names of the input data file and an output text file as in Figure C.8.

Figure C.8

 

When the file names are set, the program starts calculation. After sampling by Stan ends, posterior distribution of parameter d is presented (Figure C.9).

Figure C.9

 

Closing the form in Figure C.9, posterior distribution of τ is shown (Figure C.10).

Figure C.10

 

Closing the form in Figure C.10, posterior distribution of category boundaries are shown (Figure C.11).

Figure C.11

 

Closing the form in Figure C,11, graphs of data and model predictions are shown (Figure C.12).

Figure C.12

 

Dotted lines represent proportions of rating categories calculated from the data, and solid lines represent estimated probabilities of rating categories, which are medians of posterior distributions. Rating categories 1 to 6 correspond to the columns 2 to 7 in Figure C.7, respectively.

Closing the form of Figure C.12, the program ends. After the program ends, output text file can be opened. The content is like this:

 

Input Data File = Data6Cat.csv

 

<N, S>

[6, 23, 5, 4, 37, 25]

 

<S, N>

[42, 38, 2, 5, 11, 2]

 

Inference for Stan model: anon_model_f2c4de04275e32f42265169574310b77.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

mu_s       0.99  1.9e-3   0.12   0.76   0.91   0.99   1.06   1.22   3781    1.0

tau        0.36  1.7e-3   0.11   0.15   0.29   0.36   0.43   0.57   4000    1.0

theta[0]   0.19  7.1e-4   0.04   0.11   0.16   0.19   0.22   0.29   4000    1.0

theta[1]   1.42  1.8e-3   0.11   1.21   1.34   1.42   1.49   1.65   4000    1.0

C[0]      -1.61  1.8e-3   0.11  -1.84  -1.69  -1.61  -1.54  -1.39   4000    1.0

C[1]      -0.19  7.1e-4   0.04  -0.29  -0.22  -0.19  -0.16  -0.11   4000    1.0

C[2]        0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0   4000    nan

C[3]       0.19  7.1e-4   0.04   0.11   0.16   0.19   0.22   0.29   4000    1.0

C[4]       1.61  1.8e-3   0.11   1.39   1.54   1.61   1.69   1.84   4000    1.0

p_NS[0]    0.06  2.9e-4   0.02   0.03   0.05   0.06   0.07    0.1   3455    1.0

p_NS[1]    0.22  4.1e-4   0.03   0.17   0.21   0.22   0.24   0.28   4000    1.0

p_NS[2]    0.05  1.7e-4   0.01   0.03   0.04   0.05   0.05   0.07   4000    1.0

p_NS[3]    0.05  1.9e-4   0.01   0.03   0.04   0.05   0.06   0.08   4000    1.0

p_NS[4]    0.38  4.3e-4   0.03   0.32   0.36   0.38   0.39   0.43   3947    1.0

p_NS[5]    0.24  6.0e-4   0.04   0.17   0.22   0.24   0.27   0.32   4000    1.0

p_SN[0]    0.02  1.4e-4 8.0e-3 7.9e-3   0.01   0.02   0.02   0.04   3464    1.0

p_SN[1]    0.12  3.4e-4   0.02   0.08   0.11   0.12   0.13   0.17   3899    1.0

p_SN[2]    0.03  1.2e-4 7.6e-3   0.02   0.03   0.03   0.04   0.05   4000    1.0

p_SN[3]    0.04  1.5e-4 9.5e-3   0.02   0.03   0.04   0.04   0.06   4000    1.0

p_SN[4]    0.37  4.5e-4   0.03   0.31   0.35   0.36   0.38   0.42   4000    1.0

p_SN[5]    0.43  7.1e-4   0.04   0.34    0.4   0.43   0.46   0.51   4000    1.0

lp__     -298.1    0.03   1.42 -301.7 -298.8 -297.8 -297.1 -296.3   1956    1.0

 

Samples were drawn using NUTS at Fri Oct 25 11:12:13 2019.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

d':

Median = 0.987,   90%CI = [0.796, 1.178]

 

tau:

Median = 0.361,   90%CI = [0.182, 0.538]

 

C[1]:

Median = -1.611,   90%CI = [-1.802, -1.427]

 

C[2]:

Median = -0.189,   90%CI = [-0.272, -0.124]

 

C[3]:

Median = 0.000,   90%CI = [0.000, 0.000]

 

C[4]:

Median = 0.189,   90%CI = [0.124, 0.272]

 

C[5]:

Median = 1.611,   90%CI = [1.427, 1.802]

 

 

References

Klein, S. A. (2001). Measuring, estimating, and understanding the psychometric function: A commentary. Perception & Psychophysics, 63, 1421-1455.

 

 

Up