Up

差尺度法

　差の一対比較データに基づいて尺度構成を行うものに最尤法による差尺度法（maximum likelihood difference scaling: MLDS; Maloney & Yang, 2003）がある。統計法としては、ベイズ統計法もあるので、ベイズ統計法による差尺度法を試みた。ベイズ統計法の場合も基本となる確率モデルは最尤法の場合と同じであるが、ベイズ統計法を用いる場合は略称としてMLDSを用いることは不適切であるので、BDS（Bayesian Difference Scaling）を用いることにする。BDSのためのStanスクリプトは、リスト１に示すように簡単である。

　差尺度法では、感覚の差の比較判断を行う。刺激、、、の感覚量を、、、で表したとき、刺激と刺激の感覚の差はで表され、刺激と刺激の差はで表される。いま、との差をで表したとき、は平均、分散の正規分布に従うとする。すなわち、

である。

　刺激と刺激の差（違い）の方が刺激と刺激の差より大きいという判断（曲線の不等号記号を用いて「」と表す）はのときに行われるとすると、その確率は

で与えられる。ここで、は累積標準正規分布関数である。モデル（DS.2）において、感覚量の連続体の原点および単位が不定であるので、

とおいて原点と単位を設定する（Maloney & Yang, 2003）。ここで、Nは刺激値の総数であり、刺激は刺激値の小さいものから大きさの順に通し番号が付けられているものとする。

　刺激と刺激の差と刺激と刺激の差の比較判断が回行われ、その内「」の判断が回であったとすると、その確率は

で与えられる。したがって、比較判断の行われた４つ組の刺激、、、にわたる全体の尤度は

で与えられる。ここで、

であり、は刺激の数である。４つ組の刺激において、であれば、３つの刺激、、に対して差との比較を行うことになる。

上記のモデルのStanスクリプトを、リスト１のように用意した。

リスト１　差尺度のベイズ分析のためのStanスクリプト（ファイル名ba_diff_s.stan）

data {

int N_data;

int N_st;

int<lower = 1, upper = N_st> ID_i[N_data];

int<lower = 1, upper = N_st> ID_j[N_data];

int<lower = 1, upper = N_st> ID_s[N_data];

int<lower = 1, upper = N_st> ID_t[N_data];

int<lower = 0> N[N_data];

int<lower = 1> TN[N_data];

}

parameters {

simplex[N_st - 1] theta;

real<lower = 0.0> sgm;

}

transformed parameters {

real psi[N_st];

psi[1] = 0.0;

for (i in 2:N_st)

psi[i] = psi[i - 1] + theta[i - 1];

}

model {

for (i in 1:N_data)

N[i] ~ binomial(TN[i],

Phi((fabs(psi[ID_i[i]] - psi[ID_j[i]])

- fabs(psi[ID_s[i]] - psi[ID_t[i]])) / sgm));

}

リスト１を利用するPythonスクリプトをこのウェブサイト後ろにリスト２として示した。スクリプトファイルとサンプルデータファイルは、圧縮ファイルB_Diff_Scale.zipとしてまとめた。ダウンロード展開すれば、自由に利用できる。

リスト２のプログラム用のデータ例を図１に示す。

図１

1行目１列目のセルに、刺激数６が設定されている。

２行目に、i、j、s、tと刺激の添え字文字が書かれ、続いてNとTNが書かれている。各行には、これらの名前に相当する値が設定される。

値の設定後、CSV形式のファイルとして保存する。

一般に、互いに異なる刺激からなる４つ組における差の比較判断の場合、の条件のもとでとの差ととの差が構成されたとき（Maloney & Yan, 2003）、刺激の数が個のとき4つ組の可能な数は

で与えられる。差の比較実験においては、この4つ組の数はのときは組であるが、のときは組となり4つ組すべてに対して比較データを取得することは困難になる。Maloney & Yang (2003)のシミュレーション実験では、のときの組の4つ組すべてを用いたときのパラメータ推定値の分散と、のときの組の4つ組からランダムに10％を選んだ484組に対するデータからのパラメータ推定値の分散はほぼ等しいことが報告されている。Maloneyらは、パラメータ推定値の分散の大きさに関係するのは、刺激の総数ではなく、４つ組の比較判断が行われた試行の総数であると説明している。

