Up

弁別閾と主観的等価点のベイズ分析（CmdStanPy 対応）

２件法、３件法、４件法データの分析

２件法より３件法、４件法の方が統計学的に望ましいことを示すシミュレーション結果が報告されている（岡本、2017）。

確率モデルを用いると、弁別閾と判断の基準値を分離した分析が行えるというシミュレーションがある（岡本、1995）。

なお、psychometric関数の説明を＜岡本安晴「心理学データ分析と測定」、第11章＞で行っている。

感覚・知覚測定法の概説が、岡本（2025）にある。

弁別閾・丁度可知差異あるいは主観的等価点は、理論的には心理測定関数（psychometric function: PF）により与えられる（岡本、2025）。

Psychometric関数（PF）を求めるpsychophysics（精神物理学／心理物理学）実験としては、反応が「大きい／小さい」など（本ウェブサイトでは、物理刺激の比較判断を「大きさ」の次元を例として用いている。他の感覚、例えば「重さ」であれば、「重い／軽い」判断となる）の２カテゴリ強制選択法（２件法）が標準であるが、「わからない／同じ」という判断を認める３件法もある。３件法については、データ分析を単純な平均値に基づいて行う場合には、感覚の弁別力（閾値）の推定値が判断の基準の影響を受けるということが知られているが（岡本、1995）、この問題は確率モデルを設定して分析を行えば回避できる。設定された確率モデルに対して最尤法を用いた場合（岡本、1994、1995）、およびベイズ的方法を用いた場合（岡本、2005a）において、感覚の弁別力（閾値）と判断の基準が確率モデルによって分離されて適切に分析できることが報告されている。主観的等価点PSEを求めることを考慮して、２件法と３件法の混合型ともいえる４件法（「大きい」、「どちらかといえば大きい」、「どちらかといえば小さい」、「小さい」の４カテゴリを用いる方法）もある（Okamoto, 2012）。２件法、３件法、４件法、いずれの場合も分析の基本となるのは物理刺激と反応の関係を表すPFである。

本ウェブサイトで扱われているモデルは、標準刺激と比較刺激の心理物理学的関係についてのものであるので、恒常法や全系列法および上下法などの適応的方法の区別なく適用することができる。観察者の判断カテゴリ数により、カテゴリ数に対応したモデル構成になるので、以下に２件法、３件法および４件法の場合について順番に論じる。なお、実験デザインの立案において、比較刺激の間隔と範囲および総試行数はパラメータの推定値を求めるときに重要な問題であるが、以下に簡単に論じる。

比較刺激の範囲と間隔および総試行数について

総試行数については、上下法を用いると少ない試行数でPSEを求めることができるというのはある意味で間違っている。上下法の場合も、精度を保証するためにはそれなりの試行数が必要であることが、シミュレーションによって示されている（岡本、2019、Tables 1, 2, and 3）。各Tableは、それぞれ総試行数が40試行、100試行、400試行の場合の結果であるが、総試行数が増えるとbiasおよびRMS errorは顕著に減少している。

比較刺激値の範囲と間隔については、以下のようである。

恒常法の場合

比較刺激の範囲と刺激間間隔について少し書いておく。刺激値の範囲は、比較刺激の方が「大きい」および「小さい」の判断が出る範囲であればよい。確率モデルに基づいて分析が行われるので、全系列法のように「大きい」あるいは「小さい」の判断の確率が１になる極端な刺激値を含む必要はないが、「大きい」判断および「小さい」判断の確率の高い比較刺激値に対するデータが含まれているとパラメータの推定は安定する。刺激値の間隔については、岡本（2006a）は、JNDのと3倍とで分析に差がないことを報告している。

適応的方法の場合

試行の始めにおいては、明らかに標準刺激との差が分かる比較刺激値からスタートし、比較刺激値の増減の幅も最初は大きいものから始めて、徐々に小さくする方法がとられるが、このとき、最終的な比較刺激値の増減の幅は過度に小さくする必要はない。これは、確率モデルを用いた分析を行うからで、上の恒常法の場合と同様である。また、先行試行での反応に提示刺激が依存していることを隠すために、２系列の試行を用意してランダムに混合して実施するなどの工夫が必要である。適応的方法としては、Ψ法と呼ばれているエントロピー最小化基準によるものがあり、２件法（岡本、2006b）、３件法（岡本、2004）、４件法（Okamoto, 2012）の場合について説明されている。単純な適応的方法では、ターゲットとなる１つのパラメータに焦点を当てた提示刺激の選択が行われるが、Ψ法では、パラメータ全体の事後分布の最適化がエントロピー最適化基準によって行われる。多件法による上下法（岡本、2019）も提案されているが、単純な手続きであるので実用性が高いと思われる。紙にメモを取りながら実験を進め、実験後、データを確率モデルで分析すればよい。

