カテゴリ尺度法／評定尺度法
岡本安晴、2021.09
刺激の感覚の強さを順序カテゴリで答え、それに基づいて尺度を構成する方法は、カテゴリ尺度法（Gescheider, 1997, p. 214）、評定尺度法（難波・桑野、1998, pp. 72-86）、あるいは継次カテゴリ法／系列カテゴリ法（method of successive categories）（岡本、2006、pp. 27-36）、継次間隔法（method of successive intervals）（中谷、1969, p. 162）と呼ばれている。カテゴリ尺度法のデータは、刺激の感覚量とカテゴリ境界の関係にThurstone（1927）の判断の基本モデルを適用したカテゴリ判断の法則（The law of categorical judgment; Torgerson, 1958p. 205）を用いて分析することができる。
提示された刺激
の感覚の強さをK個のカテゴリから選んで答えるものとする。刺激の引き起こす感覚量を確率変数
で表し、平均
、標準偏差
の正規分布に従うものとする。すなわち、
である。
感覚量の連続体上におけるK個のカテゴリの境界を
で、評価カテゴリを
で表すとき、
のとき 
とする。
このとき、刺激の感覚がカテゴリ
であると判断される確率は、
である。ここで、
は累積標準正規分布関数である。
感覚連続体の原点と単位を決めるために、カテゴリ境界の両端を０と１に設定する。
、 
感覚の強さの範囲（ダイナミックレンジ）が刺激の種類に依らずほぼ同じであるという仮定（Teghtsoonian, 1971）に従えば、この制約条件は自然なものである。
感覚の標準偏差については、以下の３つの場合を想定した。
F条件： 標準偏差に特に制約条件を設定しない。
L条件： 標準偏差に制約条件
を設定する。これはStevens（1975, p. 235）がEkmanの法則と呼んでいるものを踏まえたものである。
C条件： 標準偏差は全て等しく
とする場合。これは、Fechnerの考え方に対応するものである。
上記３条件に対応してPythonプログラムを３通り作成した。これらのプログラムとデータ例のファイルは圧縮ファイルCategoryScalesFiles.zipにまとめた。
３つのプログラムの入力データの形式はすべて同じである。例を図１に示す。
第１行目は、変数名である。２行目からデータを、１列目に提示刺激、２列目から各評価カテゴリの頻度を入力する。
図１
刺激は、上の行から順番に強くなるように並べる。社会的・文化的刺激（絵画など）のように、強度の順序が不明の場合は、カテゴリに仮の数値、例えば図１の例では、左から１，２，３，４，５，６，７と数値を与えて、刺激ごとに平均値あるいは中央値を求め、それに基づいて刺激の強さの順序を推測することができる。カテゴリの数値の平均値でも大まかな刺激の数量化が出来るが、さらに本ウェブサイトの分析法を用いれば、計量心理学のモデルに基づいたより精緻な分析が可能である。社会的・文化的刺激をカテゴリ数値の平均値で順序付けたときは、図１の第１欄の刺激値には例えば順番に１から整数値を当てはめておけばよい。この数値はグラフを描くときに用いられるが、分析の結果やには影響しない。
図１のように作成したデータを、CSV形式で保存する。CSV形式のファイルは、Excelの場合、保存ファイル名のファイル拡張子として.csv選ぶとCSV形式で保存される。
以下に、F条件、L条件、C条件のプログラムについて順番に説明する。
標準偏差に制約条件を置かないときのベイズ分析のStanスクリプトをリストF.1に示すように用意した。
リストF1のStanスクリプトを利用するPythonスクリプトはリストF.2のように用意した。
リストF.2のプログラムを以下のように実行すると、入力データファイル名の設定が求められる。
(py39) PS
D:\SigmaFree> python .\CategoryScaleFree.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataL.csv
入力データファイル名（上の場合は、DataL.csv）を設定してEnterキーを押すと計算が始まる。StanスクリプトによるMCMCサンプリングが終了すると、図F.1のように、パラメータの事後分布が順番に表示される。
図F.1
図F.１の３つのグラフの提示後（フォームの右上角のX印をクリックすればフォームが閉じられる）、しばらくの計算の後、次の図F.2のグラフが表示される。
図F.2
図F.2のグラフの提示後、しばらくの計算の後、WAICの値（WAIC2）が次のように表示されて、プログラムの実行終了となる。
WAIC2 =
175.73879512540347
(py39) PS
D:\SigmaFree>
図F.1および図F.2のグラフはそれぞれファイル（拡張子.png）として保存されている。計算結果もResults.txtに保存されている。
上の分析の場合、Results.txtの内容は以下のようである。
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
標準偏差に線形制約条件を置くときのベイズ分析のStanスクリプトをリストL.1に示すように用意した。
リストL.1のStanスクリプトを利用するPythonスクリプトはリストL.2のように用意した。
リストL.2のプログラムを以下のように実行すると、入力データファイル名の設定が求められる。
(py39) PS
D:\SigmaLinear> python .\CategoryScaleL.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataL.csv
入力データファイル名（上の場合は、DataL.csv）を設定してEnterキーを押すと計算が始まる。StanスクリプトによるMCMCサンプリングが終了すると、図L.1のように、パラメータの事後分布が順番に表示される。
図L.1
図L.1の３つのグラフの提示後、しばらくの計算の後、次の図L.2のグラフが表示される。
図L.2
図L.2のグラフの提示後、しばらくの計算の後、WAICの値（WAIC2）が次のように表示されて、プログラムの実行終了となる。
WAIC2 =
165.93447882890743
(py39) PS
D:\SigmaLinear>
図L.1および図L.2のグラフは、それぞれファイル（拡張子.png）として保存されている。計算結果もResults.txtに保存されている。
上の分析の場合、Results.txtの内容は以下のようである。
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
標準偏差に等分散の制約条件を置くときのベイズ分析のStanスクリプトをリストC.1に示すように用意した。
リストC.1のStanスクリプトを利用するPythonスクリプトはリストC.2のように用意した。
リストC.2のプログラムを以下のように実行すると、入力データファイル名の設定が求められる。
(py39) PS
D:\SigmaConst> python .\CategoryScaleC.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Input file
(*.csv) = DataC.csv
入力データファイル名（上の場合は、DataC.csv）を設定してEnterキーを押すと計算が始まる。StanスクリプトによるMCMCサンプリングが終了すると、図C.1のように、パラメータの事後分布が順番に表示される。
図C.1
図C.１の３つのグラフの提示後、しばらくの計算の後、次の図C.2のグラフが表示される。
図C.2
図C.2のグラフの提示後、しばらくの計算の後、WAICの値（WAIC2）が次のように表示されて、プログラムの実行終了となる。
WAIC2 =
162.4332857781438
(py39) PS
D:\SigmaConst>
図C.1および図C.2のグラフは、それぞれファイル（拡張子.png）として保存されている。計算結果もResults.txtに保存されている。
上の分析の場合、Results.txtの内容は以下のようである。
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
Python用のStan（PyStan）の説明を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
Gescheider, G. A. (1997). Psychophysics: The fundamentals, Third edition. Mahwah: Lawrence Erlbaum Associations, Publishers.
中谷和夫（1969）．第５章 尺度構成法． 八木冕（監修）／田中良久（編）講座心理学２ 計量心理学．東京大学出版会．
難波精一郎・桑野園子（1998）．音の評価のための心理学的測定法．コロナ社．
岡本安晴（2006）．計量心理学．培風館．
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.
リストF.1. 標準偏差に制約がないモデルのStanスクリプト
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]);
}
}
リストF.2. リストF.1を用いるPythonスクリプト
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()
リストL.1. 標準偏差に線形制約を用いるモデルのStanスクリプト
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]);
}
}
リストL.2. リストL.2のStanスクリプトを用いるPythonスクリプト
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()
リストC.1. 標準偏差が等しい等分散モデルのStanスクリプト
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]);
}
}
リストC.2. リストC.1のStanスクリプトを用いるPythonスクリプト
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()