Up

３点識別試験ベイズ分析

d’の回帰分析例

岡本安晴、2021.12、2022.01

３つの刺激の内、２つが同じ刺激Aで１つが異なる刺激Bであるとき、３つの中で他の２つと異なるものとして正しくBが選ばれる（あるいは、Bが他の２つと最も異なるものとして選ばれる）確率を調べる実験（３点識別試験：triangular test, 3AFC Oddity）において、実験協力者の特性（年齢など）によりAとBの感覚の差がどのように変化しているかを調べる分析例のプログラムを作成してみた。

信号検出理論に基づいて、刺激Aおよび刺激Bの感覚が、それぞれ独立に平均および平均の分散の正規分布（ガウス分布）に従うものとする（サーストンモデル）。２つの刺激Aの感覚をおよびで表し、刺激Bの感覚をで表すと

である。

３点識別試験法において正しくBが選ばれるのは

が成り立つときである（Frijters et al., 1980, p. 177）。

いま、信号検出理論の慣例に倣って、感覚の平均値の差を標準偏差を単位として表し、

とおく。このとき、正答確率（Bが選ばれる確率）と感覚の差の関係は次式で表される（Frijters et al., 1980, eq.(1)）。

式（1）において、の正負を入れ換え（をで置き換え）ても値は変わらない。したがって、感覚量の座標の向きをBの平均の方がAの平均より大きくなる方向にとってとおいても一般性は失われない。式（１）より、感覚の差に対する正答確率が与えられる。

３点識別試験データが、実験参加者の判断の正誤と参加者の特性のペアで与えられているとする。

例えば、図１のデータ（仮想データ）では、判断の正誤が変数Rの値として１（正答）・０（誤答）で与えられ、判断者の年齢が変数Aの値として記入されている。変数CaseIDはデータの識別値で図１では通し番号が記入されている。CaseIDが１の判断者の反応は正答（１）であり、年齢は33歳である。

・

・

・

図１

感覚の差d’と年齢の関係を調べるが、d’には値が負でないという制約があるので、リンク関数（link function; Gelman et al., 2021）を用いてd’の範囲をからに変換して年齢（変数Xで表す）の1次式として表す（単回帰モデル）。すなわち、

とおく。このとき、

である。

式（１）と式（２）より、データに対する尤度はパラメータaとbの関数として与えられる。

である。上式に基づき、リスト１のMCMCサンプリングにおける事後確率の計算は以下のようになっている。Rsが判断Rを表し、Agsが独立変数Xである。

v = my_log(phi(a, 0.0, 10.0)) + my_log(phi(b, 0.0, 10.0))

for i in range(len(self.Rs)):

dprime = np.exp(a + b * self.Ags[i])

if dprime > 10.0:

dprime = 10.0

p = self.TblPDprime[f'{dprime:.2f}']

v += self.Rs[i] * my_log(p) +\

(1 - self.Rs[i]) * my_log(1- p)

独立変数Agsは標準化された値（プログラム内で自動的に標準化される）であるという前提の下で、事前分布は弱情報事前分布として平均０，標準偏差10の正規分布が設定されている。また、式（１）の計算は時間が掛かるので、予めまとめて計算しておき辞書TblPDprimeに格納している。

リスト１のプログラムでは、MCMCサンプリングはMetropolis-within-Gibbs (Robert & Casella, 2010)で行っている（岡本、2014、pp.229-232）。

MCMCサンプリングは、以下のスクリプト（リスト２内）によって４本のチェインが同時に並列実行されている。

#

# Parameter values for multiproccing, each of which runs one chain in MCMC

#

params_t = {'Rs': deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),

'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,

'n_samples': n_samples, 'chain_id': None}

#

# List of parameter values for multiprocessing

#

params = []

for i in range(4):

p_t = deepcopy(params_t)

p_t['chain_id'] = i

params.append(p_t)

print('MCMC sampling by multiprocessing started.')

with Pool() as pool:

Results = pool.map(mm.MCMCs, params) # Run the multiprocessing