Fechner問題について

　閾値と感覚尺度の関係を扱う問題はFechner問題と呼ばれているが、これについてはOkamoto (2001, 2005b)で論じている。

２件法のモデル

２件法の場合について考える。比較刺激の方が標準刺激より大きいと判断される（と表記）確率をの関数として表したものがPF（psychometric function）であるが、心理学ではPFとして累積正規分布を用いることが多い。累積標準正規分布をで表すと

と書ける。およびは、標準刺激によって決まる値であり、１つの固定したの値について考えるときは、PFはの関数となる。PFをの関数とみてをで表すと、式（A１）は

となる。２件法なので

である。

ここで、は主観的等価点を表し

である。

弁別閾（difference threshold; Gescheider, 1997）を表す丁度可知差異（JND; 岡本（2014）、pp. 168-169参照）としては、の倍を用いる。ここで、は、確率0.75に対応する標準正規分布における値である。すなわち、

である。プログラムでは

として計算している。すなわち、

である。

第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

ベイズ的分析では、式（２）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。本ウェブサイトでは、Stan用ライブラリとしてCmdStanPyを用いた。CmdStanPyの準備の簡単な説明が、このウェブサイトにある。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。

ベイズ分析では、分析結果は事後分布として与えられる。主観的等価点などの点推定値は、事後分布から求めることになる。分布の代表値としては平均値、中央値、モードがよく知られている。事後分布からの点推定値として、これらの３つの推定値をシミュレーションで比較した岡本（2022）では、モード（MAP推定値）が真値との関係で薦められている。しかし、MAPを用いるためには、モードの安定した事後分布が必要であるので、精度にこだわらない場合は、中央値が実用的と思われる。Gelmanら（2021）は、平均値より安定しているとして中央値を勧めている。

２件法の場合のベイズ分析用Stanスクリプト（pm2cat.stan）をリストA1のように用意した。ファイルはjndfiles.zipにまとめた。

リストA1のStanスクリプトを利用するPythonスクリプト（pm2cat.py）をリストA2のように用意した。リストA2のスクリプトは、CmdStanPy対応である。

リストA1のファイル（pm2cat.stan）、リストA2のファイル（pm2cat.py）、および入力データファイルを同じフォルダにおいて、そのフォルダをカレントディレクトリとして次のコマンドを実行する。

(stan) *****/pm2catfiles$ python pm2cat.py

上のコマンドを実行すると、入力データファイル名と出ファイル名を聞いてくる。

(stan) *****/pm2catfiles$ python pm2cat.py

Input Data File (*.xlsx) = Data2Cat.xlsx

Output Text File (*.txt) = Results.txt

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

・

・

・

図A1

第１行目に変数名を設定して、２行目以降にデータ値を入力する。

第１列目はデータのID用である。プログラムでは無視されるので任意の値あるいは文字列でよい。

第２列目に比較刺激値、第３列目に判断カテゴリ値を入力する。判断カテゴリ値は、「比較刺激値の方が小さい」というときは１、「比較刺激の方が大きい」というときは２を設定する。

入力データファイルとして、上の例では、図A1のファイルData2Cat.xlsxを設定している。

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

上の例では、Results.txtが設定されている。出力ファイル名を設定すると、入力データファイルが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルには多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、まずトレース図（図A2）が表示される。

図A2

図A2のWindowを閉じると、パラメータの事後分布が表示される（図A3）。

図A3

図A3のWindowを閉じると、パラメータの事後分布が表示される（図A4）。

図A4

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

図A5

青の曲線がPFである。このPFはとそれぞれの事後分布の中央値（median）を点推定値として描かれている。主観的等価点はで算出されている。

緑色の小円はデータである。ただし、小円が重ならないように、ランダムな小さい値を加えて表示してある。すなわち、プロットする点を

とするとき、

である。ここで、 はランダムな小さい正数である。

図A5のWindowを閉じると、プログラムの実行終了である。端末には以下の出力が表示されている。

Results.txt was saved.

(stan) *****/pm2catfiles$

出力ファイルResults.txtを開くと、以下のような内容である。

Input Data File = Data2Cat.xlsx

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -55.256700 0.022715 0.989831 ... 2448.74 1791.76 1.00302

