Up

上下法（階段法）データのベイズ分析

心理物理学（精神物理学、psychophysics）において主観的等価点を求める方法として上下法（Up-down method）あるいは階段法（Staircase method）と呼ばれている方法があり、優れた方法として紹介されている（相場・鳥居，2004; 大山・岩脇・宮埜，2005）。これは、提示された比較刺激が標準刺激より「強い」と判断されたときは、次の試行ではより弱い比較刺激を提示し、「弱い」と判断されたときは次の試行ではより強い比較刺激を提示するものである。すなわち、第試行における比較刺激とそれに対する判断をとで表し、刺激の変化量をで表すとき、

とおくと、

のとき,　

のとき,　,

である。

このとき、主観的等価点（Point of Subjective Equality: PSE）は提示刺激系列の折り返し点の平均値として求められる。しかし、以下に示すように、確率モデルを用いると心理測定関数（精神測定関数：岡本、2014、pp. 167-169）を求めて、主観的等価点だけでなく弁別閾（Just Noticeable Difference: JND）も求めることができる。

確率モデルを以下のように設定する（岡本、2019；岡本、2011）。

上下法では、提示刺激値が直前の判断に依存して決まるため、被験者が提示刺激値を予測する可能性が高い。これを避けるため、２系列用意して２つの系列の試行をランダムに混ぜて実施する（二重上下法；大山、2005）。上式において、２つの系列を系列１と系列２として、系列でのものであることを右上に添え字として明記すると、以下のようになる。

とおくと、

のとき,　

のとき,　,

データ分析においては、系列１のデータと系列２のデータをまとめてとおき、刺激に対する反応の確率を次式で与える。

ここで、は累積標準正規分布関数である。このとき、尤度関数を次式で与える。

、である。この尤度関数は、データが上下法以外の方法、例えば恒常法、で集められた場合にも適用できることに注意。データは、「強い」判断および「弱い」判断が確実なそれぞれ刺激値が十分に強いものと弱いものを含むと、の事後分布が拡散することを防ぐことができる。

この尤度関数に対して、事後分布を求めるStanスクリプトをリスト１に示す。の値は、−と＋を１と２で表している。

リスト１　上下法データのベイズ分析用Stanスクリプト（PM2CatPS3.stan）

data {

int N; // The number of trials

array[N] int Res; // Responses: 1 denotes -, 2 denotes +

array[N] real X; // Stimuli

real min_sgm; // The lower boundary of sgm

real max_sgm; // The upper boundary of sgm

real p75;

}

parameters {

real<lower = 0.0> v_sgm;

real v_pse;

}

transformed parameters {

real sgm;

real pse;

array[N] vector[2] p; // Probabilities for X[N]

sgm = 0.0001 + v_sgm * (max_sgm - min_sgm) * 0.5;

pse = (min_sgm + max_sgm) / 2.0 + v_pse * (max_sgm - min_sgm);

for (i in 1:N){

p[i][1] = 1 - Phi((X[i] - pse) / sgm);

p[i][2] = Phi((X[i] - pse) / sgm);

}

}

model {

v_sgm ~ exponential(0.1);

v_pse ~ normal(0.0, 10.0);

for (i in 1:N)

Res[i] ~ categorical(p[i]);

}

generated quantities {

real jnd;

jnd = sgm * p75;

}

リスト１のStanスクリプトによる分析プログラム例を付記にリスト２として示す。プログラムファイルとサンプルデータファイルは圧縮ファイルfiles2CatPS3.zipとしてまとめた。ダウンロードして自由に利用されたい。

リスト２のスクリプトをWindowsにインストールしたUbuntu/WSLで実行した。Anaconda/Ubuntu/WSL/WindowsのPython3.11仮想環境にインストールされたPyStan 3.7を利用した。Anacondaのインストールなどは、このウェブサイトを参照されたい。

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

(py311) yasuharu@PSY-PM-PC:/mnt/*****/files2CatPS3$ python MainTwoCatDiscrimPS3.py

Data file name (*.xlsx) =

データファイルは、Excel形式のファイルとして用意する。形式は図１のようである。

・

・

・

図１

第１行目に、変数名を適当に設定する。２行目からデータを設定するが、１列目はデータの識別用である。図２では通し番号が設定されているが、任意の文字列でよい。２列目に提示刺激値を設定し、３列目に反応を設定する。反応は、「弱い」判断のときは１、「強い」判断のときは２を設定する。データの行の順序は任意である。これは、判断の確率が刺激値のみによって与えられるので、例えば、図２の２行目にある０番目のデータと101行目の99番目のデータを入れ替えても尤度関数値は変わらないからである。

