Ordered Categorical Rating Scale
Psychophysical Point of View
Yasuharu Okamoto, 2022.01
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 mean 
and 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 the origin and 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 are 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 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.
For the three conditions, three Python scripts are developed, which are included in the file CategoryScalesFiles.zip.
The archived file can be downloaded and used freely on user’s responsibility, although all rights are reserved.
The Python script CategoryScaleFree.py in the folder SigmaFree is for the condition F.
The Python script CategoryScaleL.py in the folder SigmaLinear is for the condition L.
The Python script CategoryScaleC.py in the folder SigmaConst is for the condition C.
How to use the scripts is explained bellow.
All the program accept an input data file of the same format, which is shown in Figure 1.
Figure 1. File name DataL.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.
The second to the last (K+1st) columns correspond to the first to the Kth categories, values in which represent the number of rating responses for combinations of stimuli and rating categories, respectively.
The data file should be saved in CSV file format. That is, the file should be saved with the file name, which has file extension .csv.
How to use the programs are explained in the order of Condition F, Condition L, and Condition C.
Listing F1 and Listing F2 show the Stan script and Python script for Condition F.
Run the Python script in Listing F2, which uses the Stan script in Listing F1, an input file name is asked as follows:
(py39) PS
D:\SigmaFree> python .\CategoryScaleFree.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataL.csv
When the data file name (in the above example, DataL.csv) is set and the key ‘Enter’ is pressed down, calculation starts.
After MCMC sampling, graphs of posterior distributions of parametrs are displayed (Figure F1).
Figure F1
Close the window of Figure F1, the next window (Figure F2) appears.
Figure F2
Close the window of Figure F2, calculation of WAIC starts, and the result is displayed as follows.
WAIC2 =
175.73879512540347
(py39) PS
D:\SigmaFree>
Graphs in Figure F1 and Figure F2 were automatically saved in files, file extension of which is png. Some results of computation were saved in the file Results.txt.
Content of Results.txt is as follows:
Input Data File = DataL.csv
Data:
1 63 37 0 0 0 0 0
1.668 33 65 2 0 0 0 0
2.783 18 67 15 0 0 0 0
4.642 4 53 40 3 0 0 0
7.743 3 28 56 13 0 0 0
12.915 1 13 45 30 11 0 0
21.544 0 3 21 34 33 8 1
35.938 0 0 8 20 36 25 11
59.948 0 0 1 15 21 32 31
100 0 0 0 2 7 26 65
Stimulus Mode
Psi Mode Sigma
1.00
-0.02
0.08
1.67
0.04
0.09
2.78
0.10
0.11
4.64
0.20
0.12
7.74
0.29
0.14
12.91
0.41
0.17
21.54
0.59
0.18
35.94
0.75
0.21
59.95
0.89
0.23
100.00
1.09
0.21
WAIC = 175.739
Listing L1 and Listing L2 show the Stan script and Python script for
Condition L.
Run the Python script in Listing L2, which uses the Stan script in Listing L1, an input file name is asked as follows:
(py39) PS
D:\SigmaLinear> python .\CategoryScaleL.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataL.csv
When the data file name (in the above
example, DataL.csv) is set and the key ‘Enter’ is pressed down, calculation
starts.
After MCMC sampling, graphs of posterior distributions of parametrs are displayed (Figure L1).
Figure L1
Close the window of Figure L1, the next window (Figure L2) appears.
Figure L2
Close the window of Figure L2, calculation of WAIC starts, and the result is displayed as follows.
WAIC2 =
165.93447882890743
(py39) PS
D:\SigmaLinear>
Graphs in Figure L1 and Figure L2 were automatically
saved in files, file extension of which is png. Some results of computation
were saved in the file Results.txt.
Content of Results.txt is as follows:
Input Data File = DataL.csv
Data:
1 63 37 0 0 0 0 0
1.668 33 65 2 0 0 0 0
2.783 18 67 15 0 0 0 0
4.642 4 53 40 3 0 0 0
7.743 3 28 56 13 0 0 0
12.915 1 13 45 30 11 0 0
21.544 0 3 21 34 33 8 1
35.938 0 0 8 20 36 25 11
59.948 0 0 1 15 21 32 31
100 0 0 0 2 7 26 65
Stimulus Mode
Psi Mode Sigma
1.00
-0.03
0.09
1.67
0.04
0.10
2.78
0.11
0.11
4.64
0.20
0.12
7.74
0.28
0.14
12.91
0.41
0.16
21.54
0.58
0.18
35.94
0.75
0.21
59.95
0.88
0.23
100.00
1.11
0.26
WAIC = 165.934
Listing C1 and Listing C2 show the Stan script and Python script for
Condition C.
Run the Python script in Listing C2, which uses the Stan script in Listing C1, an input file name is asked as follows:
(py39) PS
D:\SigmaConst> python .\CategoryScaleC.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataC.csv
When the data file name (in the above
example, DataC.csv) is set and the key ‘Enter’ is pressed down, calculation
starts.
After MCMC sampling, graphs of posterior distributions of parametrs are displayed (Figure C1).
Figure C1
Close the window of Figure C1, the next window (Figure C2) appears.
Figure C2
Close the window of Figure C2, calculation of WAIC starts, and the result is displayed as follows.
WAIC2 =
162.4332857781438
(py39) PS
D:\SigmaConst>
Graphs in Figure C1 and Figure C2 were
automatically saved in files, file extension of which is png. Some results of
computation were saved in the file Results.txt.
Content of 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 Mode Psi
1.00
0.00
1.67
0.03
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: mode =
0.15
WAIC = 162.433
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 F1. Stan
script for the case of no constraints on standard deviations s.
data {
int NSt;
int K;
real Sts[NSt];
int Rs[NSt, K];
}
transformed data {
vector[NSt+1] a_Psis;
vector[K-2] a_C;
for (i in 1:(NSt+1)) {
a_Psis[i] = 1.0;
}
for (i in 1:(K-2)) {
a_C[i] = 1.0;
}
}
parameters {
real<lower = 0.0>
sgms[NSt];
simplex[NSt+1] prePsis;
simplex[K-2] preC;
}
transformed parameters {
real Psis[NSt];
real C[K-1];
simplex[K] theta[NSt];
real v;
v = 0.0;
for (s in 1:NSt) {
v =
v + prePsis[s];
Psis[s] = logit(v);
}
C[1] = 0.0;
C[K-1] = 1.0;
v = 0.0;
for (k in 2:(K-2)) {
v =
v + preC[k-1];
C[k]
= C[1] + v;
}
for (s in 1:NSt) {
theta[s][1] = normal_cdf(C[1], Psis[s], sgms[s]);
theta[s][K] = 1.0 - normal_cdf(C[K-1], Psis[s], sgms[s]);
for
(k in 2:(K-1)) {
theta[s][k] = normal_cdf(C[k], Psis[s], sgms[s]) -
normal_cdf(C[k-1], Psis[s], sgms[s]);
}
}
}
model {
for (s in 1:NSt) {
sgms[s] ~ exponential(0.001);
}
prePsis ~ dirichlet(a_Psis);
preC ~ dirichlet(a_C);
//a ~ exponential(0.001);
//b ~ exponential(0.001);
for (s in 1:NSt) {
Rs[s] ~ multinomial(theta[s]);
}
}
Listing F2. Python
script for the case of no constraints on standard deviations s.
import csv
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import estmode as em
import CalcWAIC2 as cw
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')
sm = pystan.StanModel(file =
'RatingFree.stan')
Data = {'NSt':NSt, 'K':K, 'Sts':Sts,
'Rs':Rs}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
Psis = fit['Psis']
sgms = fit['sgms']
Cs = fit['C']
for i in range(NSt):
sb.kdeplot(Psis.T[i], label
= '$\psi_{}$'.format(i+1))
plt.title('Posterior Ditributions of
$\psi_s$', fontsize = 20)
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig('FigPsiPost.png')
plt.show()
for i in range(NSt):
sb.kdeplot(sgms.T[i], label
= r'$\sigma_{}$'.format(i))
plt.title('Posterior Distributions of
$\sigma_s$', fontsize = 20)
plt.legend()
plt.savefig('FigSgmPost.png')
plt.show()
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()
v_psis = np.empty(NSt)
v_sgms = np.empty(NSt)
for i in range(NSt):
v_psis[i] =
em.EstMode(Psis.T[i])
v_sgms[i] =
em.EstMode(sgms.T[i])
fout.write('\n{0:>15s}{1:>15s}{2:>15s}\n'.format('Stimulus',
'Mode Psi', 'Mode Sigma'))
for vst, vpsi, vsgm in zip(Sts,
v_psis, v_sgms):
fout.write('{0:>15.2f}{1:>15.2f}{2:>15.2f}\n'.format(vst, vpsi,
vsgm))
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] =
em.EstMode(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 = '$\psi$')
plt.plot(Sts, v_psis - v_sgms, color =
'g', linestyle = '--', lw = 1,label = '$\psi - \sigma$')
plt.plot(Sts, v_psis + v_sgms, color =
'g', linestyle = '--', lw = 1, label = '$\psi+ \sigma$')
plt.xlabel('Stimulus', fontsize = 16)
plt.ylabel('$\psi$', fontsize = 16)
plt.legend(loc = 'upper left')
plt.title('Relation of Stimulus and
$\psi$', fontsize = 18)
plt.tight_layout()
plt.savefig('PsiSgmC.png')
plt.show()
Probs = np.empty((NSt, K))
for s in range(NSt):
Probs[s][0] = ss.norm(loc =
v_psis[s], scale = v_sgms[s]).cdf(v_Cs[0])
Probs[s][K-1] = 1.0 -
ss.norm(loc = v_psis[s], scale = v_sgms[s]).cdf(v_Cs[K-2])
for k in range(1, K-1):
Probs[s][k] = ss.norm(loc = v_psis[s], scale = v_sgms[s]).cdf(v_Cs[k])
-\
ss.norm(loc = v_psis[s], scale = v_sgms[s]).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()
waic2 = cw.calcwaic2(Rs, fit['theta'])
print('WAIC2 =', waic2)
fout.write('\n\nWAIC =
{0:.3f}\n'.format(waic2))
fout.close()
Listing L1. Stan
script for the case of linear constraint on standard deviations s.
data {
int NSt;
int K;
real Sts[NSt];
int Rs[NSt, K];
}
transformed data {
vector[NSt+1] a_Psis;
vector[K-2] a_C;
for (i in 1:(NSt+1)) {
a_Psis[i] = 1.0;
}
for (i in 1:(K-2)) {
a_C[i] = 1.0;
}
}
parameters {
simplex[NSt+1] prePsis;
simplex[K-2] preC;
real<lower = 0.0> a;
real<lower = 0.0> b;
}
transformed parameters {
real Psis[NSt];
real sgms[NSt];
real C[K-1];
simplex[K] theta[NSt];
real v;
v = 0.0;
for (s in 1:NSt) {
v =
v + prePsis[s];
Psis[s] = logit(v);
}
for (s in 1:NSt) {
sgms[s] = a + b * (Psis[s] - Psis[1]);
}
C[1] = 0.0;
C[K-1] = 1.0;
v = 0.0;
for (k in 2:(K-2)) {
v =
v + preC[k-1];
C[k]
= C[1] + v;
}
for (s in 1:NSt) {
theta[s][1]
= normal_cdf(C[1], Psis[s], sgms[s]);
theta[s][K] = 1.0 - normal_cdf(C[K-1], Psis[s], sgms[s]);
for
(k in 2:(K-1)) {
theta[s][k] = normal_cdf(C[k], Psis[s], sgms[s]) -
normal_cdf(C[k-1], Psis[s], sgms[s]);
}
}
}
model {
prePsis ~ dirichlet(a_Psis);
preC ~ dirichlet(a_C);
a ~ exponential(0.001);
b ~ exponential(0.001);
for (s in 1:NSt) {
Rs[s] ~ multinomial(theta[s]);
}
}
Listing L2. Python script
for the case of linear constraint on standard deviations s.
import csv
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import estmode as em
import CalcWAIC2 as cw
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')
sm = pystan.StanModel(file =
'RatingEkman.stan')
Data = {'NSt':NSt, 'K':K, 'Sts':Sts,
'Rs':Rs}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
Psis = fit['Psis']
sgms = fit['sgms']
Cs = fit['C']
a = fit['a']
b = fit['b']
for i in range(NSt):
sb.kdeplot(Psis.T[i], label
= '$\psi_{}$'.format(i+1))
plt.title('Posterior Ditributions of
$\psi_s$', fontsize = 20)
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig('FigPsiPost.png')
plt.show()
for i in range(NSt):
sb.kdeplot(sgms.T[i], label
= r'$\sigma_{}$'.format(i))
plt.title('Posterior Distributions of $\sigma_s$',
fontsize = 20)
plt.legend()
plt.savefig('FigSgmPost.png')
plt.show()
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()
mode_a = em.EstMode(a)
mode_b = em.EstMode(b)
plt.figure(figsize = (8,4))
plt.subplot(1,2,1)
sb.kdeplot(a)
plt.yticks([])
plt.xlabel('a', fontsize = 14)
plt.title('Posterior Distribution of
a' +
'\nmode = {0:.2f}'.format(mode_a),
fontsize = 14)
plt.subplot(1,2,2)
sb.kdeplot(b)
plt.title('Posterior Distribtuion of
b' +
'\nmode = {0:.2f}'.format(mode_b), fontsize = 14)
plt.yticks([])
plt.xlabel('b', fontsize = 14)
plt.tight_layout()
plt.savefig('Fig_a_b_Post.png')
plt.show()
v_psis = np.empty(NSt)
v_sgms = np.empty(NSt)
for i in range(NSt):
v_psis[i] = em.EstMode(Psis.T[i])
v_sgms[i] =
em.EstMode(sgms.T[i])
fout.write('\n{0:>15s}{1:>15s}{2:>15s}\n'.format('Stimulus',
'Mode Psi', 'Mode Sigma'))
for vst, vpsi, vsgm in zip(Sts,
v_psis, v_sgms):
fout.write('{0:>15.2f}{1:>15.2f}{2:>15.2f}\n'.format(vst, vpsi,
vsgm))
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] =
em.EstMode(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 = '$\psi$')
plt.plot(Sts, v_psis - v_sgms, color =
'g', linestyle = '--', lw = 1,label = '$\psi - \sigma$')
plt.plot(Sts, v_psis + v_sgms, color =
'g', linestyle = '--', lw = 1, label = '$\psi+ \sigma$')
plt.xlabel('Stimulus', fontsize = 16)
plt.ylabel('$\psi$', fontsize = 16)
plt.legend(loc = 'upper left')
plt.title('Relation of Stimulus and
$\psi$', fontsize = 18)
plt.tight_layout()
plt.savefig('PsiSgmC.png')
plt.show()
Probs = np.empty((NSt, K))
for s in range(NSt):
Probs[s][0] = ss.norm(loc =
v_psis[s], scale = v_sgms[s]).cdf(v_Cs[0])
Probs[s][K-1] = 1.0 -
ss.norm(loc = v_psis[s], scale = v_sgms[s]).cdf(v_Cs[K-2])
for k in range(1, K-1):
Probs[s][k] = ss.norm(loc = v_psis[s], scale = v_sgms[s]).cdf(v_Cs[k])
-\
ss.norm(loc = v_psis[s], scale = v_sgms[s]).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()
waic2 = cw.calcwaic2(Rs, fit['theta'])
print('WAIC2 =', waic2)
fout.write('\n\nWAIC =
{0:.3f}\n'.format(waic2))
fout.close()
Listing C1. Stan
script for the case of constant standard deviations s.
data {
int NSt;
int K;
real Sts[NSt];
int Rs[NSt, K];
}
transformed data {
vector[NSt+1] a_Psis;
vector[K-2] a_C;
for (i in 1:(NSt+1)) {
a_Psis[i] = 1.0;
}
for (i in 1:(K-2)) {
a_C[i] = 1.0;
}
}
parameters {
simplex[NSt+1] prePsis;
simplex[K-2] preC;
real<lower = 0.0> sgm;
}
transformed parameters {
real Psis[NSt];
real C[K-1];
simplex[K] theta[NSt];
real v;
v = 0.0;
for (s in 1:NSt) {
v =
v + prePsis[s];
Psis[s] = logit(v);
}
C[1] = 0.0;
C[K-1] = 1.0;
v = 0.0;
for (k in 2:(K-2)) {
v =
v + preC[k-1];
C[k]
= C[1] + v;
}
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 {
prePsis ~ dirichlet(a_Psis);
preC ~ dirichlet(a_C);
sgm ~ exponential(0.001);
for (s in 1:NSt) {
Rs[s] ~ multinomial(theta[s]);
}
}
Listing C2. Python
script for the case of constant standard deviations s.
import csv
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import estmode as em
import CalcWAIC2 as cw
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')
sm = pystan.StanModel(file =
'RatingFechner.stan')
Data = {'NSt':NSt, 'K':K, 'Sts':Sts,
'Rs':Rs}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
Psis = fit['Psis']
sgm = fit['sgm']
Cs = fit['C']
for i in range(NSt):
sb.kdeplot(Psis.T[i], label
= '$\psi_{}$'.format(i+1))
plt.title('Posterior Ditributions of
$\psi_s$', fontsize = 20)
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig('FigPsiPost.png')
plt.show()
for i in range(NSt):
sb.kdeplot(sgm, label =
r'$\sigma_{}$'.format(i))
plt.title('Posterior Distribution of
$\sigma$', fontsize = 20)
plt.legend()
plt.savefig('FigSgmPost.png')
plt.show()
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()
v_psis = np.empty(NSt)
for i in range(NSt):
v_psis[i] =
em.EstMode(Psis.T[i])
v_sgm = em.EstMode(sgm)
fout.write('\n{0:>15s}{1:>15s}\n'.format('Stimulus',
'Mode Psi'))
for vst, vpsi in zip(Sts, v_psis):
fout.write('{0:>15.2f}{1:>15.2f}\n'.format(vst, vpsi))
fout.write('\nSigma: mode =
{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] =
em.EstMode(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 = '$\psi$')
plt.plot(Sts, v_psis - v_sgm, color =
'g', linestyle = '--', lw = 1,label = '$\psi - \sigma$')
plt.plot(Sts, v_psis + v_sgm, color =
'g', linestyle = '--', lw = 1, label = '$\psi+ \sigma$')
plt.xlabel('Stimulus', fontsize = 16)
plt.ylabel('$\psi$', fontsize = 16)
plt.legend(loc = 'upper left')
plt.title('Relation of Stimulus and
$\psi$', fontsize = 18)
plt.tight_layout()
plt.savefig('PsiSgmC.png')
plt.show()
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()
waic2 = cw.calcwaic2(Rs, fit['theta'])
print('WAIC2 =', waic2)
fout.write('\n\nWAIC =
{0:.3f}\n'.format(waic2))
fout.close()