Up

Up-and-Down Method with Three Category Rating

Yasuharu Okamoto

 

To estimate point of subjective equality (PSE) and just noticeable difference (JND), up-down-method with three category rating (Stronger, Do not know, Weaker, e.g.) was proposed (Okamoto, 2019). This up-and-down method with three 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 , then

,

If , then

,

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, and 3, respectively.

 

Listing 1. Stan script for up-and-down method with three category rating.

 

//          Yasuharu Okamoto, 2020.04

data {

    int N;              //      The number of trials

    int Res[N];         //      Responses:  1 denotes -, 2 denotes ?, 3 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[3] p[N];             //      Probabilities for X[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 - C) / sgm);

        p[i][3] = 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 files3Cat2004.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, Do not know, and Stronger are represented by integers 1, 2, and 3, 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 two estimated psychometric functions corresponding to the rating categories and the psychometric function for no judgment criterion 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...Data3Cat.csv

 

Number of data = 100

 

    0      257.50      3

    1      237.50      3

    2      222.50      3

    3      207.50      3

    4      192.50      1

    5      207.50      2

      .

      .

      .

   95      222.50      2

   96      207.50      3

   97      192.50      1

   98      207.50      1

   99      222.50      3

 

Summary...

Inference for Stan model: anon_model_cb0f4bc94392c4f817b8c50ecc769c72.

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

C     10.12    0.04   1.68    7.2   8.95  10.01  11.15  13.76   1807    1.0

mu   203.71    0.04   1.86 199.94 202.48 203.76 204.98 207.22   2497    1.0

sgm   14.77    0.06    2.4  11.01   13.1  14.48  16.06  20.29   1865    1.0

lp__ -70.93    0.03   1.24 -74.16 -71.48 -70.61 -70.03 -69.51   1437    1.0

 

Samples were drawn using NUTS at Wed Apr 29 12:23:33 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 = 203.94    95% CI for mu = [199.94, 207.22]

 

Map for sigma = 13.42    95% CI for sigma = [11.02, 20.27]

 

MAP for C = 9.46      95% CI for C = [7.20, 13.74]

 

MAP for PSE = 203.94,    95% CI for PSE = [199.94, 207.22]

MAP for JND = 9.05,    95% CI for JND= [7.44, 13.67]

 

 

Reference

Okamoto, Y. (2019). Up-down methods with more than two reponse categories. The Japanese Journal of Psychonomic Science, 38, 90-104.

 

 

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 files3Cat2004.zip , which can be freely downloaded and used.

 

import numpy as np

import matplotlib.pyplot as plt

import pystan

import seaborn as sb

import pickle

import csv

from my_module3Cat2004 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)]

 

print('data...\n', data)

 

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]))

 

C1C2 = np.percentile(X, [20, 80])

print('C1C2 = ', C1C2)

 

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, sgm = C1C2[1] - C1C2[0],

                mu = (C1C2[1] + C1C2[0])/2)

 

 

stan_data = {'N': N, 'X':X, 'Res': Res, 'a': param_a, 'b': param_b}

#

#       Sampling by Stan

#

sm = pystan.StanModel(file = 'PM3Cat_EqC2004.stan')

 

fit = sm.sampling(data = stan_data, pars = ['C', 'mu', 'sgm'],

                  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

SmplC = fit['C']        #   Samples from 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)

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))

fout.write('\nMAP for mu = {0:<.2f}    95% CI for mu = [{1:<.2f}, {2:<.2f}]\n'.

           format(muMAP, MuL025, MuU975))

#

#       KDE plot for the samples from posterior distribution for mu

#

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, 'k--', linewidth = 3, label = 'PF')

plt.plot(xvalues, yvalues_L, 'g-', linewidth = 3, label = 'PF(>1)')

plt.plot(xvalues, yvalues_R, 'b-', linewidth = 3, label = 'PF(>2)')

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) / 2 + 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

w_2 = np.empty(len(Res))        #   Weights for response >= 3

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

 

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

    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))

       

    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_R, 'b-', linewidth = 3, label = 'PF(>2)')

    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, 'b--', linewidth = 2, label = 'Res(>2)')

    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_module3Cat2004.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] == 3:

                v += -my_log(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 - C) / 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 values for square root of sigma

    z[2] = CMed ** 0.5          #   Initial values 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):

    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

 

 

Up