Up

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

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

とおくと、

のとき,　

のとき,　,

である。

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

なお、本ウェブサイトでは、２件法の場合を扱っているが、３件法、４件法の方が精度の高い結果が得られるというシミュレーションが報告されている（岡本、2019、2025）。

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

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

とおくと、

のとき,　

のとき,　,

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

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

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

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

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

リスト１のスクリプトファイル（PM2Cat.stan）、リスト２のスクリプトファイル（main.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移動して、CmdStanPyのインストールされている環境において、次のコマンドを実行する。CmdStanPyの準備の簡単な説明を、このウェブサイトで行った。

(stan) *****/updown2catfiles$ python main.py

上のコマンドを実行すると、入力データファイル名と出力データファイル名の設定が求められる。

(stan) *****/updown2catfiles$ python main.py

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

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

入力データファイルは、図１の形式のExcelファイルとして用意する。

・

・

・

図１

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

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

入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。出力ファイル名は、任意のテキストファイル名（*.txt）である。

入力ファイル名と出力ファイル名の設定が終わると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。このコンパイルは、多少時間が掛かる。コンパイルが終わると、MCMCサンプリングが行われ、サンプリングが終了すると、MCMCサンプリングのトレース図が表示される（図２）。

図２

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

図３

点推定値として中央値が示されている。Gelmanら（2021）は、平均値と中央値を比べて、中央値の方が計算においてより安定しているとして中央値を勧めている（p. 73）。

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

図４

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

こ

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

を満たす値である。

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

図５

データ値は、次式

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

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

図６

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

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

The current value of w_k = 5.902400280245459

The next value of w_k =

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

The current value of w_k = 5.902400280245459

The next value of w_k = 9

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

図７

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

The current value of w_k = 5.902400280245459

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.

(stan) *****/updown2catfiles$

保存された出力ファイル（上の場合、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

Mean MCSE StdDev MAD ... ESS_bulk ESS_tail ESS_bulk/s R_hat

lp__ -48.864300 0.029578 1.064520 0.799333 ... 1311.68 1832.56 2474.87 1.00151

sgm 18.455100 0.112635 4.676550 3.948490 ... 1884.45 2133.41 3555.56 1.00206

pse 203.109000 0.055844 2.643130 2.518430 ... 2330.67 1874.94 4397.49 1.00196

p[1,1] 0.010657 0.000453 0.017501 0.005055 ... 1841.22 1859.70 3474.01 1.00238

p[1,2] 0.989343 0.000453 0.017501 0.005055 ... 1841.09 1859.70 3473.75 1.00238

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

p[98,2] 0.429719 0.001078 0.057183 0.058273 ... 2839.27 2499.97 5357.12 1.00180

p[99,1] 0.252625 0.001465 0.064001 0.064754 ... 2020.89 1583.46 3813.01 1.00221

p[99,2] 0.747375 0.001465 0.064001 0.064754 ... 2020.89 1583.46 3813.01 1.00221

p[100,1] 0.074306 0.001258 0.050087 0.044498 ... 1843.12 1578.06 3477.59 1.00257

p[100,2] 0.925694 0.001258 0.050087 0.044498 ... 1843.12 1578.06 3477.59 1.00257

[203 rows x 11 columns]

PSE(Med) = 203.098 Sgm(Med) = 17.502

PSE = 203.098, 95% CI = [197.87, 208.42]

JND = 11.80, 95% CI = [8.15, 20.47]

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

リスト１　上下法データのベイズ分析Stanスクリプト（PM2Cat.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

}

parameters {

real<lower = 0.0> sgm;

real pse;

}

transformed parameters {

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

for (i in 1:N){

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

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

}

}

model {

sgm ~ uniform(min_sgm, max_sgm);

for (i in 1:N)

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

}

リスト２　リスト１のStanスクリプトを用いるベイズ分析のPythonスクリプト（main.py）

from cmdstanpy import CmdStanModel

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import scipy

import arviz as az

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

