差尺度法のベイズ分析
試行単位モデル
差の一対比較データに基づいて尺度構成を行うものに最尤法による差尺度法（maximum likelihood difference scaling: MLDS; Maloney & Yang, 2003）がある。現代の統計法として、ベイズ統計法が注目されているので、ベイズ統計法による差尺度法を試みた（岡本、2025）。ベイズ統計法の場合も基本となる確率モデルは最尤法の場合と同じであるが、ベイズ統計法を用いる場合は略称としてMLDSを用いることは不適切であるので、BDS（Bayesian Difference Scaling）を用いることにする。
差尺度法では、感覚の差の比較判断を行う。刺激
、
、
、
の感覚量を
、
、
、
で表したとき、刺激と刺激の感覚の差は
で表され、刺激と刺激の差は
で表される。いま、
と
の差を
で表したとき、は平均
、分散
の正規分布に従うとする。すなわち、
である。
刺激と刺激の差（違い）の方が刺激と刺激の差より大きいという判断（曲線の不等号記号
を用いて「
」と表す）は
のときに行われるとすると、その確率は
で与えられる。ここで、
は累積標準正規分布関数である。モデル（DS.2）において、感覚量の連続体の原点および単位が未定であるので、
とおいて原点と単位を設定する（Maloney & Yang, 2003）。Nは刺激値の総数であり、刺激は刺激値の小さいものから大きさの順に通し番号が付けられているものとする。
第τ試行における４つ組の提示刺激値を
、
、
、
で表し、その時の比較判断を
で表す。
このとき、比較判断の行われた総試行にわたる尤度は
で与えられる。ここで、
であり、
は刺激の数である。４つ組の刺激において、
であれば、３つの刺激、、に対して差
と
の比較を行うことになる。
一般に、互いに異なる刺激からなる４つ組における差の比較判断の場合、
の条件のもとでとの差ととの差が構成されたとき（Maloney & Yan, 2003）、刺激の数が個のとき4つ組の可能な数は
で与えられる。差の比較実験においては、この4つ組の数
は
のときは
組であるが、
のときは
組となり4つ組すべてに対して比較データを取得することは困難になる。Maloney &
Yang (2003)のシミュレーション実験では、
のときの
組の4つ組すべてを用いたときのパラメータ推定値の分散と、のときの組の4つ組からランダムに10％を選んだ484組に対するデータからのパラメータ推定値の分散はほぼ等しいことが報告されている。Maloneyらは、パラメータ推定値の分散の大きさに関係するのは、刺激の総数ではなく、４つ組の比較判断が行われた試行の総数であると説明している。
このウェブサイトにおける実験例では、刺激として５個の刺激が用いられ、試行ごとに4つの刺激値がランダムに選ばれている。総試行数は400試行である。
モデル（DS.1）においては２つの仮説、が正規分布に従うこと、および正規分布の分散が刺激の組み合わせに関わらず一定値であることが設定されている。これらの仮説に対して差尺度法は頑健であることがシミュレーションによって確認されている（同上）。
実験例のウェブサイトの実験を行うと、データとして図１のCSV形式のファイルが得られる。
図１
1行目に、i、j、s、tと刺激を表す文字が書かれ、続いてXが書かれている。１列目はデータ（試行）の通し番号である。
i、j、s、tは提示刺激を表し、Xは比較判断を表す。X＝１は対i,jが選ばれたことを表し、X＝０は対s,tが選ばれたことを表す。
図１の形式の差尺度実験データのベイズ分析を行うStanスクリプトをリスト１のように作成した。リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２のように作成した。これらのファイルは、ba_files.zipにまとめた。
リスト１のファイル（BDS_bytrial.stan）、リスト２のファイル（BDS_bytrial.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境において、次のコマンドを実行する。CmdStanPyの説明を、このウェブサイトで行っている。
(stan) *****/ba_files$ python BDS_bytrial.py
上のコマンドを実行すると、入力データファイル名の設定が求められる。次の例では、図１のファイル名が設定されている。
(stan) *****/ba_files$
python BDS_bytrial.py
Data(*.csv) = raw_data.csv
N_st = 5
入力データファイル名を設定すると、次に刺激数の設定が求められる。上の例では、入力データraw_data.csvに合わせて、５が設定されている。
刺激数を設定すると、ファイルが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了するとトレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３の刺激の感覚強度の事後分布のグラフが表示される。
図３
は、式（DS.3）による仮定であることに注意。
図３のWindowを閉じると、σの事後分布を示すグラフが表示される（図４）。
図４
図４のフォームを閉じると、次に刺激と感覚の関係を示すグラフが表示される（図５）。
図５
感覚量の事後分布の中央値、50％確信区間、95％確信区間を求め、横軸に刺激の通し番号をとって描かれている。感覚量の値
と
は、パラメータ推定のための仮定（DS.3）であることに注意。参照線として、（１，0）から（5，1）まで直線が引かれている。感覚量がこの直線上にあるのは、観察者の比較判断が円盤の長さの特性（直径、半径？）に基づいて行われた可能性が示唆される。刺激のラベル番号と円盤の半径は線形関係にある。
図５のWindowを閉じると、比較判断Resp（X＝１：上を選択、X＝０：下を選択）と推定されたパラメータによるモデルの確率の散布図が表示される（図６）。
図６
比較判断はX＝１または０の2値しかとらないので複数のデータ点が重なり見難くなる。これを避けるためにジッターを加えてある。
データとモデルに基づく確率はよく対応していると言える。
Maloney, L. T., & Yang, J. N. (2003). Maximum likelihood difference scaling. Journal of Vision, 3, 573-585.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房
リスト１ ベイズ分析のためのStanスクリプト（BDS_bytrial.stan）
data {
int N_data;
int N_st;
/*
刺激値 i, j, s, t の値
*/
array[N_data]
int<lower = 1, upper = N_st> ID_i;
array[N_data]
int<lower = 1, upper = N_st> ID_j;
array[N_data]
int<lower = 1, upper = N_st> ID_s;
array[N_data]
int<lower = 1, upper = N_st> ID_t;
array[N_data]
int<lower=0, upper=1> X;
}
transformed data {
vector[N_st-1] a;
for (i
in 1:N_st-1)
a[i] = 1/(N_st-1.0);
}
parameters {
simplex[N_st
- 1] theta;
real<lower = 0.0, upper =
1.0> vsgm;
}
transformed parameters {
real sgm;
/*
感覚量ψ
*/
array[N_st]
real psi;
sgm
= 0.00001 + vsgm*10.0;
psi[1] = 0.0;
for (i
in 2:N_st)
psi[i] = psi[i - 1] + theta[i - 1];
}
model {
theta ~ dirichlet(a);
vsgm
~ beta(1.0, 1.0);
for (i
in 1:N_data)
X[i] ~ bernoulli(Phi((abs(psi[ID_i[i]] - psi[ID_j[i]])
- abs(psi[ID_s[i]] - psi[ID_t[i]])) / sgm));
}
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（BDS_bytrial.py）
from cmdstanpy
import CmdStanModel
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
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(input('N_st
= '))
"""
データ（i, j, s, t, X）の設定
"""
N_data = 0
ID_i = []
ID_j = []
ID_s = []
ID_t = []
X = []
for v in raw_data[1:]:
N_data
+= 1
ID_i.append(int(v[1]))
ID_j.append(int(v[2]))
ID_s.append(int(v[3]))
ID_t.append(int(v[4]))
X.append(int(v[5]))
print('N_data
= ', N_data)
for i in
range(N_data):
print("{0:>3}: {1},
{2}, {3}, {4}/ {5}".
format((i + 1), ID_i[i], ID_j[i],
ID_s[i], ID_t[i], X[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, \
'X':
X}
model = CmdStanModel(stan_file="BDS_bytrial.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
df_smpl = fit.draws_pd()
# Pandas DataFrameに変換
rng = np.random.default_rng()
print(az.summary(inf_data, var_names=['theta', 'sgm']))
az.plot_trace(inf_data, var_names=['theta', 'sgm'], figsize=(12,7))
plt.tight_layout()
plt.show()
"""
ψの事後分布
"""
Psi = [0]*N_st
for i in
range(N_st):
Psi[i]
= df_smpl[f'psi[{i+1}]'].values
plt.plot([0, 0], [0, 10], label =
r'$\psi1$') # ψ0=0 に固定
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)) # ψN_st=1 に固定
plt.legend()
plt.yticks([])
plt.xlabel(r'$\psi$', fontsize
= 16)
plt.title('Posterior Distributions', fontsize = 18)
plt.savefig('FigPsi.png')
plt.show()
"""
σの事後分布
"""
Sgm = df_smpl['sgm']
sb.kdeplot(Sgm)
plt.xlabel(r'$\sigma$', fontsize
= 14)
plt.title(r'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,
c='b', marker='o', lw=0, label = 'Median')
plt.plot(stimuli, Psi_25, c='g', marker='o', fillstyle='none', lw=0, label =
'50% CI')
plt.plot(stimuli, Psi_75, c='g', marker='o', fillstyle='none', lw=0)
plt.plot(stimuli, Psi_025, c='y', marker='o', fillstyle='none', lw=0, label =
'95% CI')
plt.plot(stimuli, Psi_975, c='y', marker='o', fillstyle='none', lw=0)
plt.plot([stimuli[0], stimuli[N_st-1]], [0,
1], c='r')
plt.title('Sensation and Stimulus', fontsize = 18)
plt.xticks(stimuli, fontsize=14)
plt.xlabel('Stimulus ID', fontsize
= 14)
plt.ylabel(r'Sensation
$\psi$', fontsize = 14)
plt.legend()
plt.tight_layout()
plt.savefig('FigSummary.png')
plt.show()
Prob_model = np.empty(N_data)
for i in
range(N_data):
Prob_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)
"""
データ値とモデル値の関係
"""
for i in
range(N_data):
plt.plot(Prob_model[i],
X[i]+rng.uniform(-0.1,0.1),
'o', fillstyle='none')
plt.title('Scattergram of (Prob_Model,
Resp+jitter)', fontsize =
14)
plt.ylabel('Resp', fontsize
= 14)
plt.yticks([0, 1], fontsize=14)
plt.xlabel('Prob.', fontsize
= 16)
plt.tight_layout()
plt.show()