３件法上下法のベイズ分析（CmdStanPy対応）

３件法による上下法（階段法）

CmdStanPy対応

サンプルプログラムとデータファイルは、圧縮ファイルupdown3catfiles.zipにまとめた。ダウンロード・解凍すれば、自由に利用することができる。

心理物理学（精神物理学、psychophysics）において主観的等価点を求める方法として上下法（Up-down method）あるいは階段法（Staircase method）と呼ばれている方法がある。上下法では、提示刺激値が直前の試行における判断に依存して決まるため、被験者が提示刺激値を予測する可能性が高い。これを避けるため、２系列用意して２つの系列の試行をランダムに混ぜて実施する（二重上下法；大山、2005）。２つの系列を系列１と系列２として、系列でのものであることを右上に添え字として明記する。なお、３件法については、全系列法においてその問題点が指摘されたことがあるが、適切な数理モデル（確率モデル）を用いればそのような問題は無い。このことについては付録　「３件法に関わる事柄など」において論じた。

３件法の場合、判断は「強い」、「同じ／わからない」、「弱い」の３通りである。系列の第試行における判断（反応）を以下のように表す。

系列の第試行における提示刺激値をで表すとき、次試行における提示刺激値を判断に応じて以下のように決める（岡本、2019）。

のとき,

判断がのときに前試行における刺激値との関係で次試行の刺激値が決められているのは、「強い」判断および「弱い」判断のデータが多い方が心理測定関数（psychometric function: PF）の傾きの推定が安定するからである。

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

ここで

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

上のモデルにおいて、

と置くと、

となる（岡本、2017）。このとき、主観的等価点PSE は

で（岡本、1994；大山、1971）、弁別閾JNDは

で（岡本、2014、p. 169）与えられる。

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

この尤度関数に対して、事後分布を求めるStanスクリプトをリスト１に示す。の値は、－、？、＋を１、２、３で表している。ファイルは圧縮ファイルupdown3catfiles.zipにまとめたので、ダウンロード展開して自由に利用できる。

リスト１のStanスクリプトによる分析プログラム例をリスト２に示す。

リスト１のスクリプトファイル（PM3Cat_EqC.stan）、リスト２のスクリプトファイル（main.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。

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

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

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

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

Output file name = Results.txt

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

・

図１

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

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

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

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

図２

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

図３

すなわちPSEの事後分布のグラフである。点推定値として中央値（median）が示されている。Gelmanら（2021）は、平均値と中央値を比べて、中央値の方が計算においてより安定しているとして中央値を勧めている。

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

図４

の事後分布のグラフである。

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

図５

の事後分布のグラフである。

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

図６

PFとデータ値のグラフが、事後分布の中央値を点推定値とするPSEとJNDの値とともに表示されている。

PFは、以下の式に基づくものである。

データ値は、１から３までの値であるが、これを０から１までの範囲に縮小して表示している。すなわち、

である。さらに、上下に少しランダムな変動を加えて、完全に重なって表示されることを防いでいる。

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

図７

図７に表示されているPF（実線）と、対応するデータ値（評定Rが２以上および３以上の比率）のKDE（Kernel Density Estimation）による平滑化されたグラフ（破線）が表示されている。図７のフォームを閉じると、次のようにw_kの値の現在値が表示され、次の値の設定が求められる。

The current value of w_k = 4.858530096670557

The nexr value of w_k =

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

The current value of w_k = 4.858530096670557

The nexr value of w_k = 9

いま、上のようにw_k=9を設定すると、図８のようになる。

図８

w_kの値として０以下の値（次の例では０）を設定すると、プログラムは終了する。

The current value of w_k = 4.858530096670557

The nexr value of w_k = 9

The current value of w_k = 9.0

The nexr value of w_k = 0

Output file Results.txt was saved.

(stan) *****/updown3catfiles$

表示されている出力ファイルResults.txtを開くと、以下のようになっている。

Data file...Data3Cat.xlsx

Number of data = 100

0 257.50 3

1 237.50 3

2 222.50 3

　　　　・

97 192.50 1

98 207.50 1