fout_nm = input('\nOutput file name (*.txt) = ')

fout = open(fout_nm, 'w')

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

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

data = data_xlsx.values

#print(data)

N = len(data) - 1

ID = []

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

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

for t in range(N):

ID.append(data[t + 1][0])

X[t] = float(data[t + 1][1])

Res[t] = int(data[t + 1][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))

min_sgm = 0.0 # lower boundary for sigma

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

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

def f_init():

"""

Initial values for pse(mu) and sigma

"""

return dict(pse = (C1C2[0] + C1C2[1]) / 2, sgm = C1C2[1] - C1C2[0])

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

# コンパイル

model = CmdStanModel(stan_file="PM2Cat.stan")

# サンプリング

fit = model.sample(data=stan_data, chains=4, inits=[f_init()]*4) # サンプリング

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

i_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

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

plt.tight_layout()

plt.show()

fit = fit.draws_pd() # 　　Pandas DataFrame型への変換

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

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

#

# Percentile points for posterior distributions

#

SgmL025 = np.percentile(SmplSgm, 2.5)

SgmMed = np.percentile(SmplSgm, 50)

SgmU975 = np.percentile(SmplSgm, 97.5)

PSEL025 = np.percentile(SmplPSE, 2.5)

PSEMed = np.percentile(SmplPSE, 50)

PSEU975 = np.percentile(SmplPSE, 97.5)

print('PSE = ', PSEMed, ' Sgm = ', SgmMed)

fout.write('\nPSE(Med) = {0:<.3f} Sgm(Med) = {1:<.3f}\n'.

format(PSEMed, SgmMed))

#

# KDE plot for the samples from posterior distribution for PSE(mu)

#

sb.kdeplot(SmplPSE)

plt.plot([PSEL025, PSEU975], [0.0, 0.0], color = 'b', linewidth = 5, label = '95% CI')

plt.plot(PSEMed, 0.0, 'go', markersize = 10, label = 'Med')

plt.yticks([])

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

plt.title(r'Posterior Distribution for PSE($\mu$)' + '\n' +

('Med = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.

format(PSEMed, PSEL025, PSEU975)),

fontsize = 16)

plt.legend()

plt.tight_layout(h_pad = 5)

plt.show()

#

# KDE plot for the samples from posterior distribution for sigma

#

sb.kdeplot(SmplSgm)

plt.plot([SgmL025, SgmU975], [0.0, 0.0], color = 'b', linewidth = 5, label = '95% CI')

plt.plot(SgmMed, 0.0, 'go', markersize = 10, label = 'Med')

plt.yticks([])

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

plt.title(r'Posterior Distribution for $\sigma$' + '\n' +

('Med = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.

format(SgmMed, SgmL025, SgmU975)),

fontsize = 16)

plt.legend()

plt.tight_layout(h_pad = 5)

plt.show()

f_phi = scipy.stats.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 bisect( f, c, Lb, Ub ):

"""

Return the value root so that f(root) = c

"""

s = 1.0 if f(Ub) > f(Lb) else -1.0 # s = 1 when f(x) is an increasing functin

while ((abs(Ub - Lb) > (1.0e-15) * (abs(Ub) + abs(Lb)))

and

(abs(Ub) + abs(Lb) > 1.0e-14)):

m = 0.5 * (Ub + Lb)

if (s * (f(m) - c)) > 0.0:

Ub = m

else:

Lb = m

return 0.5 * (Ub + Lb)

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

#

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

#

s75 = bisect(f_phi, 0.75, 0.5, 1.0) # f_phi from my_module.py

print('s75 = ', s75, ' p = ', f_phi(s75))

JND = SgmMed * s75

PSE = PSEMed

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

PSE, PSEL025, PSEU975))

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

JND, s75 * SgmL025, s75 * SgmU975))

#

# Data and the psychometric function

#

pf = PF(PSEMed, SgmMed) # Class PF defined above

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 Response >= 2

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

参考文献

相場　覚・鳥居修晃（編著）（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．

岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房

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

Home