２肢強制選択課題（信号検出理論：SDT）
評定法データのベイズ分析（CmdStanPy対応）
２肢強制選択（2 alternative-forced-choice: 2AFC）課題の場合、ノイズ刺激
とシグナル刺激
が各試行において提示されるが、空間的位置を変えて（右と左、あるいは上と下など）同時に提示される場合、あるいは提示の時間順序を変えて（第１インタバルと第２インタバル）提示される場合がある。統計モデルを考えるときは、空間順序あるいは時間順序における２つの刺激提示法をともに＜、＞あるいは＜、＞で表すことができる。＜、＞は、刺激が空間位置で提示されているときは、例えば、左側にノイズ刺激、右側にシグナル刺激が提示されることを表し、時間順序のときは第１インタバルでノイズ刺激が提示され、第２インタバルでシグナル刺激が提示されることを表す。＜、＞のときは、順序が逆になる。以後の説明では、例えば＜、＞は、は第１位置、は第２位置に提示される条件を表すということにする。刺激提示条件別に判断（反応）を分類すると表1のようになる。
|
|
シグナル刺激の提示位置として判断された位置 |
|
|
|
第１位置 |
第２位置 |
|
＜ |
誤反応 |
正反応 |
|
＜ |
正反応 |
誤反応 |
いま、1AFC課題のモデルと同じく、ノイズ刺激の感覚
が平均０、分散１の正規分布に従い、シグナル刺激の感覚
が平均
、分散
の正規分布に従うものとする。すなわち、
、 
である。
提示条件＜、＞および＜、＞において恒常誤差
および位置選択バイアス
を考える。恒常誤差は、時間誤差あるいは空間誤差を表す。位置選択バイアスは、感覚の違いの弁別とは別に第１位置あるいは第２位置の選択を好む傾向を表す。これらのバイアスをまとめて
とおき、第２位置を基準にした第１位置へのバイアス効果として表すと、感覚およびの分布は以下のようになる。
刺激提示が＜、＞のとき、
、
で与えられ、
となる。
刺激提示が＜、＞のとき、
、 
で与えられ、
となる。
表１の実験課題の場合は、を求めることができない。「
」の仮定（等分散モデル）のもとで、分析が行われる。表１のデザインの下で評定法を用いることができる（岡本、2025）。「」の仮定を置かない（不等分散モデルと呼ばれているが、等分散を仮定しないモデルということであり、データが等分散であってもよい。等分散であるかどうかわからないときは、不等分散モデルを用いればよい）ときは、表１のデザインを拡張すればよい。以下に、共通事項について説明した後、等分散モデルで評定を用いない場合、等分散モデルで奇数個の評定カテゴリを用いる場合、等分散モデルで偶数個の評定カテゴリを用いる場合、不等分散モデルの場合のベイズ分析について説明する。なお、CmdStanPyの準備の簡単な説明を、このウェブサイトで行った。ファイルは、2AFCfiles.zipにまとめた。
評定法における判断カテゴリの境界値（基準値）を
とおく。すなわち、シグナル刺激の感覚とノイズ刺激の感覚の差
が
との間にあるときカテゴリ
の判断になると考える。
のとき
観察者の判断
カテゴリ
である。
判断カテゴリの総数が４カテゴリのとき、例えば判断カテゴリは次のようになる。ここで、の値は、シグナル刺激の感覚が強いほどカテゴリの値は大きい値になるが、＜、＞と＜、＞ではシグナル刺激の提示位置が異なることに注意。
＜、＞のとき、
カテゴリ１： シグナルは第１位置である。
カテゴリ２： シグナルは多分第１位置である。
カテゴリ３： シグナルは多分第２位置である。
カテゴリ４： シグナルは第２位置である。
＜、＞のとき、
カテゴリ１： シグナルは第２位置である。
カテゴリ２： シグナルは多分第２位置である。
カテゴリ３： シグナルは多分第１位置である。
カテゴリ４： シグナルは第１位置である。
いま、判断カテゴリの総数を
とおいたとき、カテゴリ境界を
とおく。累積標準正規分布関数
に対して
、 
とおいたとき、判断カテゴリがカテゴリである確率は次のように与えられる。
＜、＞のとき、カテゴリである確率
は、式（１）より
となる。
＜、＞のとき、カテゴリである確率
は、式（２）より
となる。
表2 刺激提示条件に対するカテゴリ評定判断の頻度。
|
|
カテゴリ１ |
・・・ |
カテゴリk |
・・・ |
カテゴリK |
|
＜ |
|
・・・ |
|
・・・ |
|
|
＜ |
|
・・・ |
|
・・・ |
|
各刺激提示条件に対するカテゴリ判断の頻度が表2のようであるとき、尤度関数は次式で与えられる。
モデルのパラメータの数は
であるが、式（３）および式（４）の関数の引数の分数形から以下のようなパラメータ値の不定性が認められる。
分数の分母子に同じ数を掛けても値は変わらないので、分母の値を
とおいても一般性は失われない。これは、等分散モデル（
）あるいはより一般的には感覚の単位を
であるように設定するということである。
標準の２AFC課題について議論を続ける。
この場合は、等分散モデルで考えるが、感覚の弁別力は
と表記される。
等分散を仮定しない場合は、弁別力は
で表される。
である（岡本、2025）。
また、分子における
と
に同じ定数を加えても
であるので、パラメータ値を決めるためには制約条件が必要である。感覚の差
の判断が原点に関して対称であるとすると、これはカテゴリ境界に次の仮定をおくことである。
カテゴリ数が偶数
（
自然数）のとき
このとき、モデルのパラメータの自由度は
となる。
カテゴリ数が奇数
（ 自然数）のとき
このとき、モデルのパラメータの自由度は
となる。
表2のデータの自由度
は、表2のデータを刺激提示条件における頻度の比率データと考えれば
となる。パラメータの推定が不定とならないためには、パラメータの自由度
はデータの自由度以下でなければならない。すなわち、
であるが、この条件は満たされている。
等分散を仮定しない場合については、不等分散モデルの節で論じる。
表１のデータのベイズ分析のStanスクリプトをリストA1のように用意した。リストA1のスクリプトを用いてベイズ分析するPythonスクリプトをリストA2のように用意した。リストA1のファイル（SDT2AFC2Cat.stan）、リストA2のファイル（SDT2AFC2Cat.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移動して、CmdStanPyのインストールされた環境において、次のコマンドを実行する。
(stan)
*****/sdt_simple$ python SDT2AFC2Cat.py
実行すると、入力データファイル名と出力ファイル名の設定が求められる。
(stan) *****/sdt_simple$ python
SDT2AFC2Cat.py
入力データファイル名(*.csv) = Data2Cat.csv
出力テキストファイル名(*.txt) = Results.txt
出力データファイルは任意のテキストファイル名（ファイル拡張子.txt）であるが、入力データファイルはCSV形式のファイルとして図A１の様式で用意する。
図A１ データファイルData2Cat.csv
第１行目は変数名の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第２行目、第３行目とも、第２列目は「シグナル刺激は第１位置である（例えば、シグナル刺激は左である）」と判断された度数、第３列目は「シグナル刺激は第２位置である（例えば、シグナル刺激は右である）」と判断された度数である。データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
入力データファイル名と出力テキストファイル名を設定すると、ファイルが読み込まれ、Stanのコンパイルが始まる。コンパイルには多少の時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、まず
（）の事後分布が表示される（図A２）。
図A2
図A2のWindowを閉じると、図A3のグラフが表示される。
図A3
図A3のWindowを閉じると、図A4のグラフが表示される。
図A4
図A1の表のデータを比率に直して表示したもの（正方形）と、モデル予測値の事後分布の中央値（小円）のグラフである。正方形（データ）と小円（モデル）が重なっているので、モデルはデータによく適合していると言える。
図A4のWindowを閉じると、実行終了である。端末には、以下のメッセージが表示されている。
Results.txt was saved.
(stan) *****/sdt_simple$
ファイルResults.txtを開くと、以下のようである。
Input Data File = Data2Cat.csv
<N, S>: [30, 70]
<S, N>: [84, 16]
Mean
MCSE StdDev ...
ESS_tail ESS_bulk/s R_hat
lp__ -106.038000 0.023853 0.963004 ... 2150.88 56747.7 1.00161
mu_s 1.082450 0.002467 0.139474 ... 2555.86 106729.0 1.00049
tau
0.332082 0.002353 0.141214 ... 2626.98 120047.0 1.00110
p_NS[1] 0.299458 0.000708 0.045674 ... 2611.75 137324.0 1.00131
p_NS[2] 0.700542 0.000708 0.045674 ... 2611.75 137324.0 1.00132
p_SN[1] 0.161207 0.000665 0.035739 ... 2398.93 98051.0 1.00053
p_SN[2] 0.838793 0.000665 0.035739 ... 2398.93 98051.0 1.00053
[7 rows x 11 columns]
d':
Med. = 1.083, 90%CI = [0.852, 1.315
tau:
Med. = 0.334, 90%CI ~ [0.100, 0.562]
リストA1 ２AFCのStanスクリプト（SDT2AFC2Cat.stan）
data {
array[2] int n_NS;
array[2] int n_SN;
}
parameters {
real mu_s;
real tau;
}
transformed parameters {
simplex[2] p_NS;
simplex[2] p_SN;
p_NS[1] = Phi(-(mu_s -
tau)/sqrt(2.0));
p_NS[2] = 1 - p_NS[1];
p_SN[1] = Phi(-(mu_s +
tau)/sqrt(2.0));
p_SN[2] = 1 - p_SN[1];
}
model {
mu_s ~ normal(0.0, 100.0);
tau ~ normal(0.0, 100.0);
n_NS ~ multinomial(p_NS);
n_SN ~ multinomial(p_SN);
}
リストA2 リストA1のStanスクリプトを用いるPythonスクリプト（SDT2AFC2Cat.py）
from cmdstanpy import CmdStanModel
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import csv
import seaborn as sb
flnm = input('入力データファイル名(*.csv) = ')
with open(flnm, 'r') as f:
data = [v for v in
csv.reader(f)]
fout_nm = input('出力テキストファイル名(*.txt) = ')
fout = open(fout_nm, 'w')
fout.write('Input Data File =
{}\n'.format(flnm))
Freq_NS = []
Freq_NS.append(int(data[1][1]))
Freq_NS.append(int(data[1][2]))
Freq_SN = []
Freq_SN.append(int(data[2][1]))
Freq_SN.append(int(data[2][2]))
print(Freq_NS)
print(Freq_SN)
fout.write('<N, S>:
{}'.format(Freq_NS))
fout.write('\n<S, N>:
{}'.format(Freq_SN))
TN_NS = sum(Freq_NS)
Rating_NS = []
Rating_NS.append(Freq_NS[0])
Rating_NS.append(Freq_NS[1])
TN_SN = sum(Freq_SN)
Rating_SN = []
Rating_SN.append(TN_SN - Freq_SN[0])
Rating_SN.append(TN_SN - Freq_SN[1])
Data = {'n_NS': Rating_NS, 'n_SN':
Rating_SN}
model =
CmdStanModel(stan_file='SDT2AFC2Cat.stan')
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
fit = fit.draws_pd()
mu = fit['mu_s']
tau = fit['tau']
p_NS = np.array([fit['p_NS[1]'],
fit['p_NS[2]']]).T
p_SN = np.array([fit['p_SN[1]'],
fit['p_SN[2]']]).T
mu_L05 = np.percentile(mu, 5)
mu_med = np.percentile(mu, 50)
mu_U95 = np.percentile(mu, 95)
fout.write("\nd': \nMed. =
{0:<.3f}, 90%CI = [{1:<.3f},
{2:<.3f}".
format(mu_med, mu_L05, mu_U95))
tau_L05 = np.percentile(tau, 5)
tau_med = np.percentile(tau, 50)
tau_U95 = np.percentile(tau, 95)
fout.write('\ntau: \nMed. =
{0:<.3f}, 90%CI ~ [{1:<.3f},
{2:<.3f}]'.
format(tau_med, tau_L05, tau_U95))
p1 = np.median(p_NS, axis = 0) # [P(1|NS), P(2|NS)]
p2 = np.median(p_SN, axis = 0)
p2 = [p2[1], p2[0]]
# [P(2|SN),
P(1|SN)]
sb.kdeplot(mu)
plt.xlabel(r"d'($\mu_S$)",
fontsize = 14)
plt.title(r"Posterior Distribution
of d'($\mu_S$)" + \
'\nMed. = {0:<.3f}, 90%CI
= [{1:M<.3f}, {2:<.3f}]'.
format(mu_med, mu_L05, mu_U95), fontsize = 16)
plt.show()
sb.kdeplot(tau)
plt.xlabel(r'$\tau$', fontsize = 18)
plt.title(r'Posterior Distribution of
$\tau$' + \
'\nMed. = {0:<.3f}, 90%CI
= [{1:<.3f}, {2:<.3f}]'.
format(tau_med, tau_L05, tau_U95), fontsize = 16)
plt.show()
xcat = [1, 2]
xlabels = ['Cat-1', 'Cat-2']
y1 = [Freq_NS[0]/TN_NS,
Freq_NS[1]/TN_NS]
y2 = [Freq_SN[0]/TN_SN,
Freq_SN[1]/TN_SN]
plt.plot(xcat, y1, c='b', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<N, S>')
plt.plot(xcat, p1, c='g', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth = 3, label = 'Model/<N, S>')
plt.plot(xcat, y2, c='r', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<S, N>')
plt.plot(xcat, p2, c='m', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth=3, label = 'Model/<S, N>')
plt.xticks(xcat, xlabels, fontsize =
14)
plt.xlim(0.8, 2.2)
plt.xlabel('Rating Category', fontsize
= 14)
plt.ylabel('Probability/Proportion',
fontsize = 14)
plt.yticks(np.arange(0,1.1,0.2))
plt.title('Rating in Conditions <N,
S> and <S, N>', fontsize = 18)
plt.legend(fontsize=18, loc='center')
plt.show()
fout.close()
print('\n' + fout_nm + ' was
saved.\n')
判断のカテゴリ数Kが奇数の場合のStanスクリプトをリストB１のように用意した。リストB1のStanスクリプトを使ってベイズ分析を行うPythonスクリプトをリストB2のように用意した。リストB1のファイル（SDT2AFCOddCat.stan）、リストB2のファイル（SDT2AFCOddCat.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（ファイル）を移し、CmdStanPyのインストールされている環境において次のコマンドを実行する。なお、CmdStanPy準備の簡単な説明を、このウェブサイトで行った。
(stan)
*****/sdt_oddcat$ python SDT2AFCOddCat.py
上のコマンドを実行すると、入力データファイル名と出力データファイル名の設定が求められる。
(stan) *****/sdt_oddcat$ python
SDT2AFCOddCat.py
入力データファイル(*.csv) = Data5Cat.csv
出力テキストファイル名 = Results.txt
入力データファイルは、CSV形式のファイルとして用意する（図B1）。Excelの場合、ファイル拡張子として「.csv」を選ぶと、CSV形式のファイルとして保存される。
図B1
図B1は、判断カテゴリ数が５の場合の例である。第１行目は、変数名（カテゴリ名）の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「多分、シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第４列目に「わからない」と判断された度数、第５列目に「たぶん、シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数、第６列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。
第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「多分、シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第４列目に「わからない」と判断された度数、第５列目に「たぶん、シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数、第６列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。
データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
入力データファイル名の設定に続いて、出力ファイル名の設定を行うが、出力ファイル名はテキストファイル名（ファイル拡張子が「.txt」）であれば、任意である。
ファイル名の設定が終わると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには多少の時間が掛かる。コンパイル後、MCMCサンプリングが始まり、MCMCサンプリング後、図B2のグラフが表示される。
図B2
パラメータ
の事後分布である。
図B2のWindowを閉じると、図B3のグラフが表示される。
図B3
パラメータ
の事後分布のグラフである。
図B3のWindowを閉じると、図B4のグラフが表示される。
図B4
判断のカテゴリ境界の事後分布である。
図B4のWindowを閉じると、図B５のグラフが表示される。
図B5
図B1の表のデータを比率に直して表示したもの（正方形）と、モデル予測値の事後分布の中央値（小円）のグラフである。正方形（データ）と小円（モデル）の重なりがよいので、モデルはデータによく適合していると言える。
図B５のWindowを閉じると、実行終了である。端末には、以下のメッセージが表示されている。
Results.txt was saved.
(stan) *****/sdt_oddcat$
出力ファイルResults.txtを開くと、以下のような内容である。
Input Data File = Data5Cat.csv
<N, S>
[21, 10, 5, 18, 46]
<S, N>
[66, 10, 7, 12, 5]
Mean
MCSE StdDev ...
ESS_tail ESS_bulk/s R_hat
lp__ -270.337000 0.033033 1.459580 ... 2942.01 25956.6 1.00054
mu_s
0.931202 0.001831 0.118723 ... 3102.69 53068.0 1.00049
tau
0.385409 0.001735 0.115954 ... 3512.09 56082.4 1.00056
theta[1] 0.140015 0.000583 0.037475 ... 2969.47 50438.4 1.00009
theta[2] 0.578691 0.001213 0.073092 ... 2720.57 43876.1 1.00065
C[1]
-0.718706 0.001237 0.079069 ... 2831.07 50589.4 1.00061
C[2]
-0.140015 0.000583 0.037475 ... 2969.47 50438.4 1.00009
C[3]
0.140015 0.000583 0.037475 ... 2969.47 50438.4 1.00009
C[4]
0.718706 0.001237 0.079069 ... 2831.07 50589.4 1.00061
p_NS[1] 0.187613 0.000536 0.034582 ... 3017.60 52503.5 1.00049
p_NS[2] 0.127419 0.000259 0.016823 ... 3312.71 52788.9 1.00086
p_NS[3] 0.072758 0.000302 0.019428 ... 2859.19 50721.2 1.00072
p_NS[4] 0.160512 0.000334 0.020278 ... 2558.40 44966.5 1.00053
p_NS[5] 0.451698 0.000672 0.047737 ... 3372.57 63257.1 1.00093
p_SN[1] 0.077157 0.000330 0.020760 ... 2816.51 49380.9 1.00169
p_SN[2] 0.076247 0.000191 0.012534 ... 3283.69 53404.2 1.00102
p_SN[3] 0.051098 0.000227 0.014653 ... 2952.11 50882.3 1.00078
p_SN[4] 0.132891 0.000309 0.019267 ... 3104.66 47175.8 1.00127
p_SN[5] 0.662607 0.000663 0.044700 ... 3276.61 55393.9 1.00055
[19 rows x 11 columns]
d':
Median = 0.931, 90%CI = [0.739, 1.126]
tau:
Median = 0.388, 90%CI = [0.195, 0.579]
C[1]:
Median = -0.717, 90%CI = [-0.852, -0.592]
C[2]:
Median = -0.136, 90%CI = [-0.205, -0.085]
C[3]:
Median = 0.136, 90%CI = [0.085, 0.205]
C[4]:
Median = 0.717, 90%CI = [0.592, 0.852]
MCMCサンプリングの統計、パラメータの事後分布の中央値などが示されている。
リストB1 奇数個の評定カテゴリ数のStanスクリプト（SDT2AFCOddCat.stan）
//
// Number of rating
categories should be odd
//
//
Yasuharu Okamoto,
2019.10, 2025.10
//
functions {
real half_normal_lpdf(real
y, real sgm){
if
(y > 0.0){
return normal_lpdf(y | 0.0, sgm);
}
else {
return log(0.0);
}
}
}
data {
int K;
array[K] int n_NS;
array[K] int n_SN;
}
parameters {
real mu_s;
real tau;
array[(K - 1) / 2]
real<lower = 0.0> theta;
}
transformed parameters {
array[K - 1] real C;
simplex[K] p_NS;
simplex[K] p_SN;
C[(K - 1) / 2 + 1] =
theta[1];
C[(K - 1) / 2] = -C[(K - 1)
/ 2 + 1];
for (k in 2: ((K - 1) / 2)){
C[(K
- 1) /2 + k] = C[(K - 1) / 2 + k - 1] + theta[k];
C[(K
- 1) / 2 - k + 1] = -C[(K - 1) / 2 + k];
}
p_NS[1] = Phi((C[1] - (mu_s
- tau))/sqrt(2.0));
p_NS[K] = 1 - Phi((C[K - 1]
- (mu_s - tau))/sqrt(2.0));
for (k in 2:(K-1)) {
p_NS[k] = Phi((C[k] - (mu_s - tau))/sqrt(2.0)) -
Phi((C[k - 1] - (mu_s - tau))/sqrt(2.0));
}
p_SN[1] = Phi((C[1] - (mu_s
+ tau))/sqrt(2.0));
p_SN[K] = 1 - Phi((C[K - 1]
- (mu_s + tau))/sqrt(2.0));
for (k in 2:(K-1)) {
p_SN[k] = Phi((C[k] - (mu_s + tau))/sqrt(2.0)) -
Phi((C[k - 1] - (mu_s + tau))/sqrt(2.0));
}
}
model {
for (k in 1:((K - 1) / 2)) {
theta[k] ~ half_normal(1000.0);
}
mu_s ~ normal(0.0, 100.0);
tau ~ normal(0.0, 100.0);
n_NS ~ multinomial(p_NS);
n_SN ~ multinomial(p_SN);
}
リストB2 リストB1のStanスクリプトを用いるPythonスクリプト（SDT2AFCOddCat.py）
#
#
Yasuharu Okamoto,
2019.10, 2025.10
#
from cmdstanpy import CmdStanModel
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import csv
flnm = input('入力データファイル(*.csv) = ')
with open(flnm, 'r') as f:
data = [v for v in
csv.reader(f)]
fout_nm = input('出力テキストファイル名 = ')
fout = open(fout_nm, 'w')
fout.write('\nInput Data File =
{}\n'.format(flnm))
K = len(data[0]) - 1
print('K = ', K)
if (K % 2) == 0:
print('カテゴリ数 K は奇数でなければならない！')
import sys
sys.exit()
Freq_NS = []
for k in range(K):
Freq_NS.append(int(data[1][1
+ k]))
Freq_SN = []
for k in range(K):
Freq_SN.append(int(data[2][1
+ k]))
print(Freq_NS)
print(Freq_SN)
fout.write('\n<N,
S>\n{}\n'.format(Freq_NS))
fout.write('\n<S,
N>\n{}\n'.format(Freq_SN))
TN_NS = sum(Freq_NS)
Rating_NS = Freq_NS
TN_SN = sum(Freq_SN)
Rating_SN = []
for k in range(K):
Rating_SN.append(Freq_SN[K -
1 - k])
Data = {'K': K, 'n_NS': Rating_NS,
'n_SN': Rating_SN}
model =
CmdStanModel(stan_file='SDT2AFCOddCat.stan')
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
fit = fit.draws_pd()
mu = fit['mu_s']
tau = fit['tau']
C = []
for k in range(K-1):
C.append(fit[f'C[{k+1}]'])
C = np.array(C).T
p_NS = []
p_SN = []
for k in range(K):
p_NS.append(fit[f'p_NS[{k+1}]'])
p_SN.append(fit[f'p_SN[{k+1}]'])
p_NS = np.array(p_NS).T
p_SN = np.array(p_SN).T
mu_L05 = np.percentile(mu, 5)
mu_med = np.percentile(mu, 50)
mu_U95 = np.percentile(mu, 95)
fout.write("\nd':\nMedian =
{0:<.3f}, 90%CI =
[{1:<.3f}, {2:<.3f}]\n".
format(mu_med, mu_L05, mu_U95))
tau_L05 = np.percentile(tau, 5)
tau_med = np.percentile(tau, 50)
tau_U95 = np.percentile(tau, 95)
fout.write('\ntau:\nMedian =
{0:<.3f}, 90%CI =
[{1:<.3f}, {2:<.3f}]\n'.
format(tau_med, tau_L05, tau_U95))
C_L05 = np.zeros(K-1)
C_med = np.zeros(K-1)
C_U95 = np.zeros(K-1)
for k in range(K-1):
C_L05[k] =
np.percentile(C.T[k], 5)
C_med[k] =
np.percentile(C.T[k], 50)
C_U95[k] =
np.percentile(C.T[k], 95)
for k in range(K-1):
fout.write(('\nC[{0}]:\n' +
\
'Median = {1:<.3f},
90%CI = [{2:<.3f}, {3:<.3f}]\n').
format(k+1, C_med[k], C_L05[k], C_U95[k]))
sb.kdeplot(mu)
plt.xlabel(r"d'($\mu_S$)",
fontsize = 14)
plt.title(r"Posterior
Distribution of d'($\mu_S$)" + \
'\nMed. = {0:<.3f},
90%CI = [{1:<.3f}, {2:<.3f}]'.
format(mu_med, mu_L05, mu_U95), fontsize = 16)
plt.show()
sb.kdeplot(tau)
plt.xlabel(r'$\tau$', fontsize = 18)
plt.title(r'Posterior Distribution of
$\tau$' + \
'\nMed. = {0:<.3f},
90%CI = [{1:<.3f}, {2:<.3f}]'.
format(tau_med, tau_L05, tau_U95), fontsize = 16)
plt.show()
for k in range(K-1):
sb.kdeplot(C.T[k], label =
'C{}'.format(k+1))
plt.title('Posterior Distributions of
Category Boundaries', fontsize = 18)
plt.legend()
plt.show()
p1 = np.median(p_NS, axis = 0)
p2 = np.median(p_SN, axis = 0)
p2 = np.flip(p2)
xcat = np.arange(1, K+0.1, 1)
xlabels = ['{}'.format(int(v)) for v
in xcat]
y1 = []
y2 = []
for k in range(K):
y1.append(Freq_NS[k]/TN_NS)
y2.append(Freq_SN[k]/TN_SN)
plt.plot(xcat, y1, c='b', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<N, S>')
plt.plot(xcat, p1, c='g', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth=3, label =
'Model/<N, S>')
plt.plot(xcat, y2, c='r', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<S, N>')
plt.plot(xcat, p2, c='m', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth=3, label = 'Model/<S, N>')
plt.xticks(xcat, xlabels, fontsize =
14)
plt.xlabel('Rating Category', fontsize
= 14)
plt.ylabel('Probability/Proportion',
fontsize = 14)
plt.yticks(np.arange(0, 1.01, 0.2))
plt.title('Ratings in Conditions
<N, S> and <S, N>', fontsize = 18)
plt.legend(loc = 'upper center',
fontsize = 18)
plt.show()
fout.close()
print('\n{} was
saved.\n'.format(fout_nm))
判断のカテゴリ数Kが偶数の場合のStanスクリプトをリストC1のように用意した。リストC1のStanスクリプトを使ってベイズ分析を行うPythonスクリプトをリストC2のように用意した。リストC1のファイル（SDT2AFCEvenCat.stan）、リストC2のファイル（SDT2AFCEvenCat.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（ファイル）を移し、CmdStanPyのインストールされている環境において次のコマンドを実行する。なお、CmdStanPy準備の簡単説明を、このウェブサイトで行った。
(stan)
*****/sdt_evencat$ python SDT2AFCEvenCat.py
上のコマンドを実行すると、入力データファイル名と出力データファイル名の設定が求められる。
(stan) *****/sdt_evencat$ python
SDT2AFCEvenCat.py
入力データファイル名(*.csv) = Data4Cat.csv
出力テキストファイル名 = Results.txt
入力データファイルは、CSV形式のファイルとして用意する（図C1）。Excelの場合、ファイル拡張子として「.csv」を選ぶと、CSV形式のファイルとして保存される。
図C1
図C1は、判断カテゴリ数が４の場合の例である。第１行目は変数名（カテゴリ名）の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「多分、シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第４列目に「たぶん、シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数、第５列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。
第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「多分、シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第４列目に「たぶん、シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数、第５列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。
データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
入力データファイル名の設定に続いて、出力ファイル名の設定を行うが、出力ファイル名はテキストファイル名（ファイル拡張子が「.txt」）であれば、任意である。
ファイル名の設定が終わると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには多少の時間が掛かる。コンパイル後、MCMCサンプリングが始まり、MCMCサンプリング後、図C2のグラフが表示される。
図C2
パラメータの事後分布である。
図C2のWindowを閉じると、図C3のグラフが表示される。
図C3
パラメータの事後分布のグラフである。
図C3のWindowを閉じると、図C4のグラフが表示される。
図C4
判断のカテゴリ境界の事後分布である。中央のカテゴリ境界値、図C4の場合はC2、は
に固定されている（式（５））ので、図C4のグラフには表示されていない。
図C4のWindowを閉じると、図C５のグラフが表示される。
図C5
図C1の表のデータを比率に直して表示したもの（正方形）と、モデル予測値の事後分布の中央値（小円）のグラフである。正方形（データ）と小円（モデル）の重なりがよいので、モデルはデータによく適合していると言える。
図C５のWindowを閉じると、実行終了である。端末には、以下のメッセージが表示されている。
Results.txt was saved.
(stan) *****/2AFCfiles/sdt_evencat$
出力ファイルResults.txtを開くと、以下のような内容である。
Input Data File = Data4Cat.csv
<N, S>
[24, 9, 12, 55]
<S, N>
[74, 13, 7, 6]
Mean
MCSE StdDev ...
ESS_tail ESS_bulk/s R_hat
lp__ -210.013000 0.028870 1.235630 ... 2533.94 24819.8 1.00205
mu_s
1.050070 0.002271 0.127871 ... 3037.34 37836.9 1.00107
tau
0.433704 0.002160 0.125256 ... 2805.00 40400.1 1.00231
theta[1] 0.505888 0.001185 0.070206 ... 2548.88 41632.2 1.00079
C[1]
-0.505888 0.001185 0.070206 ... 2548.88 41632.2 1.00079
C[2]
0.000000 NaN 0.000000 ... NaN
NaN
NaN
C[3]
0.505888 0.001185 0.070206 ... 2548.88 41632.2 1.00079
p_NS[1] 0.215756 0.000644 0.038867 ... 3008.57 43619.5 1.00047
p_NS[2] 0.116810 0.000255 0.015828 ... 2472.62 45541.1 1.00071
p_NS[3] 0.136530 0.000325 0.019622 ... 2377.75 42978.1 1.00059
p_NS[4] 0.530904 0.000750 0.048415 ... 2947.36 49534.6 1.00227
p_SN[1] 0.082127 0.000455 0.022824 ... 2545.22 30228.1 1.00084
p_SN[2] 0.067085 0.000190 0.011896 ... 2381.59 45704.8 1.00128
p_SN[3] 0.097352 0.000264 0.017272 ... 2403.68 50376.1 1.00156
p_SN[4] 0.753435 0.000712 0.041785 ... 2453.56 40455.9 1.00079
[15 rows x 11 columns]
d':
Median = 1.051, 90%CI = [0.844, 1.262]
tau:
Median = 0.432, 90%CI = [0.231, 0.643]
C[1]:
Median = -0.504, 90%CI = [-0.626, -0.396]
C[2]:
Median = 0.000, 90%CI = [0.000, 0.000]
C[3]:
Median = 0.504, 90%CI = [0.396, 0.626]
MCMCサンプリングの統計、パラメータの事後分布の中央値などが示されている。
リストC1 偶数数個の評定カテゴリ数のStanスクリプト（SDT2AFCEvenCat.stan）
//
// Number of Rating
categoriees should be even
//
//
Yasuharu Okamoto, 2019.10, 2025.10
//
functions {
real half_normal_lpdf(real
y, real sgm){
if
(y > 0.0){
return
normal_lpdf(y | 0.0, sgm);
}
else {
return log(0.0);
}
}
}
data {
int K;
array[K] int n_NS;
array[K] int n_SN;
}
parameters {
real mu_s;
real tau;
array[K/2 - 1] real<lower
= 0.0> theta;
}
transformed parameters {
array[K-1] real C;
simplex[K] p_NS;
simplex[K] p_SN;
C[K/2] = 0.0;
for (k in 1: (K/2 -1)){
C[K/2 + k] = C[K/2 + k - 1] + theta[k];
C[K/2 - k] = -C[K/2 + k];
}
p_NS[1] = Phi((C[1] - (mu_s
- tau))/sqrt(2.0));
p_NS[K] = 1 - Phi((C[K - 1]
- (mu_s - tau))/sqrt(2.0));
for (k in 2:(K-1)) {
p_NS[k] = Phi((C[k] - (mu_s - tau))/sqrt(2.0)) -
Phi((C[k - 1] - (mu_s - tau))/sqrt(2.0));
}
p_SN[1] = Phi((C[1] - (mu_s
+ tau))/sqrt(2.0));
p_SN[K] = 1 - Phi((C[K - 1]
- (mu_s + tau))/sqrt(2.0));
for (k in 2:(K-1)) {
p_SN[k] = Phi((C[k] - (mu_s + tau))/sqrt(2.0)) -
Phi((C[k - 1] - (mu_s + tau))/sqrt(2.0));
}
}
model {
for (k in 1:(K/2 - 1)) {
theta[k] ~ half_normal(1000.0);
}
mu_s ~ normal(0.0, 100.0);
tau ~ normal(0.0, 100.0);
n_NS ~ multinomial(p_NS);
n_SN ~ multinomial(p_SN);
}
リストC2 リストC1のStanスクリプトを用いるPythonスクリプト（SDT2AFCEvenCat.py）
#
#
Yasuharu Okamoto, 2019.10, 2025.10
#
from cmdstanpy import CmdStanModel
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import csv
import pickle
flnm = input('入力データファイル名(*.csv) = ')
with open(flnm, 'r') as f:
data = [v for v in
csv.reader(f)]
fout_nm = input('出力テキストファイル名 = ')
fout = open(fout_nm, 'w')
fout.write('Input Data File =
{}\n'.format(flnm))
K = len(data[0]) - 1
print('K = ', K)
if (K % 2) == 1:
print('評定カテゴリ数は偶数でなければならない!')
import sys
sys.exit()
Freq_NS = []
for k in range(K):
Freq_NS.append(int(data[1][1
+ k]))
Freq_SN = []
for k in range(K):
Freq_SN.append(int(data[2][1
+ k]))
print(Freq_NS)
print(Freq_SN)
fout.write('\n<N,
S>\n{}\n'.format(Freq_NS))
fout.write('\n<S,
N>\n{}\n'.format(Freq_SN))
TN_NS = sum(Freq_NS)
Rating_NS = Freq_NS
print(Rating_NS)
TN_SN = sum(Freq_SN)
Rating_SN = []
for k in range(K):
Rating_SN.append(Freq_SN[K -
1 - k])
Data = {'K': K, 'n_NS': Rating_NS,
'n_SN': Rating_SN}
model =
CmdStanModel(stan_file='SDT2AFCEvenCat.stan')
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
fit = fit.draws_pd()
mu = fit['mu_s']
tau = fit['tau']
C = []
for k in range(K-1):
C.append(fit[f'C[{k+1}]'])
C = np.array(C).T
p_NS = []
p_SN = []
for k in range(K):
p_NS.append(fit[f'p_NS[{k+1}]'])
p_SN.append(fit[f'p_SN[{k+1}]'])
p_NS = np.array(p_NS).T
p_SN = np.array(p_SN).T
mu_L05 = np.percentile(mu, 5)
mu_med = np.percentile(mu, 50)
mu_U95 = np.percentile(mu, 95)
fout.write("\nd':\nMedian =
{0:<.3f}, 90%CI =
[{1:<.3f}, {2:<.3f}]\n".
format(mu_med, mu_L05, mu_U95))
tau_L05 = np.percentile(tau, 5)
tau_med = np.percentile(tau, 50)
tau_U95 = np.percentile(tau, 95)
fout.write('\ntau:\nMedian =
{0:<.3f}, 90%CI =
[{1:<.3f}, {2:<.3f}]\n'.
format(tau_med, tau_L05, tau_U95))
C_L05 = np.zeros(K-1)
C_med = np.zeros(K-1)
C_U95 = np.zeros(K-1)
for k in range(K-1):
C_L05[k] =
np.percentile(C.T[k], 5)
C_med[k] =
np.percentile(C.T[k], 50)
C_U95[k] =
np.percentile(C.T[k], 95)
for k in range(K-1):
fout.write(('\nC[{0}]:\n' +
\
'Median = {1:<.3f},
90%CI = [{2:<.3f}, {3:<.3f}]\n').
format(k+1, C_med[k], C_L05[k], C_U95[k]))
sb.kdeplot(mu)
plt.xlabel(r"d'($\mu_S$)",
fontsize = 14)
plt.title(r"Posterior
Distribution of d'($\mu_S$)" + \
'\nMed. = {0:<.3f},
90%CI = [{1:<.3f}, {2:<.3f}]'.
format(mu_med, mu_L05, mu_U95), fontsize = 16)
plt.show()
sb.kdeplot(tau)
plt.xlabel(r'$\tau$', fontsize = 18)
plt.title(r'Posterior Distribution of
$\tau$' + \
'\nMed. = {0:<.3f},
90%CI = [{1:<.3f}, {2:<.3f}]'.
format(tau_med, tau_L05, tau_U95), fontsize = 16)
plt.show()
for k in range(K-1):
if k + 1 != K/2:
sb.kdeplot(C.T[k], label = 'C{}'.format(k+1))
plt.title('Posterior Distributions of
Category Boundaries' + \
'\nC{} = 0 fixed'.format(K//2), fontsize = 16)
plt.legend(loc = 'lower center')
plt.show()
xcat = np.arange(1, K+0.1, 1)
xlabels = ['{}'.format(int(v)) for v
in xcat]
y1 = []
y2 = []
for k in range(K):
y1.append(Freq_NS[k]/TN_NS)
y2.append(Freq_SN[k]/TN_SN)
p1 = np.median(p_NS, axis = 0)
p2 = np.median(p_SN, axis = 0)
p2 = np.flip(p2)
plt.plot(xcat, y1, c='b', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<N, S>')
plt.plot(xcat, p1, c='g', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth=3, label = 'Model/<N, S>')
plt.plot(xcat, y2, c='r', marker='s',
markersize=20,
linewidth = 0, label = 'Data/<S, N>')
plt.plot(xcat, p2, c='m', marker='o',
markersize=30, fillstyle='none',
linewidth = 0, markeredgewidth=3, label = 'Model/<S, N>')
plt.xticks(xcat, xlabels, fontsize =
14)
plt.xlabel('Rating Category', fontsize
= 14)
plt.ylabel('Probability/Proportion',
fontsize = 14)
plt.yticks(np.arange(0, 1.01, 0.2))
plt.title('Ratings in Conditions
<N, S> and <S, N>', fontsize = 18)
plt.legend(loc = 'upper center',
fontsize = 18)
plt.show()
fout.close()
print('\n{} was
saved.\n'.format(fout_nm))
Bradyら（2023）は、記憶のパフォーマンスの研究においては２AFCを用いることを推奨しているが、さらに、不等分散モデルがROC曲線によく合うと説明している（p. 442）。2AFCを拡張して不等分散モデルが適用できるように工夫された実験法と分析モデルが提案されている（Okamoto, 2023）。不等分散が予想される場合は、この拡張2肢強制選択評定法を用いることができる。また、不等分散とは等分散を仮定しないということであり、データが等分散のものであっても不等分散モデルを適用することができる。この場合は、シグナル刺激の分散がノイズ刺激と同じ１であると推定される。
と
の値は、２AFC課題を以下のように拡張すると求めることができる。
２AFCでは、提示刺激対は＜ノイズ、シグナル＞と＜シグナル、ノイズ＞の2種類である。これに、＜ノイズ、ノイズ＞と＜シグナル、シグナル＞の２種類を加えて、４種類の刺激対に対して評定法を用いると、シグナル刺激感覚の平均値と分散の推定を行うことができる。
２AFC課題では、各試行において２つの刺激（刺激対）が与えられる。各刺激の提示位置（空間的、時間的）を＜S1, S2＞で表す。第1提示位置がS1、第２提示位置がS2である。＜n, s＞（あるいは、n-s）は、第１位置にノイズ刺激、第２位置にシグナル刺激が提示されたことを表す。
第１提示位置の刺激感覚を
、第２刺激提示位置の刺激感覚を
で表す。
第１刺激位置に対する第２刺激位置のバイアスをbで表す。このバイアスは恒常誤差のように感覚のものと反応のものの両者を含むものである。
このとき、以下のようになる。
のとき、
のとき、
のとき、
のとき、
いま、判断がK個のカテゴリによる評定で回答されるものとし、判断のカテゴリ境界を以下のようにおく。
このとき、
であれば、提示刺激対<
>に対する判断
はが与えられるとする。すなわち、
である。
例えば、判断カテゴリー数が
個のときは、評定判断は以下のようなものが考えられる。
： シグナル刺激は確かに第１位置である。
： シグナル刺激は第１位置である。
： シグナル刺激はたぶん第１刺激である。
： わからない／どちらも同じと思われる。
： シグナル刺激はたぶん第２位置である。
： シグナル刺激は第２位置である。
： シグナル刺激は確かに第２位置である。
平均
、分散
の正規分布の累積分布を
で表し、以下のように便宜的に設定する。
このとき、差
の分布は以下のように与えられる。
刺激対 
が提示されとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
いま、カテゴリ境界は、原点を中心として対称に設定されているとする。
奇数であれば、
上のモデルに基づいてベイズ分析するためのStanスクリプトを、Kが偶数の場合をリストD1に、Kが奇数の場合をリストD2に示す。これらのスクリプトファイルを用いて拡張２AFC評定法によるデータをベイズ分析するスクリプトをリストD3に示す。ファイルは、2AFCfiles.zipにまとめた。
リストD1のファイル（SDT_2AFC_KCat_even.stan）、リストD2のファイル（SDT_2AFC_KCat_odd.stan）、リストD3のファイル（AnalEx2AFCR.py）、および入力データファイルを同じフォルダに置く。このフォルダにカレントディレクトリ（フォルダ）を移して、次のコマンドをCmdStanPyのインストールされた環境で実行する。CmdStanPyの準備の簡単な説明を、このウェブサイトに用意した。
(stan)
*****/sdt_nesgm$ python AnalEx2AFCR.py
上のコマンドを実行すると、出力ファイル名と入力データファイル名の設定が求められる。
(stan) *****/sdt_nesgm$ python
AnalEx2AFCR.py
Output file (*.txt) = Results.txt
Input data file (*.xlsx) =
Data6Cat.xlsx
出力ファイル名は、任意のテキストファイル名（ファイル拡張子が「.txt」）を設定すればよい。
入力データファイルは、図D1の形式のExcelファイルとして用意する。
図D1
図D1のデータは、判断カテゴリ数Kが６個の場合である。第１行目は、実験条件を表している。第２列目が刺激対 、第３列目が刺激対 、第４列目が刺激対
、第５列目が刺激対
である。実験条件を表す刺激対の順序は、この順番で入力する。
第１列目は、判断カテゴリRを表す。カテゴリ１からカテゴリK（図D1の場合は、K=7）までの数値で順番に表す。判断カテゴリの数値Rと判断（反応）の関係は、この箇所の説明を参照されたい。例えば、カテゴリ１「シグナル刺激は確かに第１位置」、カテゴリ２「シグナル刺激は第１位置」、カテゴリ３「シグナル刺激はたぶん第１位置」、カテゴリ４「シグナル刺激はたぶん第２位置」、カテゴリ５「シグナル刺激は第２位置」、カテゴリ６「シグナル刺激は確かに第２位置」となる。
刺激提示条件と判断カテゴリの度数データを、該当するセルに書き込んでいく。
ファイル名を設定すると、ファイルが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルに多少の時間が掛かるが、コンパイル後、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、図D2のグラフが表示される。
図D2
パラメータすなわちの事後分布のグラフである。
図D2のWindowを閉じると、図D3のグラフが表示される。
図D3
パラメータ
の事後分布のグラフである。
図D3のWindowを閉じると、図D4のグラフが表示される。
図D4
判断カテゴリの境界値の事後分布のグラフである。境界値C3は、カテゴリ数が偶数の場合は、仮定（６）により中央のカテゴリ境界値は０に固定されていることにより、
で定数であるので、原点上に赤の小円盤で表示されている。
図D4のWindowを閉じると、図D5のグラフが表示される。
図D5
パラメータbの事後分布のグラフである。
図D5のWindowを閉じると、図D6のグラフが表示される。
図D6
各刺激提示条件において、カテゴリ反応がｋまでのデータの累積比率Obs(k)と、モデルの予測累積比率（事後分布の中央値）Est(k)の点
を線分で結んだものである。モデルの予測累積比率とデータの累積比率が一致すれば、この線分で結んだ折れ線は、原点と点（1，1）を通る線分（黒破線）に重なる。図Ｄ６の折れ線は、この対角線の近くにあると言える。
図D6のWindowを閉じると、実行終了である。
端末には、次のようなメッセージが表示されている。
Results.txt was saved.
(stan) *****/sdt_nesgm$
出力ファイルResults.txtを開くと、以下のような内容である。
Mean
MCSE
StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -299.271000 0.036671 1.577230e+00 ... 2502.58 8454.32 1.000630
mu_s
0.951517 0.003732 1.945770e-01 ... 2505.81 12214.60 0.999686
sgm_s 1.278720 0.003718 1.852620e-01 ... 2526.00 11063.70 1.000690
b
-0.027728 0.001971 1.157970e-01 ... 2459.44 15541.80 1.000130
preC[1] 0.626552 0.000468 2.359620e-02 ... 2664.96 11415.00 1.000610
...
...
...
... ... ...
...
...
C[5]
2.679530 0.005267 2.676970e-01 ... 2674.99 11517.60 1.000450
vsum
0.726786 0.000394 1.989950e-02 ... 2674.99 11517.60 1.000450
sgm_nn 1.414210 NaN 2.242930e-14 ... NaN
NaN
NaN
sgm_ns 1.627280 0.002934 1.462370e-01 ... 2526.00 11063.70 1.000690
sgm_ss 1.808380 0.005259 2.620000e-01 ... 2526.00 11063.70 1.000700
[68 rows x 11 columns]
リストD1 評定カテゴリ数が偶数のとき（SDT_2AFC_KCat_even.stan）
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_cond;
array[K] int ss_cond;
}
transformed data {
int hK;
vector[K/2] alpha;
hK = K / 2;
for (i in 1:hK) {
alpha[i] = 1.0;
}
}
parameters {
real mu_s;
real<lower = 0.0>
sgm_s;
real b;
simplex[hK] preC;
}
transformed parameters {
vector[K] nn_theta;
vector[K] ns_theta;
vector[K] sn_theta;
vector[K] ss_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K+1] ss_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
real sgm_ss;
C[hK] = 0.0;
vsum = 0.0;
for (i in 1:hK-1) {
vsum
+= preC[i];
C[hK+i] = vsum / (1.0 - vsum);
C[hK-i] = -C[hK+i];
}
sgm_nn = sqrt(1.0 + 1.0);
sgm_ns = sqrt(1.0 +
square(sgm_s));
sgm_ss =
sqrt(square(sgm_s)*2);
nn_cum_p[1] = 0.0;
ns_cum_p[1] = 0.0;
sn_cum_p[1] = 0.0;
ss_cum_p[1] = 0.0;
nn_cum_p[K+1] = 1.0;
ns_cum_p[K+1] = 1.0;
sn_cum_p[K+1] = 1.0;
ss_cum_p[K+1] = 1.0;
for (i in 2:K) {
nn_cum_p[i] = normal_cdf(C[i-1] | b, sgm_nn);
ns_cum_p[i] = normal_cdf(C[i-1] | mu_s + b, sgm_ns);
sn_cum_p[i] = normal_cdf(C[i-1] | -mu_s + b, sgm_ns);
ss_cum_p[i] = normal_cdf(C[i-1] | b, sgm_ss);
}
for (i in 1:K) {
nn_theta[i] = nn_cum_p[i+1] - nn_cum_p[i];
ns_theta[i] = ns_cum_p[i+1] - ns_cum_p[i];
sn_theta[i] = sn_cum_p[i+1] - sn_cum_p[i];
ss_theta[i] = ss_cum_p[i+1] - ss_cum_p[i];
}
}
model {
mu_s ~ normal(0.0, 1000.0);
sgm_s ~ exponential(0.001);
preC ~ dirichlet(alpha);
b ~ normal(0.0, 1000.0);
nn_cond ~
multinomial(nn_theta);
ns_cond ~
multinomial(ns_theta);
sn_cond ~
multinomial(sn_theta);
ss_cond ~
multinomial(ss_theta);
}
リストD2 評定カテゴリ数が奇数のとき（SDT_2AFC_KCat_odd.stan）
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_cond;
array[K] int ss_cond;
}
transformed data {
int hK;
vector[(K/2)+1] alpha;
hK = K / 2;
for (i in 1:hK+1) {
alpha[i] = 1.0;
}
}
parameters {
real mu_s;
real<lower = 0.0>
sgm_s;
real b;
simplex[hK+1] preC;
}
transformed parameters {
vector[K] nn_theta;
vector[K] ns_theta;
vector[K] sn_theta;
vector[K] ss_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K+1] ss_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
real sgm_ss;
vsum = 0.0;
for (i in 1:hK) {
vsum
+= preC[i];
C[hK+i] = vsum / (1.0 - vsum);
C[hK+1-i] =
-C[hK+i];
}
sgm_nn = sqrt(1.0 + 1.0);
sgm_ns = sqrt(1.0 +
square(sgm_s));
sgm_ss =
sqrt(square(sgm_s)*2);
nn_cum_p[1] = 0.0;
ns_cum_p[1] = 0.0;
sn_cum_p[1] = 0.0;
ss_cum_p[1] = 0.0;
nn_cum_p[K+1] = 1.0;
ns_cum_p[K+1] = 1.0;
sn_cum_p[K+1] = 1.0;
ss_cum_p[K+1] = 1.0;
for (i in 2:K) {
nn_cum_p[i] = normal_cdf(C[i-1] | b, sgm_nn);
ns_cum_p[i] = normal_cdf(C[i-1] | mu_s + b, sgm_ns);
sn_cum_p[i] = normal_cdf(C[i-1] | -mu_s + b, sgm_ns);
ss_cum_p[i] = normal_cdf(C[i-1] | b, sgm_ss);
}
for (i in 1:K) {
nn_theta[i] = nn_cum_p[i+1] - nn_cum_p[i];
ns_theta[i] = ns_cum_p[i+1] - ns_cum_p[i];
sn_theta[i] = sn_cum_p[i+1] - sn_cum_p[i];
ss_theta[i] = ss_cum_p[i+1] - ss_cum_p[i];
}
}
model {
mu_s ~ normal(0.0, 1000.0);
sgm_s ~ exponential(0.001);
preC ~ dirichlet(alpha);
b ~ normal(0.0, 1000.0);
nn_cond ~
multinomial(nn_theta);
ns_cond ~
multinomial(ns_theta);
sn_cond ~
multinomial(sn_theta);
ss_cond ~
multinomial(ss_theta);
}
リストD3 リストD1のスクリプトとリストD2のスクリプトを用いるPythonスクリプト（AnalEx2AFCR.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
outflnm = input('Output file (*.txt) =
')
fout = open(outflnm, 'w')
inflnm = input('Input data file
(*.xlsx) = ')
df_raw = pd.read_excel(inflnm)
print(df_raw)
data = np.array(df_raw.values)
print(data)
K = len(data)
print('K =', K)
n_n = data.T[1]
n_s = data.T[2]
s_n = data.T[3]
s_s = data.T[4]
print('n_n =', n_n)
print('n_s =', n_s)
print('s_n =', s_n)
print('s_s =', s_s)
model = CmdStanModel(stan_file =
'SDT_2AFC_KCat_odd.stan' if K % 2 == 1 else
'SDT_2AFC_KCat_even.stan')
Data = {'K':K, 'nn_cond':n_n,
'ns_cond':n_s, 'sn_cond':s_n, 'ss_cond':s_s}
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
fit = fit.draws_pd()
mu_s_med = np.median(fit['mu_s'])
sb.kdeplot(fit['mu_s'])
plt.xlabel(r'$\mu_s$')
plt.title(r'$\mu_s$(Med) =
{0:.3f}'.format(mu_s_med))
plt.show()
sgm_s_med = np.median(fit['sgm_s'])
sb.kdeplot(fit['sgm_s'])
plt.title(r'$\sigma_s$(Med) =
{0:.3f}'.format(sgm_s_med))
plt.xlabel(r'$\sigma_s$')
plt.show()
C = []
for k in range(K-1):
C.append(fit[f'C[{k+1}]'])
C = np.array(C).T
s_title = ''
if K % 2 == 1:
for k in range(K-1):
sb.kdeplot(C.T[k])
c_med = np.median(C.T[k])
s_title += f'C{k+1}={c_med:.2f}'
if k
< K-2:
s_title += ', '
else:
for k in range(K-1):
if k
!= (K//2) - 1:
sb.kdeplot(C.T[k])
#fit['C'].T[k])
c_med = np.median(C.T[k])
s_title +=
f'C{k+1}={c_med:.2f}'
if k < K-2:
s_title += ', '
else:
plt.plot([0], [0], marker = 'o', markersize = 15, c = 'r')
s_title += f'C{k+1}=0, '
plt.title('Med. estimates\n' +
s_title)
plt.show()
b_med = np.median(fit['b'])
sb.kdeplot(fit['b'])
plt.title('b(Med) =
{0:.3f}'.format(b_med))
plt.show()
cum_n_n = np.cumsum(n_n)
pcum_n_n = cum_n_n/cum_n_n[-1]
cum_n_s = np.cumsum(n_s)
pcum_n_s = cum_n_s / cum_n_s[-1]
cum_s_n = np.cumsum(s_n)
pcum_s_n = cum_s_n / cum_s_n[-1]
cum_s_s = np.cumsum(s_s)
pcum_s_s = cum_s_s / cum_s_s[-1]
fit_nn_cum_p = []
fit_ns_cum_p = []
fit_sn_cum_p = []
fit_ss_cum_p = []
for k in range(K+1):
fit_nn_cum_p.append(fit[f'nn_cum_p[{k+1}]'])
fit_ns_cum_p.append(fit[f'ns_cum_p[{k+1}]'])
fit_sn_cum_p.append(fit[f'sn_cum_p[{k+1}]'])
fit_ss_cum_p.append(fit[f'ss_cum_p[{k+1}]'])
fit_nn_cum_p =
np.array(fit_nn_cum_p).T
fit_ns_cum_p =
np.array(fit_ns_cum_p).T
fit_sn_cum_p =
np.array(fit_sn_cum_p).T
fit_ss_cum_p =
np.array(fit_ss_cum_p).T
est_pcum_n_n = np.median(fit_nn_cum_p,
axis = 0)[1:-1]
est_pcum_n_s = np.median(fit_ns_cum_p,
axis = 0)[1:-1]
est_pcum_s_n = np.median(fit_sn_cum_p,
axis = 0)[1:-1]
est_pcum_s_s = np.median(fit_ss_cum_p,
axis = 0)[1:-1]
plt.plot(est_pcum_n_n, pcum_n_n[:-1],
label = 'n_n')
plt.plot(est_pcum_n_s, pcum_n_s[:-1],
label = 'n_s')
plt.plot(est_pcum_s_n, pcum_s_n[:-1],
label = 's_n')
plt.plot(est_pcum_s_s, pcum_s_s[:-1],
label = 's_s')
plt.plot([0,1], [0,1], c = 'k', ls =
'--', label = 'Obs.=Est.')
plt.xlabel('Est.Cum.P')
plt.ylabel('Obs.Cum.P')
plt.legend()
plt.title('Cumulative Proportions')
plt.show()
fout.close()
print('\n', outflnm, 'was saved.\n')
Brady,T.F, Robinson,M.M., Williams,J.R., & Wixted,J.T. (2023). easuring memory is harder than you think: How to avoid problematic measurement practices in memory research. Psychonomic Bulletin & Review, 30, 421–449.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
Okamoto, Y. (2023) Extended 2AFC Rating Task. OSF. from https://osf.io/ts2eq/
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房