３件法による上下法（階段法）
2023.10.22更新
PyStanがWindowsで使い易くなった（2020.07）が、これに関する説明を
PyStan（Python用Stan）インストールの説明のウェブサイトを用意した。
サンプルプログラムとデータファイルは、圧縮ファイルfiles3Cat.zipにまとめた。ダウンロード・解凍すれば、自由に利用することができる。
PythonとPyStanの使い方については、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」、丸善出版＞で説明している。
PythonもPyStanも無料のソフトである。
心理物理学（精神物理学、psychophysics）において主観的等価点を求める方法として上下法（Up-down method）あるいは階段法（Staircase method）と呼ばれている方法がある。上下法では、提示刺激値が直前の試行における判断に依存して決まるため、被験者が提示刺激値を予測する可能性が高い。これを避けるため、２系列用意して２つの系列の試行をランダムに混ぜて実施する（二重上下法；大山、2005）。２つの系列を系列１と系列２として、系列
でのものであることを右上に添え字として明記する。
３件法の場合、判断は「強い」、「同じ／わからない」、「弱い」の３通りである。系列の第
試行における判断（反応）
を以下のように表す。
系列の第試行における提示刺激値を
で表すとき、次試行における提示刺激値を判断に応じて以下のように決める（岡本、2019）。
のとき,
,
のとき,
,
のとき,
判断がのときに前試行における刺激値との関係で次試行の刺激値が決められているのは、「強い」判断および「弱い」判断のデータが多い方が心理測定関数（psychometric function: PF）の傾きの推定が安定するからである。
データ分析においては、系列１のデータ
と系列２のデータ
をまとめて
とおき、刺激に対する反応の確率を次式で与える（岡本、2019）。
ここで
であり、
は累積標準正規分布関数である。このとき、尤度関数を次式で与える。
上のモデルにおいて、
と置くと、
となる（岡本、2017）。このとき、主観的等価点PSE は
で（岡本、2014、p. 169）与える。
この尤度関数は、データが上下法以外の方法、例えば恒常法、で集められた場合にも適用できることに注意。データは、「強い」判断および「弱い」判断がそれぞれ確実な刺激値が十分に強いものと弱いものを含むと、
の事後分布が拡散することを防ぐことができる。
この尤度関数に対して、事前分布として一様分布を設定したときの事後分布を求めるStanスクリプトをリスト１に示す。
の値は、−、？、＋を１、２、３で表している。ファイルはfiles3Cat.zip圧縮ファイルにまとめたので、ダウンロード展開して自由に利用できる。
リスト１ ３件法上下法データのベイズ分析用Stanスクリプト（PM3Cat_EqC.stan）
data {
int N; // The number of
trials
int Res[N]; // Responses:
// 1
denotes -, 2 denotes ?, 3 denotes +
real X[N]; // 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;
vector[3] p[N]; // 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スクリプトによる分析プログラム例を付記のリスト２に示す。プログラムファイルとサンプルデータファイルは圧縮ファイルfiles3Cat.zipとしてまとめた。
リスト２のスクリプトをAnaconda/WindowsのPython 3.7仮想環境にインストールされたPyStan2.17.1.0により実行した。Anacondaのインストールは簡単であるが、このウェブサイトを参照されたい。
リスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py37) PS D:\＊＊＊\files3Cat> python
MainThreeCatDiscrim.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
Data file name (*.csv) =
データファイルは、CSV形式のファイルとして用意するが（Excelで保存時に、ファイル拡張子として.csvを指定すればよい）、形式は図１のようである。
・
・
・
図１
第１行目に、変数名を適当に設定する。２行目からデータを設定するが、１列目はデータの識別用である。図２では通し番号が設定されているが、任意の文字列でよい。２列目に提示刺激値
を設定し、３列目に反応を設定する。反応は、「弱い」判断のときは１、「同じ／わからない」判断のときは２、「強い」判断のときは３を設定する。データの行の順序は任意である。これは、判断の確率が刺激値のみによって与えられるので、例えば、図２の２行目にある０番目のデータと101行目の99番目のデータを入れ替えても尤度関数値は変わらないからである。
データの設定が終われば、ファイル拡張子として.csvを選んで保存すればCSV形式で保存される。CSV形式のファイルについては、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、pp. 33-35,
101-102＞で説明している。パラメータの事後分布が発散しないように、提示刺激には「強い判断」の確率が１に近い十分強い刺激と「弱い」判断の確率が０に近い十分弱い刺激が含まれるようにする。２系列用意する上下法において、「強い」判断が容易な強い刺激値で始まる系列と「弱い」判断が容易な弱い刺激値で始まる系列の２系列を用意すればよい。
入力データファイル名を設定すると、次に出力ファイル名の設定が求められる（図３）。適当なテキストファイル名を設定する。
(py37) PS D:\＊＊＊\files3Cat> python
MainThreeCatDiscrim.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
Data file name (*.csv) = Data3Cat.csv
Output file name = Results.txt
入力ファイル名と出力ファイル名の設定が終わると、データが読み込まれ、Stanによる事後分布からのサンプリングが始まる。リスト２のプログラムでは、以下のように初期化の関数f_initとデータの辞書stan_dataが設定されて、リスト１のStanスクリプトPM3Cat_EqC.stanによるサンプリングが行われる。PythonにおけるStan（PyStan）の使用については、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、pp. 178-208＞で説明している。
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])
#
# Boundaries
and initial values
#
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}
#
# Sampling
by Stan
#
sm = pystan.StanModel(file =
'PM3Cat_EqC.stan')
fit = sm.sampling(data = stan_data,
pars = ['C', 'C1', 'C2', 'mu', 'sgm'],
init = f_init,
n_jobs = 1)
サンプリングが終了すると、パラメータmu(
)の事後分布のグラフがMAP推定値、95％ CIとともに表示される（図２）。グラフは、保存ボタンをクリックすると保存できる。点推定値は、MAP推定値が良いという報告がある（岡本、2022）。
図２
フォーム右上のX印をクリックしてフォームを閉じると、次に、sigma（）の事後分布のグラフが、MAP推定値と95％ CIとともに表示される（図3
）。
図３
図３のフォームを閉じると、次にCの事後分布のグラフがMAP推定値と95％ CIとともに表示される（図４）。
図４
図４のフォームを閉じると、次にPFとデータ値のグラフがMAP推定値に基づくPSEとJNDの推定値とともに表示される（図５）。
図５
PFは、以下の式に基づくものである。
データ値は、１から３までの値であるが、これを０から１までの範囲に縮小して表示している。すなわち、
である。さらに、上下に少しランダムな変動を加えて、完全に重なって表示されることを防いでいる。
図５のフォームを閉じると図６のフォームが表示される。
図６
図６に表示されているPF（実線）と対応するデータ値（評定Rが２以上および３以上の比率）のKDE（Kernel Density Estimation）による平滑化されたグラフ（破線）が表示されている。図６のフォームを閉じると、次のようにw_kの値の現在値が表示され、次の値の設定が求められる。
The current value of w_k = 4.677120320447509
The nexr value of w_k =
このw_kは平滑化におけるkernelの窓の幅で、この値が大きいと大まかな滑らかな曲線が得られ、小さいときは細かく変動した凸凹の曲線になる。
w_kの値として０以下の値（次の例では０）を設定すると、プログラムは終了する。
The current value of w_k =
4.677120320447509
The nexr value of w_k = 0
Output file Results.txt was saved.
(py37) PS D:\＊＊＊\files3Cat>
表示されている出力ファイルResults.txtを開くと、以下のようになっている。
Data file...Data3Cat.csv
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...
Inference for Stan model:
anon_model_ceeb5395c8b2c88c8327d0f4a1698755.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
C 10.19 0.05 1.77 7.21 8.95 10.04 11.27 14.01 1398 1.0
C1 193.54 0.07 2.73 187.7 191.79 193.72 195.46 198.37 1517 1.0
C2 213.93 0.04 2.42 209.48 212.29 213.81 215.44
219.06 4000 1.0
mu 203.73 0.03 1.88 199.9 202.52 203.74 205.02 207.25 2990 1.0
sgm 14.85 0.06 2.48 11.03 13.11 14.47 16.19 20.68 1568 1.0
lp__ -70.06 0.03 1.27 -73.25 -70.62 -69.73 -69.14
-68.61 1397 1.0
Samples were drawn using NUTS at Sun
Oct 22 16:20:36 2023.
For each parameter, n_eff is a crude
measure of effective sample size,
and Rhat is the potential scale
reduction factor on split chains (at
convergence, Rhat=1).
MAP for mu = 203.60 95% CI for mu = [199.90,
207.25]
Map for sigma = 13.87 95% CI for sigma = [11.03,
20.68]
MAP for C = 9.71 95% CI for C =
[7.22, 14.01]
MAP for PSE = 203.60, 95% CI for PSE = [199.90,
207.25]
MAP for JND = 9.35, 95% CI for JND= [7.44,
13.95]
読み込まれたデータの出力後、ベイズ分析の統計量が出力されている。データが正しく読み込まれたかどうかの確認と、ベイズ分析の統計量を出力ファイルでも読むことができる。
引用文献
岡本安晴 いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――．丸善出版
岡本安晴（2014）心理学データ分析と測定――データの見方と心の測り方――．勁草書房
岡本安晴（2019）．多件法による上下法．基礎心理学研究、38, 90-104.
岡本安晴（2022）ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――．基礎心理学研究、41, 1-7．
大山 正（1971）．感覚・知覚測定法 和田陽平・大山 正・今井省吾（編著）感覚・知覚心理学ハンドブック（pp. 32–55）誠信書房
大山正（2005）．第７章 精神物理学的測定法． 大山正・岩脇三良・宮埜嘉夫（共著）．心理学研究法――データ収集・分析から論文作成まで――．サイエンス社．
リスト１のStanスクリプトを利用して３件法による上下法のデータを分析するPythonスクリプトを以下に示す。これらのファイルは圧縮ファイルfiles3Cat.zipにまとめたので、ダウンロードして自由に利用されたい。
リスト２ リスト１のStanスクリプトによる上下法データのベイズ分析Pythonスクリプト．モジュールmy_module.pyはリスト３に示す。
import numpy as np
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
import pickle
import csv
from my_module import *
fin_nm = input('Data file name (*.csv)
= ')
fout_nm = input('Output file name = ')
fout = open(fout_nm, 'w')
fout.write('Data
file...{}\n'.format(fin_nm))
#
# Read data
from CSV file
#
with open(fin_nm, 'r') as fin:
data = [v for v in
csv.reader(fin)]
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:>5s} {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])
#
# Boundaries
and initial values
#
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}
#
# Sampling
by Stan
#
sm = pystan.StanModel(file =
'PM3Cat_EqC.stan')
fit = sm.sampling(data = stan_data,
pars = ['C', 'C1', 'C2', 'mu', 'sgm'],
init = f_init,
n_jobs = 1)
print(fit)
fout.write('\nSummary...\n{}\n'.format(fit))
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)
#
# 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(SmplMu)[0]
sgmMAP = CalcMAPEst(SmplSgm)[0]
CMAP = CalcMAPEst(SmplC)[0]
print('muMAP = {0} sgmMAP = {1} CMAP = {2}'.format(muMAP, sgmMAP,
CMAP))
fout.write('\nMAP for mu =
{0:<.2f} 95% CI for
mu = [{1:<.2f}, {2:<.2f}]\n'.
format(muMAP, 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(muMAP, 0.0, 'go', markersize
= 10, label = 'MAP')
plt.yticks([])
plt.xlabel('$\mu$', fontsize = 14)
plt.title('Posterior Distribution for
$\mu$\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(muMAP,
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('\nMap for sigma =
{0:<.2f} 95% CI for
sigma = [{1:<.2f}, {2:<.2f}]\n'.
format(sgmMAP, SgmL025, SgmU975))
sb.kdeplot(SmplSgm)
plt.plot([SgmL025, SgmU975], [0.0,
0.0], color = 'b', linewidth = 5, label = '95% CI')
plt.plot(sgmMAP, 0.0, 'go', markersize
= 10, label = 'MAP')
plt.yticks([])
plt.xlabel('$\sigma$', fontsize = 14)
plt.title('Posterior Distribution for
$\sigma$\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(sgmMAP,
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('\nMAP for C =
{0:<.2f}
95% CI for C = [{1:<.2f}, {2:<.2f}]\n'.
format(CMAP, CL025, CU975))
sb.kdeplot(SmplC)
plt.plot([CL025, CU975], [0.0, 0.0],
color = 'b', linewidth = 5, label = '95% CI')
plt.plot(CMAP, 0.0, 'go', markersize =
10, label = 'MAP')
plt.yticks([])
plt.xlabel('C', fontsize = 14)
plt.title('Posterior Distribution for
C\n' +
('MAP = {0:<.2f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(CMAP,
CL025, CU975)),
fontsize = 16)
plt.legend()
plt.tight_layout(h_pad = 5)
plt.show()
#
# 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 = sgmMAP * s75
PSE = muMAP
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 +
CMAP))
yvalues_R.append(pf.v(t -
CMAP))
plt.title('Psychometric Functions\n' +
('PSE = {0:<.2f}, JND =
{1:<.2f}, C =
{2:<.2f}'.format(PSE, JND, CMAP)),
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, CMAP)),
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))
リスト２で用いられているモジュールファイルmy_module.pyをリスト３に示す。これも、圧縮ファイルfiles3Cat.zipあるいはfiles3Cat_jn.zipにまとめられている。
リスト３ モジュールファイルmy_module.py
import numpy as np
import scipy
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