４本のチェインが独立に実行されるために必要な独立な乱数生成オブジェクトは、次のスクリプト（リスト１内）により用意している。

rn = RNp2to191()

#

# Get independent chains

# cf. http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/

#

for i in range(self.chain_id):

rn.jump(100)

この乱数生成オブジェクトの詳しい説明は、このウェブサイトを参照されたい。

リスト１のモジュール（ファイル名mymodule.py）を用いてベイズ分析を行うプログラム（main_tritestdrgrssn.py）はリスト２に示した。

プログラムファイルなどは圧縮ファイルTriTestDRgrssnFiles.zipにまとめた。ダウンロード展開すれば、自由に利用できる。

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

(py39) PS D:\xxxxx\TriTestDRgrssnFiles> python .\main_tritestdrgrssn.py

Input file (*.csv) = Data.csv

入力データは、上の図１の形式で用意する。１行目に変数名を書く。1列目はデータの識別変数を置く。変数値は、数値でも文字列でもよい。２列目は、判断の正誤を記入する。正答なら１，誤答なら０である。３列目には、回帰分析の独立変数を置く。

作成したデータは、ファイル拡張子として「.csv」を選ぶとCSV形式で保存される。入力データファイルは、CSV形式のファイルである。

入力ファイル名（上の例の場合、Data.csv）を設定してEnterキーを押すと、計算が始まる。

計算に必要な辞書の作成が終わると、MCMCサンプリングが始まるが、まずproposal分布の調整が行われ、その後、メインのサンプリングが４本のマルコフ鎖で並行して行われる。このサンプリングが終わると、まずパラメータaの４本のマルコフ鎖の分布とトレースのグラフが表示される（図２）。

図２

トレースのグラフの上には（Rhat）の値（Gelman et al., 2014, pp.284-285）が表示されている。

図２のグラフのウィンドウを閉じる（ウィンドウの右上角のX印をクリックする）と、次にパラメータbの分布とトレースが表示される（図３）。

図３

図３のグラフのウィンドウを閉じると、４本のサンプリングをまとめたパラメータaの分布とモードが表示される（図４）。

図４

図４のウィンドウを閉じると、パラメータbの４本のマルコフ鎖をまとめた分布とモードが表示される（図５）。

図５

図５のグラフのウィンドウを閉じると、正答確率の回帰曲線とデータ値のグラフが表示される（図６）。

図６

正答確率の回帰曲線は、パラメータの事後分布のモードをパラメータ値として描かれている。データ値は、独立変数を４つの区間に分割したときの各区間内の平均正答率である。データ値を回帰曲線がよく表していると言える。

図６のウィンドウを閉じると、事後分布のモードをパラメータ値とするd’の回帰曲線のグラフが表示される。

図７

図７のグラフのウィンドウを閉じると、プログラムの実行終了である。

regrslts.txt was saved.

(py39) PS D:\xxxxx\TriTestDRgrssnFiles>

スクリプトファイルのフォルダに出力されたファイルregrslts.txtを開くと、回帰モデルの式が出力されている。

Regression model of Pc on the standardized independent variable Z:

P(R = 1) = Pc(exp(0.70 + -0.36 * Z))

Regression model of Pc on the raw independent variable X:

P(R = 1) = Pc(exp(1.96 + -0.03 * X))

Regression model of d' on the standardized independent variable Z:

d' = exp(0.70 + -0.36 * Z)

Regression model of d' on the raw independent variable X:

d' = exp(1.96 + -0.03 * X)

グラフも画像ファイルとして保存されている。

参考文献

Frijters, J. E. R., Kooistra, A., & Vereijken, P. F. G. (1980). Tables od d' for the triangular method and the 3-AFC signal detection procedure. Perception & Psychophysics, 1980, 27, 176-178.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.

Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., & Rubin, D.B. (2014). Bayesian data analysis, third edition. CRC Press.

岡本安晴（2014）．心理学データ分析と測定．勁草書房

