差尺度法のベイズ分析
試行単位モデル
差の一対比較データに基づいて尺度構成を行うものに最尤法による差尺度法（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が選ばれたことを表す。
図１の形式の差尺度実験データのベイズ分析スクリプトをリスト１とリスト２に分けて作成した。Windowsで実行するときは、リスト１のスクリプトはWindowsにインストールされたUbuntu上で実行し、その後、リスト2のスクリプトを直接Windows上で実行する。リスト１ではPyStan3.10を利用しているが、これはWindows上では直接利用できないので、WindowsにWSL（Windows Subsystem for Linux）インストールしてUbuntu24をインストールした。Ubuntu24においてMinicondaのPython 3.12仮想環境を用意してリスト１のスクリプトを実行した。使用したWSLではグラフ描画ができなかったので（2025.04）、リスト１でのMCMCサンプルをファイルに保存しておいて、保存されたMCMCサンプルの分析をWindows上でPython 3.13により行いグラフを描いた。
ファイルはまとめて圧縮ファイルba_files.zipとしている。
以下の実行例では、リスト１とリスト２のスクリプトファイルと入力データファイルは同じフォルダ内に置かれているものとする。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
次の例では、図１のファイル名raw_data.csvが設定されている。
(py312) ****/ba_files$ python
BDS_bytrial_ps3.py
Data(*.csv) = raw_data.csv
入力データファイル名を設定すると、次に刺激の数N_stの設定が求められる。図１のデータの場合、刺激数は５であるので、「N_st = 5」を設定する。
N_st = 5
刺激数を設定すると計算が始まる。
計算が終了すると、以下の出力が表示されて、実行終了となる。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta[0] 0.266 0.017 0.235 0.298
0.0
0.0 23654.0 21634.0 1.0
theta[1] 0.200 0.026 0.152 0.251
0.0
0.0 20153.0 19726.0 1.0
theta[2] 0.307 0.028 0.255 0.361
0.0
0.0 22282.0 24185.0 1.0
theta[3] 0.227 0.017 0.194 0.258
0.0
0.0 28993.0 29524.0 1.0
vsgm 0.017 0.002 0.013 0.021
0.0
0.0 27430.0 27133.0 1.0
sgm 0.172 0.022 0.133 0.213
0.0
0.0 27430.0 27133.0 1.0
psi[0] 0.000 0.000 0.000 0.000
0.0
NaN 40000.0 40000.0 NaN
psi[1] 0.266 0.017 0.235 0.298
0.0
0.0 23654.0 21634.0 1.0
psi[2] 0.467 0.019 0.431 0.502
0.0
0.0 27956.0 28073.0 1.0
psi[3] 0.773 0.017 0.742 0.806
0.0
0.0 28993.0 29524.0 1.0
psi[4] 1.000 0.000 1.000 1.000
0.0
0.0 39663.0 38972.0 1.0
(py312) ****/ba_files$
psi[0]の行にNaNとあるのは、０に固定されたパラメータであることによる。
リスト１のスクリプトの終了後、リスト２をWindows上で実行する。
リスト１の出力ファイルが読み込まれ、計算が始まる。
まず、トレース図が表示される。
図２
図２のWindowを閉じると、図３の刺激の感覚強度の事後分布のグラフが表示される。
図３
は、式（DS.3）による仮定であることに注意。
図３のWindowを閉じると、σの事後分布を示すグラフが表示される（図４）。
図４
図４のフォームを閉じると、次に刺激と感覚の関係を示すグラフが表示される（図５）。
図５
感覚量の事後分布の中央値、50％確信区間、95％確信区間を求め、横軸に刺激の通し番号をとって描かれている。感覚量の値
と
は、パラメータ推定のための仮定（DS.3）であることに注意。参照線として、（１，0）から（5，1）まで直線が引かれている。感覚量がこの直線上にあるのは、観察者の比較判断が円盤の長さの特性（直径、半径？）に基づいて行われた可能性が示唆される。
刺激図形は次のコードで作成されている。
"""
N_st=5 に対して作成
大きさの異なる５枚の円盤のイメージをリストL_patchesに格納
"""
L_patches = []
for ip in range(5):
temp_s =
pygame.Surface((scl_img, scl_img))
c_r = scl_img*0.5*0.2 + 0.8
* scl_img * 0.5 * ip/4
pygame.draw.circle(temp_s,
(255,255,255), (scl_img/2, scl_img/2), c_r)
L_patches.append(temp_s)
刺激のラベル番号と円盤の半径は線形関係にある。
図５のフォームを閉じると、比較判断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章）、誠信書房
import stan
import numpy as np
import scipy.stats as ss
import arviz as az
import csv
import pickle
StanCode = """
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> N;
//array[N_data] int<lower = 1> TN;
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;
//array[N_data] real<lower=0.0, upper=1.0> 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);
//theta ~ uniform(0.0, 1.0);
vsgm
~ beta(1.0, 1.0);
/*
総試行数ＴＮにおいて（i,j）がＮ回選ばれたとする
*/
for
(i in 1:N_data)
//N[i] ~ binomial(TN[i],
//
Phi((abs(psi[ID_i[i]]
- psi[ID_j[i]])
//
- abs(psi[ID_s[i]] - psi[ID_t[i]])) / sgm));
X[i] ~ bernoulli(Phi((abs(psi[ID_i[i]] - psi[ID_j[i]])
- abs(psi[ID_s[i]] - psi[ID_t[i]])) / sgm));
}
"""
fin_nm = input('Data(*.csv) = ')
with open(fin_nm, 'r') as f:
raw_data = [v for v in
csv.reader(f)]
with open('raw_data.pkl', 'wb') as f:
pickle.dump(raw_data, f)
N_st = int(input('N_st = '))
with open('N_st.pkl', 'wb') as f:
pickle.dump(N_st, f)
"""
データ（i, j, s, t, X）の設定
"""
N_data = 0
ID_i = []
ID_j = []
ID_s = []
ID_t = []
X = []
#N = []
#TN = []
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}
sm = stan.build(StanCode, data = Data)
fit = sm.sample(num_samples=10000)
print(fit)
"""
InferenceData型としてMCMCサンプリングの保存
"""
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
with open('inf_data.pkl', 'wb') as f:
pickle.dump(inf_data, f)
print(az.summary(inf_data))
"""
DataFrame型としてMCMCサンプルの保存
"""
df_smpl = fit.to_frame()
with open('df_smpl.pkl', 'wb') as f:
pickle.dump(df_smpl, f)
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import arviz as az
import csv
import pickle
rng = np.random.default_rng()
with open('raw_data.pkl', 'rb') as f:
raw_data = pickle.load(f)
with open('N_st.pkl', 'rb') as f:
N_st = pickle.load(f)
"""
データ（i, j, s, t, N, TN）の設定
"""
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]))
with open('inf_data.pkl', 'rb') as f:
inf_data = pickle.load(f)
print(az.summary(inf_data,
var_names=['theta', 'sgm']))
az.plot_trace(inf_data, var_names=['theta', 'sgm'],
figsize=(12,7)) #, 'psi', 'sgm']) #, var_names=['psi'])
plt.tight_layout()
plt.show()
with open('df_smpl.pkl', 'rb') as f:
df_smpl = pickle.load(f)
print(df_smpl)
"""
ψの事後分布
"""
Psi = [0]*N_st
for i in range(N_st):
Psi[i] = df_smpl[f'psi.{i+1}'].values
for i in range(N_st):
print(Psi[i])
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],
#+rng.uniform(-0.1,0.1),
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()