Up

Bayesian Difference Scaling

 

MLDS (maximum likelihood difference scaling; Maloney & Young, 2003) estimates parameter values for the difference scaling model by a maximum likelihood method. A Bayesian method for the difference scaling is presented here.

 

The difference scaling model sets the following:

where,  denotes the difference of magnitudes  and , of stimuli  and , that is,

and  denotes the random variable, which corresponds to the sensation of the difference of  and .  is assumed to have a normal distribution of mean  and variance .

It is assumed that the difference  is judged to be larger than the difference , when . Hence, we have

 denotes the standard normal cumulative distribution.

To get MCMC sampling stable, the following restriction is set:

The origin and unit of the scale is set by the following:

 

When comparison of  and  is made  times, and the number of judgments that  is larger than  is , the probability is given by

The likelihood for the results of the experiment is given by

 

Stan script for the above model is shown in Listing 1.

 

Listing 1.  Stan script for Bayesian Difference Scaling (file name: ba_diff_s.stan)

data {

    int N_data;

    int N_st;

    int<lower = 1, upper = N_st> ID_i[N_data];

    int<lower = 1, upper = N_st> ID_j[N_data];

    int<lower = 1, upper = N_st> ID_s[N_data];

    int<lower = 1, upper = N_st> ID_t[N_data];

    int<lower = 0> N[N_data];

    int<lower = 1> TN[N_data];

}

 

parameters {

    simplex[N_st - 1] theta;

    real<lower = 0.0> sgm;

}

 

transformed parameters {

    real psi[N_st];

    psi[1] = 0.0;

    for (i in 2:N_st)

        psi[i] = psi[i - 1] + theta[i - 1];

}

 

model {

    for (i in 1:N_data)

        N[i] ~ binomial(TN[i],

                        Phi((fabs(psi[ID_i[i]] - psi[ID_j[i]])

                         - fabs(psi[ID_s[i]] - psi[ID_t[i]])) / sgm));

}

 

 

The Python script shown in Listing 2 uses the Stan script in the above. The script files and the sample data files are archived in the zip fileB_Diff_Scale.zip, which can be freely downloaded and used under the userfs responsibility. All rights are reserved.

 

An example of the input data for the script is shown in Figure 1.

Figure 1

 

In the first row of the first column, the number of stimuli is put, e.g., in the Figure 1, 6.

In the second row, characters denoting stimuli are set, then characters N and TN, which denote  and .

From the third row on, each data is set row by row.

Save the data with a file name of the extension g.csvh.

Run the script in Listing 2, the file name for the input data is required as is shown below:

 

(py39) PS D:\XXXXX\B_Diff_Scale> python BDS.py

INFO:numexpr.utils:NumExpr defaulting to 8 threads.

Data(*.csv) = Data.csv

 

Set the input data file name, then press down the Enter key.

Calculation starts.

After MCMC sampling, the posterior distributions of  are shown as in Figure 2.

Figure 2

 

Notice that  are the assumption to set the origin and unit (Eq. (3)).

Close the window of Figure 2, then the posterior distribution of  is displayed (Figure 3).

Figure 3

 

 

Close the window of Figure 3, lines, which show relations of the stimuli and sensations, are displayed (Figure 4).

Figure 4

 

Close the window in Figure 4, scatter gram of data  and probability  estimated by the model is shown (figure 5).B

Figure 5

 

The correspondence of the model and data shown in Figure 5 is good.

Close the windows of Figure 5, then the program completes.

 

 

Reference

Maloney, L. T., & Yang, J. N. (2003). Maximum likelihood difference scaling. Journal of Vision, 3, 573-585.

 

 

Listing 2  Python script for Bayesian Difference Scaling (file name: BDS.py)

import pystan

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import scipy.stats as ss

import csv

 

fin_nm = input('Data(*.csv) = ')

with open(fin_nm, 'r') as f:

    raw_data = [v for v in csv.reader(f)]

 

N_st = int(raw_data[0][0])

 

N_data = 0

ID_i = []

ID_j = []

