Bayesian Analysis for Triangular Test
Grid method for d’
In triangular test, three stimuli are presented, two of them, A an A’, are the same one, and the observer is required to choose the third one B, which is difficult to discriminate from A or A’.
The sensations of A, A’ and B are represented by random variables 
, 
, and 
, and assumed to have the following distributions, respectively.
The correct response is made, when the following conditions are satisfied (Frijters et al., 1980, p. 177):
Set d’ as follows
We get the following equation (Frijters et al., 1980, eq.(1)):
where 
denotes the probability
of choosing the stimulus B.
Now, the probability of k correct responses for the total number of trials N is given by the equation (2):
The script based on the above model is shown in Listing 1.
The script file is archived in the fileTriTestScriptFile.zip, which can be downloaded and used freely on the responsibility of the user. All rights are reserved.
Run the script, values for N and k are required to be set as follows:
>>>
= RESTART:
D:\yasuharu\XXXXX\TriangularTest_dprime_Grid\TriTestScriptFile\TriangularTest_dprime_grid.py
N = 50
k = 20
Set the values, then calculation starts.
After the calculation ends, the results are shown as the graph like Figure 1.
Figure 1
In Figure 1, the posterior distribution is shown with the mode and 95% HDI.
The graph is saved automatically by the code:
plt.savefig(f'FigN{N}k{k}.png')
In the case of Figure 1, the file name is FigN50k20.png.
Reference
Frijters, J. E. R., Kooistra, A., & Vereijken, P. F. G. (1980). Tables of d' for the triangular method and the 3-AFC signal detection procedure. Perception & Psychophysics, 1980, 27, 176-178.
Listing
1. Script for the Bayesian analysis
of triangular test
"""
Yasuharu Okamoto, 2021.11, 2022.01
"""
import numpy as np
import scipy.stats as ss
import scipy.integrate as si
import matplotlib.pyplot as plt
import seaborn as sb
N = int(input('N = '))
k = int(input('k = '))
class Triangular:
def __init__(self, d_prime):
self.d_prime = d_prime
def func(self, u):
v =
ss.norm.cdf(-u * (3.0**0.5) + self.d_prime * ((2/3)**0.5)) + \
ss.norm.cdf(-u * (3.0**0.5) - self.d_prime * ((2/3)**0.5))
return v * np.exp(-0.5 * (u**2)) / ((2.0 * np.pi)**0.5)
def Pc_d_prime(d_prime):
d_func =
Triangular(d_prime).func
return 2 * si.quad(d_func,
0.0, +np.inf)[0]
dprimes = np.linspace(0, 10, 1001)
Pc = np.empty(len(dprimes))
for i in range(len(dprimes)):
if i % 100 == 0:
print('{0:}/{1:}
'.format(i, len(dprimes)))
#, end = '\r')
Pc[i] =
Pc_d_prime(dprimes[i])
postPdprime = np.empty(len(dprimes))
for i in range(len(dprimes)):
postPdprime[i] = (Pc[i] **
k) * ((1 - Pc[i]) ** (N - k))
postPdprime = (postPdprime /
np.sum(postPdprime))
map_idx = np.argmax(postPdprime)
MAP_Est = map_idx / 100
print('MAP_Est =', MAP_Est)
postPcum = np.cumsum(postPdprime)
a = 0.05
Lp = 0
Rp = 0
w_t = 1000
L_t = Lp
R_t = Rp
while True:
if Lp == 0:
cumPL = 0.0
else:
cumPL = postPcum[Lp-1]
while 1.0 - postPcum[Rp]
>= a - cumPL:
Rp
+= 1
if Lp == 0:
L_t
= Lp
R_t
= Rp
w_t
= R_t - L_t
elif Rp - Lp < w_t:
L_t
= Lp
R_t
= Rp
w_t
= R_t - L_t
Lp += 1
if postPcum[Lp] >= a:
Lp
-= 1
break
Lp = 1000
Rp = 1000
w_t = 1000
L_t1 = Lp
R_t1 = Rp
while True:
if Rp == 1000:
cumPR = 0.0
else:
cumPR = 1.0 - postPcum[Rp]
while postPcum[Lp] >= a -
cumPR:
if
Lp > 0:
Lp -= 1
else:
break
if Rp == 1000:
L_t1
= Lp
R_t1
= Rp
w_t
= R_t1 - L_t1
elif Rp - Lp < w_t:
L_t1
= Lp
R_t1
= Rp
w_t
= R_t1 - L_t1
Rp -= 1
if 1.0 - postPcum[Rp] >=
a:
Rp
+= 1
break
HDI_L = ((L_t + L_t1) / 2.0) / 100.0
HDI_R = ((R_t + R_t1) / 2.0) / 100.0
plt.plot(np.linspace(0, 10, 1001),
postPdprime)
plt.plot([MAP_Est,MAP_Est], [0,
postPdprime[map_idx]], label = 'Mode = {}'.format(MAP_Est))
plt.title(f"Psterior Distribution
for d'\nN = {N}, k = {k}", fontsize = 14)
plt.xlabel("d'", fontsize =
12)
plt.plot([0, 10], [0, 0], lw = 1, c =
'k')
plt.yticks([0])
plt.plot([HDI_L, HDI_R], [0, 0], lw =
5, label = '95% HDI = [{0:.2f}, {1:.2f}]'.format(HDI_L, HDI_R))
plt.legend()
plt.savefig(f'FigN{N}k{k}.png')
plt.show()