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

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

CmdStanPy対応

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

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

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

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

のとき,

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

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

ここで

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

、

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

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

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

リスト２のスクリプトは、CmdStanPy対応である。CmdStanPyの準備の簡単な説明を、このウェブサイトに用意した。

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

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

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

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

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

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

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

グラフィカルユーザーインターフェイス, アプリケーション, テーブル, Excel

AI 生成コンテンツは誤りを含む可能性があります。

・

グラフィカルユーザーインターフェイス, アプリケーション, テーブル, Excel

AI 生成コンテンツは誤りを含む可能性があります。

図１

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

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

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

入力ファイル名と出力ファイル名の設定が終わると、データが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルには、多少時間が掛かる。コンパイルが終わると、MCMCサンプリングが始まる。

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

図２

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

図３

点推定値として中央値が示されている。Gelmanら（2021）は、平均値と中央値を比べて、中央値の方が計算においてより安定しているとして中央値を勧めている（p. 73）。

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

図４

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

図５

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

図６

PFは次式

によって、中央値を点推定値として描かれている。

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

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

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

図７

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

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

The current value of w_k = 4.670716402259204

The next value of w_k =

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

The current value of w_k = 4.670716402259204

The next value of w_k = 9

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

図８

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

The current value of w_k = 4.670716402259204

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.

(stan) *****/updwn4catfiles$

上の画面出力では、出力ファイルとしてテキストファイル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

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -94.128100 0.031063 1.298030 ... 2334.48 1243.98 1.000450

v_sgm 0.897401 0.003212 0.147883 ... 2070.81 1611.28 1.001030

v_mu -0.448059 0.002260 0.118233 ... 2471.53 1917.04 1.001270

v_C 0.519510 0.002204 0.095828 ... 2654.92 1361.73 1.001130

mu 200.443000 0.035589 1.862170 ... 2471.53 1917.04 1.001270

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

p[99,4] 0.016485 0.000274 0.014324 ... 3305.81 2048.37 1.000700

p[100,1] 0.497568 0.001293 0.068602 ... 2587.96 1957.65 1.000910

p[100,2] 0.218222 0.000565 0.036770 ... 3221.91 2959.77 1.000020

p[100,3] 0.156833 0.000480 0.026413 ... 3055.43 2139.51 0.999787

p[100,4] 0.127378 0.000744 0.042085 ... 2988.01 2343.64 1.001190

[407 rows x 11 columns]

Summary:

mean sd hdi_3% hdi_97% ... mcse_sd ess_bulk ess_tail r_hat

mu 200.443 1.862 197.146 204.075 ... 0.033 2759.0 2472.0 1.0

sgm 14.135 2.329 10.353 18.736 ... 0.046 2319.0 2071.0 1.0

C 8.182 1.509 5.672 11.122 ... 0.024 1960.0 2655.0 1.0

[3 rows x 9 columns]

Med. for PSE = 200.44, 95% CI for PSE = [196.74, 204.05]

Med. for JND = 9.34, 95% CI for JND= [7.08, 13.06]

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

リスト１　４件法データのベイズ分析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]);

}

リスト２　リスト１のスクリプトを利用するPythonスクリプト（main.py）

from cmdstanpy import CmdStanModel

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import pandas as pd

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

# コンパイル

model = CmdStanModel(stan_file="PM4CatPS3.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 = ['mu', 'sgm', 'C'])

plt.tight_layout()

plt.show()

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

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

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

CMed = np.median(dframe['C'])

# 2.5, 50 and 97.5 percentiles of posterior distributions

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

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

f_phi = ss.norm.cdf # Cumulative standard normal distribution function

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

def my_kernel(x, y, w_k):

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

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

"""

Kernal estimate

"""

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

# 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_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, CMed)),

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

参考文献

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．

岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック第三版（第2章）、誠信書房

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