Psychometric関数（PF）のベイズ推定（PyStan 3 対応）
２件法、３件法、４件法データの分析
岡本安晴、2023.07
本ウェブサイトのプログラムは、利用者の責任において、自由に御利用下さい。
２件法より３件法、４件法の方が統計学的に望ましいことを示すシミュレーション結果を報告した（岡本、2017）。
確率モデルを用いると、弁別閾と判断の基準値を分離した分析が行えるというシミュレーションがある（岡本、1995）。
なお、psychometric関数の説明を＜岡本安晴「心理学データ分析と測定」、第11章＞で行っている。
Psychometric関数（精神測定関数／心理測定関数）のベイズ推定により、閾値（丁度可知差異；Just Noticeable Difference: JND）および主観的等価点（Point of Subjective Equality: PSE）を求めるPythonプログラムを作成してみた。プログラムの使用法などは、本ウェブサイトの後半で説明することとし、まず、背景となる考え方について説明する。
Psychometric関数（PF）を求めるpsychophysics（精神物理学／心理物理学）実験としては、反応が「大きい／小さい」など（本ウェブサイトでは、物理刺激の比較判断を「大きさ」の次元を例として用いている。他の感覚、例えば「重さ」であれば、「重い／軽い」判断となる）の２カテゴリ強制選択法（２件法）が標準であるが、「わからない／同じ」という判断を認める３件法もある。３件法については、データ分析を単純な平均値に基づいて行う場合には、感覚の弁別力（閾値）の推定値が判断の基準の影響を受けるということが知られているが（岡本、1995）、この問題は確率モデルを設定して分析を行えば回避できる。設定された確率モデルに対して最尤法を用いた場合（岡本、1994、1995）、ベイズ的方法を用いた場合（岡本、2005a）において、感覚の弁別力（閾値）と判断の基準が確率モデルによって分離されて適切に分析できることが報告されている。主観的等価点PSEを求めることを考慮して、２件法と３件法の混合型ともいえる４件法（「大きい」、「どちらかといえば大きい」、「どちらかといえば小さい」、「小さい」の４カテゴリを用いる方法）もある（Okamoto, 2012）。２件法、３件法、４件法、いずれの場合も分析の基本となるのは物理刺激と反応の関係を表すPFである。
本ウェブサイトで扱われているモデルは、標準刺激と比較刺激の心理物理学的関係についてのものであり、恒常法や全系列法および上下法などの適応的方法の区別なく適用することができる。観察者の判断カテゴリ数により、カテゴリ数に対応したモデル構成になるので、以下に２件法、３件法および４件法の場合について順番に論じる。また、実験デザインの立案において、比較刺激の間隔と範囲および総試行数はパラメータの推定値を求めるときに重要な問題であるが、このことについて３つのモデルの説明の後で論じる。
まず、２件法の場合について考える。比較刺激
の方が標準刺激
より大きいと判断される（
と表記）確率
をの関数として表したものがPF（psychometric function）であるが、心理学ではPFとして累積正規分布を用いることが多い。累積標準正規分布を
で表すと
と書ける。
および
は、標準刺激によって決まる値であり、１つの固定したの値について考えるときは、PFはの関数となる。PFをの関数とみてを
で表すと、式（１）は
となる。２件法なので
である。
ここで、は主観的等価点
を表し
である。
弁別閾（difference threshold; Gescheider,
1997）を表す丁度可知差異（JND; 岡本（2014）、pp.
168-169参照）としては、の
倍を用いる。ここで、は、確率0.75に対応する標準正規分布における値である。すなわち、
である。プログラムでは
として計算している。すなわち、
である。
第
試行において提示された比較刺激を
、それに対する判断（反応）を
で表す。ただし、
とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式
で与えられる。ここで、
とおいた。
ベイズ的分析では、式（２）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。Python用のStanとしてPyStanがあるが、PyStanについては「岡本（2019）、pp. 178-210」で説明している。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。
ベイズ分析では分析結果は事後分布として与えられる。主観的等価点などの点推定値は、事後分布から求めることになる。分布の代表値としては平均値、中央値、モードがよく知られているが、事後分布からの点推定値として、これらの３つの推定値をシミュレーションで比較した岡本（2022）では、モード（MAP推定値）が真値との関係で薦められている。
３件法データの分析モデルは、２件法の拡張として考える。３件法では、判断は「標準刺激より比較刺激の方が小さい；
」、「とは同じ／どちらが大きいかわからない；
」、「よりの方が大きい；」の３つのカテゴリから選ばれる。これらを、
「」と「あるいは（
と表記）」、
および
「あるいは（
と表記）」と「」
の２通りに分けて考える。
「」と「」の組み合わせに対するＰＦを
とおく。
は、判断「」と「」を分ける判断の基準値である。
「」と「」の組み合わせに対するＰＦを
とおく。
は、判断「」と「」を分ける判断の基準値である。
判断「」の確率は、次式
で与えられる。また、
である。
JNDは、２件法の場合と同様に
で与えられる。
PSEについては、２つの判断の基準値
と
の中央に主観的等価点PSEがあると考えられるときは、次式
でPSEが推定されるが、Okamoto (2012)は４件法のデータにおいて式（３）が成り立たない例を報告している。
いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、
とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式
で与えられる。ここで、
とおいた。
ベイズ的分析では、式（4）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。Python用のStanとしてPyStanがあるが、PyStanについては「岡本（2019）、pp. 178-210」で説明している。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。
４件法では、判断は「標準刺激より比較刺激の方が小さい；
」、「どちらかといえばよりの方が小さい；」、「どちらかといえばよりの方が大きい；」、「標準刺激より比較刺激の方が大きい；
」の４カテゴリから１つを選ぶという形で行われる。これらの判断を以下の３つの２件法データの形式に分けて考える。
（A）「」と
「またはまたは（
と表記）」
（B）「または（と表記）」と
「または（と表記）」
（C）「またはまたは（
と表記）」と
「」
形式（A）によるデータに対するＰＦを
とおく。形式（B）と（C）に対しても同様に
とおく。
このとき、
である。
上記モデルにおいて、PSEはによって与えられる。すなわち、
である。ここで、
は仮定されていないことに注意。式（５）が成り立つときは、２つの判断の基準値はPSEから等距離に設定されていることになる。Okamoto (2012)では、４件法データの分析において、式（５）の成り立たない例が報告されている。
JNDは、２件法および３件法の場合と同様に
で与えられる（岡本、2014、 p. 169）。
いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、
とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式
で与えられる。ここで、
とおいた。
ベイズ的分析では、式（６）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2019）、pp. 169-177」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、ベイズ分析のライブラリStanを用いた。Python用のStanとしてPyStanがあるが、PyStanについては「岡本（2019）、pp. 178-210」で説明している。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。
総試行数については、上下法を用いると少ない試行数で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)で論じている。
Gescheider, G. A. (1997) Psychophysics: The fundamentals, Third edition. Mahwah: Lawrence Erlbaum Associates, Publishers.
岡本安晴（1994）「全系列データの最尤法による分析」金沢大学文学部論集 行動科学科篇、pp. 63-71。
岡本安晴（1995）「全系列法データの最尤法による分析について」金沢大学文学部論集 行動科学科篇、pp. 25-30。
岡本安晴（2004）「３件法に対するBayesian Up-down法」日本行動計量学会第３２回大会発表論文抄録集、pp. 4-5。
岡本安晴（2005a）「恒常法データ（３件法）のベイズ的分析」日本心理学会第６９回大会発表論文集、p. 565。
岡本安晴（2006a）「適応的恒常法の可能性について」日本心理学会第７０回大会発表論文集、p. 559。
岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館
岡本安晴（2011）「上下法のベイズ的分析」基礎心理学研究、44−55。
岡本安晴（2017）「恒常法における反応カテゴリ数の影響」基礎心理学研究、85−91．
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版
岡本安晴（2019）「多件法による上下法」基礎心理学研究、38、90-104．
岡本安晴（2022）「ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――」基礎心理学研究、41、1-7.
本ウェブサイトのスクリプトをWindowsで実行するときは、WindowsにUbuntuをインストールしてUbuntu上で実行する。WindowsにUbuntuをインストールしてのPythonの利用については、このウェブサイトで説明しているので参照されたい。Ubuntu/WSL/WindowsにインストールしたAnacondaのPython 3.11仮想環境でPyStan 3を利用するという方法だが、WindowsのプログラムとUbuntu/WSL/Winodwsのプログラムはシームレスに実行できるので使い易い。
プログラム（スクリプト）について、２件法データのプログラム、３件法データのプログラム、４件法データのプログラムと順番に説明する。
２件法の場合のベイズ分析用Stanスクリプトをリスト２Aのように用意した。
リスト２A 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]);
}
リスト２AのStanスクリプトを利用するPythonスクリプトをリスト２Bのように用意した。スクリプトファイルなどは、まとめて圧縮ファイルpm2catfiles.zipとした。ダウンロード展開すれば自由に利用できる。
リスト２BのスクリプトをAnaconda/Ubuntu/WSL/windowsのPython 3.11仮想環境において次のコマンドで実行した。
(py311) ・・・/pm2catfiles$ python
pm2cat.py
実行すると、まず入力データファイル名の設定が求められる。
(py311) ・・・/jnd_PS3/pm2catfiles$ python
pm2cat.py
Input Data
File (*.xlsx) =
入力データファイルはExcelの形式で、図２Aのように用意する。
・
・
・
図２A
第１行目に変数名を設定して、２行目以降にデータ値を入力する。
第１列目はデータのID用である。プログラムでは無視されるので任意の値あるいは文字列でよい。
第２列目に比較刺激値、第３列目に判断カテゴリ値を入力する。判断カテゴリ値は、「比較刺激値の方が小さい」というときは１、「比較刺激の方が大きい」というときは２を設定する。
入力データファイルとして、次の例では、図２AのファイルData2Cat.xlsxを設定している。
(py311) ・・・/pm2catfiles$ python
pm2cat.py
Input Data
File (*.xlsx) = Data2Cat.xlsx
Output Text
File (*.txt) = Results.txt
入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。任意のテキストファイル名を設定する。
上の例では、Results.txtが設定されている。出力ファイル名を設定すると、計算が始まる。
MCMCサンプリングが終了すると、まずトレース図（図２B）が表示される。
図２B
図２Bのフォームを閉じると、パラメータの事後分布が表示される（図２C）。
図２C
図２Cのフォームを閉じると、パラメータの事後分布が表示される（図２D）。
図２D
図２Dのフォームを閉じると、図２Eのフォームが提示される。
図２E
青の曲線がPFである。このPFはとそれぞれの事後分布のMAP推定値を点推定値として描かれている。主観的等価点は
で算出されている。
緑色の小円はデータ
である。ただし、小円が重ならないように、ランダムな小さい値を加えて表示してある。すなわち、プロットされている点を
としたとき、
ここで、
はランダムな小さい正数である。
図２Eのフォームを閉じると、プログラムの実行終了である。端末には以下の出力が表示されている。
Results.txt was saved.
(py311) ・・・/pm2catfiles$
出力ファイルResults.txtを開くと、以下のような内容である。
Input Data File = Data2Cat.xlsx
Summary:
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu 201.490 2.005 197.662 205.194 0.015 0.011 17499.0 13219.0 1.0
sgm 14.879 1.992 11.181 18.549 0.016 0.011 16531.0 13096.0 1.0
入力データファイル名と、MCMCサンプリングの統計量が出力されている。
リスト２B リスト２AのStanスクリプト（pm2cat.stan）を利用するPythonスクリプト（pm2cat.py）
import numpy as np
import pandas as pd
import stan
import arviz as az
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
""" MAP推定値とMAP推定値における推定確率密度関数の値を求める
samples: MAP推定値を求めるサンプル
a: サンプルの左右a/2を除く。デフォルト値=0.05
n_points: KDEによる確率密度関数値を求める点の数。
デフォルト値=1000
"""
import numpy as np
import scipy.stats as ss
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]
#
# MAP推定値とMAP推定値における確率密度関数値
#
return MAP_Est,
est_pdf[map_idx]
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}
with open('pm2cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_samples = 5000,
num_chains = 4,
init = [f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['mu', 'sgm'])
print('Summary =\n',smry)
fout.write(f'\nSummary:\n{smry}\n')
az.plot_trace(i_data, var_names =
['mu', 'sgm'])
plt.show()
d_fm = fit.to_frame()
map_mu = calcMAPEst(d_fm['mu'])[0]
map_sgm = calcMAPEst(d_fm['sgm'])[0]
sb.kdeplot(d_fm['mu'])
plt.title(f'Posterior distribution of
$\mu$, MAP = {map_mu:.2f}')
plt.xlabel('$\mu$', fontsize = 14)
plt.show()
sb.kdeplot(d_fm['sgm'])
plt.title(f'Posterior distribution of
$\sigma$, MAP = {map_sgm:.2f}')
plt.xlabel('$\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, map_mu,
map_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(map_mu, map_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()
３件法の場合のベイズ分析用Stanスクリプトをリスト３Aのように用意した。
リスト３A ３件法用ベイズ分析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]);
}
リスト３AのStanスクリプトを利用するPythonスクリプト（pm3cat.py）をリスト３Bのように用意した。スクリプトファイルなどは、まとめて圧縮ファイルpm3catfiles.zipとしている。ダウンロード展開して自由に利用されたい。
リスト３Bのスクリプトを、Anaconda/Ubuntu/WSL/WindowsのPython3.11仮想環境において次のコマンドにより実行した。
(py311) ・・・/pm3catfiles$ python
pm3cat.py
実行が始まると、まず入力データファイル名の設定が求められる。
Input Data
File (*.xlsx) =
入力データファイルは、図３AのようにExcelファイルの形式で用意する。
・
・
・
図３A
第１行目に変数名を書く。データ値は２行目以降に設定する。
第１列目にはデータ識別用の値を設定する。この列の値は、プログラムでは無視されるので、数値でも文字列でもよい。
第２列目に、提示刺激値、第３列目に判断カテゴリ値を設定する。判断カテゴリ値は、比較刺激の方が小さいなら１、同じなら２、大きいなら３を設定する。
次の例では、図３Aのファイル名Data3Cat.xlsxが設定されている。
(py311) ・・・/pm3catfiles$ python pm3cat.py
Input Data File (*.xlsx) =
Data3Cat.xlsx
Output Text File (*.txt) = Results.txt
入力データファイル名を設定すると、次に出力ファイル名の設定が求められるので、適当なテキストファイル名を設定する。上の例では、Data3Cat.xlsxが設定されている。
出力ファイル名を設定すると、計算が始まる。
MCMCサンプリングが終了すると、トレース図が表示される（図３B）。
図３B
図３Bのフォームを閉じると、図３Cのグラフが表示される。
図３C
パラメータC1とC2の事後分布と、次式（３）で与えられるC1とC2の中点midCsの事後分布が描かれている。２つの判断の基準値との中央に主観的等価点PSEがあると考えられるときは、次式
でPSEが推定される。
図３Cのフォームを閉じると、パラメータの事後分布が表示される（図３D）。
図３D
図３Dのフォームを閉じると、図３Eのグラフが表示される。
図３E
２つの心理測定関数
と
のグラフ、および
のグラフである。
感覚の弁別力は心理測定関数の傾きで表され、これはパラメータの逆数である。これに対して、判断の基準の厳しさは心理測定関数の平行移動の度合いで表され、心理測定関数間の横方向の距離が小さいと判断の基準は厳しいと言える。確率モデルを用いた分析により、感覚の弁別力と判断の厳しさを分離することができるが、岡本（1995）は全系列法のシミュレーションによってこのことを示している。
青の小円はデータ点を表すが、データ点
は
とランダムな小さな値を加えて、データの小円の重なりを避けている。
図３Eのフォームを閉じると、プログラムの実行終了である。
端末には、以下のメッセージが出力されている。
Results.txt was saved.
(py311) ・・・/pm3catfiles$
出力ファイルResults.txtの内容は以下のようである。
Input Data File = Data3Cat.xlsx
Data...
180
1
180
1
220
3
・
・
・
160
1
190
1
220
3
Summary:
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
sgm 14.219 1.419 11.714 16.994 0.012 0.009 13204.0 13014.0 1.0
C[0] 192.368 1.937 188.721 195.993 0.018 0.012 12105.0 11402.0 1.0
C[1] 213.753 2.073 209.823 217.511 0.014 0.010 22136.0 16050.0 1.0
読み込まれたデータの出力と、MCMCサンプリングの統計量が出力されている。
リスト３B リスト３AのStanスクリプトを利用するPythonスクリプト（pm3Cat.py）
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import stan
import arviz as az
import scipy.special as ssp
import seaborn as sb
import scipy.stats as ss
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
""" MAP推定値とMAP推定値における推定確率密度関数の値を求める
samples: MAP推定値を求めるサンプル
a: サンプルの左右a/2を除く。デフォルト値=0.05
n_points: KDEによる確率密度関数値を求める点の数。
デフォルト値=1000
"""
import numpy as np
import scipy.stats as ss
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] # MAP推定値とMAP推定値における確率密度関数値
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 #ssp.expit(p_mu * 0.9)
p_C2 = p_mu * 1.1 #ssp.expit(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}
with open('pm3cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_samples = 5000,
num_chains = 4,
init = [f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['sgm', 'C'])
print('Summary =\n',smry)
fout.write(f'\nSummary:\n{smry}\n')
az.plot_trace(i_data, var_names =
['C', 'sgm'])
plt.show()
d_fm = fit.to_frame()
C1s = d_fm['C.1']
C2s = d_fm['C.2']
midCs = (C1s + C2s) / 2.0
map_C1 = calcMAPEst(C1s)[0]
map_C2 = calcMAPEst(C2s)[0]
map_midCs = calcMAPEst(midCs)[0]
sb.kdeplot(C1s, label =
f'C1(MAP={map_C1:.2f})')
sb.kdeplot(C2s, label =
f'C2(MAP={map_C2:.2f})')
sb.kdeplot(midCs, ls = '--', label =
f'midCs(MAP={map_midCs:.2f})')
plt.xlabel('')
plt.legend()
plt.show()
map_sgm = calcMAPEst(d_fm['sgm'])[0]
sb.kdeplot(d_fm['sgm'], label =
r'$\sigma$(MAP={0:.2f})'.format(map_sgm))
plt.xlabel('$\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, map_C1,
map_sgm)
psi2 = ss.norm.cdf(phi, map_C2,
map_sgm)
psimid = ss.norm.cdf(phi, map_midCs,
map_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(map_C1, map_C2, map_sgm * 0.67449, map_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()
４件法データのベイズ分析用Stanスクリプトをリスト４Aのように用意した。
リスト４A ４件法データのベイズ分析用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]);
}
リスト４AのStanスクリプトを利用するPythonスクリプトをリスト４Bのように用意した。これらのファイルは、まとめて圧縮ファイルpm4catfiles.zipとした。ダウンロード展開すれば、自由に利用できる。
リスト４Bのスクリプトを、Anaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境で次のようにPythonコマンドで実行すると、入力ファイル名の設定が求められる。
(py311) ・・・/pm4catfiles$ python pm4cat.py
Input Data File (*.xlsx) =
入力データファイルは、Excelファイルの形式で用意する（図４A）。
・
・
・
図４A
第１行目に変数名を設定する。データ値は、２行目以降に書く。
第１列目はデータの識別用で、プログラム内では無視される。この識別用変数の値は、数値でも文字列でもよい。
第２列目は、提示刺激値である。
第３列目には観察者の判断値（反応）を設定する。比較刺激の方が小さいなら１、たぶん小さいなら２、たぶん大きいなら３、大きいなら４を設定する。
入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。
(py311) ・・・/pm4catfiles$ python pm4cat.py
Input Data File (*.xlsx) =
Data4Cat.xlsx
Output Text File (*.txt) = Results.txt
上の例では、入力データファイル名として図４Aのファイル名Data4Cat.xlsxが設定されている。出力ファイル名は、任意のテキストファイル名である。上の例では、Results.txtと設定されている。
出力ファイル名を設定すると、計算が始まる。MCMCサンプリングが終了すると、まずトレース図が表示される（図４B）。
図４B
図４Bのフォームを閉じると、判断カテゴリのパラメータの事後分布が表示される（図４C）。
図４C
図４Cのフォームを閉じると、図４Dのグラフが表示される。
図４D
「大きい」判断（カテゴリ３および４）と「小さい」判断（カテゴリ１および２）の境界（主観的等価点PSE）とカテゴリ判断の基準値C2およびC1との距離
および
の差
の事後分布のグラフである。式（Diff）の分布は正の側によっており、確率
であるので、判断の基準値C2は基準値C1より主観的等価点
から離れていると言える。
図４Dのフォームを閉じると、パラメータの事後分布が表示される（図４E）。
図４E
図４Eのフォームを閉じると、図４Fのグラフが表示される。
図４F
３本のＰＭ関数が描かれている。
心理測定関数（PM）が心理測定関数（PM1）の方に寄っていることは、図４Dと対応している。
データ点は
とランダムな小さい値を加えて、データの小円の重なりを避けている。
図４Fのフォームを閉じると、プログラムの実行終了である。
端末には、以下のようなメッセージが出力されている。
Results.txt was saved.
(py311) ・・・/jnd_PS3/pm4catfiles$
出力ファイルResults.txtを開くと以下のような内容である。
Input Data File = Data4Cat.xlsx
Data...
180
1
180
1
220
4
・
・
・
160
1
190
1
220
4
Summary:
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
C1 192.246 1.921 188.570 195.751 0.015 0.011 16108.0 14203.0 1.0
mu 200.479 1.908 196.917 204.068 0.014 0.010 19729.0 16097.0 1.0
C2 213.806 2.074 209.982 217.725 0.013 0.009 24576.0 17108.0 1.0
入力ファイル名、入力データ、MCMCサンプリングの統計量が出力されている。
リスト４B リスト４AのStanスクリプトを利用するPythonスクリプト（pm4cat.py）
import stan
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
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
""" MAP推定値とMAP推定値における推定確率密度関数の値を求める
samples: MAP推定値を求めるサンプル
a: サンプルの左右a/2を除く。デフォルト値=0.05
n_points: KDEによる確率密度関数値を求める点の数。
デフォルト値=1000
"""
import numpy as np
import scipy.stats as ss
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] # MAP推定値とMAP推定値における確率密度関数値
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}
with open('pm4cat.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = Data)
fit = sm.sample(num_samples = 5000,
num_chains = 4,
init = [f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data, var_names =
['C1', 'mu', 'C2'])
print('Summary =\n',smry)
fout.write(f'\nSummary:\n{smry}\n')
az.plot_trace(i_data, var_names =
['C1', 'mu', 'C2', 'sgm'])
plt.show()
d_fm = fit.to_frame()
map_C1 = calcMAPEst(d_fm['C1'])[0]
map_mu = calcMAPEst(d_fm['mu'])[0]
map_C2 = calcMAPEst(d_fm['C2'])[0]
sb.kdeplot(d_fm['C1'], label = f'C1(MAP={map_C1:.2f})')
sb.kdeplot(d_fm['mu'], label =
f'$\mu$(MAP={map_mu:.2f})')
sb.kdeplot(d_fm['C2'], label =
f'C2(MAP={map_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])
map_diff = calcMAPEst(diffC2muC1)[0]
sb.kdeplot(diffC2muC1)
plt.title(f'(C2-mu)-(mu-C1):
MAP={map_diff:.2f}' +
f'\nP(>0) = {prop:.3f}', fontsize = 16)
plt.show()
map_sgm = calcMAPEst(d_fm['sgm'])[0]
sb.kdeplot(d_fm['sgm'], label =
f'$\sigma$(MAP={map_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, map_C1,
map_sgm)
psiMu = ss.norm.cdf(phi, map_mu,
map_sgm)
psi2 = ss.norm.cdf(phi, map_C2,
map_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(map_sgm * 0.67449, map_C1, map_mu, map_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')