Up

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

PyStan（Python用Stan）インストールの説明ウェブサイトを用意した。

サンプルプログラムとデータファイルは圧縮ファイルfiles4Cat.zipにまとめた。Jupyter Notebook用のものも圧縮ファイルfiles4Cat_jn.zipにまとめた。

それぞれダウンロード・解凍すれば、自由に利用することができる。

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

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

心理物理学（精神物理学、psychophysics）において主観的等価点を求める方法として上下法（Up-down method）あるいは階段法（Staircase method）と呼ばれている方法がある。上下法では、提示刺激値が直前の試行における判断に依存して決まるため、被験者が提示刺激値を予測する可能性が高い。これを避けるため、２系列用意して２つの系列の試行をランダムに混ぜて実施する（二重上下法；大山、2005）。２つの系列を系列１と系列２として、系列でのものであることを右上に添え字として明記する。

４件法の場合、判断は例えば、「強い」、「どちらかといえば強い」、「どちらかといえば弱い」、「弱い」の４通りである。系列の第試行における判断（反応）を以下のように表す。

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

のとき,

のとき,

のとき,

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

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

ここで

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

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

この尤度関数に対して、事前分布として一様分布を設定したときの事後分布を求めるStanスクリプトをリスト１に示す（岡本、早期公開2019.6.21）。の値は、−、−？、＋？、＋を１、２、３、４で表している。

リスト１　４件法上下法データのベイズ分析用Stanスクリプト

data {

int N;

int Res[N];

real X[N];

real min_sgm;

real max_sgm;

}

parameters {

real<lower = 0.0> sgm;

real mu;

real<lower = 0.0> C;

}

transformed parameters {

vector[4] p[N];

for (i in 1:N){

p[i][1] = 1 - Phi((X[i] - mu + C) / sgm);

p[i][2] = Phi((X[i] - mu + C) / sgm) - Phi((X[i] - mu) / sgm);

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

p[i][4] = Phi((X[i] - mu - C) / sgm);

}

}

model {

sgm ~ uniform(min_sgm, max_sgm);

for (i in 1:N)

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

}

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

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

図１

データファイルは、CSV形式のファイルとして用意するが（Excelで保存時に、ファイル拡張子として.csvを指定すればよい）、形式は図２のようである。

・

・

・

図２

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

データの設定が終われば、ファイル拡張子として.csvを選んで保存すればよい。CSV形式のファイルについては、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、pp. 33-35, 101-102＞で説明している。パラメータの事後分布が発散しないように、提示刺激には「強い判断」の確率が１に近い十分強い刺激と「弱い」判断の確率が０に近い十分弱い刺激が含まれるようにする。２系列用意する上下法において、「強い」判断が容易な強い刺激値で始まる系列と「弱い」判断が容易な弱い刺激値で始まる系列の２系列を用意すればよい。

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

図３

入力ファイル名と出力ファイル名の設定が終わると、データ読み込まれ、Stanによる事後分布からのサンプリングが始まる。リスト２のプログラムでは、以下のように初期化の関数f_initとデータの辞書stan_dataが設定されて、リスト１のStanスクリプトPM4Cat.stanによるサンプリングが行われる。PythonにおけるStan（PyStan）の使用については、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、pp. 178-208＞で説明している。

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

min_sgm = 0.0

max_sgm = np.amax(X) - np.amin(X)

def f_init():

"""

Initial values for C, mu, and sigma

"""