Robert, C. P., & Casella, G. (2010). Monte Carlo statistical methods. Springer

リスト１　MCMCサンプリングなどのモジュール（ファイル名mymodule.py）

"""

Yasuharu Okamoto, 2021.12

"""

import numpy as np

import scipy.stats as ss

import scipy.integrate as si

from RNGen import * # Random number generator

CumPhi = ss.norm.cdf

phi = ss.norm.pdf

def f_phi(z):

return phi(z)

class k_Pc:

def __init__(self, d):

self.d = d

def f(self,u):

"""

Frijters, Kooistra & Vereijken.(1980).

Perception & Psychophysics, 1980, 27, 176-178.

"""

v = (CumPhi(-u*(3**0.5) + self.d * ((2/3)**0.5))

+ CumPhi(-u*(3**0.5) - self.d * ((2/3)**0.5))) \

* phi(u)

return v

def f_Pc(d):

int_k_Pc = k_Pc(d)

v = si.quad(int_k_Pc.f, 0.0, +np.inf)

return v[0] * 2.0

def calcPDprimes():

"""

Dictionary: key = d'

value = probability of the correct response

"""

pdprimes = {}

dprimes = np.linspace(0.0, 10.0, 1001)

for i, d in enumerate(dprimes):

if i % 100 == 0:

print(i, end = '\r')

p = f_Pc(d)

pdprimes[f'{d:.2f}'] = p

return pdprimes

def scale_sgm(acpt_r, ck):

"""

Check the size of acceptance proportion

"""

v = 1.0

if (acpt_r > 0.5):

ck += 1

v = 1.0 + 9.0 * (acpt_r - 0.5) / 0.5

elif (acpt_r < 0.3):

ck += 1

v = 1.0 / (1.0 + 9.0 * (0.3 - acpt_r) / 0.3)

return v, ck

def my_log(v):

if v <= 0.0:

return -750.0

else:

return math.log(v)

class MCMC:

"""

Metropolis-within-Gibbs algorithm

cf. Robert,C.P. & Casella, G. (2010)

Monte Carlo Statistical Methods. p.393

岡本安晴 (2014)．心理学データ分析と測定．pp.229-232.

"""

def __init__(self, Rs, Ags, TblPDprime, init_a, init_b,

sgm_a, sgm_b, n_samples, chain_id):

self.Rs = Rs

self.Ags = Ags

self.TblPDprime = TblPDprime

self.n_samples = n_samples

self.chain_id = chain_id

self.a_coeff = np.empty(self.n_samples + 1)

self.b_coeff = np.empty(self.n_samples + 1)

self.a_coeff[0] = init_a

self.b_coeff[0] = init_b

self.sgm_a = sgm_a

self.sgm_b = sgm_b

self.naccpt_a = 0

self.naccpt_b = 0

def LL(self, a, b):

"""

Log likelihood of the parameters a and b

"""

v = my_log(phi(a, 0.0, 10.0)) + my_log(phi(b, 0.0, 10.0))

for i in range(len(self.Rs)):

dprime = np.exp(a + b * self.Ags[i])

if dprime > 10.0:

dprime = 10.0

p = self.TblPDprime[f'{dprime:.2f}']

v += self.Rs[i] * my_log(p) +\

(1 - self.Rs[i]) * my_log(1- p)

return v

def do_mcmc(self):

"""

Doing MCMC sampling

"""

import tkinter as tk # For message window

if self.chain_id < 0:

msg0 = 'Adjusting sgm...step-{}'.format(-self.chain_id)

else:

msg0 = 'Proc-{}'.format(self.chain_id)

root = tk.Tk()

root.title(msg0)

cnvs = tk.Canvas(width = 700, height = 300, background = '#00ffff')

cnvs.pack()

cnvs.create_text(350, 150, text = '', font = ('', 90), fill = 'blue',

tags = 'MyText')

rn = RNp2to191()

#

# Get independent chains

# cf. http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/

#

for i in range(self.chain_id):

rn.jump(100)