v_mu 2.015200 0.000355 0.020207 ... 2568.80 3124.25 1.00138

v_sgm 0.148276 0.000390 0.019777 ... 2126.67 2610.05 1.00212

theta[1,1] 0.924744 0.000640 0.032212 ... 2222.13 2592.27 1.00076

theta[1,2] 0.075256 0.000640 0.032212 ... 2222.13 2592.27 1.00076

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

theta[179,2] 0.217999 0.000927 0.048984 ... 2322.19 2752.74 1.00113

theta[180,1] 0.107539 0.000719 0.039091 ... 2226.27 3058.18 1.00241

theta[180,2] 0.892461 0.000719 0.039091 ... 2226.27 3058.18 1.00241

mu 201.520000 0.035491 2.020710 ... 2568.80 3124.25 1.00138

sgm 14.827700 0.039032 1.977680 ... 2126.67 2610.05 1.00211

[365 rows x 11 columns]

入力データファイル名と、MCMCサンプリングの統計が出力されている。

リストA1　2件法用のStanスクリプト（pm2cat.stan）

data {

int N;

array[N] real St;

array[N] int Res;

}

parameters {

real v_mu;

real<lower = .0> v_sgm;

}

transformed parameters {

array[N] simplex[2] theta;

real mu;

real sgm;

/*

パラメータの変域をデータに合わせる

*/

mu = v_mu * 100;

sgm = v_sgm * 100.0 + 0.0001;

for (i in 1:N) {

theta[i][2] = normal_cdf(St[i] | mu, sgm);

theta[i][1] = 1.0 - theta[i][2];

}

}

model {

v_mu ~ normal(0.0, 10.0);

v_sgm ~ exponential(0.1);

for (i in 1:N)

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

}

リストA2　リストA1のStanスクリプトを利用するPythonスクリプト（pm2cat.py）

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import arviz as az

import scipy.stats as ss

import matplotlib.pyplot as plt

import seaborn as sb

finnm = input('Input Data File (*.xlsx) = ')

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl, header = None)

data = data_xlsx.values

foutnm = input('Output Text File (*.txt) = ')

fout = open(foutnm, 'w')

fout.write('Input Data File = {}\n'.format(finnm))

print('data[0] =\n', data[0])

St = data[1:,1]

Res = data[1:,2]

for v in zip(St,Res):

print(v)

N = len(St)

def f_init():

return {'v_mu':np.median(St)/100.0,

'v_sgm': np.std(St) / 100.0}

Data = {'N':N, 'St':St, 'Res':Res}

model = CmdStanModel(stan_file='pm2cat.stan')

fit = model.sample(data=Data, chains=4, inits=[f_init()]*4)

print(fit.summary())

fout.write(f'\n{fit.summary()}\n')

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

az.plot_trace(inf_data, var_names = ['mu', 'sgm']) # トレースプロット

plt.tight_layout()

plt.show()

d_fm = fit.draws_pd() # PandasのDataFrame型に変換

med_mu = np.median(d_fm['mu'])

med_sgm = np.median(d_fm['sgm'])

sb.kdeplot(d_fm['mu'])

plt.title(fr'Posterior distribution of $\mu$, med = {med_mu:.2f}')

plt.xlabel(r'$\mu$', fontsize = 14)

plt.show()

sb.kdeplot(d_fm['sgm'])

plt.title(fr'Posterior distribution of $\sigma$, med = {med_sgm:.2f}')

plt.xlabel(r'$\sigma$', fontsize = 14)

plt.show()

min_phi = np.min(St)

max_phi = np.max(St)

phis = np.linspace(min_phi, max_phi, 1000)

psis = ss.norm.cdf(phis, med_mu, med_sgm)

p_st = []

p_res = []

for i in range(0, N):

p_st.append(St[i])

if Res[i] == 1:

p_res.append(0.1 * ss.uniform.rvs())

else:

p_res.append(1.0 - 0.1 * ss.uniform.rvs())

plt.title("Psychometric function PSE = {0:.5} JND = {1:.5}".

format(med_mu, med_sgm * 0.67449))

plt.plot(phis, psis, 'b', linewidth = 5)

plt.plot(p_st, p_res, 'go')

plt.xlabel('$S_{comp}$', fontsize = 14)

plt.legend(['PM', 'Data'], loc = 2)

plt.show()

fout.close()

print()

print(foutnm + ' was saved.')

print()

３件法のモデル