データの設定が終われば、保存する。パラメータの事後分布が発散しないように、提示刺激には「強い判断」の確率が１に近い十分強い刺激と「弱い」判断の確率が０に近い十分弱い刺激が含まれるようにする（岡本、2017）。２系列用意する上下法において、「強い」判断が容易な強い刺激値で始まる系列と「弱い」判断が容易な弱い刺激値で始まる系列の２系列を用意すればよい。

入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。適当なテキストファイル名を設定する。

Data file name (*.xlsx) = Data2Cat.xlsx

Output file name (*.txt) = Results.txt

入力ファイル名と出力ファイル名の設定が終わると、データが読み込まれ、Stanによる事後分布からのサンプリングが始まる。サンプリングが終了すると、MCMCサンプリングのトレース図が表示される（図２）。

図２

図２のフォームを閉じると、PSEの事後分布が表示される（図３）。

図３

点推定値として中央値が示されている。Gelmanら（2021）は、平均値と中央値を比べて、中央値の方が計算においてより安定しているとして中央値を勧めている（p. 73）。また、モードについては、データと事前分布に対して全体的に最良の適合を与えると評価している（Gelmanら、2021、p. 107）。岡本（2022）は、点推定値としてモードが最も良かったという報告をしている。点推定値に基づいて心理測定関数を描くときは、モードを用いることにする。モードは、ベイズ分析の事後分布においてはMAP（maximum a posteriori；岡本、2019）推定値と呼ばれている。

図３のフォームを閉じると、JNDの事後分布のグラフが表示される（図４）。

図４

JNDは、に対して次式で与えられる（岡本、2014、p. 169）。

ここで、は累積標準正規分布関数に対して

を満たす値である。

図４のフォームを閉じると、次に心理測定関数PFとデータ値のグラフがMAP推定値に基づくPSEとJNDの推定値とともに表示される（図５）。PFはこれらのMAP推定値に基づいて描かれている。

図５

データ値は、次式

により、０あるいは１に変換した後、上下に少しランダムな変動を加えて、完全に重なって表示されることを防いでいる。

図５のフォームを閉じると、図６のフォームが表示される。

図６

データ（判断R=2の比率）をKDEで滑らかにした曲線（緑の破線）とパラメータPSEとJNDの推定値に基づくPF曲線（青の曲線）が描かれている。

図６のフォームを閉じると、次のw_k値の設定が求められる。

The current value of w_k = 5.456718057006637

The next value of w_k =

このw_kは平滑化におけるkernelの窓の幅で、この値が大きいと大まかな滑らかな曲線が得られ、小さいときは細かく変動した凸凹の曲線になる。

The current value of w_k = 5.456718057006637

The next value of w_k = 9

上のように9を設定すると図７のようなグラフになる。

図７

w_kに０以下の値を設定すると、プログラムの実行終了である。

The current value of w_k = 5.456718057006637

The next value of w_k = 9

The current value of w_k = 9.0

The next value of w_k = 0

Output file Results.txt was saved.

(py311) yasuharu@PSY-PM-PC:/mnt/*****/files2CatPS3$

保存された出力ファイル（上の場合、Results.tx）を開くと以下のようになっている。

Data file...Data2Cat.xlsx

Number of data = 100

0 250.00 2

1 230.00 2

2 215.00 2

　　　・

　　　・

　　　・

97 200.00 1

98 215.00 1

99 230.00 2

Summary:

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat

v_sgm 0.373 0.110 0.221 0.551 0.003 0.002 2380.0 1566.0 1.0

v_pse 1.532 0.027 1.481 1.582 0.001 0.000 2019.0 1384.0 1.0

sgm 18.645 5.520 11.065 27.563 0.154 0.116 2380.0 1566.0 1.0

pse 203.177 2.707 198.132 208.203 0.063 0.045 2019.0 1384.0 1.0

p[0, 0] 0.011 0.023 0.000 0.042 0.001 0.001 2347.0 1466.0 1.0

... ... ... ... ... ... ... ... ... ...

p[98, 0] 0.255 0.065 0.139 0.376 0.001 0.001 2308.0 1321.0 1.0

p[98, 1] 0.745 0.065 0.624 0.861 0.001 0.001 2308.0 1321.0 1.0

p[99, 0] 0.076 0.052 0.004 0.168 0.001 0.001 2328.0 1489.0 1.0

p[99, 1] 0.924 0.052 0.832 0.996 0.001 0.001 2328.0 1489.0 1.0

jnd 12.576 3.723 7.463 18.591 0.104 0.078 2380.0 1566.0 1.0