99 222.50 3

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -70.066500 0.032344 1.291820 ... 2173.76 1549.01 1.00017

sgm 14.742600 0.060024 2.458680 ... 2050.74 1630.66 1.00197

C1 193.712000 0.065285 2.684840 ... 1555.76 1602.00 1.00118

w_C 20.149900 0.085995 3.425470 ... 1830.14 1481.14 1.00029

C2 213.862000 0.041800 2.506100 ... 2964.05 3315.75 1.00085

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

p[100,1] 0.027728 0.000389 0.018354 ... 2226.33 2247.63 1.00218

p[100,2] 0.247939 0.001179 0.059961 ... 2812.64 2315.68 1.00111

p[100,3] 0.724334 0.001437 0.070558 ... 2766.94 2192.44 1.00171

mu 203.787000 0.034438 1.952170 ... 2918.90 2900.90 1.00072

C 10.075000 0.042998 1.712740 ... 1830.14 1481.14 1.00029

[307 rows x 11 columns]

Med for mu = 203.80 95% CI for mu = [199.84, 207.56]

Med for sigma = 14.41 95% CI for sigma = [10.95, 20.52]

Med for C = 9.92 95% CI for C = [7.18, 13.90]

Med for PSE = 203.80, 95% CI for PSE = [199.84, 207.56]

Med for JND = 9.72, 95% CI for JND= [7.39, 13.84]

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

附記　３件法に関わる事柄など

３件法の用いられた全系列法において、単純な四則演算で分析が行われると感覚と判断が分離できない（岡本、1995）という問題があり、３件法の使用が控えられてきたという経緯があると思う。しかし、３件法においても感覚と判断を区別する適切な数理モデル（確率モデル）により、感覚と判断を分離した分析が行えることが報告されている（岡本、1994）。

刺激値の間隔については、Ψ法による適応的恒常法において刺激値の間隔の影響を調べたものが報告されている（岡本、2006）。

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

data {

int N; // The number of trials

array[N] int Res; // Responses:

1 denotes -, 2 denotes ?, 3 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 C1;

real<lower = 0.0> w_C;

}

transformed parameters {

real C2;

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

C2 = C1 + w_C;

for (i in 1:N){

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

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

p[i][3] = Phi((X[i] - C2) / sgm);

}

model {

sgm ~ uniform(min_sgm, max_sgm);

for (i in 1:N)

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

}

generated quantities {

real mu;

real C;

mu = (C1 + C2) / 2;

C = C2 - mu;

}

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

from cmdstanpy import CmdStanModel

import pandas as pd

import arviz as az

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

import scipy

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

fout_nm = input('Output file name = ')

fout = open(fout_nm, 'w')

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

xlsx_fl = pd.ExcelFile(fin_nm)

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

data = data_xlsx.values

print('data...\n', 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\n'.format(N))

for i in range(N):

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

fout.write('{0:>5} {1:>10.2f} {2:>5d}\n'.format(ID[i], X[i], Res[i]))

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

print('C1C2 = ', C1C2)

def f_init():

"""

Initial values for C, mu and sigma

"""

return dict(sgm = C1C2[1] - C1C2[0],

C1 = C1C2[0],

w_C = C1C2[1] - C1C2[0])

min_sgm = 0.0 # lower boundary for sigma

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

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

# コンパイル

model = CmdStanModel(stan_file="PM3Cat_EqC.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', 'C1', 'C2'])

plt.tight_layout()

plt.show()

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

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

SmplC = fit['C'] # Samples from the posterior distribution for C

SmplMu = fit['mu'] # Samples from the posterior distribution for mu

# Percentile points for posterior distributions

SgmL025 = np.percentile(SmplSgm, 2.5)

SgmMed = np.percentile(SmplSgm, 50)

SgmU975 = np.percentile(SmplSgm, 97.5)

CL025 = np.percentile(SmplC, 2.5)

CMed = np.percentile(SmplC, 50)

CU975 = np.percentile(SmplC, 97.5)

C1Med = np.percentile(fit['C1'], 50)

C2Med = np.percentile(fit['C2'], 50)

MuL025 = np.percentile(SmplMu, 2.5)

MuMed = np.percentile(SmplMu, 50)

MuU975 = np.percentile(SmplMu, 97.5)

print('muMed = {0} sgmMed = {1} CMed = {2}'.format(MuMed, SgmMed, CMed))

fout.write('\nMed for mu = {0:<.2f} 95% CI for mu = [{1:<.2f}, {2:<.2f}]\n'.

format(MuMed, MuL025, MuU975))

# KDE plot for the samples from posterior distribution for mu

sb.kdeplot(SmplMu)

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

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

plt.yticks([])

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

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

('Med = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(MuMed, MuL025, MuU975)),

fontsize = 16)

plt.legend()

plt.tight_layout(h_pad = 5)

plt.show()

# KDE plot for the samples from posterior distribution for sigma

fout.write('\nMed for sigma = {0:<.2f} 95% CI for sigma = [{1:<.2f}, {2:<.2f}]\n'.

format(SgmMed, SgmL025, SgmU975))

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

# KDE plot for the samples from posterior distribution for C

fout.write('\nMed for C = {0:<.2f} 95% CI for C = [{1:<.2f}, {2:<.2f}]\n'.

format(CMed, CL025, CU975))

sb.kdeplot(SmplC)

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

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

plt.yticks([])

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

plt.title('Posterior Distribution for C\n' +

('Med = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(CMed, CL025, CU975)),

fontsize = 16)

plt.legend()

plt.tight_layout(h_pad = 5)

plt.show()

f_phi = scipy.stats.norm.cdf # Cumulative standard normal distribution function

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)

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

s75 = bisect(f_phi, 0.75, 0.5, 1.0)

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

JND = SgmMed * s75

PSE = MuMed

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

PSE, MuL025, MuU975))

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

JND, s75 * SgmL025, s75 * SgmU975))

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)