　モデル（DS.1）においては２つの仮説、が正規分布に従うこと、および正規分布の分散が刺激の組み合わせに関わらず一定値であることが設定されている。これらの仮説に対して差尺度法は頑健であることがシミュレーションによって確認されている（同上）。

図１のデータは、感覚尺度を構成する刺激値の範囲から６個の刺激を選び、４つ組による差の比較判断実験を行ったところ、図１に示す結果を得たとするものである。刺激は１から６までの通し番号で表されていて、刺激の強さの順であることに注意。各４つ組に対して３０回の判断が行われているが、これは１人の個人からのデータと考えてもよいし、当該の感覚に関して同じ感覚特性をもつ複数の観察者による判断をまとめたものと考えてもよい。

リスト２のプログラムを実行すると、入力データファイル名の入力が求められる。

次の例では、データファイルとして図１のファイルData.csvが設定されている。

(py39) PS D:\XXXXX\B_Diff_Scale> python BDS.py

INFO:numexpr.utils:NumExpr defaulting to 8 threads.

Data(*.csv) = Data.csv

データファイル名を設定してEnterキーを押すと、計算が始まる。

差の比較判断の判断数と各４つ組に対する総判断数に対してモデル（DS.2）が適用され、事後分布が求められる。

まず、刺激の感覚強度の事後分布が、図２のように表示される。

図２

は、式（DS.3）による仮定であることに注意。

図２の表示されているフォームを閉じると、σの事後分布を示すグラフが表示される（図３）。

図３

図３のフォームを閉じると、次に刺激と感覚の関係を示す折れ線グラフが表示される（図４）。

図４

感覚量の事後分布の中央値、50％確信区間、95％確信区間を求め、横軸に刺激の通し番号をとって描かれている。感覚量の値とは、パラメータ推定のための仮定（DS.3）であることに注意。

図４のフォームを閉じると、データにおける比率N/TNと推定されたパラメータによるモデルに基づく確率の散布図が表示される（図５）。

図５

データにおける比率とモデルに基づく確率はよく対応していると言える。

引用文献

Maloney, L. T., & Yang, J. N. (2003). Maximum likelihood difference scaling. Journal of Vision, 3, 573-585.

リスト２　差尺度法のPythonスクリプト（ファイル名BDS.py）

import pystan

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import scipy.stats as ss

import csv

import pickle

fin_nm = input('Data(*.csv) = ')

with open(fin_nm, 'r') as f:

raw_data = [v for v in csv.reader(f)]

#print(raw_data)

N_st = int(raw_data[0][0])

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]))

#sys.exit()

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 = pystan.StanModel(file = 'ba_diff_s.stan')

with open('sm.pkl', 'wb') as f:

pickle.dump(sm, f)

with open('sm.pkl', 'rb') as f:

sm = pickle.load(f)

fit = sm.sampling(data = Data, n_jobs = 1)

with open('fit.pkl', 'wb') as f:

pickle.dump(fit, f)

with open('fit.pkl', 'rb') as f:

fit = pickle.load(f)

print(fit)

Psi = fit['psi'].T

plt.plot([0, 0], [0, 10], label = '$\psi1$')

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))

plt.legend()

plt.yticks([])

plt.xlabel('$\psi$', fontsize = 16)

plt.title('Posterior Distributions', fontsize = 18)

plt.savefig('FigPsi.png')

plt.show()

Sgm = fit['sgm']

sb.kdeplot(Sgm)

plt.xlabel('$\psi$', fontsize = 14)

plt.title('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, 'b-', linewidth = 2, label = 'Median')

plt.plot(stimuli, Psi_25, 'g-.', linewidth = 1, label = '50% CI')

plt.plot(stimuli, Psi_75, 'g-.', linewidth = 1)

plt.plot(stimuli, Psi_025, 'y--', linewidth = 1, label = '95% CI')

plt.plot(stimuli, Psi_975, 'y--', linewidth = 1)

plt.legend()

plt.title('Sensation and Stimulus', fontsize = 18)

plt.xlabel('Stimulus', fontsize = 14)

plt.ylabel('Sensation', fontsize = 14)

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 (Prop_Data, Prob_Model)', fontsize = 16)

plt.ylabel('N/TN', fontsize = 14)

plt.xlabel('Prob.', fontsize = 16)

plt.tight_layout()

plt.savefig('FigScatter.png')

plt.show()

なお、PyStan (Python用Stan) の説明を、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。

Up