Triangular Test
A logistic regression model
Yasuharu Okamoto, 2021.12
In triangular test, three stimuli are presented, where two of them are the same. The same two stimuli are denoted by A and A’, and the other by B. The observer is required to choose the one, which the observer judged as the B. This is the standard method of the triangular test, and discussed in the first half of the website.
We also consider the following situation:
For the same two stimuli, A and A’, the third stimuli to be discriminated/detected, three stimuli, B, C, and BC, are prepared. The stimulus BC has the properties of B and C. The triangular tests consist of the three sets of stimuli, (A, A’, B), (A, A’, C), and (A, A’, BC). We analyze the interaction of B and C in BC by a logistic regression model, which is presented in the second half of the website.
The programs in the website are archived in the file TriTestLogisticFiles.zip, which can be downloaded and used freely on the user’s responsibility. All rights are reserved.
First, consider the standard case of the triangular test.
Denote the probability of detection by 
. If the observer detects the stimulus B with probability Pd, the
observer responds correctly.
If the observer cannot detect the stimulus B, he/she chooses randomly among the three. In this case, the probability of correct response is 1/3.
Hence, the probability of correct response Pc is given as follows:
When the total number of trials of the
triangular test is N, using the binomial distribution, we have the probability 
that the number of
correct responses is k as follows:
Stan script based on the equation (2) is shown in Listing 1.
Python script, which uses the Stan script in Listing 1, is shown in Listing 2.
Run the script in Listing 2, values of k and N are required to be set as follows:
(py39) PS
D:\xxxxx\TriTestLogisticFiles> python .\TriTest.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
k = 50
N = 100
In the above example, 50 and 100 are set for k and N.
After the value of N is set, calculation starts.
After the calculation ends, the posterior distribution of Pd is shown as in Figure 1.
Figure 1
In Figure 1, it is shown that the mode is 0.24, and the 95% HDI is [0.11, 0.40].
The graph is saved automatically with the name FigPd.png.
The logistic transformation of probability p is given by
Hence, we have
Logistic regression model is constructed by this transformation (Gelman et al., 2014, sec. 3.7).
Let the detection/discrimination probabilities PdB, PdC, and PdBC for stimuli B, C, and BC, respectively, be given as follows
fB, fC, and fBC are factors, which represent the main effects of stimuli B and C, and the interaction by B and C.
Stan script based on the above model is shown in Listing 3.
The Python script in Listing 4 uses the Stan script in Listing 3.
Run the script in Linsting 4, values for k and N for the three sets of triangular tests, (A, A’, B), (A, A’, C), and (A, A’, BC), are required to be set as follows.
(py39) PS
D:\xxxxx\TriTestLogisticFiles> python .\TriTestLogistic.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
kB = 50
NB = 100
kC = 34
NC = 100
kBC = 80
NBC = 100
In the above example, ks and Ns are set to be 50 and 100 for (A, A’, B), 34 and 100 for (A, A’, C), and 80 and 100 for (A, A’, BC).
After the input of NBC = 100, calculation starts, and the posterior distributions of factors fB, fC, and fBC are drawn as in Figure 2.
Figure 2
In the case where values for effects fB or fC be negative with large absolute values, probabilities of PdB or PdC be small ones (eqs.(3) and (4)). The interaction fBC is 0 if the model (5) is additive one with no interaction between B and C.
In Figure 2, the interaction fBC is clearly positive, which means that the effect of combination of B and C is strong.
When the form in Figure 2 is closed, calculation for the posterior distributions of PdB, PdC, and PdBC starts.
The posterior distributions of PdB, PdC, and PdBC are displayed as in Figure 3.
Figure 3
Detection probability PdC is negligible, but when C is combined with B, C boosts the detection probability PdBC far above that of PdB.
Graphs in Figures 2 and 3 are automatically saved in files FigFctrs.png and FigPd.png.
As another example, the following values are set to ks and Ns.
(py39) PS
D:\xxxxx\TriTestLogisticFiles> python .\TriTestLogistic.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
kB = 60
NB = 100
kC = 60
NC = 100
kBC = 40
NBC = 100
In this case, the results are as follows.
Figure 4 shows the posterior distributions of fB, fC, and fBC.
Figure 4
The interaction fBC is distributed in the negative region.
The posterior distributions of PdB, PdC, and PdBC are shown in Figure 5.
Figure 5
Detection probabilities PdB and PdC are about 0.4, but that of PdBC is negligible, which goes along with the negative value of fBC.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, third edition. CRC Press.
Listing 1. Stan script for the model (2). (file name: tri_test.stan)
data {
int N;
int k;
}
parameters {
real<lower = 0.0, upper =
1.0> Pd; // Discrimination Probability
}
transformed parameters {
real Pc;
// Correct Response
Probability
Pc = Pd + (1.0 - Pd) *
(1/3.0);
}
model {
Pd ~ uniform(0.0, 1.0);
k ~ binomial(N, Pc);
}
Listing 2. Python script, which uses the Stan script in listig 1. (file name: TriTest.py)
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
import pickle
def CalcMAPEst(samples, a = 0.05,
n_points = 10000):
"""
Calculation of the 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
"""
print('Calculating the
mode.')
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]
def calcHDI(sample, a = 0.05):
"""
Calculation of 100(1-a)% HDI for sample
Chen-Shao HPD Estimation algorithm
"""
print('Calculating HDI.')
smpl = np.sort(sample)
n = len(smpl)
ckW = max(smpl) - min(smpl)
cki = 0
RB = int(n * (1-a))
L0 = 0
U0 = n-1
while (RB + cki < n):
if
smpl[RB + cki] - smpl[cki] < ckW:
ckW = smpl[RB + cki] - smpl[cki]
L0 = cki
U0 = RB + cki
cki
+= 1
ckW = max(smpl) - min(smpl)
cki = 1
LB = n - int(n * (1-a))
L1 = 0
U1 = n-1
while (LB - cki >= 0):
if
smpl[n - cki] - smpl[LB - cki] < ckW:
ckW = smpl[n - cki] - smpl[LB - cki]
U1 = n - cki
L1 = LB - cki
cki
+= 1
return smpl[int((L1+L0)/2)],
smpl[int((U1+U0)/2)]
k = int(input('k = '))
N = int(input('N = '))
Data = {'N': N, 'k': k}
sm = pystan.StanModel(file =
'tri_test.stan')
fit = sm.sampling(data = Data, iter =
10_000, n_jobs = 1)
print(fit)
sb.kdeplot(fit['Pd'])
plt.hist(fit['Pd'], density = True,
color = 'g', alpha = 0.5)
MAPEst = CalcMAPEst(fit['Pd'])[0]
plt.scatter([MAPEst], [0], marker =
'o', s = 200,
c = 'magenta', label = f'mode = {MAPEst:.2f}')
L_bndry, R_bndry = calcHDI(fit['Pd'],
0.05)
plt.plot([L_bndry, R_bndry], [0, 0], c
= 'r', lw = 5,
label = '95% HDI [{0:.2f}, {1:.2f}]'.format(L_bndry, R_bndry))
plt.yticks([])
plt.xlabel('Pd', fontsize = 12)
plt.title(f'Posterior distribution\nk
= {k}, N = {N}', fontsize = 16)
plt.legend()
plt.savefig('FigPd.png')
plt.show()
Listing 3. Stan script of the logistic regression model (file name: tri_test_logistic.stan)
//
// Yasuharu
Okamoto, 2021.12
//
data {
int NB;
int kB;
int NC;
int kC;
int NBC;
int kBC;
}
parameters {
real fB; // factor for B
real fC; // factor for C
real fBC; // interaction factor of B and C
}
transformed parameters {
real<lower = 0.0, upper =
1.0> PdB;
real<lower = 0.0, upper =
1.0> PdC;
real<lower = 0.0, upper =
1.0> PdBC;
real<lower = 0.0, upper =
1.0> PcB;
real<lower = 0.0, upper =
1.0> PcC;
real<lower = 0.0, upper =
1.0> PcBC;
//
// Logistic regression
model
//
PdB = 1.0 / (1.0 +
exp(-fB));
PdC = 1.0 / (1.0 +
exp(-fC));
PdBC = 1.0 / (1.0 + exp(-(fB
+ fC + fBC)));
PcB = PdB + (1.0 - PdB) *
(1.0 / 3.0);
PcC = PdC + (1.0 - PdC) *
(1.0 / 3.0);
PcBC = PdBC + (1.0 - PdBC) *
(1.0 / 3.0);
}
model {
//
// Weakly informative priors
//
fB ~ normal(0.0, 5.0);
fC ~ normal(0.0, 5.0);
fBC ~ normal(0.0, 5.0);
//
// Distributions of the correct
responses
//
kB ~ binomial(NB, PcB);
kC ~ binomial(NC, PcC);
kBC ~ binomial(NBC, PcBC);
}
Listing 4. Python script, which uses the Stan script in Listing 3 (file name: TriTestLogistic.py)
import numpy as np
import scipy.stats as ss
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
def CalcMAPEst(samples, a = 0.05,
n_points = 10000):
"""
Calculation of the 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
"""
print('Calculating the
mode.')
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]
def calcHDI(sample, a = 0.05):
"""
Calculation of 100(1-a)% HDI for sample
Chen-Shao HPD Estimation algorithm
"""
print('Calculating the
HDI.')
smpl = np.sort(sample)
n = len(smpl)
ckW = max(smpl) - min(smpl)
cki = 0
RB = int(n * (1-a))
L0 = 0
U0 = n-1
while (RB + cki < n):
if
smpl[RB + cki] - smpl[cki] < ckW:
ckW = smpl[RB + cki] - smpl[cki]
L0 = cki
U0 = RB + cki
cki
+= 1
ckW = max(smpl) - min(smpl)
cki = 1
LB = n - int(n * (1-a))
L1 = 0
U1 = n-1
while (LB - cki >= 0):
if
smpl[n - cki] - smpl[LB - cki] < ckW:
ckW =
smpl[n - cki] - smpl[LB - cki]
U1 = n - cki
L1 = LB - cki
cki
+= 1
return smpl[int((L1+L0)/2)],
smpl[int((U1+U0)/2)]
kB = int(input('kB = '))
NB = int(input('NB = '))
kC = int(input('kC = '))
NC = int(input('NC = '))
kBC = int(input('kBC = '))
NBC = int(input('NBC = '))
sm = pystan.StanModel(file =
'tri_test_logistic.stan')
Data = {'kB': kB, 'NB': NB, 'kC': kC,
'NC': NC, 'kBC': kBC, 'NBC': NBC}
fit = sm.sampling(data = Data, iter =
10_000, n_jobs = 1)
print(fit)
#
# Posterior
distributions of factors fB, fC, and fBC
#
MAPfB = CalcMAPEst(fit['fB'])[0]
MAPfC = CalcMAPEst(fit['fC'])[0]
MAPfBC = CalcMAPEst(fit['fBC'])[0]
L_HDIfB, R_HDIfB = calcHDI(fit['fB'])
L_HDIfC, R_HDIfC = calcHDI(fit['fC'])
L_HDIfBC, R_HDIfBC =
calcHDI(fit['fBC'])
sb.kdeplot(fit['fB'], label = 'fB,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfB, L_HDIfB, R_HDIfB))
sb.kdeplot(fit['fC'], label = 'fC,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfC, L_HDIfC, R_HDIfC))
sb.kdeplot(fit['fBC'], label = 'fBC,
mode = {0:.2f}, \n
95% HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfBC, L_HDIfBC, R_HDIfBC))
plt.xlabel('Factor', fontsize = 14)
plt.yticks([])
plt.title('Posterior distribtion\n' +
\
('kB = {0:}, NB = {1:}, kC = {2:}, NC = {3:}\n'.format(kB, NB, kC, NC))
+\
('kBC = {0:}, NBC = {1:}'.format(kBC, NBC)), fontsize = 16)
plt.legend()
plt.tight_layout()
plt.savefig('FigFctrs.png')
plt.show()
#
# Posterior
distributions of PdB, PdC, and PdBC
#
MAPPB = CalcMAPEst(fit['PdB'])[0]
MAPPC = CalcMAPEst(fit['PdC'])[0]
MAPPBC = CalcMAPEst(fit['PdBC'])[0]
L_HDIPB, R_HDIPB = calcHDI(fit['PdB'])
L_HDIPC, R_HDIPC = calcHDI(fit['PdC'])
L_HDIPBC, R_HDIPBC = calcHDI(fit['PdBC'])
sb.kdeplot(fit['PdB'], label = 'PdB,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPB, L_HDIPB, R_HDIPB))
sb.kdeplot(fit['PdC'], label = 'PdC,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPC, L_HDIPC, R_HDIPC))
sb.kdeplot(fit['PdBC'], label = 'PdBC,
mode = {0:.2f}, \n
95% HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPBC, L_HDIPBC, R_HDIPBC))
plt.xlabel('Prob. of detection',
fontsize = 14)
plt.yticks([])
plt.title('Posterior distribtion\n' +
\
('kB = {0:}, NB = {1:}, kC = {2:}, NC = {3:}\n'.format(kB, NB, kC, NC))
+\
('kBC = {0:}, NBC = {1:}'.format(kBC, NBC)), fontsize = 16)
plt.legend()
plt.tight_layout()
plt.savefig('FigPd.png')
plt.show()