ベイズ差尺度法のデータ前処理

ベイズ差尺度法のデータ前処理と分析

差尺度法（Maloney & Yang, 2003）は、2対の刺激について各対の感覚の差の比較であるので、観察者には容易な課題と思われるが、実験者には差尺度の分析ができるようにデータをまとめるのは厄介と思われる。差尺度法のベイズ分析のためのデータのフォーマット例がこのウェブサイトの図１に示されているが、実験の素データをこの形式にまとめるスクリプトを作成してみた。

差尺度法では、刺激iとjの感覚の差と、刺激ｓとｔの感覚の差の比較が求められ、どちらの対の感覚の差がより強いかの判断が求められる。

実験の素データの例（仮想データ）を図１に示す。

図１　素データ例のＣＳＶファイル（raw_data.csv）

図１は、素データをCSV形式のファイルとしたものである。１行目に変数ラベルが書かれている。１列目はデータの識別値であり、分析には用いられないので何でもよい。図１の場合は、０から始まる試行数である。データは２行目から置かれている。１行目は、刺激値iが６、刺激値ｊが５，刺激値ｓが４の場合に、対iとｊの差の方が対ｓとｔの差より大きいと判断されたことを示している。Ｘの値は、対iとｊの方が対ｓとｔより差が強い判断されたＸ＝１であり、対iとｊの方が弱いと判断されたときはＸ＝０である。Ｘは感覚量を表すものではなく、どちらの対の感覚の差の方が大きいと判断されたかを示す指標であることに注意。３行目のデータは、刺激iが１，刺激ｊが２，刺激ｓが６、刺激ｔが３のときに、刺激ｓとｔの感覚の差の方が大きいと判断されたことを表している（Ｘ＝０）。

刺激値は１以上の整数値で表されていることに注意。K個の刺激値があるときは、１，２，・・・、Kで表される。差尺度法の場合、刺激値は順序尺度であればよい。物理値との関係を求めるときは、刺激値を表す整数値を物理値に戻してグラフを描くなどすればよい。差尺度法における刺激値は、単に刺激値を区別するだけのものであるが、Stanスクリプト内において、刺激値を表す整数値が感覚の強さと単調関係にあることが仮定されている。

図１の形式の素データをまとめるスクリプトをリスト１のように用意した。ファイルは圧縮ファイルsmplfiles.zipにまとめられている。

リスト１と図１のファイルを同じフォルダに置いておいてリスト１を実行する。

入力データファイル名の設定が求められる。次の例では、図１のファイルraw_data.csvが設定されている。

Input Data File (*.csv) = raw_data.csv

入力データファイル名を設定すると計算が始まり、計算が進むと刺激値の総数の入力と出力ファイル名の設定が求められる。次の例では、刺激値の総数として６，出力ファイル名としてAggre_data.csvが設定されている。

Number of stimuli = 6

Output File Name (*.csv) = Aggre_data.csv

出力ファイルはCSVファイルの形式で出力されるので、ファイル拡張子は「.csv」である。

出力ファイル名を設定すると、まとめられたデータは保存されてスクリプトの実行終了である。

出力ファイルを開くと、図２のようになっている。

図２　図１のデータをまとめたもの

得られた図２のファイルに対して差尺度法を適用してみる。スクリプトは、このウェブサイトで説明されているものにコメントを追加するなどの多少の修正を加えたものであり、リスト２に示した。

リスト２は、Python 3.7においてpystan2.17.1.0を利用して実行した。また、pystanの実行においてC++コンパイラーが用いられるが、使用PCにはVisual Studio C++がインストールされている。Python 3.7の環境は、Minicondaを利用すると簡単にできる。Miniconda/Anacondaの説明のウェブサイトも用意している。

実行すると、入力データファイル名の設定が求められる。次の例では、図２のファイル名が設定されている。

Data(*.csv) = Aggre_data.csv

ファイル名を入力すると、ファイルを読み込んで計算が始まる。

MCMCサンプリングが終わると、まずψの事後分布が表示される（図３）。

図３

図３のWindowを閉じると、σの事後分布が表示される（図４）。

図４

図４のWindowを閉じると、図５のグラフが表示される。

図５

各刺激値に対する感覚量の事後分布の中央値、50％CIの端点、95％CIの端点が示されている。

図５のWindowを閉じると、図６のグラフが表示される。

図６

データ値N/TNとモデル値Prob(X=1|刺激値対)の散布グラフである。点（0，0）から点（1，1）へ引かれた対角線の周りにデータ点が集まっていることがわかる。

図１の素データを図２のようにまとめるとき、刺激対iとｊ、および刺激対ｓとｔの取り方により、図６の散布図は異なってくる。これは刺激対iとｊの組と刺激対ｓとｔの組を入れ換えることによって生じる。本ウェブサイトの図６とこのウェブサイトの図５が異なるのは、このことによる。本ウェブサイトの図２のデータのまとめ方と、このウェブサイトの図１のまとめ方を比較すると分かる。どちらも素データは同じである。

リスト１　素データをまとめるスクリプト例（Windows上でも実行可：aggregate_data.py）

import numpy as np

import pandas as pd

rng = np.random.default_rng()

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)

print(df_data_val)

"""

素データをまとめるための辞書型変数の作成

"""

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 rng.random() < 0.5:

a, b, c, d = c, d, a, b

X = 0 if X ==1 else 1

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

リスト２　ベイズ差尺度法スクリプト（BDS.py）：リスト３のStanスクリプトを使用

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

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

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

print(fit)

"""

ψの事後分布

"""

Psi = fit['psi'].T

plt.plot([0, 0], [0, 10], label = '$\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('$\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, 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.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 (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()

リスト３　差尺度用Stanスクリプト例（ba_diff_s.stan）

data {

int N_data;

int N_st;

刺激値 i, j, s, t の値

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 {

総試行数ＴＮにおいて（i,j）がＮ回選ばれたとする

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

}

参考文献

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

Home