ID_s = []

ID_t = []

N = []

TN = []

for v in raw_data[2:]:

    N_data += 1

    ID_i.append(int(v[0]))

    ID_j.append(int(v[1]))

    ID_s.append(int(v[2]))

    ID_t.append(int(v[3]))

    N.append(int(v[4]))

    TN.append(int(v[5]))

 

print('N_data = ', N_data)

 

for i in range(N_data):

    print("{0:>3}: {1}, {2}, {3}, {4}/  {5}, {6}".

          format((i + 1), ID_i[i], ID_j[i], ID_s[i], ID_t[i], N[i], TN[i]))

 

Data = {'N_data': N_data, 'N_st': N_st, \

        'ID_i': ID_i, 'ID_j': ID_j, 'ID_s': ID_s, 'ID_t': ID_t, \

        'N': N, 'TN': TN}

 

 

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

 

 

fit = sm.sampling(data = Data, n_jobs = 1)

 

print(fit)

 

Psi = fit['psi'].T

plt.plot([0, 0], [0, 10], label = '$\psi1$')

for i in range(1, N_st-1):

    sb.kdeplot(Psi[i], label = r'$\psi${}'.format(i+1))

plt.plot([1, 1], [0, 10], label = r'$\psi${}'.format(N_st))

plt.legend()

plt.yticks([])

plt.xlabel('$\psi$', fontsize = 16)

plt.title('Posterior Distributions', fontsize = 18)

plt.savefig('FigPsi.png')

plt.show()

 

Sgm = fit['sgm']

sb.kdeplot(Sgm)

plt.xlabel('$\psi$', fontsize = 14)

plt.title('Posterir distribution of $\sigma$', fontsize = 18)

plt.yticks([])

plt.savefig('FigSgm.png')

plt.show()

 

Sgm_med = np.percentile(Sgm, 50)

print('Sgm_med =', Sgm_med)

 

 

Psi_med = []

Psi_025 = []

Psi_25 = []

Psi_75 = []

Psi_975 = []

 

for i in range(N_st):

    q025, q25, med, q75, q975 = np.percentile(Psi[i], [2.5, 25, 50, 75, 97.5])

    Psi_med.append(med)

    Psi_025.append(q025)

    Psi_25.append(q25)

    Psi_75.append(q75)

    Psi_975.append(q975)

   

   

stimuli = []

for i in range(N_st):

    stimuli.append(i + 1)

 

 

plt.plot(stimuli, Psi_med, 'b-', linewidth = 2, label = 'Median')

plt.plot(stimuli, Psi_25, 'g-.', linewidth = 1, label = '50% CI')

plt.plot(stimuli, Psi_75, 'g-.', linewidth = 1)

plt.plot(stimuli, Psi_025, 'y--', linewidth = 1, label = '95% CI')

plt.plot(stimuli, Psi_975, 'y--', linewidth = 1)

plt.legend() 

plt.title('Sensation and Stimulus', fontsize = 18)

plt.xlabel('Stimulus', fontsize = 14)

plt.ylabel('Sensation', fontsize = 14)

plt.tight_layout()

plt.savefig('FigSummary.png')

plt.show()

 

Prop_data = np.empty(N_data)

Prop_model = np.empty(N_data)

for i in range(N_data):

    Prop_data[i] = N[i] / TN[i]

    Prop_model[i] = ss.norm(0, 1).cdf((abs(Psi_med[ID_i[i]-1] - Psi_med[ID_j[i]-1]) -

                                   abs(Psi_med[ID_s[i]-1] - Psi_med[ID_t[i]-1])) / Sgm_med)

 

plt.figure(figsize = (5,5))

plt.scatter(Prop_model, Prop_data)

plt.plot([0,1], [0,1])

plt.title('Scattergram of (Prop_Data, Prob_Model)', fontsize = 16)

plt.ylabel('N/TN', fontsize = 14)

plt.xlabel('Prob.', fontsize = 16)

plt.tight_layout()

plt.savefig('FigScatter.png')

plt.show()

 

 

Up