A Generalized Simple Regression Model of d’
Bayesian analysis
Yasuharu Okamoto, 2022.01
For the triangular test, i.e., 3AFC oddity
task, a generalized simple regression model of d’ was
proposed. The three stimuli, A, A’, B, are assumed to
invoke sensations, 
, 
, , and 
, respectively, where A and A’ are the same
physical stimuli, which are different from B. and have normal
distributions of mean 
and varinac 
, and has a normal
distribution of mean 
and variance , that is,
The observer chooses the correct one, if and only if the following conditions are satisfied (Frijters et al., 1980, p. 177):
According to signal detection theory (SDT: Macmillan & Creelman, 2005; Green
& Swets, 1966), define 
as follows:
Then the probability 
of the correct response
is given by the following equation (1) (Frijters et al. 1980,
eq (1)):
When 
is replaced by 
in equation (1), we get
the same value. Hence, we can assume
Consider the case, where the data is given by a pair, person’s response R and attribute X.
R is 1, if the response is correct and 0 otherwise. X is, e.g., the person’s age.
The data might be summarized as like in Figure 1. In the first row, variable names are written. In the first column, data identification value, numerical value or string, are set. In the second column, response Rs are set , 1 if success, 0 otherwise. In the third column, person’s attribute are set, e.g., ages.
The data with CaseID 1 in Figure 1 shows that the response is correct (R = 1), and the age is 33 years old.
.
.
.
Figure 1. An artificial input data file (file name: Data.csv)
To construct a regression model of d’, transform d’ by the following link function (Gelman et al., 2021):
Hence, we have
By equations (1) and (2), the likelihood function can be given as a function of parameters a and b, that is,
Script in Listing 1 calculates the log of the posterior distribution of a and b as follows:
v = my_log(phi(a, 0.0, 10.0)) + my_log(phi(b,
0.0, 10.0))
for
i in range(len(self.Rs)):
dprime = np.exp(a + b * self.Ags[i])
if dprime > 10.0:
dprime = 10.0
p = self.TblPDprime[f'{dprime:.2f}']
v += self.Rs[i] * my_log(p) +\
(1 - self.Rs[i]) * my_log(1- p)
Rs represents R, and Ags represents X.
Independent variable Ags is automatically standardized by the program. The prior distributions of a and b are set to be normal distributions with mean 0 and standard deviation 10 as weak informative ones.
Calculation of integral (1) takes rather long time, so values of d’ and corresponding correct probability Pc are calculated before the MCMC sampling and stored in a dictionary by the following code in Listing 2.
TblPDprime = mm.calcPDprimes()
MCMC sampling is executed by Metropolis-within-Gibbs (Robert & Casella, 2010) as shown in Listing 1.
mymodule.py in Listing 1 contains code for MCMC sampling and so on. The main script, which uses the module mymodule, is shown in Listing 2.
The files of Listing 1 and Listing 2 and so on are archived in TriTestDRgrssnFiles.zip, which can be freely downloaded and used on the user’s responsibility. All rights are reserved.
Run the script main_tritestdrgrssn.py in Listing 2, the name of the input data file is asked.
(py39) PS
D:\xxxxx\TriTestDRgrssnFiles> python .\main_tritestdrgrssn.py
Input file (*.csv) = Data.csv
Enter the file name of the input data (in the above example, Data.csv), then calculation begins.
Before the main MCMC sampling, adjustment of the proposal distributions is done. The main MCMC sampling is composed of four Markov chains, The four chains are executed simultaneously by multiprocessing.
The script for the four Markov chains is as follows (Listing 2):
#
#
Parameter values for multiproccing, each of which runs one chain in MCMC
#
params_t = {'Rs':
deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),
'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,
'n_samples': n_samples, 'chain_id': None}
#
# List of parameter values for
multiprocessing
#
params = []
for i in range(4):
p_t
= deepcopy(params_t)
p_t['chain_id'] = i
params.append(p_t)
print('MCMC sampling by
multiprocessing started.')
with Pool() as pool:
Results = pool.map(mm.MCMCs, params) # Run the multiprocessing
MCMC sampling is executed by four sanpling chains using independent random generators, which is set up by the following code (Listing 1):
rn =
RNp2to191()
#
# Get independent chains
# cf.
http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/
#
for
i in range(self.chain_id):
rn.jump(100)
For more information about the random generator, see this website.
When the sampling ends, sampling distributions and traces of parameter a from the four chains are displayed (Figure 2).
Figure 2
Above the graph of traces, 
(Rhat) (Gelman
et al., 2014, pp.284-285) is shown.
Close the window in which the graphs in Figure 2 are displayed, then graphs of sampling distributions and traces of parameter b for four chains are displayed (Figure 3).
Figure 3
Close the window of Figure 3, then the posterior distribution of parameter a, gathered from the four chains, is displayed with the value of the mode (Figure 4).
Figure 4
Close the window of Figure 4, then the posterior distribution of parameter b, gathered from the four chains, is displayed with the value of the mode (Figure 5).
Figure 5
Close the window of Figure 5, the graph of the regression curve and the data points are displayed (Figure 6).
The regression curve is drawn with the parameter values a and b set to be the modes.
The data points are proportions of correct responses for the four categories of the ages, i.e., independent variables.
It seems that the regression curve is well fitted to the data points.
Figure 6
Close the window of Figure 6, the regression curve of d’, which is based on the modes of the parameter distributions, is shown (Figure 7).
Figure 7
Close the window of Figure 7, then the program ends as follows:
regrslts.txt was saved.
(py39) PS D:\xxxxx\TriTstDRgrssn\TriTestDRgrssnFiles>
The file regrslts.txt contains the following test.
Regression model of Pc on the
standardized independent variable Z:
P(R = 1) = Pc(exp(0.70 + -0.36 * Z))
Regression model of Pc on the raw
independent variable X:
P(R = 1) = Pc(exp(1.96 + -0.03 * X))
Regression model of d' on the
standardized independent variable Z:
d' = exp(0.70 + -0.36 * Z)
Regression model of d' on the raw
independent variable X:
d' = exp(1.96 + -0.03 * X)
The Figures of the graphs are also saved in files in the same folder of the script files automatically.
Frijters, J. E. R., Kooistra, A., & Vereijken, P. F. G. (1980). Tables od d' for the triangular method and the 3-AFC signal detection procedure. Perception & Psychophysics, 1980, 27, 176-178.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., & Rubin, D.B. (2014). Bayesian data analysis, third edition. CRC Press.
Green, D. M., & Swets, J. A. (1966). Signal detection theory and psychophysics. John Wiley & Sons, Inc.
Macmillan, N. A., & Creelman, C. D. (2005). Detection theory: A user’s guide, second edition. Lawrence Erlbaum Associates, Publishers.
Robert, C. P., & Casella, G. (2010). Monte Carlo statistical methods. Springer
Listing 1. Module for MCMC sampling and so on. (file name: mymodule.py; All the files are in , which can be freely downloaded and used under user’s responsibility).
"""
Yasuharu Okamoto, 2021.12
"""
import numpy as np
import scipy.stats as ss
import scipy.integrate as si
from RNGen import *
# Random number
generator
CumPhi = ss.norm.cdf
phi = ss.norm.pdf
def f_phi(z):
return phi(z)
class k_Pc:
def __init__(self, d):
self.d = d
def f(self,u):
"""
Frijters, Kooistra & Vereijken.(1980).
Perception & Psychophysics, 1980, 27, 176-178.
"""
v =
(CumPhi(-u*(3**0.5) + self.d * ((2/3)**0.5))
+ CumPhi(-u*(3**0.5) - self.d * ((2/3)**0.5))) \
* phi(u)
return v
def f_Pc(d):
int_k_Pc = k_Pc(d)
v = si.quad(int_k_Pc.f, 0.0,
+np.inf)
return v[0] * 2.0
def calcPDprimes():
"""
Dictionary: key = d'
value = probability of the correct response
"""
pdprimes = {}
dprimes = np.linspace(0.0,
10.0, 1001)
for i, d in
enumerate(dprimes):
if i
% 100 == 0:
print(i, end = '\r')
p =
f_Pc(d)
pdprimes[f'{d:.2f}'] = p
return pdprimes
def scale_sgm(acpt_r, ck):
"""
Check the size of acceptance proportion
"""
v = 1.0
if (acpt_r > 0.5):
ck
+= 1
v =
1.0 + 9.0 * (acpt_r - 0.5) / 0.5
elif (acpt_r < 0.3):
ck
+= 1
v =
1.0 / (1.0 + 9.0 * (0.3 - acpt_r) / 0.3)
return v, ck
def my_log(v):
if v <= 0.0:
return -750.0
else:
return math.log(v)
class MCMC:
"""
Metropolis-within-Gibbs algorithm
cf. Robert,C.P. & Casella, G. (2010)
Monte Carlo Statistical Methods. p.393
岡本安晴 (2014)．心理学データ分析と測定．pp.229-232.
"""
def __init__(self, Rs, Ags,
TblPDprime, init_a, init_b,
sgm_a, sgm_b, n_samples, chain_id):
self.Rs = Rs
self.Ags = Ags
self.TblPDprime = TblPDprime
self.n_samples = n_samples
self.chain_id = chain_id
self.a_coeff = np.empty(self.n_samples + 1)
self.b_coeff = np.empty(self.n_samples + 1)
self.a_coeff[0] =
init_a
self.b_coeff[0] = init_b
self.sgm_a = sgm_a
self.sgm_b = sgm_b
self.naccpt_a = 0
self.naccpt_b = 0
def LL(self, a, b):
"""
Log likelihood of the parameters a and b
"""
v =
my_log(phi(a, 0.0, 10.0)) + my_log(phi(b, 0.0, 10.0))
for
i in range(len(self.Rs)):
dprime = np.exp(a + b * self.Ags[i])
if dprime > 10.0:
dprime = 10.0
p = self.TblPDprime[f'{dprime:.2f}']
v += self.Rs[i] * my_log(p) +\
(1 - self.Rs[i]) * my_log(1- p)
return v
def do_mcmc(self):
"""
Doing MCMC sampling
"""
import tkinter as tk
# For message window
if
self.chain_id < 0:
msg0 = 'Adjusting sgm...step-{}'.format(-self.chain_id)
else:
msg0 = 'Proc-{}'.format(self.chain_id)
root
= tk.Tk()
root.title(msg0)
cnvs
= tk.Canvas(width = 700, height = 300, background = '#00ffff')
cnvs.pack()
cnvs.create_text(350, 150, text = '', font = ('', 90), fill = 'blue',
tags = 'MyText')
rn = RNp2to191()
#
# Get independent chains
# cf.
http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/
#
for
i in range(self.chain_id):
rn.jump(100)
prev_LL = self.LL(self.a_coeff[0], self.b_coeff[0]) # The initial values
den
= prev_LL
# The initial value
for
t in range(self.n_samples):
if (t % 100) == 0:
cnvs.itemconfig('MyText', text = f'{t}/{self.n_samples}')
cnvs.update()
#
# The next value for
the parameter a
#
y = rn.normalMS(self.a_coeff[t], self.sgm_a)
num = self.LL(y, self.b_coeff[t])
den = prev_LL
accpt = math.exp(num - den)
if rn.uni() < accpt:
self.a_coeff[t + 1] = y
self.naccpt_a += 1
prev_LL = num
else:
self.a_coeff[t + 1] = self.a_coeff[t]
prev_LL = den
#
# The next value for
the parameter b
#
y = rn.normalMS(self.b_coeff[t], self.sgm_b)
num = self.LL(self.a_coeff[t+1], y)
den = prev_LL
accpt = math.exp(num - den)
if rn.uni() < accpt:
self.b_coeff[t + 1] = y
self.naccpt_b += 1
prev_LL = num
else:
self.b_coeff[t + 1] = self.b_coeff[t]
prev_LL = den
root.destroy()
return self.naccpt_a, self.naccpt_b, self.n_samples, \
self.a_coeff[1:], self.b_coeff[1:]
def MCMCs(params):
"""
The
function for muliprocessing
"""
my_mcmc = MCMC(params['Rs'],
params['Ags'], params['TblPDprime'],
params['init_a'],
params['init_b'], params['sgm_a'], params['sgm_b'],
params['n_samples'], params['chain_id'])
r = my_mcmc.do_mcmc()
return r
def calcRhat(psi):
"""
Gelman,A., Carlin,J.B.,Stern,H.S., Dunson,D.B.,
Vehtari,A, & Rubin,D.B. (2014)
Bayesian Data Analysis, third edition. pp.284-285.
"""
psi = np.array(psi)
psidj = psi.mean(axis = 1)
psidd = psi.mean()
n = len(psi[0])
m = len(psi)
B = (((psidj - psidd) **
2).sum()) * n / (m - 1)
sj2 = []
for j in range(m):
sj2.append(((psi[j] - psidj[j]) ** 2).sum() / (n - 1))
sj2 = np.array(sj2)
W = (sj2.sum()) / m
var_hat = W * (n - 1) / n +
B / n
Rhat = (var_hat / W) ** 0.5
return Rhat
def EstMode(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
"""
import numpy as np
import scipy.stats as ss
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
mode_idx =
np.argmax(est_pdf)
Mode_Est =
coord[mode_idx]
return Mode_Est
Listing 2. A program of analyzing for triangular test data, which uses the script in Listing 1 (file name: main_tritestdrgrssn.py).
"""
Yasuharu Okamoto, 2021.12, 2022.01
"""
import csv
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
from multiprocessing.pool import Pool
import mymodule as mm # a personal module, file name:
mymodule.py
if __name__ == '__main__':
fin_nm = input('Input file
(*.csv) = ')
with open(fin_nm, 'r') as f:
Data_in =
[v for v in csv.reader(f)]
Data_in = np.array(Data_in)
Rs = [int(v) for v in
Data_in.T[1][1:]]
Ags = [float(v) for v in
Data_in.T[2][1:]]
print(len(Rs), len(Ags))
Rs = np.array(Rs)
Ags = np.array(Ags)
Ags_raw = deepcopy(Ags)
Ags_mean = Ags.mean()
Ags_std = Ags.std()
Ags = (Ags - Ags_mean) /
Ags_std
# Construction of a dictionary for
d's and probabilities of correct response
print("Preparing the
dictionary of d's and probabilties")
TblPDprime =
mm.calcPDprimes()
#
# Adjusting the
proposal distributions
#
i_step = -1
sgm_a = 0.5
sgm_b = 0.5
init_a = 0.0
init_b = 0.0
print('Adjusting the
proposals')
while True:
print('step-{0:} started.'.format(-i_step))
my_mcmc = mm.MCMC(deepcopy(Rs), deepcopy(Ags), deepcopy(TblPDprime),
init_a, init_b, sgm_a, sgm_b, 2000, i_step)
n_accpt_a, n_accpt_b, L_chain, chain_a, chain_b = my_mcmc.do_mcmc()
acpt_r_a = n_accpt_a / L_chain
acpt_r_b = n_accpt_b / L_chain
ck_a
= 0
v,
ck_a = mm.scale_sgm(acpt_r_a, ck_a)
# Check the acceptance
proportions
init_a = chain_a.mean()
if
ck_a > 0:
sgm_a *= v
# Adjust the sigma of
the normal distribution for proposal
ck_b
= 0
v,
ck_b = mm.scale_sgm(acpt_r_b, ck_b)
init_b = chain_b.mean()
if
ck_b > 0:
sgm_b *= v
if
ck_a + ck_b == 0: # The sizes of the sigmas OK?
break
# Adjustment ends,
i_step -= 1
n_samples = 5000 # Sampling size for each chain
in MCMC sampling
#
#
Parameter values for multiproccing, each of which runs one chain in MCMC
#
params_t = {'Rs':
deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),
'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,
'n_samples': n_samples, 'chain_id': None}
#
# List of parameter values for
multiprocessing
#
params = []
for i in range(4):
p_t
= deepcopy(params_t)
p_t['chain_id'] = i
params.append(p_t)
print('MCMC sampling by
multiprocessing started.')
with Pool() as pool:
Results = pool.map(mm.MCMCs, params) # Run the multiprocessing
print('MCMC sampling
ended.')
fout_nm = 'regrslts.txt'
fout = open(fout_nm,
'w') # Output text file
#
# Posterior
distributon and trace of the parameter a
#
plt.figure(figsize = (12,5))
plt.subplot(1,2,1)
param_a = [[]]*4
for i in range(4):
param_a[i] = Results[i][3]
for i in range(4):
sb.kdeplot(param_a[i]) # KDE posteriors for the parameter a
plt.xlabel('a', fontsize =
16)
plt.yticks([])
plt.title('Posterior
distributions of a', fontsize = 14)
Rhat_a =
mm.calcRhat(param_a)
plt.subplot(1,2,2)
for i in range(4):
plt.plot(range(len(param_a[0])), param_a[i]) # Traces of sampling of the
parameter a
plt.xlabel('trial', fontsize
= 14)
plt.title(f'Traces of
a, Rhat = {Rhat_a:.2f}', fontsize =
14)
plt.savefig('FigDistTrace_a.png')
plt.show()
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
param_b = [[]] * 4
for i in range(4):
param_b[i] = Results[i][4]
for i in range(4):
sb.kdeplot(param_b[i])
plt.xlabel('b', fontsize =
16)
plt.yticks([])
plt.title('Posterior
distributions of b', fontsize = 14)
Rhat_b =
mm.calcRhat(param_b)
plt.subplot(1,2,2)
for i in range(4):
plt.plot(range(len(param_b[i])), param_b[i])
plt.xlabel('trial', fontsize
= 14)
plt.title(f'Traces of
b, Rhat = {Rhat_b:.2f}', fontsize =
14)
plt.savefig('FigDistTrace_b.png')
plt.show()
a_whole = []
b_whole = []
for i in range(4):
a_whole += [v for
v in param_a[i]]
b_whole += [v for v in param_b[i]]
a_mode = mm.EstMode(a_whole)
b_mode = mm.EstMode(b_whole)
print(a_mode, b_mode)
# Posterior
distribution of a
sb.kdeplot(a_whole)
plt.title(f'Posterior
distribution of a, mode =
{a_mode:.2f}', fontsize = 16)
plt.xlabel('a', fontsize =
16)
plt.yticks([])
plt.savefig('Fig_a.png')
plt.show()
# Posterior
distribution of b
sb.kdeplot(b_whole)
plt.title(f'Posterior
distribution of b, mode =
{b_mode:.2f}', fontsize = 16)
plt.xlabel('b', fontsize =
16)
plt.yticks([])
plt.savefig('Fig_b.png')
plt.show()
ags_min = np.min(Ags)
ags_max = np.max(Ags)
def catAgs(v):
"""
Categorization into 4 categories
"""
tv =
int((v - ags_min) * 4.0 / (ags_max - ags_min))
if
tv > 3:
tv = 3
return tv
x_coord_z =
np.linspace(ags_min, ags_max, 101)
# Grid for x-axis
y_coord_z =
np.empty(len(x_coord_z))
for i, v in
enumerate(x_coord_z):
d =
np.exp(a_mode + b_mode * v)
# d' predicted by the
regression model
y_coord_z[i] = mm.f_Pc(d)
# Predicted probability
plt.plot(x_coord_z,
y_coord_z, label = 'Model')
ags_count = np.zeros(4)
r_counts = np.zeros(4) # Number of responses in each of 4
categories
rc_counts = np.zeros(4) # Number of crrect responses in each
of 4 categories
for rv, av in zip(Rs, Ags):
ic =
catAgs(av)
r_counts[ic] += 1
if
rv == 1:
rc_counts[ic] += 1
ave_rc = rc_counts /
r_counts
# proportions of crrect
responses in each category
ave_ags = np.empty(4)
# means of the independent variable in
each category
step_ags = (ags_max -
ags_min) / 8.0 # half of the category width
for i in range(4):
ave_ags[i] = ags_min + step_ags + i * step_ags * 2 # mid points of the categories
plt.plot(ave_ags, ave_rc,
'o', label = 'data')
#
# Conversion of
standardized scores back to raw scores
#
x_ags = [(ags_min +
i*step_ags*2) for i in range(5)]
# positions of the raw
score labels
x_ags = np.array(x_ags)
x_ags_raw = x_ags * Ags_std
+ Ags_mean
x_ags_lbl =
[f"{int(v)}" for v in x_ags_raw]
# labels for raw
scores
plt.xticks(x_ags, x_ags_lbl)
plt.xlabel('Independent
variable', fontsize = 14)
plt.ylabel('Probability',
fontsize = 14)
plt.title('Regression Curve
of Probability', fontsize = 16)
plt.legend()
plt.savefig('Fig_regression_p.png')
plt.show()
fout.write('\nRegression
model of Pc on the standardized independent variable Z:\n')
fout.write('P(R = 1) =
Pc(exp({0:.2f} + {1:.2f} * Z))\n'.format(a_mode, b_mode))
a_raw = a_mode - b_mode *
Ags_mean / Ags_std
b_raw = b_mode / Ags_std
fout.write('\nRegression
model of Pc on the raw independent variable X:\n')
fout.write('P(R = 1) =
Pc(exp({0:.2f} + {1:.2f} * X))\n'.format(a_raw, b_raw))
x_coord_z =
np.linspace(ags_min, ags_max, 101)
# Grid for x-axis
y_coord_d =
np.empty(len(x_coord_z))
for i, v in
enumerate(x_coord_z):
y_coord_d[i] = np.exp(a_mode + b_mode * v) # d' predicted by the regression model
plt.plot(x_coord_z,
y_coord_d)
plt.xticks(x_ags, x_ags_lbl)
plt.xlabel('Independent
variable', fontsize = 14)
plt.ylabel("d'",
fontsize = 16)
plt.title("Regression
curve of d'", fontsize = 16)
plt.savefig('Fig_regression_d.png')
plt.show()
fout.write("\nRegression model of d' on the standardized
independent variable Z:\n")
fout.write("d' =
exp({0:.2f} + {1:.2f} * Z)\n".format(a_mode, b_mode))
fout.write("\nRegression model of d' on the raw independent
variable X:\n")
fout.write("d' =
exp({0:.2f} + {1:.2f} * X)\n".format(a_raw, b_raw))
fout.close()
print('\n', fout_nm, 'was
saved.\n')