prev_LL = self.LL(self.a_coeff[0], self.b_coeff[0]) # The initial values

den = prev_LL # The initial value

for t in range(self.n_samples):

if (t % 100) == 0:

cnvs.itemconfig('MyText', text = f'{t}/{self.n_samples}')

cnvs.update()

#

# The next value for the parameter a

#

y = rn.normalMS(self.a_coeff[t], self.sgm_a)

num = self.LL(y, self.b_coeff[t])

den = prev_LL

accpt = math.exp(num - den)

if rn.uni() < accpt:

self.a_coeff[t + 1] = y

self.naccpt_a += 1

prev_LL = num

else:

self.a_coeff[t + 1] = self.a_coeff[t]

prev_LL = den

#

# The next value for the parameter b

#

y = rn.normalMS(self.b_coeff[t], self.sgm_b)

num = self.LL(self.a_coeff[t+1], y)

den = prev_LL

accpt = math.exp(num - den)

if rn.uni() < accpt:

self.b_coeff[t + 1] = y

self.naccpt_b += 1

prev_LL = num

else:

self.b_coeff[t + 1] = self.b_coeff[t]

prev_LL = den

root.destroy()

return self.naccpt_a, self.naccpt_b, self.n_samples, \

self.a_coeff[1:], self.b_coeff[1:]

def MCMCs(params):

"""

The function for muliprocessing

"""

my_mcmc = MCMC(params['Rs'], params['Ags'], params['TblPDprime'],

params['init_a'],

params['init_b'], params['sgm_a'], params['sgm_b'],

params['n_samples'], params['chain_id'])

r = my_mcmc.do_mcmc()

return r

def calcRhat(psi):

"""

Gelman,A., Carlin,J.B.,Stern,H.S., Dunson,D.B.,

Vehtari,A, & Rubin,D.B. (2014)

Bayesian Data Analysis, third edition. pp.284-285.

"""

psi = np.array(psi)

psidj = psi.mean(axis = 1)

psidd = psi.mean()

n = len(psi[0])

m = len(psi)

B = (((psidj - psidd) ** 2).sum()) * n / (m - 1)

sj2 = []

for j in range(m):

sj2.append(((psi[j] - psidj[j]) ** 2).sum() / (n - 1))

sj2 = np.array(sj2)

W = (sj2.sum()) / m

var_hat = W * (n - 1) / n + B / n

Rhat = (var_hat / W) ** 0.5

return Rhat

def EstMode(samples, a = 0.05, n_points = 10000):

"""

Calculatte a MAP estimate from a KDE graph on [Lp, Up]

Lp and Up are 100*a/2 and 100(1-a/2) percentile points of samples

"""

import numpy as np

import scipy.stats as ss

Lp, Up = np.percentile(samples, [100 * a/2, 100 * (1 - a/2)]) # import numpy as np

coord = np.linspace(Lp, Up, n_points)

est_pdf = ss.gaussian_kde(samples).pdf(coord) # import scipy.stats as ss

mode_idx = np.argmax(est_pdf)

Mode_Est = coord[mode_idx]

return Mode_Est

リスト２　リスト１のモジュールを用いて３点識別試験データを分析するプログラム（ファイル名main_tritestdrgrssn.py）

"""

Yasuharu Okamoto, 2021.12, 2022.01

"""

import csv

from copy import deepcopy

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

from multiprocessing.pool import Pool

import mymodule as mm # a personal module, file name: mymodule.py

if __name__ == '__main__':

fin_nm = input('Input file (*.csv) = ')

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

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

Data_in = np.array(Data_in)

Rs = [int(v) for v in Data_in.T[1][1:]]

Ags = [float(v) for v in Data_in.T[2][1:]]

print(len(Rs), len(Ags))

Rs = np.array(Rs)

Ags = np.array(Ags)

Ags_raw = deepcopy(Ags)

Ags_mean = Ags.mean()

Ags_std = Ags.std()

Ags = (Ags - Ags_mean) / Ags_std

