上下法（階段法）データのベイズ分析

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

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

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

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

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

心理物理学（精神物理学、psychophysics）において主観的等価点を求める方法として上下法（Up-down method）あるいは階段法（Staircase method）と呼ばれている方法がある。これは、提示された比較刺激が標準刺激より「強い」と判断されたときは、次の試行ではより弱い比較刺激を提示し、「弱い」と判断されたときは次の試行ではより強い比較刺激を提示するものである。すなわち、第試行における比較刺激とそれに対する判断をとで表し、刺激の変化量をで表すとき、

とおくと、

のとき,

である。

上下法では、提示刺激値が直前の判断に依存して決まるため、被験者が提示刺激値を予測する可能性が高い。これを避けるため、２系列用意して２つの系列の試行をランダムに混ぜて実施する（二重上下法；大山、2005）。上式において、２つの系列を系列１と系列２として、系列でのものであることを右上に添え字として明記すると、以下のようになる。

とおくと、

のとき,

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

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

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

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

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

data {

int N; // The number of trials

int Res[N]; // Responses: 1 denotes -, 2 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 pse;

}

transformed parameters {

vector[2] p[N]; // Probabilities for X[N]

for (i in 1:N){

p[i][1] = 1 - Phi((X[i] - pse) / sgm);

p[i][2] = Phi((X[i] - pse) / sgm);

}

model {

sgm ~ uniform(min_sgm, max_sgm);

for (i in 1:N)

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

}

リスト１のStanスクリプトによる分析プログラム例を付記にリスト２として示す。プログラムファイルとサンプルデータファイルは圧縮ファイルfiles2Cat.zipとしてまとめた。ダウンロードして自由に利用されたい。Jupyter Notebook用も圧縮ファイルfiles2Cat_jn.zipとして用意している。

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

図１

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

・

図２

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

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

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

図３

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

min_sgm = 0.0 # lower boundary for sigma

max_sgm = np.amax(X) - np.amin(X) # upper boundary for sigma

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

def f_init():

"""

Initial values for pse(mu) and sigma

"""

return dict(pse = (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 = 'PM2Cat.stan')

fit = sm.sampling(data = stan_data, pars = ['sgm', 'pse'],

init = f_init, n_jobs = 1)

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

図４

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

図５

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

図６

データ値は、次式

により、０あるいは１に変換した後、上下に少しランダムな変動を加えて、完全に重なって表示されることを防いでいる。

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

図７

データ（判断R=2の比率）をKDEで滑らかにした曲線（緑の破線）とパラメータPSEとJNDの推定値に基づくPF曲線（青の曲線）が描かれている。

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

図８

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

図９

プログラムの終了後、表示されている出力ファイルを開くと以下のようになっている。

Data file...Data2Cat.csv

Number of data = 100

0 250.00 2

1 230.00 2

2 215.00 2

3 200.00 2

4 185.00 1

5 200.00 1

　　　・

95 200.00 1

96 215.00 2

97 200.00 1

98 215.00 1

99 230.00 2

Summary...

Inference for Stan model: anon_model_c563ee5b69f9c9fbc61d408eed6fd018.

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 18.57 0.12 4.9 12.14 15.29 17.61 20.74 30.68 1704 1.0

pse 203.05 0.06 2.71 197.63 201.33 203.03 204.77 208.27 2001 1.0

lp__ -48.87 0.03 1.14 -52.17 -49.26 -48.54 -48.1 -47.8 1135 1.0

Samples were drawn using NUTS at Wed Jul 10 13:35:13 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).

mu(MAP) = 203.094 sigma(MAP) = 15.718

PSE = 203.094, 95% CI = [197.63, 208.27]

JND = 10.60, 95% CI = [8.19, 20.67]

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

引用文献

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

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

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

付記　Pythonスクリプト

上下法データの分析を行うPythonスクリプトをリスト２およびリスト３として用意した。メインスクリプトがリスト２であり、リスト３のモジュールがimportされている。これらのファイルは、圧縮ファイルfiles2Cat.zipとしてまとめられているので、ダウンロードすれば解凍して自由に利用できる。Jupyter Notebook用も圧縮ファイルfiles2Cat_jn.zipとしてまとめた。

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

import numpy as np

import matplotlib.pyplot as plt

import pystan

import pickle

import seaborn as sb

import csv

import scipy.optimize as sciopt

from my_module import *

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

fout_nm = input('\nOutput file name (*.txt) = ')

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'.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:>5s}{1:>7.2f}{2:>5d}\n'.format(i, s, r))

