順序カテゴリ評定尺度
PyStan 3 対応
刺激の感覚の強さを順序カテゴリで答えたデータに基づく尺度構成を、Thurstone（1927）の判断の基本モデルを適用したカテゴリ判断の法則（The law of categorical judgment; Torgerson, 1958, p. 205）を用いたベイズ分析により行った。確率モデルを、サーストンに基づいて、以下のように構成した。
提示された刺激
の感覚の強さをK個のカテゴリから選んで答えるものとする。
刺激の引き起こす感覚量を確率変数
で表し、平均
、標準偏差
の正規分布に従うものとする。すなわち、
である。
感覚量の連続体上におけるK個のカテゴリの境界を
で、刺激の評価カテゴリを
で表すとき、
のとき 
とする。ただし、
である。
評価カテゴリは、通常、言葉で表される。例えば、「非常に弱い」、「弱い」、「どちらかと言えば弱い」、「どちらでもない」、「どちらかと言えば強い」、「強い」、「非常に強い」である。この場合、カテゴリの総数は
であり、各カテゴリに整数値を順に、「非常に弱い」には１，「弱い」には２，・・・、「非常に強い」には７というように、１からKまでの整数値を順番にカテゴリに割り当てて考える。
このとき、刺激の感覚がカテゴリ
であると判断される確率は、
である。ここで、
は累積標準正規分布関数である。
感覚連続体の原点と単位を決めるために、カテゴリ境界の両端を０と１に設定する。
、 
感覚の強さの範囲（ダイナミックレンジ）が刺激の種類に依らずほぼ同じであるという仮定（Teghtsoonian, 1971）に従えば、この制約条件は自然なものである。
感覚の標準偏差については、以下の制約を置く。
これは、Stevens（1975, p. 235）がEkmanの法則と呼んでいるものを踏まえたものである。
上式において、
のとき、
となる。これは、Fechnerの考え方に対応するものである。
ベイズ分析を行うためのStanスクリプトをリスト１のように用意した。PyStan 3対応である。
リスト１ 順序カテゴリ尺度構成Stanスクリプト例（RatingEkman.stan）
data {
int NSt;
int K;
array[NSt] real Sts;
array[NSt, K] int Rs;
}
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 {
array[NSt] real Psis;
array[NSt] real sgms;
array[K-1] real C;
array[NSt] simplex[K] theta;
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.1);
b ~ exponential(0.1);
for (s in 1:NSt) {
Rs[s] ~ multinomial(theta[s]);
}
}
リスト１のStanスクリプトを用いてデータを分析するPythonスクリプトMainCatScal.pyをリスト２のように作成した。サンプルデータファイルなどとともに圧縮ファイルCatScalePS3.zipにまとめたので、ダウンロード展開して自由に利用されたい、PythonでのStan（PyStan）の利用については、岡本（2019）などを参照されたい。
Windowsでリスト２のスクリプトを実行するときは、WindowsにUbuntuをインストールして、WindowsにインストールされたUbuntu上で実行する。WindowsへのUbuntuのインストールはWSL (Windows Subsystem for Linux)により簡単に行えるが、このウェブサイトで説明している。
リスト２のスクリプトは、UbuntuにインストールされたAnaconda上のPython 3.11仮想環境において、次のコマンドで実行できる。
(py311) ・・・/CatScalePS3$ python MainCatScal.py
上のコマンドを打ち込んで実行が始まると、入力データファイル名の設定が求められる。
(py311) ・・・/CatScalePS3$ python MainCatScal.py
Input file (*.xlsx) =
入力データファイルはExcelファイルの形式で用意する（図１）。
図１
第１行目は、変数名である。２行目からデータを、１列目に提示刺激、２列目に各評価カテゴリの頻度を入力する。刺激は、上の行から順番に強くなるように並べる。社会的・文化的刺激（絵画など）のように、強度の順序が不明の場合は、カテゴリに仮の数値、例えば図１の例では、左から１，２，３，４，５，６，７と数値を与えて、刺激ごとに平均値あるいは中央値を求め、それに基づいて刺激の強さの順序を推測することができる。カテゴリの数値の平均値でも大まかな刺激の数量化が出来るが、さらに本ウェブサイトの分析法を用いれば、計量心理学のモデルに基づいたより精緻な分析が可能である。社会的・文化的刺激をカテゴリ数値の平均値で順序付けたときは、図１の第１欄の刺激値には例えば順番に１から整数値を当てはめておけばよい。この数値はグラフを描くときに用いられるが、分析の結果やには影響しない。
入力データファイル名（上の場合は、Data.xlsx）を設定してEnterキーを押すと、次に出力ファイル名の設定が求められる。
Input file (*.xlsx) = Data.xlsx
Output file (*.txt) = Results.txt
適当なテキストファイル名を設定すればよい。上の例では、Results.txtと設定されている。出力ファイル名を設定すると、データが読み込まれ、計算が始まる。
StanスクリプトによるMCMCサンプリングが終了すると、MCMCサンプリングのトレース図が表示される（図２）。
図２
図２のフォームを閉じると、の事後分布が表示される（図３）。
図３
サーストンの比較判断モデルの場合、感覚量の原点と単位は不定であるので、カテゴリ境界の制約条件
、によって原点と単位が決められていることに注意。感覚量
の分布が負の領域にあるのは、このためである。
図３のフォームを閉じると、図４のグラフが表示される。
図４
感覚量の強いものに対する標準偏差の分布が右にシフトしていることが分かる（Ekmanの法則）。
図４のフォームを閉じると、の点推定値との点推定値の散布図が表示される（図５）。点推定値は事後分布の中央値としているが、これはGelmanら（2021）に倣ったものである。
図５
図５のフォームを閉じると、判断の境界値の事後分布が表示される（図６）。
図６
境界値と
は原点と単位を決めるための条件として設定されているので、縦線で位置が示されている。
図６のフォームを閉じると、図７のグラフが表示される。
図７
感覚量が青線で、その上下にの幅の位置が緑色点線で描かれている。判断の境界値が黄色の横線で描かれている。図７のグラフより、これらの関係を視覚的に理解することができる。
図７を閉じると、図８の累積比率を表す曲線が表示される。
図８
各刺激値ごとの評定カテゴリ値に対する累積比率を描いたものである。データに対するものを点線で、モデルに基づくものを実線で示している。
図８のフォームを閉じると、プログラムの実行終了となり、次のような出力が端末に表示される。
Results.txt
was saved.
(py311) ・・・/CatScalePS3$
出力ファイルResults.txtが保存されたことが表示されている。
上の例の場合、出力ファイルの内容は以下のようになっている。
Input Data File = Data.xlsx
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 Med.
Psi Med. Sigma
1.00
-0.04
0.13
1.67 0.05
0.14
2.78
0.13
0.14
4.64
0.25
0.15
7.74
0.34
0.16
12.91
0.47
0.16
21.54
0.63
0.17
35.94
0.78
0.18
59.95
0.90
0.19
100.00
1.08
0.20
感覚量と標準偏差の値が事後分布の中央値で示されているが、中央値を点推定値としているのはGelmanら（2021）に倣ったものである。
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
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.
リスト２ リスト１のStanスクリプトを用いる順序カテゴリ評定尺度構成Pythonスクリプト（MainCatScal.py）
import pandas as pd
import numpy as np
import stan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import arviz as az
pd.set_option('display.max_rows',
10000)
pd.set_option('display.max_columns',
10000)
fin_nm = input('Input file (*.xlsx) =
')
fout_nm = input('Output file (*.txt) =
')
fout = open(fout_nm, 'w')
fout.write('Input Data File =
{}\n'.format(fin_nm))
xlsx_fl = pd.ExcelFile(fin_nm)
data_xlsx =
pd.read_excel(xlsx_fl, header =
None)
data = data_xlsx.values
var_name = data[0]
values = data[1:]
NSt = len(data) - 1
K = len(data[0]) - 1
Sts = np.empty(NSt)
Rs = np.empty((NSt, K), dtype = 'int')
for j in range(NSt):
Sts[j] = data[j+1][0]
for k in range(K):
Rs[j][k] = data[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}
def f_init():
return {'a':1.0, 'b':0.1,
'prePsis':[1.0/(NSt+1)]*(NSt+1),
'preC':[1.0/(K-2)]*(K-2)}
with open('RatingEkman.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_chains=4, init=[f_init()]*4)
print(fit)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['Psis', 'sgms'])
print('Summary =\n',smry)
az.plot_trace(i_data, var_names =
['Psis', 'sgms'])
plt.show()
d_fm = fit.to_frame()
Psis = []
for i in range(NSt):
Psis.append(d_fm[f'Psis.{i+1}'])
Psis = np.array(Psis).T
sgms = []
for i in range(NSt):
sgms.append(d_fm[f'sgms.{i+1}'])
sgms = np.array(sgms).T
Cs = []
for i in range(K-1):
Cs.append(d_fm[f'C.{i+1}'])
Cs = np.array(Cs).T
a = d_fm['a']
b = d_fm['b']
for i in range(NSt):
sb.kdeplot(Psis.T[i], label
= '$\psi_{}$'.format(i+1))
plt.title('Posterior Ditributions of
$\psi_i$', 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_i$', fontsize = 20)
plt.legend()
plt.savefig('FigSgmPost.png')
plt.show()
Est_Psis = np.empty(NSt)
Est_sgms = np.empty(NSt)
for i in range(NSt):
Est_Psis[i] =
np.median(Psis.T[i])
Est_sgms[i] =
np.median(sgms.T[i])
plt.plot(Est_Psis, Est_sgms, 'o')
plt.xlabel('$\psi_i$', fontsize = 14)
plt.ylabel('$\sigma_i$', fontsize =
14)
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 $C_k$',
fontsize = 20)
plt.yticks([])
plt.legend()
plt.savefig('FigCsPost.png')
plt.show()
med_a = np.median(a)
med_b = np.median(b)
v_psis = np.empty(NSt)
v_sgms = np.empty(NSt)
for i in range(NSt):
v_psis[i] =
np.median(Psis.T[i])
v_sgms[i] =
np.median(sgms.T[i])
fout.write('\n{0:>15s}{1:>15s}{2:>15s}\n'.format\
('Stimulus', 'Med. Psi', 'Med. 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] =
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 = '$\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()
fout.close()
print('\n' + fout_nm + ' was
saved.\n')