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

PyStan 3対応

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

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

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

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

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

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

のとき,

判断がのときに前試行における刺激値との関係で次試行の刺激値が決められているのは、「強い」判断および「弱い」判断のデータの多い方が心理測定関数（psychometric function: PF）の傾きの推定が安定するからである。精神測定関数（心理測定関数）については、岡本（2014）, pp. 167-172で説明している。

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

ここで

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

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

この尤度関数に対して、事後分布を求めるStanスクリプトをリスト１に示す（岡本、2019）。の値は、－、－？、＋？、＋を１、２、３、４で表している。

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

data {

int N;

array[N] int Res;

array[N] real X;

real mid_x;

real max_x;

}

parameters {

real<lower = 0.0> v_sgm;

real v_mu;

real<lower = 0.0> v_C;

}

transformed parameters {

real mu;

real C;

real sgm;

array[N] vector[4] p;

mu = mid_x + v_mu * (max_x - mid_x);

sgm = 0.001 + v_sgm * (max_x - mid_x);

C = v_C * (max_x - mid_x);

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 {

v_mu ~ normal(0.0, 10.0);

v_C ~ exponential(0.1);

v_sgm ~ exponential(0.1); // uniform(min_sgm, max_sgm);

for (i in 1:N)

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

}

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

リスト２のスクリプトはWindowsにインストールされたUbuntu上で実行した。WindowsにWSL（Windows Subsystem for Linux）によりUbuntuをインストールして、UbuntuにおいてAnacondaをインストールした後、Python 3.11仮想環境を構築した。このPython 3.11仮想環境においてリスト２を実行した。WindowsにWSLによりUbuntuをインストールしてPyStan 3を使う方法については、このウェブサイトを参照されたい。

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

(py311) ******/files4CatPS3$ python MainFourCatDiscrimPS3.py

Data file name (*.xlsx) =

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

・

図１

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

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

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

(py311) ******/files4CatPS3$ python MainFourCatDiscrimPS3.py

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

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

入力ファイル名と出力ファイル名の設定が終わると、データ読み込まれ、Stanによる事後分布からのサンプリングが始まる。

サンプリングが終了すると、サンプリングのトレース図が表示される（図２）。

図２

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

図３

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

図３のフォームを閉じると、パラメータsgmの事後分布のグラフが表示される（図４）。

図４

図４のフォームを閉じると、パラメータCの事後分布のグラフが表示される（図５）。

図５

図５のフォームを閉じると、PF（心理測定関数）とデータの散布図がPSE、JNDおよびCのMAP推定値とともに表示される（図６）。

図６

PFは次式

によってMAP推定値を使って描かれている。

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

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

図７

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

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

The current value of w_k = 4.422319255494346

The next value of w_k =

W_kの値として大きい値を設定すると、KDEの窓の幅が広がり、より滑らかな曲線となる。

The current value of w_k = 4.422319255494346

The next value of w_k = 9

上は、w_k=9を設定した場合であるが。下図（図８）のようなグラフになる。

図８

値として０以下のものを設定すると、プログラムの実行終了である。端末画面には以下のような出力が表示されている。

The current value of w_k = 4.422319255494346

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.

(py311) ******/files4CatPS3$

上の画面出力では、出力ファイルとしてテキストファイルResults.txtが保存されたことが表示されている。

テキストファイルResults.txtの内容は、以下のようである。

Data file...Data4Cat.xlsx

Number of data = 100

0 257.50 4

1 237.50 4

2 222.50 4

　　　　　・

97 192.50 2

98 177.50 1

99 192.50 1

Summary:

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \

mu 200.404 1.865 196.922 204.000 0.038 0.027 2401.0 2031.0

sgm 14.143 2.313 10.314 18.713 0.049 0.035 2394.0 2362.0

C 8.140 1.458 5.696 11.119 0.030 0.021 2450.0 2696.0

r_hat

mu 1.0

sgm 1.0

C 1.0

MAP for PSE = 200.43, 95% CI for PSE = [196.67, 204.05]

MAP for JND = 8.84, 95% CI for JND= [7.12, 13.21]

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

引用文献

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

岡本安晴（2011）上下法データのベイズ的分析 : 確率モデルによるデータの活用．基礎心理学研究、30、44-55．

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

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

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

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

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

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

付記　Pythonスクリプト

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

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

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import stan

from my_module import *

import pandas as pd

import arviz as az

pd.set_option('display.max_columns', 10000)

pd.set_option('display.max_rows', 10000)

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

# Values for rescaling in MCMC sampling

mid_x, max_x = np.percentile(X, [50, 95])

def f_init():

"""

Initial values for v_C, v_mu, and v_sigma

"""

return dict(v_C = 0.5, v_mu = 0.0, v_sgm = 1.0)

stan_data = {'N': N, 'X':X, 'Res': Res, 'mid_x': mid_x, 'max_x': max_x}

# Sampling by Stan

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

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

fit = sm.sample(num_chains = 4, init = [f_init()]*4)

dframe = fit.to_frame()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

print(az.summary(i_data, var_names = ['mu', 'sgm', 'C']))

fout.write('\nSummary:\n')

fout.write("{0}\n".format(az.summary(i_data, var_names = ['mu', 'sgm', 'C'])))

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

plt.tight_layout()

plt.show()

ck = az.plot_posterior(i_data, var_names = 'mu', point_estimate = 'median')

plt.show()

ck = az.plot_posterior(i_data, var_names = 'sgm', point_estimate = 'median')

plt.show()

ck = az.plot_posterior(i_data, var_names = 'C', point_estimate = 'median')

plt.show()

# MAP estimates of parameters

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

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

CMAP = CalcMAPEst(dframe['C'])[0]

# 2.5 and 97.5 percentiles of posterior distributions

MuL025, MuU975 = np.percentile(dframe['mu'], [2.5, 97.5])

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

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

s75 = ss.norm.ppf(0.75, 0, 1)

print('s75 = ', s75, ' p = ', f_phi(s75)) # f_phi from my_module.py

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をリスト３に示す。これらは、圧縮ファイルfiles4CatPS3.zipにまとめられている。

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

import numpy as np

import scipy.stats as ss

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

"""

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]

f_phi = ss.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 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