Ordered Categorical Rating Scale
Psychophysical Point of View
Yasuharu Okamoto, 2022.01, 2026.01
The scripts have been revised for CmdStanPy. A simple installation of CmdStanPy is shown at this website.
In the method of ordered categorical scaling, an observer chooses the appropriate category from the ordered categories corresponding to the strength of the sensation caused by the presented stimulus. The method has been discussed by many researchers (e.g., Gescheider, 1997; Torgerson, 1958; Thurstone, 1927).
Suppose the stimulus 
invokes sensation 
, which has a normal distribution of a mean 
and a variance 
. That is,
Suppose that there are 
boundaries 
s of categories, and observer’s rating 
on the stimulus is category 
when the following
condition holds:
when 
where
Then we have
where 
is the cumulative
distribution function of the standard normal distribution.
Set an origin and a unit by the following constraint:
, 
This condition corresponds to the dynamic range, which is supposed to be approximately the same for all modalities (Teghtsoonian, 1971).
As to standard deviation 
of sensation , the following three cases would be considered.
Condition F:
No constraints on , that is, free condition.
Condition L:
Linear relation between the standard
deviation and the mean of sensation of stimulus is assumed. That is, the
following relation
is set.
This is a generalization of what Stevens (1975, p. 235) calls Ekman’s law, from which Stevens’ power law can be derived.
Condition C:
All s are the same, that is constant,
This condition corresponds to Fechner’s thinking, and leads to Fechner’s logarithmic law.
In the following, only the condition C is considered. The script files are archived in the file SigmaConst.zip. The archived file can be downloaded and used freely.
The script accepts an input data file of the format shown in Figure 1.
Figure 1. File name DataC.csv
In the first row, variable names are set. The values of the variables are set in rows under the first row.
The first column contains values of the stimuli. Values of the weakest to the strongest stimuli are set from the second to the last (bottom) row.
Columns from the second to the last (K+1st) ones correspond to categories from the first to the Kth ones. The number of rating responses for combinations of stimuli and rating categories are set in the corresponding cells.
The data file should be saved in CSV file format. That is, the file should be saved with the file name, which has the file extension .csv.
Listing 1 shows a Stan script (RatingFechner.stan)
of Bayesian analysis for the Condition C, Listing 2 shows
a Python script (CategoryScaleC.py), which uses the script in Listing 1. An
input data file (e.g., DataC.csv) is put in the folder, in which the scripts in
Listings 1 and 2 are put.
Execute
the script file CategoryScaleC.py in an environment with cmdstanpy
installed as shown bellow.
(stan) ****/SigmaConst$ ls
CategoryScaleC.py DataC.csv RatingFechner.stan
(stan) ****/SigmaConst$ python CategoryScaleC.py
Input file (*.csv) =
DataC.csv
After the
input data file name is set, the Stan script is compiled and MCMC sampling is
executed.
A trace
plot is shown like Figure 2.
Figure 2
After MCMC sampling, graphs of posterior distributions of parameters are displayed (Figure 3).
Figure 3
Close the window of Figure 3, the next graphs (Figure 4) appears.
Figure 4
Close the window of Figure 4, the script ends.
Content of the output text file Results.txt is as follows:
Input Data File = DataC.csv
Data:
1 51 36 13 0 0 0 0
1.668 43 46 10 1 0 0 0
2.783 23 48 27 2 0 0 0
4.642 10 44 36 10 0 0 0
7.743 1 21 54 22 2 0 0
12.915 0 11 37 42 9 1 0
21.544 0 0 23 49 25 3 0
35.938 0 0 0 19 49 28 4
59.948 0 0 0 1 32 37 30
100 0 0 0 0 3 16 81
Stimulus Med. Psi
1.00
0.00
1.67
0.02
2.78
0.11
4.64
0.19
7.74
0.31
12.91
0.42
21.54
0.53
35.94
0.76
59.95
0.92
100.00
1.13
Sigma: med =
0.15
Gescheider, G. A. (1997). Psychophysics: The fundamentals, Third edition. Mahwah: Lawrence Erlbaum Associations, Publishers.
Stevens, S. S. (1975). Psychophysics: Introduction to its perceptual, neural, and social prospects. New York: Wiley.
Tesghtsoonian, R. (1971). On the exponents in Stevens’ law and the constant in Ekman’s law. Psychological Review, 78, 71-80.
Thurstone, L. L. (1927). A law of comparative judgment. Psychological Review, 34, 273-286.
Torgerson, W. S. (1958). Theory and methods of scaling. New York: John Wiley & Sons, Inc.
Listing 1 Stan script of Bayesian analysis for the Condition C (RatingFechner.stan)
data {
int NSt;
int K;
array[NSt] real Sts;
array[NSt, K] int Rs;
}
transformed data {
vector[K-2] a_C;
for (i
in 1:(K-2)) {
a_C[i] = 1.0;
}
}
parameters {
//simplex[NSt+1]
prePsis;
array[NSt] real Psis;
simplex[K-2] preC;
real<lower = 0.0001> sgm;
}
transformed parameters {
array[K-1] real C;
array[NSt] simplex[K] theta;
C[1]
= 0.0;
C[K-1] = 1.0;
for (k in 2:(K-2)) {
C[k]
= C[k-1] + preC[k-1];
}
for (s in 1:NSt)
{
theta[s][1] = normal_cdf(C[1]
| Psis[s], sgm);
theta[s][K] = 1.0 - normal_cdf(C[K-1] | Psis[s], sgm);
for
(k in 2:(K-1)) {
theta[s][k] = normal_cdf(C[k] | Psis[s], sgm) -
normal_cdf(C[k-1] | Psis[s],
sgm);
}
}
}
model {
preC
~ dirichlet(a_C);
sgm
~ uniform(0.0001, 1000);
for (s in 1:NSt)
{
Psis[s] ~ normal(0.0, 10.0);
Rs[s] ~ multinomial(theta[s]);
}
}
Listing 2 Python script, which uses the script in Listing1 (CategoryScaleC.py).
import csv
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import arviz
as az
fout = open('Results.txt',
'w')
fin_nm = input('Input
file (*.csv) = ')
with open(fin_nm, 'r') as f:
Data_in
= [v for v in csv.reader(f)]
fout.write('Input Data File = {}\n'.format(fin_nm))
NSt = len(Data_in)
- 1
K = len(Data_in[0]) - 1
Sts = np.empty(NSt)
Rs = np.empty((NSt, K), dtype = 'int')
for j in range(NSt):
Sts[j]
= float(Data_in[j+1][0])
for k
in range(K):
Rs[j][k] = int(Data_in[j+1][k+1])
fout.write('\nData:\n')
for j in range(NSt):
print(f'{Sts[j]:>10g}',
end = '')
fout.write(f'{Sts[j]:>10g}')
for k
in range(K):
print(f'{Rs[j][k]:>5d}', end = '')
fout.write(f'{Rs[j][k]:>5d}')
print()
fout.write('\n')
Data = {'NSt':NSt, 'K':K,
'Sts':Sts,
'Rs':Rs}
model = CmdStanModel(stan_file='RatingFechner.stan')
fit = model.sample(data=Data)
print(fit.diagnose())
print(fit.summary())
infdata = az.from_cmdstanpy(fit) # Transforming to Arviz InfereceData
az.plot_trace(infdata, var_names=['Psis', 'C'])
plt.tight_layout()
plt.savefig('Fig_trace.png')
plt.show()
fit = fit.draws_pd()
#
Transforming to Pandas DataFrame
print(fit.keys())
Psis = []
for s in range(NSt):
Psis.append(fit[f'Psis[{s+1}]'])
Psis = np.array(Psis).T
sgm = fit['sgm']
Cs = []
for k in range(K-1):
Cs.append(fit[f'C[{k+1}]'])
Cs = np.array(Cs).T
plt.figure(figsize=(14,4))
plt.subplot(131)
for i in range(NSt):
sb.kdeplot(Psis.T[i], label = r'$\psi_{}$'.format(i+1))
plt.title(r'Posterior
Ditributions of $\psi_s$', fontsize = 20)
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.subplot(132)
sb.kdeplot(sgm)
plt.title(r'Posterior
Distribution of $\sigma$', fontsize = 20)
plt.legend()
plt.subplot(133)
plt.plot([0, 0], [0, 5], label = r'$C_0$')
for i in range(1, K-2):
sb.kdeplot(Cs.T[i], label = r'$C_{}$'.format(i+1))
plt.plot([1,1], [0, 5], label = r'$C_{}$'.format(K-1))
plt.title('Posterior Distributions of Cs', fontsize = 20)
plt.yticks([])
plt.legend()
plt.savefig('FigCsPost.png')
plt.show()
plt.figure(figsize=(9,4))
plt.subplot(121)
v_psis = np.empty(NSt)
for i in range(NSt):
v_psis[i] = np.median(Psis.T[i])
v_sgm = np.median(sgm)
fout.write('\n{0:>15s}{1:>15s}\n'.format('Stimulus',
'Med. Psi'))
for vst, vpsi in zip(Sts,
v_psis):
fout.write('{0:>15.2f}{1:>15.2f}\n'.format(vst, vpsi))
fout.write('\nSigma:
med = {0:>15.2f}\n'.format(v_sgm))
v_Cs = np.empty(K-1)
v_Cs[0] = 0.0
v_Cs[K-2] = 1.0
for k in range(1,
K-2):
v_Cs[k]
= np.median(Cs.T[k])
for k in range(K-1):
plt.plot([Sts[0], Sts[-1]], [v_Cs[k], v_Cs[k]], color = 'y')
plt.plot([], 'y-', label = 'C')
plt.plot(Sts, v_psis, color = 'b', lw = 3,
label = r'$\psi$')
plt.plot(Sts, v_psis - v_sgm, color = 'g', linestyle = '--', lw = 1,
label = r'$\psi - \sigma$')
plt.plot(Sts, v_psis + v_sgm, color = 'g', linestyle = '--', lw = 1,
label = r'$\psi+ \sigma$')
plt.xlabel('Stimulus', fontsize
= 16)
plt.ylabel(r'$\psi$', fontsize
= 16)
plt.legend(loc = 'upper left')
plt.title(r'Relation
of Stimulus and $\psi$', fontsize = 18)
plt.tight_layout()
plt.subplot(122)
Probs = np.empty((NSt, K))
for s in range(NSt):
Probs[s][0] = ss.norm(loc = v_psis[s],
scale = v_sgm).cdf(v_Cs[0])
Probs[s][K-1] = 1.0 - ss.norm(loc = v_psis[s],
scale = v_sgm).cdf(v_Cs[K-2])
for k
in range(1, K-1):
Probs[s][k] = ss.norm(loc
= v_psis[s], scale = v_sgm).cdf(v_Cs[k]) -\
ss.norm(loc = v_psis[s], scale = v_sgm).cdf(v_Cs[k-1])
print('Predicted Prob(Rating|Stimulus):')
for s in range(NSt):
for k
in range(K):
print(f' {Probs[s][k]:7.3f}', end = '')
print()
Pred_Rs = np.empty((NSt, K))
N_Rs = np.sum(Rs, axis = 1)
for i in range(NSt):
Pred_Rs[i] = N_Rs[i] * Probs[i]
print('Predicted Frequencies of Rating for
Stimulus:')
for s in range(NSt):
for k
in range(K):
print(f' {Pred_Rs[s][k]:7.1f}',
end = '')
print()
cum_data_Rs = np.cumsum(Rs, axis = 1)
cum_pred_Rs = np.cumsum(Pred_Rs, axis =
1)
print(cum_data_Rs)
for s in range(NSt):
if s == 0:
plt.plot(range(K), cum_data_Rs[s]/cum_data_Rs[s][K-1],
'g:', lw = 3,
label = 'Data')
else:
plt.plot(range(K), cum_data_Rs[s]/cum_data_Rs[s][K-1],
'g:', lw = 3 )
if s == 0:
plt.plot(range(K), cum_pred_Rs[s]/cum_pred_Rs[s][K-1],
'b-',
label = 'Model')
else:
plt.plot(range(K), cum_pred_Rs[s]/cum_pred_Rs[s][K-1],
'b-')
plt.xticks(range(K), range(1,
K+1))
plt.xlabel('Rating', fontsize
= 14)
plt.ylabel('Cum. Prop.', fontsize
= 14)
plt.legend()
plt.title('Data and Model Prediction', fontsize = 18)
plt.tight_layout()
plt.savefig('FigDataModel.png')
plt.show()
fout.close()
print('Results.txt was saved.')