# Boundaries and initial values

min_sgm = 0.0 # lower boundary for sigma

max_sgm = np.amax(X) - np.amin(X) # upper boundary for sigma

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

def f_init():

"""

Initial values for pse(mu) and sigma

"""

return dict(pse = (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 = 'PM2Cat.stan')

fit = sm.sampling(data = stan_data, pars = ['sgm', 'pse'],

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

SmplPSE = fit['pse'] # Samples from the posterior distribution for pse(mu)

# Percentile points for posterior distributions

SgmL025 = np.percentile(SmplSgm, 2.5)

SgmMed = np.percentile(SmplSgm, 50)

SgmU975 = np.percentile(SmplSgm, 97.5)

PSEL025 = np.percentile(SmplPSE, 2.5)

PSEMed = np.percentile(SmplPSE, 50)

PSEU975 = np.percentile(SmplPSE, 97.5)

# Calculate MAP estimates for mu(pse) and sigma

muMAP, sgmMAP = calcMAP(X, Res, PSEMed, SgmMed) # calcMAP from my_module.py

print('mu = ', muMAP, ' sgm = ', sgmMAP)

fout.write('\nmu(MAP) = {0:<.3f} sigma(MAP) = {1:<.3f}\n'.format(muMAP, sgmMAP))

# KDE plot for the samples from posterior distribution for PSE(mu)

sb.kdeplot(SmplPSE)

plt.plot([PSEL025, PSEU975], [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, PSEL025, PSEU975)),

fontsize = 16)

plt.legend()

plt.tight_layout(h_pad = 5)

plt.show()

# KDE plot for the samples from posterior distribution for sigma

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

# 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('\nPSE = {0:<.3f}, 95% CI = [{1:<.2f}, {2:<.2f}]'.format(

PSE, PSEL025, PSEU975))

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

for i in range(1001):

t = minX + i * (maxX - minX) / 1000.0

xvalues.append(t)

yvalues.append(pf.v(t)) # Values of the psychometric function

plt.title('Psychometric Function\n' +

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

fontsize = 16)

plt.plot(xvalues, yvalues, linewidth = 3, label = 'PF')

plt.yticks([0.0, 0.5, 0.75, 1.0])

# Resposes are scaled to values from 0 to 1, and added small fluctuations

Res_rn = (Res - 1) + 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

for t in range(len(Res)):

w_1[t] = 0 if Res[t] < 2 else 1

X_range = np.arange(minX, maxX, (maxX - minX) / 100)

w_k = JND / 2

while True:

pm_1 = [] # KDE values for Response >= 2

for v in X_range:

pm_1.append(k_ave(v, X, w_1, w_k))

plt.title('Observed and Estimated PFs\n' +

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

fontsize = 16)

plt.plot(xvalues, yvalues, color = 'b', linewidth = 3, label = 'PF')

plt.yticks([0.0, 0.5, 1.0])

plt.plot(X_range, pm_1, 'g--', linewidth = 2, label = 'Res(>1)')

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

リスト２で用いられているモジュールファイルmy_module.pyをリスト３に示す。これらは、圧縮ファイルfiles2Cat.zipにまとめられている。Jupyter Notebook用は圧縮ファイルfiles2Cat_jn.zipとして用意した。

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

import numpy as np

import scipy

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

v = 0

for t in range(self.N):

if self.Res[t] == 2:

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

else:

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

return v

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

"""

Calculation of MAP estimates

"""

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

z = np.empty(2)

z[0] = PSEMed # Initial value for pse(mu)

z[1] = SgmMed ** 0.5 # Initila value for square root of sigma

import scipy.optimize

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

return muMAP, sgmMAP

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