Up-and-Down Method with Four Category Rating
Yasuharu Okamoto
To estimate point of subjective equality (PSE) and just noticeable difference (JND), up-down-method with four category rating (‘Stronger’, ‘Probably Stronger’, ‘Probably Weaker’, ‘Weaker’, e.g.) was proposed (Okamoto, 2019). This up-and-down method with four category rating produced more accurate estimates than those by the traditional one using two responses (‘Stronger’ and ‘Weaker’).
The proposed procedure is as follows:
To prevent the observer from anticipating the next stimulus, two series of trials, series 1 and series 2, are prepared. Series 1 and series 2 are randomly interleaved.
Denote the response at trial 
of series 
by 
, and put as this
Denote the stimulus presented at trial of series by 
, and choose the stimulus to be presented at trial 
as follows:
If 
,
If 
,
If 
,
It is recommended that start values of
series 1 and 2 should be sufficiently large or small, so that for example, 
and 
.
Step size 
can be set to be large
value at trial 1, e.g. 
, and is decreasing, e.g.,
Size of is kept to be larger
than or equal to some value 
so that too small value
to discriminate is avoided. Size of may be about one or two
JNDs.
In analyses, the data from the two series, 
and 
are merged and denoted
as 
.
The following psychometric model is assumed.
where,
and 
is the cumulative
standard normal distribution function.
Then, we have the following likelihood function
where 
、
.
Stan script for this likelihood function is shown in Listing 1. In the following, responses -, -?, +?, and + are denoted by 1, 2, 3, and 4, respectively.
Listing 1. Stan script for up-and-down method with four category rating.
data {
int N;
//
The number of trials
int Res[N]; // Response: 1 denotes -, 2 denotes -?, 3 denotes +?,
4 denotes +
real X[N]; // Stimuli
real a;
//
Parameter alpha of a gamma distribution for the prior distribution
real b ;
//
Parameter beta of a gamma distribution for the prior distribution
}
parameters {
real<lower = 0.0> sgm;
real mu;
real<lower = 0.0> C;
}
transformed parameters {
vector[4] p[N];
for (i in 1:N){
p[i][1] = 1 - Phi((X[i] - mu + C) / sgm);
p[i][2]
= Phi((X[i] - mu + C) / sgm) - Phi((X[i] - mu) / sgm);
p[i][3] = Phi((X[i] - mu) / sgm) - Phi((X[i] - mu - C) / sgm);
p[i][4] = Phi((X[i] - mu - C) / sgm);
}
}
model {
sgm ~ gamma(a, b);
mu ~ gamma(a, b);
C ~ gamma(a, b);
for (i in 1:N)
Res[i] ~ categorical(p[i]);
}
Python script, which uses the Stan script of Listing 1, is shown in Listing 2 at the end of this website. Files (Python script file, Stan script file, Sample data file) are archived in this file files4Cat2004.zip , which can be used freely.
When the Python script of Listing 2 are run, input data file name is asked (Figure 1).
Figure 1
The data file should be prepared as CSV format, which can be easily prepared if saved by selecting the file name extension as .csv. The data should be arranged as shown in Figure 2.
・
・
・
Figure 2
The first row shows names of variables.
Data values are set from the second row. The first column represents trial
numbers (any string is OK), the second column represents the stimulus values
presented, and the third column represents the responses 
s. Rating responses, ‘Weaker’, ‘Probably Weaker’,
‘Probably Stronger’, and ‘Stronger’ are represented by integers 1, 2,
3, and 4, respectively.
After the input file name is set, output file name is asked (Figure 3). The name of the output file should be text file name (*.txt).
Figure 3
After the output file name is set, the data file is read in and sampling by Stan starts.
When sampling ends, posterior distribution
of parameter mu(
) , its MAP estimate, and 95% CI are displayed (Figure 4). The graph
can be saved by clicking the save button.
Figure 4
When the window in Figure 4 is closed by
clicking the icon X at the upper right corner, the next window where the
posterior distribution of sigma（
） is displayed is presented (Figure 5).
Figure 5
When the window in Figure 5 is closed, the window where posterior distribution of C is displayed is presented (Figure 6).
Figure 6
When the window in Figure 6 is closed, the window where the three estimated psychometric functions corresponding to the four rating categories are drawn is displayed (Figure 7).
Figure 7
The three curves are drawn by the following
equations:
In the graph of Figure 7, rating values s of data points 
are transformed into range
of 0 to 1, i.e.
so that they are displayed with the
psychometric functions in the same graph. Small random values are added to the transformed
values to keep the points from overlapping
When the window in Figure 7 is closed, the window in Figure 8 is displayed.
Figure 8
In Figure 8, proportions of the responses, 
, 
, and 
, are drawn by kernel density estimations as curves of broken lines.
When the window in Figure 8 is closed, the next value of w_k, which is the width of kernel density estimation, is asked (Figure 9). A larger value of w_k produces smoother curves, and a smaller value variable ones.
Figure 9
When w_k is set to be smaller than or
equal to 0, the program ends (Figure 10).
Figure 10
After the program ends, the output file can be opened. The content is as follows.
Data file ...Data4Cat.csv
Number of data = 100
0 257.50 4
1 237.50 4
2 222.50 4
3 207.50 4
4 192.50 3
5 177.50 1
.
.
.
95 192.50 3
96 207.50 4
97 192.50 2
98 177.50 1
99 192.50 1
Summary...
Inference for Stan model:
anon_model_3fb9b3a405c246ad59b9839157950b06.
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
sgm 14.13 0.06 2.33 10.43 12.51 13.84 15.47 19.28 1783 1.0
C 8.14 0.03 1.47 5.61 7.15 7.99 9.01 11.55 2338 1.0
mu 200.46 0.03 1.85 196.73 199.26 200.45 201.64
204.12 2817 1.0
lp__ -88.69 0.03 1.33 -92.03 -89.27 -88.34
-87.73 -87.2 1612 1.0
Samples were drawn using NUTS at Wed
Apr 29 15:15:59 2020.
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).
MAP for mu = 200.55 95% CI for mu = [196.74,
204.11]
Map for sigma = 12.83 95% CI for sigma = [10.43,
19.26]
MAP for C = 7.60 95% CI for C = [5.61,
11.55]
MAP for PSE = 200.55, 95% CI for PSE = [196.74,
204.11]
MAP for JND = 8.65, 95% CI for JND= [7.04,
12.99]
Listing 2 Python script which uses the Stan script in Linsting1. Files of Listings 1, 2, and 3, and the sample data are archived in this file files4Cat2004.zip , which can be freely downloaded and used.
import numpy as np
import matplotlib.pyplot as plt
import pystan
import pickle
import seaborn as sb
import csv
from my_module4Cat2004 import *
fin_nm = input('Data file name (*.csv)
= ')
fout_nm = input('Output file name = ')
fout = open(fout_nm, 'w')
fout.write('Data file
...{}\n'.format(fin_nm))
#
# Read data
from CSV file
#
with open(fin_nm, 'r') as fin:
data = [v for v in csv.reader(fin)]
N = len(data) - 1
ID = []
X = np.empty(N, dtype = float)
Res = np.empty(N, dtype = int)
for t in range(N):
ID.append(data[t + 1][0])
X[t] = float(data[t + 1][1])
Res[t] = int(data[t + 1][2])
print('Number of data = ', N)
fout.write('\nNumber of data =
{}\n\n'.format(N))
for i in range(N):
print('{0} {1}
{2}'.format(ID[i], X[i], Res[i]))
fout.write('{0:>5s} {1:>10.2f} {2:>5d}\n'.format(ID[i], X[i],
Res[i]))
#
# Boundaries
and initial values
#
C1C2 = np.percentile(X, [20, 80])
param_a = 1.0 # Parameter
a of a gamma distribution for the prior distribution
param_b = 0.001 # Parameter
b of a gamma distribution for the prior distribution
def f_init():
"""
Initial values for C, mu, and sigma
"""
return dict(C = (C1C2[1] -
C1C2[0]) / 2.0, mu = (C1C2[0] + C1C2[1]) / 2,
sgm = C1C2[1] - C1C2[0])
stan_data = {'N': N, 'X':X, 'Res':
Res, 'a': param_a, 'b': param_b}
#
# Sampling
by Stan
#
sm = pystan.StanModel(file =
'PM4Cat2004.stan')
fit = sm.sampling(data = stan_data,
init = f_init,
pars = ['sgm', 'C', 'mu'], n_jobs = 1)
print(fit)
fout.write('\nSummary...\n{}\n'.format(fit))
SmplSgm = fit['sgm']
# Samples from the
posterior distribution for sigma
SmplC = fit['C']
# Samples fron the
posterior distribution for C
SmplMu = fit['mu']
# Samples from the
posterior distribution for mu
#
# Percentile
points for posterior distributions
#
SgmL025 = np.percentile(SmplSgm, 2.5)
SgmMed = np.percentile(SmplSgm, 50)
SgmU975 = np.percentile(SmplSgm, 97.5)
print('SgmL25 = {0:.5f} SgmU975 =
{1:.5f}'.format(SgmL025, SgmU975))
CL025 = np.percentile(SmplC, 2.5)
CMed = np.percentile(SmplC, 50)
CU975 = np.percentile(SmplC, 97.5)
MuL025 = np.percentile(SmplMu, 2.5)
MuMed = np.percentile(SmplMu, 50)
MuU975 = np.percentile(SmplMu, 97.5)
#
# MAP
estimates for mu, sigma, and C
#
muMAP, sgmMAP, CMAP = calcMAP(X, Res, MuMed,
SgmMed, CMed, param_a, param_b)
print('muMAP = {0} sgmMAP = {1} CMAP = {2}'.format(muMAP, sgmMAP,
CMAP))
#
# KDE plot
for the samples from posterior distribution for mu
#
fout.write('\nMAP for mu =
{0:<.2f} 95% CI for
mu = [{1:<.2f}, {2:<.2f}]\n'.
format(muMAP, MuL025, MuU975))
sb.kdeplot(SmplMu)
plt.plot([MuL025, MuU975], [0.0, 0.0],
color = 'b', linewidth = 5, label = '95% CI')
plt.plot(muMAP, 0.0, 'go', markersize
= 10, label = 'MAP')
plt.yticks([])
plt.xlabel('$\mu$', fontsize = 14)
plt.title('Posterior Distribution for
$\mu$\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(muMAP,
MuL025, MuU975)),
fontsize = 16)
plt.legend()
plt.tight_layout(h_pad = 5)
plt.show()
#
# KDE plot
for the samples from posterior distribution for sigma
#
fout.write('\nMap for sigma =
{0:<.2f} 95% CI for
sigma = [{1:<.2f}, {2:<.2f}]\n'.
format(sgmMAP, SgmL025, SgmU975))
sb.kdeplot(SmplSgm)
plt.plot([SgmL025, SgmU975], [0.0, 0.0],
color = 'b', linewidth = 5, label = '95% CI')
plt.plot(sgmMAP, 0.0, 'go', markersize
= 10, label = 'MAP')
plt.yticks([])
plt.xlabel('$\sigma$', fontsize = 14)
plt.title('Posterior Distribution for
$\sigma$\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(sgmMAP,
SgmL025, SgmU975)),
fontsize = 16)
plt.legend()
plt.tight_layout(h_pad = 5)
plt.show()
#
# KDE plot
for the samples from posterior distribution for C
#
fout.write('\nMAP for C =
{0:<.2f}
95% CI for C = [{1:<.2f}, {2:<.2f}]\n'.
format(CMAP, CL025, CU975))
sb.kdeplot(SmplC)
plt.plot([CL025, CU975], [0.0, 0.0],
color = 'b', linewidth = 5, label = '95% CI')
plt.plot(CMAP, 0.0, 'go', markersize =
10, label = 'MAP')
plt.yticks([])
plt.xlabel('C', fontsize = 14)
plt.title('Posterior Distribution for
C\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(CMAP,
CL025, CU975)),
fontsize = 16)
plt.legend()
plt.tight_layout(h_pad = 5)
plt.show()
#
# Prob(z
< s75) = 0.75 for the standard normal distribution
#
s75 = bisect(f_phi, 0.75, 0.5,
1.0)
# f_phi from
my_module.py
print('s75 = ', s75, ' p = ', f_phi(s75))
JND = sgmMAP * s75
PSE = muMAP
fout.write('\nMAP for PSE =
{0:<.2f}, 95% CI for
PSE = [{1:<.2f}, {2:<.2f}]'.format(
PSE, MuL025, MuU975))
fout.write('\nMAP for JND =
{0:<.2f}, 95% CI for
JND= [{1:<.2f}, {2:<.2f}]'.format(
JND, s75 * SgmL025, s75 * SgmU975))
#
# Data and
the psychometric function
#
pf = PF(muMAP, sgmMAP)
# Class PF from
my_module.py
minX = X[np.argmin(X)]
maxX = X[np.argmax(X)]
xvalues = []
yvalues = []
yvalues_L = []
yvalues_R = []
for i in range(1001):
t = minX + i * (maxX - minX)
/ 1000.0
xvalues.append(t)
yvalues.append(pf.v(t))
yvalues_L.append(pf.v(t +
CMAP))
yvalues_R.append(pf.v(t -
CMAP))
plt.title('Psychometric Functions\n' +
('PSE = {0:<.2f}, JND =
{1:<.2f}, C =
{2:<.2f}'.format(PSE, JND, CMAP)),
fontsize = 16)
plt.plot(xvalues, yvalues_L, 'g-',
linewidth = 3, label = 'PF(>1)')
plt.plot(xvalues, yvalues, 'k-',
linewidth = 3, label = 'PF(>2)')
plt.plot(xvalues, yvalues_R, 'b-',
linewidth = 3, label = 'PF(>3)')
plt.yticks([0.0, 0.5, 0.75, 1.0])
#
# Responses are
scaled to values from 0 to 1, and added small fluctuations
#
Res_rn = (Res - 1) / 3 +
np.random.uniform(-0.03, 0.03, len(Res))
plt.scatter(X, Res_rn, alpha = 0.5,
color = 'g', label = 'Data')
plt.plot([minX, maxX], [0.0, 0.0],
color = 'k')
plt.plot([minX, maxX], [1.0, 1.0],
color = 'k')
plt.plot([minX, PSE, PSE], [0.5, 0.5,
0.0], color = 'k')
plt.plot([minX, PSE + JND, PSE + JND],
[0.75, 0.75, 0.0], color = 'k')
plt.legend(loc = 7)
plt.xlabel('Stimulus value', fontsize
= 12)
plt.ylabel('Probability', fontsize =
12)
plt.show()
#
# Observed
responses smoothed by Kernel density estimation
#
w_1 = np.empty(len(Res))
# Weights for Response
>= 2
w_2 = np.empty(len(Res))
# Weights for response
>= 3
w_3 = np.empty(len(Res))
# Weights for response
>= 4
for t in range(len(Res)):
w_1[t] = 0 if Res[t] < 2
else 1
w_2[t] = 0 if Res[t] < 3
else 1
w_3[t] = 0 if Res[t] < 4
else 1
X_range = np.arange(minX, maxX, (maxX
- minX) / 100)
w_k = JND / 2
while True:
pm_1 = [] # KDE values for Response >= 2
pm_2 = []
# KDE values for
Response >= 3
pm_3 = []
# KDE values for
Response >= 4
for v in X_range:
pm_1.append(k_ave(v, X, w_1, w_k))
pm_2.append(k_ave(v, X, w_2, w_k))
pm_3.append(k_ave(v, X, w_3, w_k))
plt.title('Observed and
Estimated PFs\n' +
('PSE = {0:<.2f}, JND =
{1:<.2f}, C = {2:<.2f}'.format(PSE,
JND, CMAP)),
fontsize = 16)
plt.plot(xvalues, yvalues_L,
'g-', linewidth = 3, label = 'PF(>1)')
plt.plot(xvalues, yvalues,
'k-', linewidth = 3, label = 'PF(>2)')
plt.plot(xvalues, yvalues_R,
'b-', linewidth = 3, label = 'PF(>3)')
plt.yticks([0.0, 0.5, 1.0])
plt.plot(X_range, pm_1,
'g--', linewidth = 2, label = 'Res(>1)')
plt.plot(X_range, pm_2,
'k--', linewidth = 2, label = 'Res(>2)')
plt.plot(X_range, pm_3,
'b--', linewidth = 2, label = 'Res(>3)')
plt.xlabel('Stimulus value',
fontsize = 12)
plt.legend(loc = 'upper
left')
plt.show()
print('\nThe current value
of w_k = {}'.format(w_k))
w_k = float(input('The next
value of w_k = '))
if w_k <= 0.0:
break
fout.close()
print('\nOutput file {0} was
saved.\n'.format(fout_nm))
Script in Listing 2 uses the module in Listing 3.
Listing 3 module my_module4Cat2004.py, which is used by the script in Listing 2.
import numpy as np
import scipy.optimize
import scipy.stats as ss
f_phi = scipy.stats.norm.cdf # Cumulative standard normal
distribution function
def my_log(x):
"""
To avoid the error message with 0 values for log function
"""
if x < 1.0e-323:
return -744.0
else:
return np.log(x)
class LLog:
"""
Negative log likelihood
"""
def __init__(self, X, Res,
param_a, param_b):
self.N = len(X)
self.X = X
self.Res = Res
self.a = param_a
self.b = param_b
def v(self, z):
mu =
z[0]
sgm
= z[1] ** 2
# Squared to avoid
negative values
C =
z[2] ** 2
# Squared to avoid
negative values
v =
-ss.gamma.logpdf(mu, self.a, scale = 1/self.b) \
-ss.gamma.logpdf(sgm, self.a, scale = 1/self.b) \
-ss.gamma.logpdf(C, self.a, scale = 1/self.b)
for
t in range(self.N):
if self.Res[t] == 4:
v += -my_log(f_phi((self.X[t] - mu - C) / sgm))
elif self.Res[t] == 3:
v += -my_log(f_phi((self.X[t] - mu) / sgm) -
f_phi((self.X[t] - mu - C) / sgm))
elif self.Res[t] == 2:
v += -my_log(f_phi((self.X[t] - mu + C) / sgm) -
f_phi((self.X[t] - mu) / sgm))
else:
v
+= -my_log(1.0 - f_phi((self.X[t] - mu + C) / sgm))
return v
def calcMAP(X, Res, PSEMed, SgmMed,
CMed, param_a, param_b):
"""
Calculate MAP estimates
"""
llog = LLog(X, Res, param_a,
param_b)
# Negative log
likelihood
z = np.empty(3)
z[0] = PSEMed
# Initial value for mu
z[1] = SgmMed ** 0.5
# Initial value for
square root of sigma
z[2] = CMed ** 0.5
# Initial value for
square root of C
res =
scipy.optimize.minimize(llog.v, z, method = 'Nelder-Mead',
options = {'xatol': 1.0e-7, 'fatol': 1.0e-12})
muMAP = res.x[0]
sgmMAP = res.x[1] ** 2
CMAP = res.x[2] ** 2
return muMAP, sgmMAP, CMAP
class PF:
"""
Psychometric function with the parameters mu and sigma
"""
def __init__(self, mu, sgm):
self.mu = mu
self.sgm = sgm
def v(self, x):
return f_phi((x - self.mu) / self.sgm)
def bisect( f, c, Lb, Ub ):
"""
Return the value root so that f(root) = c
"""
s = 1.0 if f(Ub) > f(Lb)
else -1.0
# s = 1 when f(x) is
an increasing functin
while ((abs(Ub - Lb) >
(1.0e-15) * (abs(Ub) + abs(Lb)))
and
(abs(Ub) + abs(Lb) > 1.0e-14)):
m = 0.5 * (Ub + Lb)
if (s * (f(m) - c)) > 0.0:
Ub = m
else:
Lb = m
return 0.5 * (Ub + Lb)
def my_kernel(x, y, w_k):
"""
Kernel function with width w_k
"""
return np.exp(-((x - y) /
w_k) ** 2)
def k_ave(v, X, W, w_k):
"""
Average by kernel function my_kernel
"""
n = len(X)
num = 0.0
den = 0.0
for i in range(n):
k_v =
my_kernel(v, X[i], w_k)
num
+= W[i] * k_v
den
+= k_v
return num / den