３件法データの分析モデルは、２件法の拡張として考える。３件法では、判断は「標準刺激より比較刺激の方が小さい；」、「とは同じ／どちらが大きいかわからない；」、「よりの方が大きい；」の３つのカテゴリから選ばれる。これらを、

「」と「あるいは（と表記）」、

および

「あるいは（と表記）」と「」

の２通りに分けて考える。

「」と「」の組み合わせに対するＰＦを

とおく。は、判断「」と「」を分ける判断の基準値である。

「」と「」の組み合わせに対するＰＦを

とおく。は、判断「」と「」を分ける判断の基準値である。

判断「」の確率は、次式

で与えられる。また、

である。

JNDは、２件法の場合と同様に

で与えられる。

PSEについては、２つの判断の基準値との中央に主観的等価点PSEがあると考えられるときは、次式

でPSEが推定されるが、Okamoto (2012)は４件法のデータにおいて式（B1）が成り立たない例を報告している。

いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

ベイズ的分析では、式（B2）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。本ウェブサイトでは、CmdStanPyによりStanを用いている。CmdStanPyの準備についての簡単な説明がこのウェブサイトにある。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。

３件法の場合のベイズ分析用StanスクリプトをリストB1のように用意した。ファイルはjndfiles.zipにまとめている。

リストB1のStanスクリプトを利用するPythonスクリプトをリストB2のように用意した。

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

(stan) *****/pm3catfiles$ python pm3cat.py

実行が始まると、入力データファイル名と出力ファイル名の設定が求められる。

(stan) *****/pm3catfiles$ python pm3cat.py

Input Data File (*.xlsx) = Data3Cat.xlsx

Output Text File (*.txt) = Results.txt

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

・

・

・

図B1

第１行目に変数名を書く。データ値は２行目以降に設定する。

第１列目にはデータ識別用の値を設定する。この列の値は、プログラムでは無視されるので、数値でも文字列でもよい。

第２列目に、提示刺激値、第３列目に判断カテゴリ値を設定する。判断カテゴリ値は、比較刺激の方が小さいなら１、同じなら２、大きいなら３を設定する。

入力データファイル名を設定すると、次に出力ファイル名の設定が求められるので、適当なテキストファイル名を設定する。上の例では、Results.txtが設定されている。

出力ファイル名を設定すると、計算が始まり、Stanスクリプトがコンパイルされる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了すると、トレース図が表示される（図B2）。

図B2

図B２のWindowを閉じると、図B3のグラフが表示される。

図B3

パラメータC1とC2の事後分布と、次式（B３）で与えられるC1とC2の中点midCsの事後分布が描かれている。２つの判断の基準値との中央に主観的等価点PSEがあると考えられるときは、次式

でPSEが推定される。

図B3のWindowを閉じると、パラメータの事後分布が表示される（図B4）。

図B4

図B4のWindowを閉じると、図B5のグラフが表示される。

図B5

２つの心理測定関数

と

のグラフ、および

のグラフである。

感覚の弁別力は心理測定関数の傾きで表され、これはパラメータの逆数である。これに対して、判断の基準の厳しさは心理測定関数の平行移動の度合いで表され、心理測定関数間の横方向の距離が小さいと判断の基準は厳しいと言える。確率モデルを用いた分析により、感覚の弁別力と判断の厳しさを分離することができるが、岡本（1995）は全系列法のシミュレーションによってこのことを示している。

青の小円はデータ点を表すが、データ点は

とランダムな小さな値を加えて、データの小円の重なりを避けている。

図B5のWindowを閉じると、プログラムの実行終了である。

端末には、以下のメッセージが出力されている。

Results.txt was saved.

(stan) *****/pm3catfiles$

出力ファイルResults.txtの内容は以下のようである。

Input Data File = Data3Cat.xlsx

Data...

180 1

180 1

220 3

　　　　・

　　　　・

　　　　・

160 1

190 1

220 3

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -100.043000 0.033617 1.299490 ... 2398.79 766.587 1.00206

v_sgm 0.142091 0.000304 0.014444 ... 2569.76 1199.090 1.00121

v_C1 1.924390 0.000419 0.019249 ... 2638.59 1064.680 1.00102

v_C2 0.213148 0.000547 0.025312 ... 2002.59 1082.950 1.00104

theta[1,1] 0.808709 0.000916 0.044743 ... 2332.26 1226.820 1.00075

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

theta[180,2] 0.301735 0.000795 0.052519 ... 3002.74 2200.100 1.00139

theta[180,3] 0.670133 0.000859 0.057039 ... 2948.79 2225.850 1.00124