# Data and the psychometric function

pf = PF(MuMed, SgmMed)

minX = X[np.argmin(X)]

maxX = X[np.argmax(X)]

xvalues = []

yvalues = []

yvalues_L = []

yvalues_R = []

for i in range(1001):

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

xvalues.append(t)

yvalues.append(pf.v(t))

yvalues_L.append(pf.v(t + CMed))

yvalues_R.append(pf.v(t - CMed))

plt.title('Psychometric Functions\n' +

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

fontsize = 16)

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

plt.plot(xvalues, yvalues_L, 'g-', linewidth = 3, label = 'PF(>1)')

plt.plot(xvalues, yvalues_R, 'b-', linewidth = 3, label = 'PF(>2)')

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

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

Res_rn = (Res - 1) / 2 + 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

w_2 = np.empty(len(Res)) # Weights for response >= 3

for t in range(len(Res)):

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

w_2[t] = 0 if Res[t] < 3 else 1

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

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

w_k = JND / 2

while True:

pm_1 = [] # KDE values for Response >= 2

pm_2 = [] # KDE values for Response >= 3

for v in X_range:

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

pm_2.append(k_ave(v, X, w_2, w_k))

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

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

fontsize = 16)

plt.plot(xvalues, yvalues_L, 'g-', linewidth = 3, label = 'PF(>1)')

plt.plot(xvalues, yvalues_R, 'b-', linewidth = 3, label = 'PF(>2)')

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

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

plt.plot(X_range, pm_2, 'b--', linewidth = 2, label = 'Res(>2)')

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 nexr value of w_k = '))

if w_k <= 0.0:

break

fout.close()

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

引用文献

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

岡本安晴（1994）全系列法データの最尤法による分析例．金沢大学文学部論集. 行動科学科篇、14、63ー71．

岡本安晴（1995）全系列法データの最尤法による分析について : シミュレーションによる検討．金沢大学文学部論集. 行動科学科篇、15、25ー30．

岡本安晴（2006）適応的恒常法の可能性について．日本心理学会第70回大会発表．

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

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

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

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

大山　正（1971）．感覚・知覚測定法　和田陽平・大山　正・今井省吾（編著）感覚・知覚心理学ハンドブック（pp. 32–55）誠信書房

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

Home