return dict(C = (C1C2[1] - C1C2[0]) / 2.0, mu = (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}

#

# Sampling by Stan

#

sm = pystan.StanModel(file = 'PM4Cat.stan')

fit = sm.sampling(data = stan_data, init = f_init,

pars = ['sgm', 'C', 'mu'], n_jobs = 1)

サンプリングが終了すると、パラメータmu()の事後分布のグラフがMAP推定値、95％ CIとともに表示される（図４）。グラフは、保存ボタンをクリックすると保存できる。

図４

フォーム右上のX印をクリックしてフォームを閉じると、次に、sigma（）の事後分布のグラフが、MAP推定値と95％ CIとともに表示される（図５）。

図５

図５のフォームを閉じると、次にCの事後分布のグラフがMAP推定値と95％ CIとともに表示される（図６）。

図６

図６のフォームを閉じると、次にPFとデータ値のグラフがMAP推定値に基づくPSEとJNDの推定値とともに表示される（図７）。

図７

PFは次式

によって描かれている。

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

による変換を行っている。さらに、上下に少しランダムな変動を加えて、完全に重なって表示されることを防いでいる。図７のフォームを閉じると図８のフォームが提示される。

図８

データ（評定値R）がR>1、R>2、,R>3である比率をKDEにより滑らかに均した曲線と、推定されたパラメータ値PSE、JND、CによるPF曲線が描かれている。

図８のフォームを閉じると、w_kの現在値の表示と次の値の設定が求められる。

図９

w_kはKDEの窓の大きさであり、この値が大きいと滑らかな大まかな曲線が得られ、小さいときは凹凸の激しい曲線となる。w_kの値として０以下の値（図９の場合）を設定すると、プログラムは終了する。

図１０

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

Data file ...Data4Cat.csv

Number of data = 100

0 257.50 4

1 237.50 4

2 222.50 4

3 207.50 4

4 192.50 3

5 177.50 1

　　　　　　　・

　　　　　　　・

　　　　　　　・

95 192.50 3

96 207.50 4

97 192.50 2

98 177.50 1

99 192.50 1

Summary...

Inference for Stan model: anon_model_82d30ceeafb737a852ee33cf7874bf36.

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

sgm 14.19 0.05 2.4 10.48 12.59 13.82 15.45 20.13 1944 1.0

C 8.2 0.03 1.52 5.62 7.12 8.08 9.08 11.62 2122 1.0

mu 200.48 0.04 1.86 196.74 199.3 200.5 201.68 204.08 2305 1.0

lp__ -88.47 0.04 1.35 -91.98 -89.05 -88.1 -87.48 -86.98 1385 1.0

Samples were drawn using NUTS at Wed Jul 10 11:14:43 2019.

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 = 200.55 95% CI for mu = [196.74, 204.08]

Map for sigma = 12.84 95% CI for sigma = [10.49, 20.12]

MAP for C = 7.60 95% CI for C = [5.62, 11.61]

MAP for PSE = 200.55, 95% CI for PSE = [196.74, 204.08]

MAP for JND = 8.66, 95% CI for JND= [7.07, 13.57]

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

引用文献

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

岡本安晴（早期公開2019.6.21）．多件法による上下法．基礎心理学研究、38巻、１号.

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

付記　Pythonスクリプト

４件法による上下法のデータの分析のためのPythonスクリプトをリスト２に、リスト２でimportされているモジュールのPythonスクリプトをリスト３に示す。これらのファイルは圧縮ファイルfiles4Cat.zipにまとめたので、ダウンロードして回答すれば自由に利用できる。Jupyter Notebook用は圧縮ファイルfiles4Cat_jn.zipにまとめた。

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

import numpy as np

import matplotlib.pyplot as plt

import pystan

import pickle

import seaborn as sb

import csv

from my_module import *

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

fout_nm = input('\nOutput 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)]

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

#

# Boundaries and initial values

#

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

min_sgm = 0.0

max_sgm = np.amax(X) - np.amin(X)

def f_init():

"""

Initial values for C, mu, and sigma

"""

return dict(C = (C1C2[1] - C1C2[0]) / 2.0, mu = (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}

#

# Sampling by Stan

#

sm = pystan.StanModel(file = 'PM4Cat.stan')

fit = sm.sampling(data = stan_data, init = f_init,

pars = ['sgm', 'C', 'mu'], 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 fron 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)

print('SgmL25 = {0:.5f} SgmU975 = {1:.5f}'.format(SgmL025, SgmU975))

CL025 = np.percentile(SmplC, 2.5)

CMed = np.percentile(SmplC, 50)

CU975 = np.percentile(SmplC, 97.5)

MuL025 = np.percentile(SmplMu, 2.5)

MuMed = np.percentile(SmplMu, 50)

MuU975 = np.percentile(SmplMu, 97.5)

#

# MAP estimates for mu, sigma, and C

#

muMAP, sgmMAP, CMAP = calcMAP(X, Res, MuMed, SgmMed, CMed)

print('muMAP = {0} sgmMAP = {1} CMAP = {2}'.format(muMAP, sgmMAP, CMAP))

#

# KDE plot for the samples from posterior distribution for mu

#

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

format(muMAP, MuL025, MuU975))

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_L, 'g-', linewidth = 3, label = 'PF(>1)')

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

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

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) / 3 + np.random.uniform(-0.03, 0.03, 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

w_3 = np.empty(len(Res)) # Weights for response >= 4

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

w_3[t] = 0 if Res[t] < 4 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

pm_3 = [] # KDE values for Response >= 4

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

pm_3.append(k_ave(v, X, w_3, 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, 'k-', linewidth = 3, label = 'PF(>2)')

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

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, 'k--', linewidth = 2, label = 'Res(>2)')

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

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

リスト２で用いられているモジュールファイルmy_module.pyをリスト３に示す。これらは、圧縮ファイルfiles4Cat.zipにまとめられている。Jupyter Notebook用は圧縮ファイルfiles4Cat_jn.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] == 4:

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

elif self.Res[t] == 3:

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

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) / 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 value for square root of sigma

z[2] = CMed ** 0.5 # Initial value 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):

"""

Kernel function with width w_k

"""

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

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

"""

Average by kernel function my_kernel

"""

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