sgm 14.209200 0.030374 1.444360 ... 2569.76 1199.090 1.00121

C[1] 192.439000 0.041904 1.924910 ... 2638.59 1064.680 1.00101

C[2] 213.753000 0.029645 2.103990 ... 2926.77 2553.720 1.00115

[547 rows x 11 columns]

読み込まれたデータの出力と、MCMCサンプリングの統計が出力されている。

リストB1　３件法用ベイズ分析Stanスクリプト（pm3cat.stan）

data {

int N;

array[N] real St;

array[N] int Res;

}

parameters {

real<lower = .0> v_sgm;

real v_C1;

real<lower = 0.0> v_C2;

}

transformed parameters {

array[N] simplex[3] theta;

real sgm;

array[2] real C;

sgm = v_sgm * 100.0 + 0.0001;

C[1] = v_C1 * 100.0;

C[2] = C[1] + v_C2 * 100.0 + 0.0001;

for (i in 1:N) {

theta[i][2] = normal_cdf((St[i] - C[1])/sgm | 0.0, 1.0);

theta[i][3] = normal_cdf((St[i] - C[2])/sgm | 0.0, 1.0);

theta[i][1] = 1.0 - theta[i][2];

theta[i][2] = theta[i][2] - theta[i][3];

}

}

model {

v_sgm ~ exponential(0.1);

v_C1 ~ normal(0.0, 10.0);

v_C2 ~ exponential(0.1);

for (i in 1:N)

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

}

リストB2　リストB1のStanスクリプトを利用するPythonスクリプト（pm3Cat.py）

from cmdstanpy import CmdStanModel

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

import arviz as az

import scipy.special as ssp

import seaborn as sb

import scipy.stats as ss

finnm = input('Input Data File (*.xlsx) = ')

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl, header = None)

data = data_xlsx.values

foutnm = input('Output Text File (*.txt) = ')

fout = open(foutnm, 'w')

fout.write('Input Data File = {}\n'.format(finnm))

print('data[0] =\n', data[0])

St = data[1:,1]

Res = data[1:,2]

fout.write('\nData...\n')

for v in zip(St,Res):

print(v)

fout.write(f'{v[0]:>10}{v[1]:>10}\n')

N = len(St)

def f_init():

p_mu = np.median(St)

p_sgm = np.std(St)

p_C1 = p_mu * 0.9

p_C2 = p_mu * 1.1

return {'v_sgm':p_sgm/100.0,

'v_C1':p_C1/100.0, 'v_C2':p_C2/100.0}

Data = {'N':N, 'St':St, 'Res':Res}

model = CmdStanModel(stan_file='pm3cat.stan')

fit = model.sample(data=Data, chains=4, inits=[f_init()]*4)

print(fit.summary())

fout.write(f'\n{fit.summary()}\n')

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

az.plot_trace(inf_data, var_names = ['C', 'sgm']) # トレースプロット

plt.tight_layout()

plt.show()

d_fm = fit.draws_pd() # Pandas DataFrame型

C1s = d_fm['C[1]']

C2s = d_fm['C[2]']

midCs = (C1s + C2s) / 2.0

med_C1 = np.median(C1s)

med_C2 = np.median(C2s)

med_midCs = np.median(midCs)

sb.kdeplot(C1s, label = f'C1(med={med_C1:.2f})')

sb.kdeplot(C2s, label = f'C2(med={med_C2:.2f})')

sb.kdeplot(midCs, ls = '--', label = f'midCs(med={med_midCs:.2f})')

plt.xlabel('')

plt.legend()

plt.show()

med_sgm = np.median(d_fm['sgm'])

sb.kdeplot(d_fm['sgm'], label = r'$\sigma$(med={0:.2f})'.format(med_sgm))

plt.xlabel(r'$\sigma$', fontsize = 14)

plt.legend()

plt.show()

min_St = np.min(St)

max_St = np.max(St)

min_x = min_St - 0.05 * (max_St - min_St)

max_x = max_St + 0.05 * (max_St - min_St)

plt.axis([min_x, max_x, 0.0, 1.0])

phi = np.linspace(min_x, max_x, 1000)

psi1 = ss.norm.cdf(phi, med_C1, med_sgm)

psi2 = ss.norm.cdf(phi, med_C2, med_sgm)

psimid = ss.norm.cdf(phi, med_midCs, med_sgm)

p_st = []

p_res = []

for i in range(0, N):

p_st.append(St[i])

if Res[i] == 1:

p_res.append(0.1 * np.random.uniform())

elif Res[i] == 2:

p_res.append(0.45 + 0.1 * np.random.uniform())

else:

p_res.append(1.0 - 0.1 * np.random.uniform())

plt.title(("Psychometric functions\nC1 = {0:.5}, midC = {3:.5}" +

", C2 = {1:.5}, JND = {2:.5}").

format(med_C1, med_C2, med_sgm * 0.67449, med_midCs))

plt.plot(phi, psi1, 'b--', linewidth = 5, label = 'PM1')

plt.plot(phi, psi2, 'b-.', linewidth = 5, label = 'PM2')

plt.plot(phi, psimid, 'b:', linewidth = 4, label = 'midPM')

plt.plot(p_st, p_res, 'go', label = 'Data')

plt.xlabel('$S_{comp}$', fontsize = 14)

plt.legend(handlelength = 7, loc = 2)

plt.show()

fout.close()

print()

print(foutnm + ' was saved.')

print()

４件法のモデル

４件法では、判断は「標準刺激より比較刺激の方が小さい；」、「どちらかといえばよりの方が小さい；」、「どちらかといえばよりの方が大きい；」、「標準刺激より比較刺激の方が大きい；」の４カテゴリから１つを選ぶという形で行われる。これらの判断を以下の３つの２件法データの形式に分けて考える。

（A）「」と

「またはまたは（と表記）」

（B）「または（と表記）」と

「または（と表記）」

（C）「またはまたは（と表記）」と

「」

形式（A）によるデータに対するＰＦを

とおく。形式（B）と（C）に対しても同様に

とおく。

このとき、

である。

上記モデルにおいて、PSEはによって与えられる。すなわち、

である。ここで、

は仮定されていないことに注意。式（C1）が成り立つときは、２つの判断の基準値はPSEから等距離に設定されていることになる。Okamoto (2012)では、４件法データの分析において、式（C1）の成り立たない例が報告されている。

JNDは、２件法および３件法の場合と同様に

で与えられる（岡本、2014、 p. 169）。

いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

ベイズ的分析では、式（C2）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。CmdStanPyによってStanを用いている。CmdStanPyの準備の簡単な説明が、このウェブサイトにある。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。

４件法データのベイズ分析用StanスクリプトをリストC1のように用意した。

リストC1のStanスクリプトを利用するPythonスクリプトをリストC2のように用意した。これらのファイルは、jndfiles.zipにまとめた。

リストC1のファイル（pm4cat.stan）、リストC2のファイル（pm4cat.py）、および入力データファイルを同じフォルダに置き、そのフォルダにカレントディレクトリ（カレントフォルダ）を移動して、CmdStanPyのインストールされている環境で次のコマンドを実行する。CmdStanPyの準備についての簡単な説明が、このウェブサイトにある。

(stan) *****/pm4catfiles$ python pm4cat.py

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

(stan) *****/pm4catfiles$ python pm4cat.py

Input Data File (*.xlsx) = Data4Cat.xlsx

Output Text File (*.txt) = Results.txt

入力データファイルは、Excelファイルの形式で用意する（図C1）。

・

・

・

図C1

第１行目に変数名を設定する。データ値は、２行目以降に書く。

第１列目はデータの識別用で、プログラム内では無視される。この識別用変数の値は、数値でも文字列でもよい。

第２列目は、提示刺激値である。

第３列目には観察者の判断値（反応）を設定する。比較刺激の方が小さいなら１、たぶん小さいなら２、たぶん大きいなら３、大きいなら４を設定する。

入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。

出力ファイルは、任意のテキストファイル名である。上の例では、Results.txtが設定されている。

出力ファイル名を設定すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まり、MCMCサンプリング後、トレース図が表示される（図C2）。

図C3

図C3のWindowを閉じると、判断カテゴリのパラメータの事後分布が表示される（図C4）。

図C4

図C4のフォームを閉じると、図C5のグラフが表示される。

図C5

「大きい」判断（カテゴリ３および４）と「小さい」判断（カテゴリ１および２）の境界（主観的等価点PSE）とカテゴリ判断の基準値C2およびC1との距離およびの差

の事後分布のグラフである。式（Diff）の分布は正の側によっており、確率であるので、判断の基準値C2は基準値C1より主観的等価点から離れていると言える。

図C5のWindowを閉じると、パラメータの事後分布が表示される（図C6）。

図C6

図C6のWindowを閉じると、図C7のグラフが表示される。

図C7

３本のＰＭ関数が描かれている。

