Up-and-Down Method with Two Category Rating
Yasuharu Okamoto
To estimate point of subjective equality (PSE) and just noticeable difference (JND), it is proposed that data of up-and-down method should be analyzed by Bayesian method (Okamoto, 2019). Estimation by Bayesian method is more accurate than that by the traditional method of arithmetic averaging. The model used by the Bayesian method 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 
, then
If 
, then
,
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, 
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 and 2, respectively.
Listing 1 Stan script for up-and-down method with two category rating
data {
int N; // The number of
trials
int Res[N]; // Responses: 1 denotes -, 2 denotes +
real X[N]; // Stimuli
real a; // Parameter alpha
of gamma distribution
real b;
//
Parameter beta of gamma distribution
}
parameters {
real<lower = 0.0> sgm;
real pse;
}
transformed parameters {
vector[2] p[N]; // Probabilities for
X[N]
for (i in 1:N){
p[i][1] = 1 - Phi((X[i] - pse) / sgm);
p[i][2] = Phi((X[i] - pse) / sgm);
}
}
model {
sgm ~ gamma(a, b);
pse ~ 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 files2Cat2004.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’ and ‘Stronger’ are
represented by integers 1 and 2, 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 the psychometric function is drawn is displayed (Figure 6).
Figure 6
In the graph of Figure 6, 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 6 is closed, the window in Figure 7 is displayed.
Figure 7
In Figure 7, proportions of the responses 
are drawn by kernel
density estimation as a broken line curve.
When the window in Figure 7 is closed, the next value of w_k, which is the width of kernel density estimation, is asked (Figure 8). A larger value of w_k produces a smoother curve, and a smaller value variable one.
Figure 8
When w_k is set to be smaller than or equal to 0, the program ends (Figure 9).
Figure 9
After the program ends, the output file can be opened. The content is as follows.
Data file...Data2Cat.csv
Number of data = 100
0 250.00 2
1 230.00 2
2 215.00 2
3 200.00 2
4 185.00 1
5 200.00 1
.
.
.
95 200.00 1
96 215.00 2
97 200.00 1
98 215.00 1
99 230.00 2
Summary...
Inference for Stan model:
anon_model_cd52dc34f0f20fa0cecbd0262b2a71cf.
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 18.47 0.12 4.73 12.22 15.3 17.69 20.6 29.34 1535 1.0
pse 203.09 0.06 2.71 197.51 201.41 203.07 204.75
208.66 1921 1.0
lp__ -49.07 0.03 1.11 -52.02 -49.47 -48.73 -48.3 -48.02 1212 1.0
Samples were drawn using NUTS at Wed
Apr 29 10:54:12 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).
mu(MAP) = 203.089 sigma(MAP) = 15.708
PSE = 203.089, 95% CI = [197.52, 208.63]
JND = 10.60, 95% CI = [8.25, 19.78]
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 files2Cat2004.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
import scipy.optimize as sciopt
from my_module2004 import *
fin_nm = input('Data file name (*.csv)
= ')
fout_nm = input('Output file name
(*.txt) = ')
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'.format(N))
for i in range(N):
print('{0} {1}
{2}'.format(ID[i], X[i], Res[i]))
for i, s, r in zip(ID, X, Res):
fout.write('{0:>5s}{1:>7.2f}{2:>5d}\n'.format(i, s, r))
#
# 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 pse(mu) and sigma
"""
return dict(pse = (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 =
'PM2Cat2004.stan')
fit = sm.sampling(data = stan_data,
pars = ['sgm', 'pse'], init = f_init, n_jobs = 1)
print(fit)
fout.write('\nSummary...\n{}\n'.format(fit))
SmplSgm = fit['sgm']
# Samples from the
posterior distribution for sigma
SmplPSE = fit['pse']
# Samples from the
posterior distribution for pse(mu)
#
# Percentile
points for posterior distributions
#
SgmL025 = np.percentile(SmplSgm, 2.5)
SgmMed = np.percentile(SmplSgm, 50)
SgmU975 = np.percentile(SmplSgm, 97.5)
PSEL025 = np.percentile(SmplPSE, 2.5)
PSEMed = np.percentile(SmplPSE, 50)
PSEU975 = np.percentile(SmplPSE, 97.5)
#
# Calculate
MAP estimates for mu(pse) and sigma
#
#
calcMAP from
my_module.py
#
muMAP, sgmMAP = calcMAP(X, Res,
PSEMed, SgmMed, param_a, param_b)
print('mu = ', muMAP, ' sgm = ', sgmMAP)
fout.write('\nmu(MAP) =
{0:<.3f} sigma(MAP)
= {1:<.3f}\n'.format(muMAP, sgmMAP))
#
# KDE plot
for the samples from posterior distribution for PSE(mu)
#
sb.kdeplot(SmplPSE)
plt.plot([PSEL025, PSEU975], [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,
PSEL025, PSEU975)),
fontsize = 16)
plt.legend()
plt.tight_layout(h_pad = 5)
plt.show()
#
# KDE plot for
the samples from posterior distribution for sigma
#
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()
#
# 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('\nPSE = {0:<.3f}, 95% CI = [{1:<.2f},
{2:<.2f}]'.format(
PSE, PSEL025, PSEU975))
fout.write('\nJND = {0:<.2f}, 95% CI = [{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 = []
for i in range(1001):
t = minX + i * (maxX - minX)
/ 1000.0
xvalues.append(t)
yvalues.append(pf.v(t)) # Values of the psychometric
function
plt.title('Psychometric Function\n' +
('PSE = {0:<.2f}, JND =
{1:<.2f}'.format(PSE, JND)),
fontsize = 16)
plt.plot(xvalues, yvalues, linewidth =
3, label = 'PF')
plt.yticks([0.0, 0.5, 0.75, 1.0])
#
# Resposes
are scaled to values from 0 to 1, and added small fluctuations
#
Res_rn = (Res - 1) +
np.random.uniform(-0.05, 0.05, 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
for t in range(len(Res)):
w_1[t] = 0 if Res[t] < 2
else 1
X_range = np.arange(minX, maxX, (maxX -
minX) / 100)
w_k = JND / 2
while True:
pm_1 = []
# KDE values for
Response >= 2
for v in X_range:
pm_1.append(k_ave(v, X, w_1, w_k))
plt.title('Observed and
Estimated PFs\n' +
('PSE = {0:<.2f}, JND =
{1:<.2f}'.format(PSE, JND)),
fontsize = 16)
plt.plot(xvalues, yvalues,
color = 'b', linewidth = 3, label = 'PF')
plt.yticks([0.0, 0.5, 1.0])
plt.plot(X_range, pm_1,
'g--', linewidth = 2, label = 'Res(>1)')
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 {} was
saved.'.format(fout_nm))
Script in Listing 2 uses the module in Listing 3.
Listing 3 module my_module2004.py, which is used by the script in Listing 2.
import numpy as np
import scipy
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
v =
-ss.gamma.logpdf(mu, self.a, scale = 1/self.b) \
-ss.gamma.logpdf(sgm, self.a, scale = 1/self.b)
for
t in range(self.N):
if self.Res[t] == 2:
v
+= -my_log(f_phi((self.X[t] - mu) / sgm))
else:
v += -my_log(1.0 - f_phi((self.X[t] - mu) / sgm))
return v
def calcMAP(X, Res, PSEMed, SgmMed,
param_a, param_b):
"""
Calculation of MAP estimates
"""
llog = LLog(X, Res, param_a,
param_b)
# Negative log
likelihood
z = np.empty(2)
z[0] = PSEMed
# Initial value for pse(mu)
z[1] = SgmMed ** 0.5 #
Initila value for square root of sigma
import scipy.optimize
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
return muMAP, sgmMAP
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):
return np.exp(-((x - y) /
w_k) ** 2)
def k_ave(v, X, W, w_k):
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