# Construction of a dictionary for d's and probabilities of correct response

print("Preparing the dictionary of d's and probabilties")

TblPDprime = mm.calcPDprimes()

#

# Adjusting the proposal distributions

#

i_step = -1

sgm_a = 0.5

sgm_b = 0.5

init_a = 0.0

init_b = 0.0

print('Adjusting the proposals')

while True:

print('step-{0:} started.'.format(-i_step))

my_mcmc = mm.MCMC(deepcopy(Rs), deepcopy(Ags), deepcopy(TblPDprime), init_a, init_b, sgm_a, sgm_b, 2000, i_step)

n_accpt_a, n_accpt_b, L_chain, chain_a, chain_b = my_mcmc.do_mcmc()

acpt_r_a = n_accpt_a / L_chain

acpt_r_b = n_accpt_b / L_chain

ck_a = 0

v, ck_a = mm.scale_sgm(acpt_r_a, ck_a) # Check the acceptance proportions

init_a = chain_a.mean()

if ck_a > 0:

sgm_a *= v # Adjust the sigma of the normal distribution for proposal

ck_b = 0

v, ck_b = mm.scale_sgm(acpt_r_b, ck_b)

init_b = chain_b.mean()

if ck_b > 0:

sgm_b *= v

if ck_a + ck_b == 0: # The sizes of the sigmas OK?

break # Adjustment ends,

i_step -= 1

n_samples = 5000 # Sampling size for each chain in MCMC sampling

#

# Parameter values for multiproccing, each of which runs one chain in MCMC

#

params_t = {'Rs': deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),

'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,

'n_samples': n_samples, 'chain_id': None}

#

# List of parameter values for multiprocessing

#

params = []

for i in range(4):

p_t = deepcopy(params_t)

p_t['chain_id'] = i

params.append(p_t)

print('MCMC sampling by multiprocessing started.')

with Pool() as pool:

Results = pool.map(mm.MCMCs, params) # Run the multiprocessing

print('MCMC sampling ended.')

fout_nm = 'regrslts.txt'

fout = open(fout_nm, 'w') # Output text file

#

# Posterior distributon and trace of the parameter a

#

plt.figure(figsize = (12,5))

plt.subplot(1,2,1)

param_a = [[]]*4

for i in range(4):

param_a[i] = Results[i][3]

for i in range(4):

sb.kdeplot(param_a[i]) # KDE posteriors for the parameter a

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

plt.yticks([])

plt.title('Posterior distributions of a', fontsize = 14)

Rhat_a = mm.calcRhat(param_a)

plt.subplot(1,2,2)

for i in range(4):

plt.plot(range(len(param_a[0])), param_a[i]) # Traces of sampling of the parameter a

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

plt.title(f'Traces of a, Rhat = {Rhat_a:.2f}', fontsize = 14)

plt.savefig('FigDistTrace_a.png')

plt.show()

plt.figure(figsize=(12,5))

plt.subplot(1,2,1)

param_b = [[]] * 4

for i in range(4):

param_b[i] = Results[i][4]

for i in range(4):

sb.kdeplot(param_b[i])

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

plt.yticks([])

plt.title('Posterior distributions of b', fontsize = 14)

Rhat_b = mm.calcRhat(param_b)

plt.subplot(1,2,2)

for i in range(4):

plt.plot(range(len(param_b[i])), param_b[i])

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

plt.title(f'Traces of b, Rhat = {Rhat_b:.2f}', fontsize = 14)

plt.savefig('FigDistTrace_b.png')

plt.show()

a_whole = []

b_whole = []

for i in range(4):

a_whole += [v for v in param_a[i]]

b_whole += [v for v in param_b[i]]

a_mode = mm.EstMode(a_whole)

b_mode = mm.EstMode(b_whole)

print(a_mode, b_mode)

# Posterior distribution of a

sb.kdeplot(a_whole)

plt.title(f'Posterior distribution of a, mode = {a_mode:.2f}', fontsize = 16)

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