心理測定関数（PM）が心理測定関数（PM1）の方に寄っていることは、図C4と対応している。

データ点は

とランダムな小さい値を加えて、データの小円の重なりを避けている。

図C7のWindowを閉じると、プログラムの実行終了である。

端末には、以下のようなメッセージが出力されている。

Results.txt was saved.

(stan) *****/pm4catfiles$

出力ファイルResults.txtを開くと、以下のような内容である。

Input Data File = Data4Cat.xlsx

Data...

180 1

180 1

220 4

・

・

・

160 1

190 1

220 4

Mean MCSE StdDev ... ESS_tail ESS_bulk/s R_hat

lp__ -129.747000 0.032340 1.421640 ... 2832.49 737.711 0.999813

v_sgm 0.144078 0.000228 0.013578 ... 2669.63 1360.110 1.000680

v_C1 1.922540 0.000347 0.019251 ... 2561.07 1172.710 1.000330

v_mu 0.082176 0.000321 0.017392 ... 2891.58 1122.670 1.000530

v_C2 0.133954 0.000356 0.021735 ... 3205.57 1411.190 0.999475

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

theta[180,4] 0.665212 0.000805 0.055996 ... 3741.32 1852.600 1.000400

sgm 14.407900 0.022801 1.357770 ... 2669.63 1360.110 1.000680

C1 192.254000 0.034715 1.925110 ... 2561.07 1172.710 1.000330

mu 200.471000 0.028651 1.913700 ... 3477.88 1699.250 1.000960

C2 213.867000 0.028595 2.064490 ... 3398.54 2001.070 1.000180

[729 rows x 11 columns]

入力ファイル名、入力データ、MCMCサンプリングの統計が出力されている。

リストC1　4件法用Stanスクリプト（pm4cat.stan）

data {

int N;

array[N] real St;

array[N] int Res;

}

parameters {

real<lower = .0> v_sgm;

real v_C1;

real<lower = 0.0> v_mu;

real<lower = 0.0> v_C2;

}

transformed parameters {

array[N] simplex[4] theta;

real sgm;

real C1;

real mu;

real C2;

sgm = v_sgm * 100.0 + 0.0001;

C1 = v_C1 * 100.0;

mu = C1 + v_mu * 100.0 + 0.0001;

C2 = mu + v_C2 * 100.0 + 0.0001;

for (i in 1:N) {

theta[i][2] = normal_cdf(St[i] | C1, sgm);

theta[i][3] = normal_cdf(St[i] | mu, sgm);

theta[i][4] = normal_cdf(St[i] | C2, sgm);

theta[i][1] = 1.0 - theta[i][2];

theta[i][2] = theta[i][2] - theta[i][3];

theta[i][3] = theta[i][3] - theta[i][4];

}

}

model {

v_sgm ~ exponential(0.1);

v_C1 ~ normal(0.0, 10.0);

v_mu ~ exponential(0.1);

v_C2 ~ exponential(0.1);

for (i in 1:N)

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

}

リストC2　リストC1のStanスクリプト（pm4cat.stan）を用いるPythonスクリプト（pm4cat.py）

from cmdstanpy import CmdStanModel

import pandas as pd

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

finnm = input('Input Data File (*.xlsx) = ')

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl, header = None)

data = data_xlsx.values

foutnm = input('Output Text File (*.txt) = ')

fout = open(foutnm, 'w')

fout.write('Input Data File = {}\n'.format(finnm))

print('data[0] =\n', data[0])

St = data[1:,1]

Res = data[1:,2]

fout.write('\nData...\n')

for v in zip(St,Res):

print(v)

fout.write(f'{v[0]:>10}{v[1]:>10}\n')

N = len(St)

def f_init():

p_mu = np.median(St)

p_sgm = np.std(St)

p_C1 = p_mu * 0.9

d_v = p_mu * 0.1

return {'v_sgm':p_sgm/100.0,

'v_C1':p_C1/100.0, 'v_mu': d_v/100.0,

'v_C2':d_v/100.0}

Data = {'N':N, 'St':St, 'Res':Res}

model = CmdStanModel(stan_file='pm4cat.stan')

fit = model.sample(data=Data, chains=4, inits=[f_init()]*4)

print(fit.summary())

fout.write(f'\n{fit.summary()}\n')

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

az.plot_trace(inf_data,

var_names = ['C1','mu','C2', 'sgm']) # トレースプロット

plt.tight_layout()

plt.show()

d_fm = fit.draws_pd() # Pandas DataFrame

