３点識別試験のベイズ分析
グリッド法
３つの刺激の内、２つが同じ刺激Aで１つが異なる刺激Bであるとき、３つの中で他の２つと異なるものとして正しくBが選ばれる頻度を調べる実験（３点識別試験：triangular test, 3AFC Oddity）のデータのベイズ分析によってAとBの感覚量の差
の推定を行うプログラムを作成してみた。プログラムのスクリプトファイルは圧縮ファイルTriTestScriptFiles.zipとしたので、ダウンロード展開して自由に御利用頂きたい。
信号検出理論に基づいて、刺激Aおよび刺激Bの感覚量が、それぞれ独立に平均
および平均
の分散
の正規分布（ガウス分布）に従うものとする（サーストンモデル）。２つの刺激Aの感覚量を
および
で表し、刺激Bの感覚量を
で表すと
である。
３点試験法において正しくBが選ばれるのは
が成り立つときである（Frijters et al., 1980, p. 177）。
いま、信号検出理論の慣例に倣って、感覚量の平均値の差を標準偏差
を単位として表し、
とおく。このとき、正答確率（Bが選ばれる確率）
と感覚量の差
の関係は次式で表される（Frijters et al., 1980, eq.(1)）。
式（1）において、
の正負を入れ換え（を
で置き換え）ても値は変わらない。したがって、感覚量の座標の向きをBの平均の方がAの平均より大きくなる方向にとって
とおいても一般性は失われない。式（１）より、感覚量の差に対する正答確率が与えられる。
総試行数
の実験において、Bが選ばれた試行数をｋとおくと、ｋは試行数Nの２項分布における成功試行数と考えて、確率は
で与えられる。この確率モデル（２）に基づいての事後分布をグリッド法に依って求めるスクリプトをリスト１に示す。
リスト１のスクリプトを実行すると、まずNの入力が求められる。
>>>
= RESTART:
D:\xxx\TriTestScriptFiles\TriangularTest_grid.py
N =
Nの値を入力すると、次にｋの値の入力が求められる。
N = 50
k =
ｋの値を入力する。次の例では、Ｎの値として50を入力した後、ｋの値20が設定されている。
N = 50
k = 20
計算が終了すると、d’の事後分布が図のように表示される。
図１
図中キャプションに、事後分布のモード（d’の推定値）と95％HDIが表示されている。
図は、次のコードによって保存される。
plt.savefig(f'FigN{N}k{k}.png')
図１の場合は、FigN50k20.pngである。
以下に、Ｎ=50に対するＫ＝10、30、40、50の場合の分析結果を示す。
N=50、k＝10の場合は、正答数ｋがランダムな回答の確率1/3より小さいので、事後分布は左端により、モードはmode =0である。
これに対して、Ｎ＝50、k=50は全試行正答であるので、d’の値はいくらでも大きい値であってもデータと矛盾しないが、d’の事前分布が10までの範囲に設定されているので、事後分布は右端に寄った分布となっている。
引用文献
Frijters, J. E. R., Kooistra, A., & Vereijken, P. F. G. (1980). Tables of d' for the triangular method and the 3-AFC signal detection procedure. Perception & Psychophysics, 1980, 27, 176-178.
"""
Yasuharu Okamoto, 2021.11, 2022.01
"""
import numpy as np
import scipy.stats as ss
import scipy.integrate as si
import matplotlib.pyplot as plt
import seaborn as sb
N = int(input('N = '))
k = int(input('k = '))
class Triangular:
def __init__(self, d_prime):
self.d_prime = d_prime
def func(self, u):
v = ss.norm.cdf(-u
* (3.0**0.5) + self.d_prime * ((2/3)**0.5)) + \
ss.norm.cdf(-u * (3.0**0.5) - self.d_prime * ((2/3)**0.5))
return v * np.exp(-0.5 * (u**2)) / ((2.0 * np.pi)**0.5)
def Pc_d_prime(d_prime):
d_func = Triangular(d_prime).func
return 2 * si.quad(d_func,
0.0, +np.inf)[0]
dprimes = np.linspace(0, 10, 1001)
Pc = np.empty(len(dprimes))
for i in range(len(dprimes)):
if i % 100 == 0:
print('{0:}/{1:}
'.format(i, len(dprimes)))
#, end = '\r')
Pc[i] =
Pc_d_prime(dprimes[i])
postPdprime = np.empty(len(dprimes))
for i in range(len(dprimes)):
postPdprime[i] = (Pc[i] **
k) * ((1 - Pc[i]) ** (N - k))
postPdprime = (postPdprime /
np.sum(postPdprime))
map_idx = np.argmax(postPdprime)
MAP_Est = map_idx / 100
print('MAP_Est =', MAP_Est)
postPcum = np.cumsum(postPdprime)
a = 0.05
Lp = 0
Rp = 0
w_t = 1000
L_t = Lp
R_t = Rp
while True:
if Lp == 0:
cumPL = 0.0
else:
cumPL = postPcum[Lp-1]
while 1.0 - postPcum[Rp]
>= a - cumPL:
Rp
+= 1
if Lp == 0:
L_t
= Lp
R_t
= Rp
w_t
= R_t - L_t
elif Rp - Lp < w_t:
L_t
= Lp
R_t
= Rp
w_t
= R_t - L_t
Lp += 1
if postPcum[Lp] >= a:
Lp
-= 1
break
Lp = 1000
Rp = 1000
w_t = 1000
L_t1 = Lp
R_t1 = Rp
while True:
if Rp == 1000:
cumPR = 0.0
else:
cumPR = 1.0 - postPcum[Rp]
while postPcum[Lp] >= a -
cumPR:
if
Lp > 0:
Lp -= 1
else:
break
if Rp == 1000:
L_t1
= Lp
R_t1
= Rp
w_t
= R_t1 - L_t1
elif Rp - Lp < w_t:
L_t1
= Lp
R_t1
= Rp
w_t
= R_t1 - L_t1
Rp -= 1
if 1.0 - postPcum[Rp] >=
a:
Rp
+= 1
break
HDI_L = ((L_t + L_t1) / 2.0) / 100.0
HDI_R = ((R_t + R_t1) / 2.0) / 100.0
plt.plot(np.linspace(0, 10, 1001),
postPdprime)
plt.plot([MAP_Est,MAP_Est], [0,
postPdprime[map_idx]], label = 'Mode = {}'.format(MAP_Est))
plt.title(f"Psterior Distribution
for d'\nN = {N}, k = {k}", fontsize = 14)
plt.xlabel("d'", fontsize =
12)
plt.plot([0, 10], [0, 0], lw = 1, c =
'k')
plt.yticks([0])
plt.plot([HDI_L, HDI_R], [0, 0], lw =
5, label = '95% HDI = [{0:.2f}, {1:.2f}]'.format(HDI_L, HDI_R))
plt.legend()
plt.savefig(f'FigN{N}k{k}.png')
plt.show()