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

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

PythonとPyStanの使い方については、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」、丸善出版＞で説明している。

PythonもPyStanも無料のソフトである。

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

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

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

のとき,

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

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

ここで

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

上のモデルにおいて、

と置くと、

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

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

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

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

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

リスト１　３件法上下法データのベイズ分析用Stanスクリプト（PM3Cat_EqCPS3.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 d_mid; // Midpoint of the stimuli

real d_sd; // Rough dispersion measure of the stimuli

}

parameters {

real<lower = 0.0, upper = 1.0> v_sgm;

real v_C1;

real<lower = 0.0> v_C2;

}

transformed parameters {

real<lower = 0.0> sgm;

real C1;

real C2;

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

sgm = v_sgm * d_sd * 10.0;

C1 = v_C1 * d_sd * 10.0 + d_mid;

C2 = C1 + v_C2 * d_sd * 10.0;

for (i in 1:N){

p[i][1] = 1 - normal_cdf((X[i] - C1) / sgm | 0, 1);

p[i][2] = normal_cdf((X[i] - C1) / sgm | 0, 1)

- normal_cdf((X[i] - C2) / sgm | 0, 1);

p[i][3] = normal_cdf((X[i] - C2) / sgm | 0, 1);

}

model {

v_sgm ~ exponential(0.1);

v_C1 ~ normal(0.0, 10.0);

v_C2 ~ exponential(0.1);

for (i in 1:N)

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

}

generated quantities {

real pse;

real jnd;

pse = (C1 + C2) / 2;

jnd = 0.6745 * sgm;

}

リスト１のStanスクリプトによる分析プログラム例を付記Aのリスト２に示す。プログラムファイルとサンプルデータファイルは圧縮ファイルfiles3CatPS3.zipとしてまとめた。

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

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

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

Data file name (*.xlsx) =

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

・

図１

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

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

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

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

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

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

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

図２

図２を閉じると、パラメータpseの事後分布が表示される（図３）。

図３

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

図３のフォームを閉じると、次にjndのグラフが表示される（図４）。

図４

図４のフォームを閉じると、MAP推定値を求めるためのしばらくの計算時間の後にPFとデータ値のグラフがMAP推定値に基づくPSEとJNDの推定値とともに表示される（図５）。

図５

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

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

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

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

図６

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

The current value of w_k = 4.723328567390731

The nexr value of w_k =

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

The current value of w_k = 4.723328567390731

The nexr value of w_k = 9

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

図７

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

The current value of w_k = 4.723328567390731

The nexr value of w_k = 9

The current value of w_k = 9.0

The nexr value of w_k = 0

Results.txt was saved.

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

表示されている出力ファイル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

Summary:

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat

v_sgm 0.099 0.016 0.072 0.128 0.000 0.000 1429.0 1834.0 1.0

v_C1 -0.092 0.018 -0.127 -0.060 0.000 0.000 1581.0 1960.0 1.0

v_C2 0.135 0.023 0.094 0.179 0.001 0.000 1351.0 2018.0 1.0

sgm 14.792 2.363 10.829 19.193 0.063 0.045 1429.0 1834.0 1.0

C1 193.667 2.676 188.517 198.448 0.068 0.048 1581.0 1960.0 1.0

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

p[99, 0] 0.028 0.017 0.003 0.060 0.000 0.000 2099.0 2221.0 1.0

p[99, 1] 0.251 0.058 0.138 0.357 0.001 0.001 2088.0 2779.0 1.0

p[99, 2] 0.722 0.068 0.592 0.846 0.002 0.001 1967.0 2441.0 1.0

pse 203.803 1.880 200.051 207.048 0.033 0.023 3239.0 3053.0 1.0

jnd 9.977 1.594 7.304 12.946 0.043 0.030 1429.0 1834.0 1.0

[308 rows x 9 columns]

MAP for PSE = 204.03, 95% CI for PSE = [200.04, 207.39]

MAP for JND = 9.45, 95% CI for JND= [7.40, 13.52]

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

引用文献

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

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

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

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

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

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

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

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

付記A　Pythonスクリプト

リスト１のStanスクリプトを利用して３件法による上下法のデータを分析するPythonスクリプトを以下に示す。これらのファイルは圧縮ファイルfiles3CatPS3.zipにまとめたので、ダウンロードして自由に利用されたい。

リスト２　リスト１のStanスクリプトによる上下法データのベイズ分析Pythonスクリプト（MainThreeCatDiscrimPS3.py）。モジュールmy_module.pyはリスト３に示す。

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import stan

import seaborn as sb

import pickle

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

C1, X_mid, C2 = np.percentile(X, [20, 50, 80])

print('C1,C2 = ', C1, C2)

X_sd = (C2 - C1) / 2

print('X_sd =', X_sd, ' X_mid =', X_mid)

stan_data = {'N': N, 'X': X, 'Res': Res, 'd_mid': X_mid, 'd_sd' : X_sd}

# Sampling by Stan

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

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

fit = sm.sample()

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', 'C1', 'C2'])

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 C

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

"""

import numpy as np

import scipy.stats as ss

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]

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

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

C1MAP = CalcMAPEst(dframe['C1'])[0]

C2MAP = CalcMAPEst(dframe['C2'])[0]

SgmL025 = np.percentile(dframe['sgm'], 2.5)

SgmU975 = np.percentile(dframe['sgm'], 97.5)

MuL025 = np.percentile(dframe['pse'], 2.5)

MuU975 = np.percentile(dframe['pse'], 97.5)

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

s75 = 0.6745

JND = sgmMAP * s75

PSE = muMAP

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

PSE, MuL025, MuU975))

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

JND, s75 * SgmL025, s75 * SgmU975))

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

PSE, MuL025, MuU975))

fout.write('\nMAP for JND = {0:<.2f}, 95% CI for JND= [{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 = []

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 + (C2MAP-C1MAP)/2))

yvalues_R.append(pf.v(t - (C2MAP-C1MAP)/2))

plt.title('Psychometric Functions\n' +

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

(PSE, JND, (C2MAP-C1MAP)/2)),

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

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, (C2MAP-C1MAP)/2)),

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(f'\n{fout_nm} was saved.\n')

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

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

import numpy as np

import scipy.optimize

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

def my_log(x):

"""

To avoid the error message with 0 values for log function

"""

if x < 1.0e-323:

return -744.0

else:

return np.log(x)

class LLog:

"""

Negative log likelihood

"""

def __init__(self, X, Res):

self.N = len(X)

self.X = X

self.Res = Res

def v(self, z):

mu = z[0]

sgm = z[1] ** 2 # Squared to avoid negative values

C = z[2] ** 2 # Squared to avoid negative values

v = 0

for t in range(self.N):

if self.Res[t] == 3:

v += -my_log(f_phi((self.X[t] - mu - C) / sgm))

elif self.Res[t] == 2:

v += -my_log(f_phi((self.X[t] - mu + C) / sgm) -

f_phi((self.X[t] - mu - C) / sgm))

else:

v += -my_log(1.0 - f_phi((self.X[t] - mu + C) / sgm))

return v

def calcMAP(X, Res, PSEMed, SgmMed, CMed):

"""

Calculate MAP estimates

"""

llog = LLog(X, Res) # Negative log likelihood

z = np.empty(3)

z[0] = PSEMed # Initial value for mu

z[1] = SgmMed ** 0.5 # Initial values for square root of sigma

z[2] = CMed ** 0.5 # Initial values for square root of C

res = scipy.optimize.minimize(llog.v, z, method = 'Nelder-Mead',

options = {'xatol': 1.0e-7, 'fatol': 1.0e-12})

muMAP = res.x[0]

sgmMAP = res.x[1] ** 2

CMAP = res.x[2] ** 2

return muMAP, sgmMAP, CMAP

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

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

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

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

引用文献

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

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

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