med_C1 = np.median(d_fm['C1'])

med_mu = np.median(d_fm['mu'])

med_C2 = np.median(d_fm['C2'])

sb.kdeplot(d_fm['C1'], label = fr'C1(med={med_C1:.2f})')

sb.kdeplot(d_fm['mu'], label = fr'$\mu$(med={med_mu:.2f})')

sb.kdeplot(d_fm['C2'], label = fr'C2(med={med_C2:.2f})')

plt.xlabel('')

plt.legend()

plt.show()

diffC2muC1 = d_fm['C2'] + d_fm['C1'] - (d_fm['mu'] * 2)

prop = np.mean([1 if v > 0.0 else 0 for v in diffC2muC1])

med_diff = np.median(diffC2muC1)

sb.kdeplot(diffC2muC1)

plt.title(f'(C2-mu)-(mu-C1): med={med_diff:.2f}' +

f'\nP(>0) = {prop:.3f}', fontsize = 16)

plt.show()

med_sgm = np.median(d_fm['sgm'])

sb.kdeplot(d_fm['sgm'], label = fr'$\sigma$(med={med_sgm:.2f})')

plt.xlabel('')

plt.legend()

plt.show()

min_St = np.min(St)

max_St = np.max(St)

min_x = min_St - 0.05 * (max_St - min_St)

max_x = max_St + 0.05 * (max_St - min_St)

plt.axis([min_x, max_x, 0.0, 1.0])

phi = np.linspace(min_x, max_x, 1000)

psi1 = ss.norm.cdf(phi, med_C1, med_sgm)

psiMu = ss.norm.cdf(phi, med_mu, med_sgm)

psi2 = ss.norm.cdf(phi, med_C2, med_sgm)

p_st = []

p_res = []

for i in range(0, N):

p_st.append(St[i])

if Res[i] == 1:

p_res.append(0.1 * np.random.uniform())

elif Res[i] == 2:

p_res.append(1.0 / 3.0 - 0.05 + 0.1 * np.random.uniform())

elif Res[i] == 3:

p_res.append(2.0 / 3.0 - 0.05 + 0.1 * np.random.uniform())

else:

p_res.append(1.0 - 0.1 * np.random.uniform())

plt.title(("Psychometric functions JND = {0:.5}" +

"\nC1 = {1:.5} PSE = {2:.5} C2 = {3:.5}").

format(med_sgm * 0.67449, med_C1, med_mu, med_C2))

plt.plot(phi, psi1, 'b--', linewidth = 5)

plt.plot(phi, psiMu, 'g', linewidth = 5)

plt.plot(phi, psi2, 'b-.', linewidth = 5)

plt.plot(p_st, p_res, 'go')

plt.xlabel('$S_{comp}$', fontsize = 14)

plt.legend(['PM1', 'PM', 'PM2', 'Data'], loc = 2)

plt.show()

fout.close()

print()

print(foutnm + ' was saved.\n')

参考文献

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

Gescheider, G. A. (1997) Psychophysics: The fundamentals, Third edition. Mahwah: Lawrence Erlbaum Associates, Publishers.

岡本安晴（1994）「全系列データの最尤法による分析」金沢大学文学部論集 行動科学科篇、pp. 63-71。

岡本安晴（1995）「全系列法データの最尤法による分析について」金沢大学文学部論集 行動科学科篇、pp. 25-30。

Okamoto, Y. (2001). On Fechner’s problem. Studies and Essays, Behavioral Sciences and Philosophy: Kanazawa University, 1-25.

岡本安晴（2004）「３件法に対するBayesian Up-down法」日本行動計量学会第３２回大会発表論文抄録集、pp. 4-5。

岡本安晴（2005a）「恒常法データ（３件法）のベイズ的分析」日本心理学会第６９回大会発表論文集、p. 565。

Okamoto, Y. (2005b). Generalization of Fechner’s approach. Japan Women’s University Journal: Faculty of Integrated Arts and Social Sciences, 107-121.

岡本安晴（2006a）「適応的恒常法の可能性について」日本心理学会第７０回大会発表論文集、p. 559。

岡本安晴（2006b）「計量心理学」培風館

岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館

岡本安晴（2011）「上下法のベイズ的分析」基礎心理学研究、44−55。

Okamoto, Y. (2012). An experimental analysis of psychometric functions in a threshold discrimination task with four response categories. Japanese Psychological Research, 368-377.

岡本安晴（2014）「心理学データ分析と測定」勁草書房

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

岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版

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

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

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

Home