信号検出理論
２肢強制選択課題とその拡張版
評定法データのベイズ分析
岡本安晴、2019.10、2023.05
２肢強制選択（2 alternative-forced-choice: 3AFC）課題の場合、ノイズ刺激
とシグナル刺激
が各試行において提示されるが、空間的位置を変えて（右と左、あるいは上と下など）同時に提示される場合、あるいは提示の時間順序を変えて（第１インタバルと第２インタバル）提示される場合がある。統計モデルを考えるときは、空間順序あるいは時間順序における２つの刺激提示法をともに＜、＞あるいは＜、＞で表してよい。＜、＞は、刺激が空間位置で提示されているときは、例えば、左側にノイズ刺激、右側にシグナル刺激が提示されることを表し、時間順序のときは第１インタバルでノイズ刺激が提示され、第２インタバルでシグナル刺激が提示されることを表す。＜、＞のときは、順序が逆になる。以後の説明では、例えば＜、＞は、は第１位置、は第２位置に提示される条件を表すということにする。刺激提示条件別に判断（反応）を分類すると表1のようになる。
表1 2AFC課題における刺激提示条件と反応（判断）。
|
|
シグナル刺激の提示位置として判断された位置 |
|
|
|
第１位置 |
第２位置 |
|
＜ |
誤反応 |
正反応 |
|
＜ |
正反応 |
誤反応 |
いま、1AFC課題のモデルと同じく、ノイズ刺激の感覚
が平均０、分散１の正規分布に従い、シグナル刺激の感覚
が平均
、分散
の正規分布に従うものとする（等分散ガウス型モデル）。すなわち、
、 
である。
提示条件＜、＞および＜、＞において恒常誤差
および位置選択バイアス
を考える。恒常誤差は、時間誤差あるいは空間誤差を表す。位置選択バイアスは、感覚の違いの弁別とは別に第１位置あるいは第２位置の選択を好む傾向を表す。これらのバイアスをまとめて
とおき、第２位置を基準にした第１位置へのバイアス効果として表すと、感覚およびの分布は以下のようになる。
刺激提示が＜、＞のとき、
、 
で与えられ、
となる。
刺激提示が＜、＞のとき、
、 
で与えられ、
となる。
評定法における判断カテゴリの境界値（基準値）を
とおく。すなわち、シグナル刺激の感覚とノイズ刺激の感覚の差
が
との間にあるときカテゴリ
の判断になると考える。
のとき
観察者の判断
カテゴリ
である。
判断カテゴリの総数が４カテゴリのとき、例えば判断カテゴリは次のようになる。ここで、の値は、シグナル刺激の感覚が強いほどカテゴリの値は大きい値になるが、＜、＞と＜、＞ではシグナル刺激の提示位置が異なることに注意。
＜、＞のとき、
カテゴリ１： シグナルは第１位置である。
カテゴリ２： シグナルは多分第１位置である。
カテゴリ３： シグナルは多分第２位置である。
カテゴリ４： シグナルは第２位置である。
＜、＞のとき、
カテゴリ１： シグナルは第２位置である。
カテゴリ２： シグナルは多分第２位置である。
カテゴリ３： シグナルは多分第１位置である。
カテゴリ４： シグナルは第１位置である。
いま、判断カテゴリの総数を
とおいたとき、カテゴリ境界を
とおく。累積標準正規分布関数
に対して
、 
とおいたとき、判断カテゴリがカテゴリである確率は次のように与えられる。
＜、＞のとき、カテゴリである確率
は、式（１）より
となる。
＜、＞のとき、カテゴリである確率
は、式（２）より
となる。
表2 刺激提示条件に対するカテゴリ評定判断の頻度。
|
|
カテゴリ１ |
・・・ |
カテゴリk |
・・・ |
カテゴリK |
|
＜ |
|
・・・ |
|
・・・ |
|
|
＜ |
|
・・・ |
|
・・・ |
|
各刺激提示条件に対するカテゴリ判断の頻度が表2のようであるとき、尤度関数は次式で与えられる。
モデルのパラメータの数は
であるが、式（３）および式（４）の関数の引数の分数形から以下のようなパラメータ値の不定性が認められる。
分数の分母子に同じ数を掛けても値は変わらないので、分母の値を
とおいても一般性は失われない。これは、等分散モデル（
）あるいはより一般的には感覚の単位を
であるように設定するということである。
しかし、Brady et al. (2023; https://link.springer.com/article/10.3758/s13423-022-02179-w)は、記憶のパフォーマンスの研究においては２AFCを用いることを推奨しているが、さらに、不等分散モデルがROC曲線によく合うと説明している（p. 442）。2AFCを拡張して不等分散モデルが適用できるように工夫された実験法と分析モデルが提案されている（https://osf.io/ts2eq/）。不等分散が予想される場合は、この拡張2肢強制選択評定法を用いるべきであろう。拡張版については、本ウェブサイトの後半で説明した。
標準の２AFC課題について議論を続ける。
この場合は、等分散モデルで考えるが、感覚の弁別力は
と表記される。
また、分子におけるとに同じ定数を加えても
であるので、パラメータ値を決めるためには制約条件が必要である。感覚の差の判断が原点に関して対称であるとすると、これはカテゴリ境界に次の仮定をおくことである。
カテゴリ数が偶数
（
自然数）のとき
このとき、モデルのパラメータの自由度は
となる。
カテゴリ数が奇数
（ 自然数）のとき
このとき、モデルのパラメータの自由度は
となる。
表2のデータの自由度
は、表2のデータを刺激提示条件における頻度の比率データと考えれば
となる。パラメータの推定が不定とならないためには、パラメータの自由度
はデータの自由度以下でなければならない。すなわち、
であるが、この条件は満たされている。
上記のモデルに基づくベイズ分析のためのStanスクリプトとそれを用いるPythonスクリプトを、K＝２（標準的な２AFC課題）、Kが奇数、Kが４以上の偶数の３つの場合に分けて作成してみた。なお、Python用のStan（PyStan）については、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。信号検出理論の基礎的な説明を拙著＜岡本安晴「心理学データ分析と測定」勁草書房＞あるいは＜岡本安晴「計量心理学」培風館＞で行っている。
K=2の場合（標準的な２AFC課題）
StanスクリプトリストA.1のように用意した。
リストA.1 K=2の場合のStanスクリプト（ファイル名SDT2AFC2Cat.stan）
data {
int n_NS[2];
int n_SN[2];
}
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);
}
リストA.1のStanスクリプトファイルSDT2AFC2Cat.stanを利用するPythonスクリプトを、リストA.2のように用意した。
リストA.2 StanスクリプトファイルSDT2AFC2Cat.stanを利用するPythonスクリプト（ファイル名SDT2AFC2Cat.py）
import numpy as np
import pystan
import matplotlib.pyplot as plt
import csv
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}
sm = pystan.StanModel(file =
'SDT2AFC2Cat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fout.write('\n{}\n'.format(fit))
mu = fit['mu_s']
tau = fit['tau']
p_NS = fit['p_NS']
p_SN = fit['p_SN']
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)
p2 = np.median(p_SN, axis = 0)
p2 = [p2[1], p2[0]]
plt.hist(mu)
plt.xlabel("d'($\mu_S$)",
fontsize = 14)
plt.title("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()
plt.hist(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, 'b:', linewidth = 10,
label = 'Data/<N, S>')
plt.plot(xcat, p1, 'g-', linewidth =
3, label = 'Model/<N, S>')
plt.plot(xcat, y2, 'r:', linewidth =
10, label = 'Data/<S, N>')
plt.plot(xcat, p2, 'm-', linewidth =
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.title('Rating in Conditions <N,
S> and <S, N>', fontsize = 18)
plt.legend(loc = 'upper center',
fontsize = 10)
plt.show()
fout.close()
print('\n' + fout_nm + ' was
saved.\n')
リストA.1とリストA.2のファイルおよびサンプルデータファイルは、まとめて圧縮ファイルRating2Cat.zipとした。ダウンロード解凍すれば、自由に利用できる。
リストA.2のファイルSDT2AFC2Cat.pyを実行すると、入力データファイル名と出力データファイル名の設定が求められる（図A.1）。
図A.1
出力データファイルは任意のテキストファイル名（ファイル拡張子.txt）であるが、入力データファイルはCSV形式のファイルとして図A.2の様式で用意する。
図A.2 データファイルData2Cat.csv
第１行目は変数名の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第２行目、第３行目とも、第２列目は「シグナル刺激は第１位置である（例えば、シグナル刺激は左である）」と判断された度数、第３列目は「シグナル刺激は第２位置である（例えば、シグナル刺激は右である）」と判断された度数である。データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
図A.1のように入力データファイル名と出力テキストファイル名を設定すると、計算が始まる。Stanによるサンプリングが終了すると、まず
（）の事後分布のヒストグラムが表示される（図A.3）。
図A.3
事後分布の中央値がMed.=1.075、90％CIが[0.850, 1.316]であることが示されている。
図A.3のふぉむを閉じると、τの事後分布のヒストグラムが表示される（図A.4）。
図A.4
事後分布の中央値がMed.=0.335、90％CIが[0.106, 0.579]であることが示されている。
図A.4のフォームを閉じると、各判断カテゴリの確率／比率を示すグラフが折れ線グラフで表示される（図A.5）。
図A.5
データが破線で比率を表し、モデルが実線で各判断カテゴリの確率の中央値である。「シグナル刺激は第１位置である」判断が「Cat-1」で、「シグナル刺激は第２位置である」判断が「Cat-2」で表されている。
図A.5のフォームを閉じると、プログラムの実行終了となる。実行終了後、出力ファイルを適当なエディタで開くと、以下のようになっている。
Input Data File = Data2Cat.csv
<N, S>: [30, 70]
<S, N>: [84, 16]
Inference for Stan model:
anon_model_76fee9688e040eaa9d47a70bb9cf69e2.
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
mu_s 1.08 2.6e-3 0.14 0.8 0.99 1.08 1.18 1.37 3058 1.0
tau 0.34 2.3e-3 0.14 0.06 0.24 0.34 0.43 0.62 3648 1.0
p_NS[0] 0.3 7.2e-4 0.05 0.22 0.27 0.3 0.33 0.39 4000 1.0
p_NS[1] 0.7 7.2e-4 0.05 0.61 0.67 0.7 0.73 0.78 4000 1.0
p_SN[0] 0.16 6.8e-4 0.04 0.1 0.13 0.16 0.18 0.24 2877 1.0
p_SN[1] 0.84 6.8e-4 0.04 0.76 0.82 0.84 0.87 0.9 2877 1.0
lp__ -106.0 0.02 1.04 -108.7 -106.4 -105.7 -105.3
-105.0 1800 1.0
Samples were drawn using NUTS at Tue
Oct 22 10:30:19 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).
d':
Med. = 1.082, 90%CI = [0.847, 1.320
tau:
Med. = 0.339, 90%CI ~ [0.105, 0.577]
Kが奇数の場合
判断のカテゴリ数Kが奇数の場合のStanスクリプトをリストB.1のように用意した。
リストB.1 Kが奇数の場合のStanスクリプト（ファイル名SDT2AFCOddCat.stan）
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;
int n_NS[K];
int n_SN[K];
}
parameters {
real mu_s;
real tau;
real<lower = 0.0>
theta[(K - 1) / 2];
}
transformed parameters {
real C[K - 1];
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);
}
リストB.1のStanスクリプトSDT2AFCOddCat.stanを利用するPythonスクリプトを、リストB.2のように用意した。
リストB.2 StanスクリプトSDT2AFCOddCat.stanを利用するPythonスクリプト（ファイル名SDT2AFCOddCat.py）
import numpy as np
import pystan
import matplotlib.pyplot as plt
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}
sm = pystan.StanModel(file =
'SDT2AFCOddCat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fout.write('\n{}\n'.format(fit))
mu = fit['mu_s']
tau = fit['tau']
C = fit['C']
p_NS = fit['p_NS']
p_SN = fit['p_SN']
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]))
plt.hist(mu)
plt.xlabel("d'($\mu_S$)",
fontsize = 14)
plt.title("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()
plt.hist(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):
plt.hist(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, 'b--', linewidth =
3, alpha = 0.7, label = 'Data/<N, S>')
plt.plot(xcat, p1, 'g-', linewidth =
3, alpha = 0.7, label = 'Model/<N, S>')
plt.plot(xcat, y2, 'r--', linewidth =
3, alpha = 0.7, label = 'Data/<S, N>')
plt.plot(xcat, p2, 'm-', linewidth =
3, alpha = 0.7, label = 'Model/<S, N>')
plt.xticks(xcat, xlabels, fontsize =
14)
plt.xlabel('Rating Category', fontsize
= 14)
plt.ylabel('Probability/Proportion',
fontsize = 14)
plt.title('Ratings in Conditions
<N, S> and <S, N>', fontsize = 18)
plt.legend(loc = 'upper center',
fontsize = 10)
plt.show()
fout.close()
print('\n{} was
saved.\n'.format(fout_nm))
リストB.1とリストB.2のファイルおよびサンプルデータファイルは、まとめて圧縮ファイルRatingOddCat.zipとした。ダウンロード解凍すれば、自由に利用できる。
リストB.2のファイルSDT2AFCOddCat.pyを実行すると、入力データファイル名と出力データファイル名の設定が求められる（図B.1）。
図B.1
出力データファイルは任意のテキストファイル名（ファイル拡張子.txt）であるが、入力データファイルはCSV形式のファイルとして図B.2の様式で用意する。
図B.2 データファイルData3Cat.csv
第１行目は変数名の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「わからない」と判断された度数、第４列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第3列目に「わからない」、第４列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第２行目、第３行目とも、第２列目は「シグナル刺激は第１位置である（例えば、シグナル刺激は左である）」と判断された度数、第３列目は「わからない」、第４列目は「シグナル刺激は第２位置である（例えば、シグナル刺激は右である）」と判断された度数である。データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
図B.1のように入力データファイル名と出力テキストファイル名を設定すると、計算が始まる。Stanによるサンプリングが終了すると、まず（）の事後分布のヒストグラムが表示される（図B.3）。
図B.3
事後分布の中央値がMed.=1.055、90％CIが[0.834, 1.262]であることが示されている。
図B.3のふぉむを閉じると、τの事後分布のヒストグラムが表示される（図B.4）。
図B.4
事後分布の中央値がMed.=0.403、90％CIが[0.194, 0.614]であることが示されている。
図B.4のフォームを閉じると、カテゴリ境界の事後分布が示される（図B.5）。
図B.5
図B.5のフォームを閉じると、各判断カテゴリの確率／比率を示すグラフが折れ線グラフで表示される（図B.6）。
図B.6
データが破線で比率を表し、モデルが実線で各判断カテゴリの確率の中央値である。「シグナルは第１位置である」判断が「カテゴリ１」で、「わからない」判断が「カテゴリ２」、「シグナルは第２位置である」判断が「カテゴリ３」で表されている。
図B.6のフォームを閉じると、プログラムの実行終了となる。実行終了後、出力ファイルを適当なエディタで開くと、以下のようになっている。
Input Data File = Data3Cat.csv
<N, S>
[20, 26, 54]
<S, N>
[77, 12, 11]
Inference for Stan model:
anon_model_4f4c02dce251e874ae19399d5b4bfa08.
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
mu_s 1.05 2.0e-3 0.13 0.8 0.97 1.05 1.14 1.3 4000 1.0
tau
0.4 2.1e-3 0.13 0.15 0.32 0.4 0.49 0.65 3685 1.0
theta[0] 0.46 1.2e-3 0.07 0.33 0.41 0.46 0.5 0.6 3380 1.0
C[0] -0.46 1.2e-3 0.07 -0.6 -0.5 -0.46 -0.41 -0.33 3380 1.0
C[1] 0.46 1.2e-3 0.07 0.33 0.41 0.46 0.5 0.6 3380 1.0
p_NS[0] 0.22 6.2e-4 0.04 0.15 0.19 0.22 0.24 0.3 4000 1.0
p_NS[1] 0.23 5.3e-4 0.03 0.17 0.2 0.23 0.25 0.3 4000 1.0
p_NS[2] 0.55 7.7e-4 0.05 0.46 0.52 0.55 0.59 0.65 4000 1.0
p_SN[0] 0.09 4.9e-4 0.02 0.05 0.07 0.09 0.11 0.15 2583 1.0
p_SN[1] 0.15 4.4e-4 0.03 0.1 0.13 0.15 0.17 0.21 4000 1.0
p_SN[2] 0.76 6.9e-4 0.04 0.67 0.73 0.76 0.79 0.83 3707 1.0
lp__ -181.4 0.03 1.25 -184.6 -181.9 -181.0 -180.4
-179.9 1942 1.0
Samples were drawn using NUTS at Wed
Oct 23 11:33:50 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).
d':
Median = 1.055, 90%CI = [0.834, 1.262]
tau:
Median = 0.403, 90%CI = [0.194, 0.614]
C[1]:
Median = -0.455, 90%CI = [-0.575, -0.353]
C[2]:
Median = 0.455, 90%CI = [0.353, 0.575]
次に、判断カテゴリ数が５個の場合（K=5）の例を示す。
データ例を図B.7に示す。
図B.7 データファイル名Data5Cat.csv
プログラムを起動して、図B.8のように入力データファイル名と出力データファイル名を設定する。
図B.8
入力データファイル名と出力データファイル名を設定すると計算が始まり、Stanによるサンプリングが終了すると、まずd’の事後分布が表示される（図B.9）。
図B.9
図B.9のフォームを閉じると、τの事後分布が表示される（図B.10）。
図B.10
図B.10のフォームを閉じると、カテゴリ境界の事後分布が表示される（図B.11）。
図B.11
図B.11のフォームを閉じると、判断カテゴリの比率と確率が表示される（図B.12）。
図B.12
データが破線で比率を表し、モデルが実線で各判断カテゴリの確率の中央値である。例えば、「シグナルは第１位置である」判断が「カテゴリ１」、「たぶんシグナルは第１位置である」判断が「カテゴリ２」、「わからない」判断が「カテゴリ３」、「たぶんシグナルは第２位置である」判断が「カテゴリ４」、「シグナルは第２位置である」判断が「カテゴリ５」で表されている。
図B.12のフォームを閉じると、プログラムの実行終了となる。実行終了後、出力ファイルをテキストエディタで開くと以下のようになっている。
Input Data File = Data5Cat.csv
<N, S>
[21, 10, 5, 18, 46]
<S, N>
[66, 10, 7, 12, 5]
Inference for Stan model:
anon_model_4f4c02dce251e874ae19399d5b4bfa08.
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
mu_s 0.93 1.9e-3 0.12 0.68 0.85 0.93 1.01 1.17 4000 1.0
tau
0.38 1.9e-3 0.12 0.16 0.3 0.38 0.46 0.62 4000 1.0
theta[0] 0.14 5.9e-4 0.04 0.08 0.11 0.14 0.16 0.22 3924 1.0
theta[1] 0.58 1.1e-3 0.07 0.44 0.53 0.58 0.62 0.72 4000 1.0
C[0] -0.72 1.2e-3 0.08 -0.87 -0.77 -0.72 -0.66 -0.57 4000 1.0
C[1] -0.14 5.9e-4 0.04 -0.22 -0.16 -0.14 -0.11 -0.08 3924 1.0
C[2] 0.14 5.9e-4 0.04 0.08 0.11 0.14 0.16 0.22 3924 1.0
C[3] 0.72 1.2e-3 0.08 0.57 0.66 0.72 0.77 0.87 4000 1.0
p_NS[0] 0.19 5.6e-4 0.04 0.12 0.16 0.19 0.21 0.26 4000 1.0
p_NS[1] 0.13 2.6e-4 0.02 0.1 0.12 0.13 0.14 0.16 4000 1.0
p_NS[2] 0.07 3.1e-4 0.02 0.04 0.06 0.07 0.09 0.11 3901 1.0
p_NS[3] 0.16 3.2e-4 0.02 0.12 0.15 0.16 0.17 0.2 4000 1.0
p_NS[4] 0.45 7.6e-4 0.05 0.36 0.42 0.45 0.48 0.55 4000 1.0
p_SN[0] 0.08 3.4e-4 0.02 0.04 0.06 0.08 0.09 0.12 3837 1.0
p_SN[1] 0.08 2.0e-4 0.01 0.05 0.07 0.08 0.08 0.1 4000 1.0
p_SN[2] 0.05 2.4e-4 0.01 0.03 0.04 0.05 0.06 0.08 3703 1.0
p_SN[3] 0.13 3.0e-4 0.02 0.1 0.12 0.13 0.14 0.17 4000 1.0
p_SN[4] 0.66 7.2e-4 0.05 0.57 0.63 0.66 0.69 0.75 4000 1.0
lp__ -270.3 0.03 1.41 -273.8 -271.0 -270.0 -269.3
-268.5 2250 1.0
Samples were drawn using NUTS at Wed
Oct 23 12:01:07 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).
d':
Median = 0.927, 90%CI = [0.727, 1.131]
tau:
Median = 0.383, 90%CI = [0.190, 0.579]
C[1]:
Median = -0.715, 90%CI = [-0.845, -0.592]
C[2]:
Median = -0.136, 90%CI = [-0.206, -0.086]
C[3]:
Median = 0.136, 90%CI = [0.086, 0.206]
C[4]:
Median = 0.715, 90%CI = [0.592, 0.845]
Kが４以上の偶数の場合
判断のカテゴリ数Kが４以上の偶数である場合のStanスクリプトをリストC.1のように用意した。
リストC.1 Kが４以上の偶数の場合のStanスクリプト（ファイル名SDT2AFCEvenCat.stan）
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;
int n_NS[K];
int n_SN[K];
}
parameters {
real mu_s;
real tau;
real<lower = 0.0>
theta[K/2 - 1];
}
transformed parameters {
real C[K - 1];
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);
}
リストC.1のStanスクリプトSDT2AFCEvenCat.stanを利用するPythonスクリプトを、リストC.2のように用意した。
リストC.2 StanスクリプトSDT2AFCEvenCat.stanを利用するPythonスクリプト（ファイル名SDT2AFCEvenCat.py）
import numpy as np
import pystan
import matplotlib.pyplot as plt
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}
sm = pystan.StanModel(file =
'SDT2AFCEvenCat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fout.write('\n{}\n'.format(fit))
mu = fit['mu_s']
tau = fit['tau']
C = fit['C']
p_NS = fit['p_NS']
p_SN = fit['p_SN']
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]))
plt.hist(mu)
plt.xlabel("d'($\mu_S$)",
fontsize = 14)
plt.title("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()
plt.hist(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:
plt.hist(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, 'b--', linewidth =
3, alpha = 0.7, label = 'Data/<N, S>')
plt.plot(xcat, p1, 'g-', linewidth =
3, alpha = 0.7, label = 'Model/<N, S>')
plt.plot(xcat, y2, 'r--', linewidth =
3, alpha = 0.7, label = 'Data/<S, N>')
plt.plot(xcat, p2, 'm-', linewidth =
3, alpha = 0.7, label = 'Model/<S, N>')
plt.xticks(xcat, xlabels, fontsize =
14)
plt.xlabel('Rating Category', fontsize
= 14)
plt.ylabel('Probability/Proportion',
fontsize = 14)
plt.title('Ratings in Conditions
<N, S> and <S, N>', fontsize = 18)
plt.legend(loc = 'upper center',
fontsize = 10)
plt.show()
fout.close()
print('\n{} was
saved.\n'.format(fout_nm))
リストC.1とリストC.2のファイルおよびサンプルデータファイルは、まとめて圧縮ファイルRatingEvenCat.zipとした。ダウンロード解凍すれば、自由に利用できる。
リストC.2のファイルSDT2AFCEvenCat.pyを実行すると、入力データファイル名と出力データファイル名の設定が求められる（図C.1）。
図C.1
出力データファイルは任意のテキストファイル名（ファイル拡張子.txt）であるが、入力データファイルはCSV形式のファイルとして図C.2の様式で用意する。
図C.2 データファイルData4Cat.csv
第１行目は変数名の名前を設定するが、プログラムでは用いないので任意の文字列でよい。ただし、第１行目第１列目に文字列IDを設定すると、Excelの場合、読み込み時にエラーメッセージが表示されることがあるが、これは無視する。第２行目に刺激提示条件＜N、S>（例えば、左にノイズ刺激、右にシグナル刺激）におけるデータを設定する。第２列目に「シグナルは第１位置である（シグナル刺激は左である）」と判断された度数、第３列目に「多分シグナル刺激は第１位置である（多分シグナル刺激は左である）」と判断された度数、第４列目に「多分シグナル刺激は第２位置である（多分シグナル刺激は右である）」と判断された度数、第５列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第３行目に刺激提示条件＜S、N＞（例えば、左にシグナル刺激、右にノイズ刺激）におけるデータを設定する。第２列目に「シグナル刺激は第１位置である（シグナル刺激は左である）」と判断された度数、第3列目に「多分シグナル刺激は左である」と判断された度数、第４列目に「多分シグナル刺激は第２位置である（多分シグナル刺激は右である）」と判断された度数、第５列目に「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数を設定する。第２行目、第３行目とも、第２列目は「シグナル刺激は第１位置である（例えば、シグナル刺激は左である）」と判断された度数、第３列目は「多分シグナル刺激は第１位置である（多分シグナル刺激は左である）」と判断された度数、第４列目は「多分シグナル刺激は第２位置である（多分シグナル刺激は右である）」と判断された度数、第５列目は「シグナル刺激は第２位置である（シグナル刺激は右である）」と判断された度数である。データの設定が終われば、ファイル拡張子として「.csv」を選んで保存すると、CSV形式のファイルとして保存される。
図C.1のように入力データファイル名と出力テキストファイル名を設定すると、計算が始まる。Stanによるサンプリングが終了すると、まず（）の事後分布のヒストグラムが表示される（図C.3）。
図C.3
事後分布の中央値がMed.=1.052、90％CIが[0.849, 1.264]であることが示されている。
図C.3のふぉむを閉じると、τの事後分布のヒストグラムが表示される（図C.4）。
図C.4
事後分布の中央値がMed.=0.433、90％CIが[0.227, 0.631]であることが示されている。
図C.4のフォームを閉じると、カテゴリ境界の事後分布が示される（図C.5）。
図C.5
カテゴリ境界C2=0は制約条件として固定されているので、C2のヒストグラムは表示されていない、
図C.5のフォームを閉じると、各判断カテゴリの確率／比率を示すグラフが折れ線グラフで表示される（図C.6）。
図C.6
データが破線で比率を表し、モデルが実線で各判断カテゴリの確率の中央値である。「シグナル刺激は第１位置である」判断が「カテゴリ１」、「多分シグナル刺激は第１位置である」判断が「カテゴリ２」、「多分シグナル刺激は第２位置である」判断が「カテゴリ３」、「シグナル刺激は第２位置である」判断が「カテゴリ４」で表されている。
図C.6のフォームを閉じると、プログラムの実行終了となる。実行終了後、出力ファイルを適当なエディタで開くと、以下のようになっている。
Input Data File = Data4Cat.csv
<N, S>
[24, 9, 12, 55]
<S, N>
[74, 13, 7, 6]
Inference for Stan model:
anon_model_f2c4de04275e32f42265169574310b77.
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
mu_s 1.05 2.0e-3 0.12 0.81 0.97 1.05 1.14 1.3 4000 1.0
tau
0.43 2.0e-3 0.12 0.19 0.35 0.43 0.52 0.68 4000 1.0
theta[0] 0.51 1.1e-3 0.07 0.38 0.46 0.5 0.55 0.65 4000 1.0
C[0] -0.51 1.1e-3 0.07 -0.65 -0.55 -0.5 -0.46 -0.38 4000 1.0
C[1]
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 4000 nan
C[2] 0.51 1.1e-3 0.07 0.38 0.46 0.5 0.55 0.65 4000 1.0
p_NS[0] 0.22 6.0e-4 0.04 0.15 0.19 0.21 0.24 0.29 4000 1.0
p_NS[1] 0.12 2.5e-4 0.02 0.09 0.11 0.12 0.13 0.15 4000 1.0
p_NS[2] 0.14 3.1e-4 0.02 0.1 0.12 0.14 0.15 0.18 4000 1.0
p_NS[3] 0.53 7.6e-4 0.05 0.44 0.5 0.53 0.56 0.63 4000 1.0
p_SN[0] 0.08 3.8e-4 0.02 0.04 0.07 0.08 0.1 0.13 3404 1.0
p_SN[1] 0.07 1.9e-4 0.01 0.05 0.06 0.07 0.07 0.09 4000 1.0
p_SN[2] 0.1 2.7e-4 0.02 0.07 0.09 0.1 0.11 0.13 4000 1.0
p_SN[3] 0.75 6.5e-4 0.04 0.67 0.73 0.76 0.78 0.83 4000 1.0
lp__ -209.9 0.03 1.22 -213.0 -210.4 -209.6 -209.1
-208.5 1661 1.0
Samples were drawn using NUTS at Wed
Oct 23 15:01:23 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).
d':
Median = 1.052, 90%CI = [0.849, 1.264]
tau:
Median = 0.433, 90%CI = [0.227, 0.631]
C[1]:
Median = -0.503, 90%CI = [-0.623, -0.396]
C[2]:
Median = 0.000, 90%CI = [0.000, 0.000]
C[3]:
Median = 0.503, 90%CI = [0.396, 0.623]
次に、判断カテゴリ数が６個の場合（K=６）の例を示す。
データ例を図C.7に示す。
図C.7 データファイル名Data6Cat.csv
プログラムを起動して、図C.8のように入力データファイル名と出力データファイル名を設定する。
図C.8
入力データファイル名と出力データファイル名を設定すると計算が始まり、Stanによるサンプリングが終了すると、まずd’の事後分布が表示される（図C.9）。
図C.9
図C.9のフォームを閉じると、τの事後分布が表示される（図C.10）。
図C.10
図C.10のフォームを閉じると、カテゴリ境界の事後分布が表示される（図C.11）。
図C.11
図C.11のフォームを閉じると、判断カテゴリの比率と確率が表示される（図C.12）。
図C.12
データが破線で比率を表し、モデルが実線で各判断カテゴリの確率の中央値である。例えば、「確かにシグナル刺激は第１位置である」判断が「カテゴリ１」、「シグナル刺激は第１位置である」判断が「カテゴリ２」、「多分シグナル刺激は第１位置である」判断が「カテゴリ３」、「多分シグナルは第２位置である」判断が「カテゴリ４」、「シグナルは第２位置である」判断が「カテゴリ５」、「確かにシグナル刺激は第２位置である」判断が「カテゴリ６」で表されている。
図C.12のフォームを閉じると、プログラムの実行終了となる。実行終了後、出力ファイルをテキストエディタで開くと以下のようになっている。
Input Data File = Data6Cat.csv
<N, S>
[6, 23, 5, 4, 37, 25]
<S, N>
[42, 38, 2, 5, 11, 2]
Inference for Stan model:
anon_model_f2c4de04275e32f42265169574310b77.
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
mu_s 0.98 1.9e-3 0.12 0.75 0.9 0.98 1.06 1.21 3659 1.0
tau
0.36 1.8e-3 0.11 0.14 0.29 0.36 0.44 0.58 3867 1.0
theta[0] 0.19 7.1e-4 0.04 0.11 0.16 0.19 0.22 0.29 4000 1.0
theta[1] 1.42 1.9e-3 0.11 1.2 1.34 1.41 1.49 1.64 3720 1.0
C[0] -1.61 2.0e-3 0.12 -1.83 -1.68 -1.6 -1.53 -1.39 3409 1.0
C[1] -0.19 7.1e-4 0.04 -0.29 -0.22 -0.19 -0.16 -0.11 4000 1.0
C[2]
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 4000 nan
C[3] 0.19 7.1e-4 0.04 0.11 0.16 0.19 0.22 0.29 4000 1.0
C[4] 1.61 2.0e-3 0.12 1.39 1.53 1.6 1.68 1.83 3409 1.0
p_NS[0] 0.06 2.9e-4 0.02 0.03 0.05 0.06 0.07 0.1 3717 1.0
p_NS[1] 0.22 4.2e-4 0.03 0.17 0.21 0.22 0.24 0.28 4000 1.0
p_NS[2] 0.05 1.7e-4 0.01 0.03 0.04 0.05 0.05 0.07 4000 1.0
p_NS[3] 0.05 1.9e-4 0.01 0.03 0.04 0.05 0.06 0.08 4000 1.0
p_NS[4] 0.37 4.7e-4 0.03 0.32 0.36 0.37 0.39 0.43 3591 1.0
p_NS[5] 0.24 6.1e-4 0.04 0.17 0.22 0.24 0.27 0.33 4000 1.0
p_SN[0] 0.02 1.5e-4 8.1e-3 7.8e-3 0.01 0.02 0.02 0.04 2755 1.0
p_SN[1] 0.12 3.6e-4 0.02 0.08 0.11 0.12 0.14 0.16 3664 1.0
p_SN[2] 0.03 1.2e-4 7.5e-3 0.02 0.03 0.03 0.04 0.05 4000 1.0
p_SN[3] 0.04 1.5e-4 9.3e-3 0.02 0.03 0.04 0.04 0.06 4000 1.0
p_SN[4] 0.36 4.7e-4 0.03 0.31 0.34 0.36 0.38 0.42 4000 1.0
p_SN[5] 0.43 7.3e-4 0.05 0.34 0.4 0.43 0.46 0.52 4000 1.0
lp__ -298.2 0.03 1.48 -301.9 -298.8 -297.8 -297.1
-296.3 1982 1.0
Samples were drawn using NUTS at Wed
Oct 23 15:26:01 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).
d':
Median = 0.981, 90%CI = [0.789, 1.178]
tau:
Median = 0.363, 90%CI = [0.177, 0.547]
C[1]:
Median = -1.605, 90%CI = [-1.799, -1.424]
C[2]:
Median = -0.187, 90%CI = [-0.271, -0.124]
C[3]:
Median = 0.000, 90%CI = [0.000, 0.000]
C[4]:
Median = 0.187, 90%CI = [0.124, 0.271]
C[5]:
Median = 1.605, 90%CI = [1.424, 1.799]
２AFC課題では、シグナル刺激感覚の平均値と分散を個別に推定することはできない。このことをグラフで見てみる。
いま、ノイズ刺激感覚
とシグナル刺激感覚
が次の正規分布に従うとする。
このとき
となるので、正答確率は
で与えられる。
正答確率
ののグラフは、図Ex.1のようである。
図Ex.1
図Ex.1のグラフはリストEx.1のスクリプトで描いている。スクリプトファイルは、圧縮ファイルEx2AFCR.zipに纏めたので、ダウンロード展開して自由に利用されたい。
シグナル感覚の平均値
が大きくなると、正答確率も上がるが、標準偏差
が大きくなると正答確率は下がる。
同じ等高線上のとの値の組は同じ正答確率Pcを与えるので、正答確率Pcからとの値の組を決めることはできない。との値の組は、２AFC課題を以下のように拡張すると求めることができる。
２AFCでは、提示刺激対は＜ノイズ、シグナル＞と＜シグナル、ノイズ＞の2種類である。これに、＜ノイズ、ノイズ＞と＜シグナル、シグナル＞の２種類を加えて、４種類の刺激対に対して評定法を用いると、シグナル刺激感覚の平均値と分散の推定を行うことができる。
２AFC課題では、各試行において２つの刺激（刺激対）が与えられる。各刺激の提示位置（空間的、時間的）を＜S1, S2＞で表す。第1提示位置がS1、第２提示位置がS2である。＜n, s＞（あるいは、n-s）は、第１位置にノイズ刺激、第２位置にシグナル刺激が提示されたことを表す。
第１提示位置の刺激感覚を
、第２刺激提示位置の刺激感覚を
で表す。
第１刺激位置に対する第２刺激位置のバイアスをbで表す。このバイアスは恒常誤差のように感覚のものと反応のものの両者を含むものである。
このとき、以下のようになる。
のとき、
のとき、
のとき、
のとき、
いま、判断がK個のカテゴリによる評定で回答されるものとし、判断のカテゴリ境界を以下のようにおく。
このとき、
であれば、提示刺激対<
>に対する判断
はが与えられるとする。すなわち、
である。
例えば、判断カテゴリー数が
個のときは、評定判断は以下のようなものが考えられる。
： シグナル刺激は確かに第１位置である。
： シグナル刺激は第１位置である。
： シグナル刺激はたぶん第１刺激である。
： わからない／どちらも同じと思われる。
： シグナル刺激はたぶん第２位置である。
： シグナル刺激は第２位置である。
： シグナル刺激は確かに第２位置である。
平均
、分散
の正規分布の累積分布を
で表し、以下のように便宜的に設定する。
このとき、差
の分布は以下のように与えられる。
刺激対 
が提示されとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
いま、カテゴリ境界は、原点を中心として対称に設定されているとする。
したがって、Kが偶数であれば、
奇数であれば、
上のモデルに基づいてベイズ分析するためのStanスクリプトを、Kが偶数の場合をリストEx.２に、Kが奇数の場合をリストEx.3に示す。これらのスクリプトファイルを用いて拡張２AFC評定法によるデータをベイズ分析するスクリプトをリストEx.4に示す。これらのファイルは圧縮ファイルEx2AFCR.zipにまとめた。ダウンロード展開して自由に利用されたい。
なお、PyStan 3用のスクリプトをこのウェブサイトに用意している。
リストEx.4の実行例を以下に示す。スクリプトの実行は、WindowsにインストールされたAnacondaのPython3.7仮想環境においてPyStan2.17.1.0を用いて行った。
Stanスクリプトの実行時に、コンパイルに関わるエラーが出たときは、C++コンパイラをインストールしてみるとよい。使用PCでは、Visual StudioのC++がインストールされている。
Ubuntuの場合も、Stanスクリプトの実行時にコンパイルに関わるエラーが表示されたときは、C++コンパイラをインストールしてみるとよい。C++を含めた複数の言語に対応するものとしてg++がある。Ubuntuの場合、コマンド「g++」を実行すると、「g++は見つかりません。次のコマンドでインストールして下さい。」というようなメッセージが表示されたら、指示通り、コマンドを入力実行するとインストールされる。g++のインストール後は、Stanスクリプトはコンパイル実行されるはずである。
まず、カテゴリ数が６個の偶数の場合について示し、続いてカテゴリ数が７個の奇数の場合を示す。
表Ex.1のようなデータが与えられたとする。
|
表Ex.1 カテゴリ数が６個の場合 |
||||
|
判断カテゴリR |
n-n |
n-s |
s-n |
s-s |
|
1 |
2 |
1 |
6 |
6 |
|
2 |
5 |
1 |
10 |
3 |
|
3 |
15 |
13 |
24 |
12 |
|
4 |
21 |
21 |
7 |
21 |
|
5 |
6 |
7 |
3 |
6 |
|
6 |
1 |
7 |
0 |
2 |
リストEx.4のファイルAnalEx2AFCR.pyを実行する。コマンドで実行する場合は、以下のようである。
(py37)…> python AnalEx2AFCR.py
プロンプト “n_n =” が表示されたら、カテゴリ１から６までの提示刺激対条件<noise, noise>に対する判断回答数を入力する。表Ex.1のデータであれば、次のように数値を半角空白文字で区切って並べて最後にEnterキーを押して入力する。
n_n = 2 5 15
21 6 1
Enterキーを押すと、次の入力プロンプトが表示されるので、同様に入力していく。
n_n = 2 5 15
21 6 1
n_s =
表Ex.1のデータをすべて入力すると、端末画面は以下のようになっている。
n_n = 2 5 15 21 6 1
n_s = 1 1 13 21 7 7
s_n = 6 10 24 7 3 0
s_s = 6 3 12 21 6 2
４つの刺激対に対するデータの入力が終わると、MCMCサンプリングが始まる。サンプリングが終了すると、まず、の事後分布とMAP推定値が表示される（図Ex.1）。
図Ex.1
MAP推定値が良いことは、岡本（2022；https://doi.org/10.14947/psychono.41.1）に報告されている。
図Ex.1のフォームを閉じると、の事後分布とMAP推定値が表示される（図Ex.2）。
図Ex.2
図Ex.2のフォームを閉じると、カテゴリ境界値C1, C2, C4, C5の事後分布とMAP推定値のグラフが表示される（図Ex.3）。
図Ex.3
カテゴリ数が偶数であるので、C3は０に固定されているが、点（０，０）の赤い半円で表されている。
図Ex.3のフォームを閉じると、パラメータｂの事後分布とMAP推定値のグラフが表示される（図Ex.4）。
図Ex.4
図Ex.4のフォームを閉じると、図Ex.5のグラフが表示される（図Ex.5）。
図Ex.5
各提示刺激対に対して、評定カテゴリ判断の累積確率のデータ値（縦軸）とモデルによる予測値（横軸）のグラフが表示されている。破線の対角線は、データ値とモデル予測値が完全に一致した場合のグラフである。データ値と予測値の関係は、この対角線にほぼ沿っているので、モデルのデータに対する当てはまりは良いといえる。
図Ex.5のフォームを閉じると、プログラムの終了である。
グラフの図は、スクリプトファイルと同じフォルダにファイルとして保存されている。
リストEx.4のスクリプトAnalEx2AFCR.pyを次のように実行する。
(py37) …> python AnalEx2AFCR.py
まず、提示刺激対＜ノイズ、ノイズ＞に対するデータの入力が求められる。
表Ex.2の場合のデータを以下のように入力する。
|
表Ex.2 カテゴリ数７個の場合のデータ |
||||
|
評定カテゴリ |
n-n |
n-s |
s-n |
s-s |
|
1 |
2 |
0 |
5 |
0 |
|
2 |
7 |
3 |
9 |
7 |
|
3 |
9 |
10 |
14 |
13 |
|
4 |
17 |
14 |
11 |
7 |
|
5 |
11 |
9 |
10 |
13 |
|
6 |
1 |
10 |
1 |
9 |
|
7 |
3 |
4 |
0 |
1 |
各値を半角空白文字（２個以上でも可）で区切って並べ、最後にEnterキーを押す。
n_n = 2 7 9
17 11 1 3
Enterキーを押すと、次の入力プロンプトが表示されるので、入力していく。
n_n = 2 7 9 17 11 1 3
n_s = 0 3 10 14 9 10 4
s_n = 5 9 14 11 10 1 0
s_s = 0 7 13 7 13 9 1
＜シグナル、シグナル＞に対する入力を終えると、ＭＣＭＣサンプリングが始まる。
サンプリングが終了すると、まず、の事後分布とＭＡＰ推定値が表示される（図Ex.6）。
図Ex.6
MAP推定値が良いという報告がある（岡本、2022；https://doi.org/10.14947/psychono.41.1）
図Ex.6のフォームを閉じると、の事後分布とMAP推定値が表示される（図Ex.7）。
図Ex.7
図Ex.7を閉じると、評定境界値C1, C2, C3, C4, C5, C6の事後分布とMAP推定値が表示される（図Ex.8）。
図Ex.8
図Ex.8を閉じると、パラメータｂの事後分布とＭＡＰ推定値が表示される（図Ex.9）。
図Ex.9
図Ex.9のフォームを閉じると、図Ex.10のグラフが表示される。
図Ex.10
各提示刺激対に対して、評定カテゴリ判断の累積確率のデータ値（縦軸）とモデルによる予測値（横軸）のグラフが表示されている。破線の対角線は、データ値とモデル予測値が完全に一致した場合のグラフである。データ値と予測値の関係は、この対角線にほぼ沿っているので、モデルのデータに対する当てはまりは良いといえる。
図Ex.10のフォームを閉じると、プログラムの終了である。
グラフの図は、スクリプトファイルと同じフォルダにファイルとして保存されている。
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import scipy.stats as ss
X = np.linspace(0.1, 2.0, 1000) # mu_s
Y = np.linspace(0.4, 2.0, 100) # sigma_s
Xg, Yg = np.meshgrid(X, Y)
def Pc(mu, sgm):
return ss.norm.cdf(mu / ((1 + sgm**2)**0.5))
Z = np.empty((len(Y), len(X)))
for i in range(len(X)):
for j in range(len(Y)):
Z[j][i] = Pc(X[i], Y[j])
fig = plt.figure(figsize = (7,6))
ax = fig.add_subplot(111)
C = ax.contourf(Xg, Yg, Z, levels =
20)
ax.set_xlabel('$\mu_s$', fontsize =
16);
ax.set_ylabel('$\sigma_s$', fontsize =
16)
ax.set_title('Pc')
fig.colorbar(C, ax=ax, fraction=0.02,
pad=0.1)
plt.tight_layout()
plt.savefig('FigPcCont.png')
plt.show()
data {
int K;
int nn_cond[K];
int ns_cond[K];
int sn_cond[K];
int ss_cond[K];
}
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);
}
data {
int K;
int nn_cond[K];
int ns_cond[K];
int sn_cond[K];
int ss_cond[K];
}
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);
}
リストEx.4 拡張２肢強制選択評定法データのベイズ分析スクリプト（AnalEx2AFCR.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import pystan
import sys
def CalcMAPEst(samples, a = 0.05,
n_points = 10000):
"""
Calculatte
a MAP estimate from a KDE graph on [Lp, Up]
Lp
and Up are 100*a/2 and 100(1-a/2) percentile points of samples
"""
Lp, Up =
np.percentile(samples, [100 * a/2, 100 * (1 - a/2)]) # import numpy as np
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
# import scipy.stats as ss
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
n_n = [int(v) for v in input('n_n =
').split()]
K = len(n_n)
if K < 3:
print('Number of categories
should be larger than 2.')
sys.exit()
n_s = [int(v) for v in input('n_s =
').split()]
if len(n_s) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
s_n = [int(v) for v in input('s_n =
').split()]
if len(s_n) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
s_s = [int(v) for v in input('s_s =
').split()]
if len(s_s) != K:
print('Numbers of categories
are not consistent.')
sys.exit()
print('n_n =', n_n)
print('n_s =', n_s)
print('s_n =', s_n)
print('s_s =', s_s)
sm = pystan.StanModel(file =
'SDT_2AFC_KCat_odd.stan' if K % 2 == 1 else 'SDT_2AFC_KCat_even.stan')
print('Compiling ended.')
Data = {'K':K, 'nn_cond':n_n,
'ns_cond':n_s, 'sn_cond':s_n, 'ss_cond':s_s}
fit = sm.sampling(data = Data, iter =
10000, n_jobs = 1)
print(fit)
mu_s_map = CalcMAPEst(fit['mu_s'])[0]
sb.kdeplot(fit['mu_s'])
plt.xlabel('$\mu_s$')
plt.title('$\mu_s$(MAP) =
{0:.3f}'.format(mu_s_map))
plt.savefig('Fig_mu.png')
plt.show()
sgm_s_map =
CalcMAPEst(fit['sgm_s'])[0]
sb.kdeplot(fit['sgm_s'])
plt.title('$\sigma_s$(MAP) =
{0:.3f}'.format(sgm_s_map))
plt.xlabel('$\sigma_s$')
plt.savefig('Fig_sgm.png')
plt.show()
s_title = ''
if K % 2 == 1:
for k in range(K-1):
sb.kdeplot(fit['C'].T[k])
c_map = CalcMAPEst(fit['C'].T[k])[0]
s_title += f'C{k+1}={c_map:.2f}'
if k
< K-2:
s_title += ', '
else:
for k in range(K-1):
if k
!= (K//2) - 1:
sb.kdeplot(fit['C'].T[k])
c_map = CalcMAPEst(fit['C'].T[k])[0]
s_title += f'C{k+1}={c_map:.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('MAP estimates\n' + s_title)
plt.savefig('FigC.png')
plt.show()
b_map = CalcMAPEst(fit['b'])[0]
sb.kdeplot(fit['b'])
plt.title('b(MAP) =
{0:.3f}'.format(b_map))
plt.savefig('Fig_b.png')
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]
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.savefig('FigCumEstObs.png')
plt.show()
データ分析用のPython解説書＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を用意した。