[205 rows x 9 columns]

PSE(MAP) = 203.051, 95% CI = [197.83, 208.41]

JND(MAP) = 10.91, 95% CI = [8.05, 20.54](

読み込まれたデータの出力後、ベイズ分析の統計量が出力されている。データが正しく読み込まれたかどうかの確認と、ベイズ分析の統計量を出力ファイルでも読むことができる。

引用文献

相場　覚・鳥居修晃（編著）（2004）．改訂版知覚心理学、放送大学教育振興会

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

岡本安晴（2011）上下法データのベイズ的分析 : 確率モデルによるデータの活用、基礎心理学研究、30巻、p. 44-55

岡本安晴（2014）「心理学データ分析と測定――データの見方と心の測り方――」勁草書房

岡本安晴（2017）恒常法における反応カテゴリー数の影響、基礎心理学研究、36、p. 85-91

岡本安晴（2019）「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――」丸善出版

岡本安晴（2019）．多件法による上下法．基礎心理学研究、38巻、90-104.

岡本安晴（2022）ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――．基礎心理学研究、41, 1-7．

大山正・岩脇三良・宮埜嘉夫（共著）（2005）．心理学研究法――データ収集・分析から論文作成まで――．サイエンス社．

付記　Pythonスクリプト

上下法データの分析を行うPythonスクリプトをリスト２およびリスト３として用意した。メインスクリプトがリスト２であり、リスト３のモジュールがimportされている。これらのファイルは、圧縮ファイルfiles2CatPS3.zipとしてまとめられているので、ダウンロードすれば解凍して自由に利用できる。

リスト２　リスト１のStanスクリプトによる上下法データのベイズ分析Pythonスクリプト（MainTwoCatDiscrimPS3.py）．

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import pandas as pd

import stan

import pickle

import seaborn as sb

import scipy.optimize as sciopt

from my_module import *

import arviz as az

fin_nm = input('Data file name (*.xlsx) = ')

fout_nm = input('Output file name (*.txt) = ')

fout = open(fout_nm, 'w')

fout.write('Data file...{}\n'.format(fin_nm))

#

# Read data from Excel file

#

xlsx_fl = pd.ExcelFile(fin_nm)

data_xlsx = pd.read_excel(xlsx_fl, header = None)

data = data_xlsx.values

var_name = data[0]

values = data[1:]

N = len(values)

ID = []

X = np.empty(N, dtype = float)

Res = np.empty(N, dtype = int)

for t in range(N):

ID.append(values[t][0])

X[t] = float(values[t][1])

Res[t] = int(values[t][2])

print('Number of data = ', N)

fout.write('\nNumber of data = {}\n'.format(N))

for i in range(N):

print('{0} {1} {2}'.format(ID[i], X[i], Res[i]))

for i, s, r in zip(ID, X, Res):

fout.write('{0:>5}{1:>7.2f}{2:>5d}\n'.format(i, s, r))

#

#

# Values for rescaling in MCMC sampling

#

min_sgm = 0.0 # lower boundary

max_sgm = np.amax(X) - np.amin(X) # upper boundary

C1C2 = np.percentile(X, [20, 80])

def f_init():

"""

Initial values

"""

return dict(v_pse = 0.0, v_sgm = 1.0)

stan_data = {'N': N, 'X': X, 'Res': Res, 'min_sgm': min_sgm, 'max_sgm': max_sgm,

'p75': ss.norm.ppf(0.75, 0, 1) }

#

# Sampling by Stan

#

with open('PM2CatPS3.stan', 'r') as f_stan:

sm = stan.build(f_stan.read(), stan_data)

fit = sm.sample(num_chains = 4, init = [f_init()]*4)

dframe = fit.to_frame()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

print(az.summary(i_data))

fout.write('\nSummary:\n')

fout.write(f'{az.summary(i_data)}\n')

az.plot_trace(i_data, var_names = ['pse', 'jnd'])

plt.tight_layout()

plt.show()

ck = az.plot_posterior(i_data, var_names = 'pse', point_estimate = 'median')

plt.show()

print(ck)

ck = az.plot_posterior(i_data, var_names = 'jnd', point_estimate = 'median')

plt.show()

#

# MAP estimates for mu, sigma and JND

#

#

muMAP = CalcMAPEst(dframe['pse'])[0]

sgmMAP = CalcMAPEst(dframe['sgm'])[0]

jndMAP = CalcMAPEst(dframe['jnd'])[0]

SmplSgm = dframe['sgm'] # Samples from the posterior distribution for sigma

SmplPSE = dframe['pse'] # Samples from the posterior distribution for pse(mu)

#

# Percentile points for posterior distributions

#

SgmL025 = np.percentile(SmplSgm, 2.5)

SgmU975 = np.percentile(SmplSgm, 97.5)

PSEL025 = np.percentile(SmplPSE, 2.5)

PSEU975 = np.percentile(SmplPSE, 97.5)

#

# Prob(z < s75) = 0.75 for the standard normal distribution

#

s75 = ss.norm.ppf(0.75, 0, 1)

JND = jndMAP

PSE = muMAP

fout.write('\nPSE(MAP) = {0:<.3f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(

PSE, PSEL025, PSEU975))