plt.yticks([])

plt.savefig('Fig_a.png')

plt.show()

# Posterior distribution of b

sb.kdeplot(b_whole)

plt.title(f'Posterior distribution of b, mode = {b_mode:.2f}', fontsize = 16)

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

plt.yticks([])

plt.savefig('Fig_b.png')

plt.show()

ags_min = np.min(Ags)

ags_max = np.max(Ags)

def catAgs(v):

"""

Categorization into 4 categories

"""

tv = int((v - ags_min) * 4.0 / (ags_max - ags_min))

if tv > 3:

tv = 3

return tv

x_coord_z = np.linspace(ags_min, ags_max, 101) # Grid for x-axis

y_coord_z = np.empty(len(x_coord_z))

for i, v in enumerate(x_coord_z):

d = np.exp(a_mode + b_mode * v) # d' predicted by the regression model

y_coord_z[i] = mm.f_Pc(d) # Predicted probability

plt.plot(x_coord_z, y_coord_z, label = 'Model')

ags_count = np.zeros(4)

r_counts = np.zeros(4) # Number of responses in each of 4 categories

rc_counts = np.zeros(4) # Number of crrect responses in each of 4 categories

for rv, av in zip(Rs, Ags):

ic = catAgs(av)

r_counts[ic] += 1

if rv == 1:

rc_counts[ic] += 1

ave_rc = rc_counts / r_counts # proportions of crrect responses in each category

ave_ags = np.empty(4) # means of the independent variable in each category

step_ags = (ags_max - ags_min) / 8.0 # half of the category width

for i in range(4):

ave_ags[i] = ags_min + step_ags + i * step_ags * 2 # mid points of the categories

plt.plot(ave_ags, ave_rc, 'o', label = 'data')

#

# Conversion of standardized scores back to raw scores

#

x_ags = [(ags_min + i*step_ags*2) for i in range(5)] # positions of the raw score labels

x_ags = np.array(x_ags)

x_ags_raw = x_ags * Ags_std + Ags_mean

x_ags_lbl = [f"{int(v)}" for v in x_ags_raw] # labels for raw scores

plt.xticks(x_ags, x_ags_lbl)

plt.xlabel('Independent variable', fontsize = 14)

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

plt.title('Regression Curve of Probability', fontsize = 16)

plt.legend()

plt.savefig('Fig_regression_p.png')

plt.show()

fout.write('\nRegression model of Pc on the standardized independent variable Z:\n')

fout.write('P(R = 1) = Pc(exp({0:.2f} + {1:.2f} * Z))\n'.format(a_mode, b_mode))

a_raw = a_mode - b_mode * Ags_mean / Ags_std

b_raw = b_mode / Ags_std

fout.write('\nRegression model of Pc on the raw independent variable X:\n')

fout.write('P(R = 1) = Pc(exp({0:.2f} + {1:.2f} * X))\n'.format(a_raw, b_raw))

x_coord_z = np.linspace(ags_min, ags_max, 101) # Grid for x-axis

y_coord_d = np.empty(len(x_coord_z))

for i, v in enumerate(x_coord_z):

y_coord_d[i] = np.exp(a_mode + b_mode * v) # d' predicted by the regression model

plt.plot(x_coord_z, y_coord_d)

plt.xticks(x_ags, x_ags_lbl)

plt.xlabel('Independent variable', fontsize = 14)

plt.ylabel("d'", fontsize = 16)

plt.title("Regression curve of d'", fontsize = 16)

plt.savefig('Fig_regression_d.png')

plt.show()

fout.write("\nRegression model of d' on the standardized independent variable Z:\n")

fout.write("d' = exp({0:.2f} + {1:.2f} * Z)\n".format(a_mode, b_mode))

fout.write("\nRegression model of d' on the raw independent variable X:\n")

fout.write("d' = exp({0:.2f} + {1:.2f} * X)\n".format(a_raw, b_raw))

fout.close()

print('\n', fout_nm, 'was saved.\n')

Up