差尺度法のベイズ分析
差の一対比較データに基づいて尺度構成を行うものに最尤法による差尺度法（maximum likelihood difference scaling: MLDS; Maloney & Yang, 2003）がある。現代の統計法として、ベイズ統計法が注目されているので、ベイズ統計法による差尺度法を試みた（岡本、2025、pp.51-53）。ベイズ統計法の場合も基本となる確率モデルは最尤法の場合と同じであるが、ベイズ統計法を用いる場合は略称としてMLDSを用いることは不適切であるので、BDS（Bayesian Difference Scaling）を用いることにする。
なお、ファイルはまとめて圧縮ファイルba_files.zipとした。
差尺度法では、感覚の差の比較判断を行う。刺激
、
、
、
の感覚量を
、
、
、
で表したとき、刺激と刺激の感覚の差は
で表され、刺激と刺激の差は
で表される。いま、
と
の差を
で表したとき、は平均
、分散
の正規分布に従うとする。すなわち、
である。
刺激と刺激の差（違い）の方が刺激と刺激の差より大きいという判断（曲線の不等号記号
を用いて「
」と表す）は
のときに行われるとすると、その確率は
で与えられる。ここで、
は累積標準正規分布関数である。モデル（DS.2）において、感覚量の連続体の原点および単位が未定であるので、
とおいて原点と単位を設定する（Maloney & Yang, 2003）。ここで、Nは刺激値の総数であり、刺激は刺激値の小さいものから大きさの順に通し番号が付けられているものとする。
刺激と刺激の差と刺激と刺激の差の比較判断が
回行われ、その内「」の判断が
回であったとすると、その確率は2項分布
で与えられる。したがって、比較判断の行われた４つ組の刺激、、、にわたる全体の尤度は
で与えられる。ここで、
であり、
は刺激の数である。４つ組の刺激において、
であれば、３つの刺激、、に対して差
と
の比較を行うことになる。
一般に、互いに異なる刺激からなる４つ組における差の比較判断の場合、
の条件のもとでとの差ととの差が構成されたとき（Maloney & Yan, 2003）、刺激の数が個のとき4つ組の可能な数は
で与えられる。差の比較実験においては、この4つ組の数
は
のときは
組であるが、
のときは
組となり4つ組すべてに対して比較データを取得することは困難になる。Maloney & Yang (2003)のシミュレーション実験では、
のときの
組の4つ組すべてを用いたときのパラメータ推定値の分散と、のときの組の4つ組からランダムに10％を選んだ484組に対するデータからのパラメータ推定値の分散はほぼ等しいことが報告されている。Maloneyらは、パラメータ推定値の分散の大きさに関係するのは、刺激の総数ではなく、４つ組の比較判断が行われた試行の総数であると説明している。このウェブサイトにおける実験例では、刺激として５個の刺激が用いられ、4つ組の数は９組である。５個のものから４個を取り出す方法は４通りであるが、さらに４個の組（i,j,s,t）から差尺度の実験用に（(i,j),(s,t)）、（(i,s),(j,t)）、（(i,t),(s,j)）の組み合わせがあることに注意。
モデル（DS.1）においては２つの仮説、が正規分布に従うこと、および正規分布の分散が刺激の組み合わせに関わらず一定値であることが設定されている。これらの仮説に対して差尺度法は頑健であることがシミュレーションによって確認されている（同上）。
差尺度法のための実験例をこのウェブサイトに用意した。実験結果は、図１のように得られる。
図１
本ウェブサイトに用意したスクリプトは、実験結果が提示刺激の組ごとにまとめられた形式のデータを入力ファイルとするので、データをまとめるスクリプトをリスト１のように作成した。
リスト１のスクリプトを実行すると、入力データファイル名を聞いてくる。入力データファイルは、図１の形式のCSVファイルを用意する。次の例では、図１のファイルraw_disc.csvが設定されている。実験結果のファイルraw_disc.csvは、実験を実施したスクリプトの置かれているフォルダ内にあるので、コピーしてリスト１のスクリプトの置かれているフォルダ内に貼り付けてあるものとする。
Input Data
File (*.csv) = raw_disc.csv
入力データファイル名が設定されると、計算が始まり、刺激刺激の数の入力が求められる。次の例では、図１のデータに合わせて５が設定されている。
Number of stimuli = 5
Output File Name (*.csv) =
aggre_disc.csv
刺激数を設定すると、次に出力ファイル名の設定が求められる。出力ファイルはCSV形式のファイルなので、ファイル拡張子は「.csv」とする。上の例では、aggre_disc.csvと設定されている。
CSV形式のファイルaggre_disc.csvは、図２のようである。
図２
1行目１列目のセルに、刺激数５が設定されている。
２行目に、i、j、s、tと刺激の添え字文字が書かれ、続いてNとTNが書かれている。各行には、これらの名前に相当する値が設定される。
図２の形式のCSVファイルが用意できたら、リスト２のスクリプトでMCMCサンプリングを行う。
リスト２を実行すると、入力データファイル名の設定が求められる。次の例では、上のリスト１の実行例で作成されたaggre_disc.csvが設定されている。
(py312) ****/ba_files$ python
BDS_ps3.py
Data(*.csv) = aggre_disc.csv
入力データファイル名が設定されると、計算が始まる。
計算の終了後、以下のように端末に出力されてスクリプトの実行終了である。
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta[0] 0.349 0.021 0.309 0.388
0.0
0.0 13355.0 16800.0 1.0
theta[1] 0.233 0.033 0.171 0.295
0.0
0.0 13027.0 15834.0 1.0
theta[2] 0.295 0.031 0.234 0.352
0.0 0.0 13796.0 16803.0 1.0
theta[3] 0.123 0.023 0.079 0.166
0.0
0.0 15993.0 16316.0 1.0
vsgm 0.013 0.002 0.009 0.017
0.0
0.0 17264.0 19612.0 1.0
sgm 0.132 0.022 0.093 0.174 0.0 0.0 17264.0 19612.0 1.0
psi[0] 0.000 0.000 0.000 0.000
0.0
NaN 40000.0 40000.0 NaN
psi[1] 0.349 0.021 0.309 0.388
0.0
0.0 13355.0 16800.0 1.0
psi[2] 0.581 0.017 0.551 0.614 0.0 0.0 18978.0 23290.0 1.0
psi[3] 0.877 0.023 0.834 0.921
0.0
0.0 15993.0 16316.0 1.0
psi[4] 1.000 0.000 1.000 1.000
0.0
0.0 39639.0 39038.0 1.0
(py312) ****/ba_files$
パラメータpsi[0]の行にNaNとあるのは、このパラメータの値が０に固定されているからである。
上のリスト２のスクリプトの実行で得られたMCMCサンプリングはファイルに保存されている。このファイルは、次のリスト３の実行において読み込まれ、分析される。
リスト３を実行する。
まずトレース図が表示される（図３）。
図３
図３のWindowを閉じると、ψの事後分布のグラフが表示される（図４）。
図４
は、式（DS.3）による仮定であることに注意。
図４のWindowを閉じると、σの事後分布を示すグラフが表示される（図５）。
図５
図５のWindowを閉じると、次に刺激と感覚の関係を示すグラフが表示される（図６）。
図６
感覚量の事後分布の中央値、50％確信区間、95％確信区間を求め、横軸に刺激の通し番号をとって描かれている。感覚量の値
と
は、パラメータ推定のための仮定（DS.3）であることに注意。
提示刺激値の円盤は、次のコードで設定されている。
"""
大きさの異なる５枚の円盤のイメージをリスト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)
刺激の番号は、円盤の半径と1次式の関係にある。図６に示されている感覚の強さは、半径の上に凸な関数である。
図６のフォームを閉じると、データにおける比率N/TNと推定されたパラメータによるモデルに基づく確率の散布図が表示される（図７）。
図７
データにおける比率とモデルに基づく確率はよく対応していると言える。
Maloney, L. T., & Yang, J. N. (2003). Maximum likelihood difference scaling. Journal of Vision, 3, 573-585.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房
import numpy as np
import pandas as pd
fin_nm = input('Input Data File
(*.csv) = ')
df_data = pd.read_csv(fin_nm)
print(df_data)
df_data_val = np.array(df_data.values,
dtype=int)
"""
素データをまとめるための辞書型変数の作成
"""
dict_data = {}
for v in df_data_val:
a = v[1]
b = v[2]
if a > b:
a, b = b, a
c = v[3]
d = v[4]
if c > d:
c, d
= d, c
X = v[5]
a = int(a)
b = int(b)
c = int(c)
d = int(d)
if (a, b, c, d) in
dict_data:
dict_data[(a, b, c, d)]['TN'] += 1
if X
== 1:
dict_data[(a, b,
c, d)]['N'] += 1
else:
if X
== 1:
dict_data[(a, b, c, d)] = {'TN':1, 'N':1}
else:
dict_data[(a, b, c, d)] = {'TN':1, 'N':0}
print(dict_data)
"""
(i,j)と(s,t)の比較と、(s,t)と(i,j)の比較データをまとめる
"""
ck_dict = {}
for k in dict_data.keys():
ki, kj, ks, kt = k
N = dict_data[k]['N']
TN = dict_data[k]['TN']
if not k in ck_dict.keys():
if
(ks, kt, ki, kj) in ck_dict.keys():
ck_dict[(ks, kt, ki, kj)]['N'] += TN - N
ck_dict[(ks, kt, ki, kj)]['TN'] += TN
else:
ck_dict[k] = {'N':N, 'TN':TN}
dict_data = ck_dict
"""
まとめたデータをCSV形式のファイルとして出力
"""
for k in dict_data.keys():
print(k, dict_data[k])
n_st = int(input('Number of stimuli =
'))
fout_nm = input('Output File Name
(*.csv) = ')
with open(fout_nm, 'w') as f:
f.write(f'{n_st},,,,,\n')
f.write('i,j,s,t,N,TN\n')
for k in dict_data.keys():
i =
k[0]
j =
k[1]
s = k[2]
t =
k[3]
TN =
dict_data[k]['TN']
N =
dict_data[k]['N']
f.write(f'{i},{j},{s},{t},{N},{TN}\n')
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;
}
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);
/*
総試行数ＴＮにおいて（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));
}
"""
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(raw_data[0][0]) # 刺激値の総数
"""
データ（i, j, s, t, N, TN）の設定
"""
N_data = 0
ID_i = []
ID_j = []
ID_s = []
ID_t = []
N = []
TN = []
for v in raw_data[2:]:
N_data += 1
ID_i.append(int(v[0]))
ID_j.append(int(v[1]))
ID_s.append(int(v[2]))
ID_t.append(int(v[3]))
N.append(int(v[4]))
TN.append(int(v[5]))
print('N_data = ', N_data)
for i in range(N_data):
print("{0:>3}: {1},
{2}, {3}, {4}/ {5}, {6}".
format((i + 1), ID_i[i], ID_j[i], ID_s[i], ID_t[i], N[i], TN[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, \
'N':
N, 'TN': TN}
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
with open('raw_data.pkl', 'rb') as f:
raw_data = pickle.load(f)
N_st = int(raw_data[0][0]) # 刺激値の総数
"""
データ（i, j, s, t, N, TN）の設定
"""
N_data = 0
ID_i = []
ID_j = []
ID_s = []
ID_t = []
N = []
TN = []
for v in raw_data[2:]:
N_data += 1
ID_i.append(int(v[0]))
ID_j.append(int(v[1]))
ID_s.append(int(v[2]))
ID_t.append(int(v[3]))
N.append(int(v[4]))
TN.append(int(v[5]))
print('N_data = ', N_data)
for i in range(N_data):
print("{0:>3}: {1},
{2}, {3}, {4}/ {5}, {6}".
format((i + 1), ID_i[i], ID_j[i], ID_s[i], ID_t[i], N[i], TN[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))
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(np.arange(N_st)+1,
fontsize=12)
plt.xlabel('Stimulus ID', fontsize =
14)
plt.ylabel(r'Sensation $\psi$',
fontsize = 14)
plt.legend()
plt.tight_layout()
plt.savefig('FigSummary.png')
plt.show()
Prop_data = np.empty(N_data)
Prop_model = np.empty(N_data)
for i in range(N_data):
Prop_data[i] = N[i] / TN[i]
Prop_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)
"""
データ値とモデル値の関係
"""
plt.figure(figsize = (5,5))
plt.scatter(Prop_model, Prop_data)
plt.plot([0,1], [0,1])
plt.title('Scattergram of (Prob_Model,
Prop_Data)', fontsize = 14)
plt.ylabel('N/TN', fontsize = 14)
plt.xlabel('Prob.', fontsize = 16)
plt.tight_layout()
plt.savefig('FigScatter.png')
plt.show()