fout.write('\nJND(MAP) = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(

JND, s75 * SgmL025, s75 * SgmU975))

#

# Data and the psychometric function

#

pf = PF(muMAP, sgmMAP) # Class PF from my_module.py

minX = X[np.argmin(X)]

maxX = X[np.argmax(X)]

xvalues = []

yvalues = []

for i in range(1001):

t = minX + i * (maxX - minX) / 1000.0

xvalues.append(t)

yvalues.append(pf.v(t)) # Values of the psychometric function

plt.title('Psychometric Function\n' +

('PSE = {0:<.2f}, JND = {1:<.2f}'.format(PSE, JND)),

fontsize = 16)

plt.plot(xvalues, yvalues, linewidth = 3, label = 'PF')

plt.yticks([0.0, 0.5, 0.75, 1.0])

#

# Resposes are scaled to values from 0 to 1, and added small fluctuations

#

Res_rn = (Res - 1) + np.random.uniform(-0.05, 0.05, len(Res))

plt.scatter(X, Res_rn, alpha = 0.5, color = 'g', label = 'Data')

plt.plot([minX, maxX], [0.0, 0.0], color = 'k')

plt.plot([minX, maxX], [1.0, 1.0], color = 'k')

plt.plot([minX, PSE, PSE], [0.5, 0.5, 0.0], color = 'k')

plt.plot([minX, PSE + JND, PSE + JND], [0.75, 0.75, 0.0], color = 'k')

plt.legend(loc = 7)

plt.xlabel('Stimulus value', fontsize = 12)

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

plt.show()

#

# Observed responses smoothed by Kernel density estimation

#

w_1 = np.empty(len(Res)) # Weights for Response = 2

for t in range(len(Res)):

w_1[t] = 0 if Res[t] < 2 else 1

X_range = np.arange(minX, maxX, (maxX - minX) / 100)

w_k = JND / 2

while True:

pm_1 = [] # KDE values for w_1

for v in X_range:

pm_1.append(k_ave(v, X, w_1, w_k))

plt.title('Observed and Estimated PFs\n' +

('PSE = {0:<.2f}, JND = {1:<.2f}'.format(PSE, JND)),

fontsize = 16)

plt.plot(xvalues, yvalues, color = 'b', linewidth = 3, label = 'PF')

plt.yticks([0.0, 0.5, 1.0])

plt.plot(X_range, pm_1, 'g--', linewidth = 2, label = 'Res(>1)')

plt.xlabel('Stimulus value', fontsize = 12)

plt.legend(loc = 'upper left')

plt.show()

print('\nThe current value of w_k = {}'.format(w_k))

w_k = float(input('The next value of w_k = '))

if w_k <= 0.0:

break

fout.close()

print('\nOutput file {} was saved.'.format(fout_nm))

リスト２で用いられているモジュールファイルmy_module.pyをリスト３に示す。これらは、圧縮ファイルfiles2CatPS3.zipにまとめられている。

リスト３　モジュールファイルmy_module.py

import numpy as np

import scipy.stats as ss

def CalcMAPEst(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

"""

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

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

est_pdf = ss.gaussian_kde(samples).pdf(coord)

map_idx = np.argmax(est_pdf)

MAP_Est = coord[map_idx]

return MAP_Est, est_pdf[map_idx]

f_phi = ss.norm.cdf # Cumulative standard normal distribution function

class PF:

"""

Psychometric function with the parameters mu and sigma

"""

def __init__(self, mu, sgm):

self.mu = mu

self.sgm = sgm

def v(self, x):

return f_phi((x - self.mu) / self.sgm)

def my_kernel(x, y, w_k):

return np.exp(-((x - y) / w_k) ** 2)

def k_ave(v, X, W, w_k):

n = len(X)

num = 0.0

den = 0.0

for i in range(n):

k_v = my_kernel(v, X[i], w_k)

num += W[i] * k_v

den += k_v

return num / den

Up