３点識別試験
弁別確率のベイズ推定
３つの刺激の内、２つが同じ刺激Aで１つが異なる刺激Bであるとき、３つの中で他の２つと異なるものとして正しくBが選ばれる頻度を調べる実験（３点識別試験：triangular test）を想定する。なお、刺激Aに対して、2種類の刺激特性BとCを考え、BとCの組み合わせに対する交互作用をロジスティック回帰モデル（Gelmanら、2014、3.7節；岡本、2014、第４章）により分析するモデルを、本ウェブサイトの後半で扱う。
本ウェブサイトのプログラムファイルは圧縮ファイルTriTestPdFiles.zipとして纏めたので、自由にダウンロードして利用されたい。
いま、正しくBが異なると判る確率を
と置く。このPdはproportion of distinguishersを表すとされているが、信号検出理論（signal detection theory）に倣えばprobability of
detection/discriminationを表すと考えた方が使い易いと思われる。同じ観察者に、配置を変えて複数回提示する場合は、確率と考えればよい。確率PdでBが異なることが検出できるとき、確率
でどれが異なるか判らないということになる。このとき、３つの刺激からランダムに１つが他の２つと異なるものとして選ばれるが、このランダムに選ばれた刺激がBである確率は
である。したがって、正しくBが異なるものとして検出される場合と、ランダムにBが選ばれる場合とを合わせると、Bが選ばれる確率は
となる。確率
は見かけの正答確率Pc（probability of correct response）である。
総試行数
の実験において、Bが選ばれた試行数をｋとおくと、ｋは試行数Nの２項分布における成功試行数と考えて、確率は
で与えられる。
確率モデル（２）に基づいてPdの事後分布を求めるStanスクリプトを本ウェブサイトの後半に示すリスト１のように用意した。
リスト１のStanスクリプトを実行するPythonスクリプトをリスト２のように用意した。
リスト２のスクリプトを実行すると、総試行数Nと正答試行数ｋの入力が求められる。
(py39) PS D:\XXXXX\TriTestPdFiles>
python .\TriTest.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
N = 50
k = 25
上の例では、Nの値として50、ｋの値として25が設定されている。
ｋの値を設定すると、計算が始まる。
計算が終了すると、Pdの事後分布が図１のように表示される。
図１
図１では、弁別確率Pdのモードが0.27であること、95％HDIが[0.06, 0.44]であることが、キャプションに表示されている。
図１のグラフは、ファイルFigPd.pngとして保存されている。
まず、刺激Aに対して、2つの刺激特性BとCについて、Bのみの刺激B、Cのみの刺激C、およびBとCの２つの特性をもつ刺激BCについて、それぞれ３点識別試験を行った場合を考える。特性BあるいはCが単独の場合に比べて、BとCが同時の場合、その効果を調べるためにロジスティック回帰モデルを用いる。
続いて、識別確率が年齢とともにどう変わるかをロジスティック回帰モデルで調べる場合を、この後半で扱う。
ロジスティック回帰モデルでは、確率ｐのロジット変換値が線形モデルで表されるとする（Gelmanら、2014，3.7節；岡本、2014、第４章）。
ロジット変換は次式で表される。
これより、次式を得る。
ロジスティック回帰モデルでは、
が回帰モデルで表される。
刺激B、C、BCの弁別確率PdB、PdC、PdBCを、刺激Bの主効果fB、刺激Cの主効果fC、交互採用fBCを設定して次式のように表す。
上式のモデルのStanスクリプトをリスト３に示した。スクリプトファイルは、圧縮ファイルTriTestPdFiles.zipにまとめた。ダウンロード・展開すれば、自由に利用できる。
リスト３のStanスクリプトを用いるPythonスクリプトはリスト４に示した。
リスト４のスクリプトを実行すると、まず刺激B、刺激C、刺激BCに対する正答試行数ｋおよび総試行数の入力が求められる。次の例では、刺激Ａ、A’、Ｂの３点の組み合わせの場合、成功試行数kBが25試行、総試行数NBが50試行のデータであることを示している。刺激ＣおよびＢＣに対しても同様である。
(py39) PS D:\xxxxx\TriTestPdFiles>
python .\TriTestLogistic.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
kB = 25
NB = 50
kC = 17
NC = 50
kBC = 40
NBC = 50
最後のNBC = 50を入力すると、計算が始まる。
計算が終了すると、効果fB、fC、fBCの事後分布のグラフが表示される（図２）。
図２
主効果fBとfCは、負の値の方向に絶対値が大きいほど効果がない（Pdが０に近付く）ことを表すが、交互作用fBCは０であれば交互作用は無い、０からはなれた値であれば交互作用ありとなる。図５ではfBCの95％HDIは［2.75, 10.09］と原点０から右方向に分布が明瞭にシフトしている。
図２のフォームを閉じると、次に弁別確率Pdの事後分布の統計量の計算が始まり、終了すると図３のグラフが表示される。
図３
刺激BのPdはモードが0.23であるが、刺激CのPdはモードが0.01と０に近い。すなわち、刺激BはAと確率0.23で弁別できているが、刺激Cは弁別できていない。刺激BとCの混合では、弁別確率PdBCはモードが0.70でよく弁別されていると言える。刺激Cは隠し味みたいなもので、刺激Bと組み合わされるとその効果が明瞭になると言えるかもしれない。
図２と図３のグラフは、自動的にファイル名Figfctrs.pngとFigPd.pngに保存される。
データの例として、上とは異なる次の場合を分析してみた。
(py39) PS D:\xxxxd\TriTestPdFiles>
python .\TriTestLogistic.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
kB = 30
NB = 50
kC = 30
NC = 50
kBC = 20
NBC = 50
各効果のグラフは図４のようになった。
図４
交互作用fBCの事後分布が０から負の方向にかなり寄っている。
図５
刺激Bと刺激Cの弁別確率PdBとPdCはモードが0.38でかなり弁別されていると言えるが、刺激BCの場合は、刺激Bおよび刺激Cの複合であるのに関わらずPdBCのモードは0.00であり、弁別されていない。２つの刺激BとCの混合刺激は弁別が難しくなり、これは交互作用fBCの分布が負の領域に偏っていることに表されている（図４）。
上で説明したように、ロジット変換は次式で表される。
これより、次式を得る。
ロジスティック回帰モデルでは、が回帰モデルで表される。いまの場合、独立変数である年齢をAで表せば、
である。このとき、識別確率Pdは次式で与えられる。
したがって、正答確率Pcは
となる。
上のモデルをStanスクリプトで表したものをリストRA1に示す。
リストRA1のStanスクリプトを利用するPythonスクリプトをリストRA2に示す。
リストRA2のプログラム（TriLogRegAge.py）では、入力データは図RA1に示す様式のCSVファイルとして用意する。
・
・
・
図RA1
１行目に変数名を書く。１列目がでーたの識別変数で数値でも文字列でもよい。２列目に判断の正誤をかく。正答であれば１，誤答であれば０である。３列目に回帰モデルの独立変数の値を書く。図RA1のデータでは、データの識別変数の値が１である先頭のデータは、判断が正答（１で表されている）で、独立変数（図RA1の場合の場合、年齢）が33歳である。
図RA1のように作成されたデータは、保存のときにファイル拡張子として「.csv」を選ぶと、CSV形式で保存される。
リストRA2のスクリプト（TriLogRegAge.py）を実行すると入力ファイル名の設定が求められる。
(py39) PS D:\xxxxx\TriTestPdFiles>
python TriLogRegAge.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
Input file (*.csv) = Data.csv
上の例では、入力ファイル名Data.csvが設定されいる。
入力ファイル名を設定して、Enterキーを押すと計算が始まる。MCMCサンプリングが終了すると、まずパラメータaとbの事後分布のグラフが表示される（図RA2）。
図RA２
図RA２の表示されているフォームを閉じると、図RA３のフォームが表示される。
図RA３
図RA３より、独立変数Aの係数bの事後分布は負の領域にあり、mode≒-1.09（この値は、MCMCサンプリングに基づくものであるので、プログラムの実行ごとに多少変動がある）で、独立変数Aの値が大きくなる程、識別確率Pdと正答確率Pcは小さくなると考えられる。
図RA３のフォームを閉じると、図RA４のフォームが表示される。
図RA４
識別確率Pd（赤曲線）と正答確率Pc（青曲線）のグラフは、事後分布のmodeをパラメータ値とするグラフである。
緑の小円は、年齢の範囲を４等分した区間に分け、各区間内の年齢の正答比率である。正答確率Pcのグラフとデータの小円はよく合っていると言える。
なお、Python用のStanであるPyStanの説明を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, third edition. CRC Press.
リスト１ ３点試験法のStanスクリプト（ファイル名tri_test.stan）
data {
int N;
int k;
}
parameters {
real<lower = 0.0, upper =
1.0> Pd; // Discrimination Probability
}
transformed parameters {
real Pc;
// Correct Response Probability
Pc = Pd + (1.0 - Pd) *
(1/3.0);
}
model {
Pd ~ uniform(0.0, 1.0);
k ~ binomial(N, Pc);
}
リスト２ リスト１のStanスクリプトを利用するPythonスクリプト（ファイル名TriTest.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
import pickle
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
"""
print('Calculating the
mode.')
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]
def calcHDI(sample, a = 0.05):
"""
Calculation of 100(1-a)% HDI for sample
Chen-Shao HPD Estimation algorithm
"""
print('Calculating HDI.')
smpl = np.sort(sample)
n = len(smpl)
ckW = max(smpl) - min(smpl)
cki = 0
RB = int(n * (1-a))
L0 = 0
U0 = n-1
while (RB + cki < n):
if
smpl[RB + cki] - smpl[cki] < ckW:
ckW = smpl[RB + cki] - smpl[cki]
L0 = cki
U0 = RB + cki
cki
+= 1
ckW = max(smpl) - min(smpl)
cki = 1
LB = n - int(n * (1-a))
L1 = 0
U1 = n-1
while (LB - cki >= 0):
if
smpl[n - cki] - smpl[LB - cki] < ckW:
ckW = smpl[n - cki] - smpl[LB - cki]
U1 = n - cki
L1 = LB - cki
cki
+= 1
return smpl[int((L1+L0)/2)],
smpl[int((U1+U0)/2)]
N = int(input('N = '))
k = int(input('k = '))
Data = {'N': N, 'k': k}
sm = pystan.StanModel(file =
'tri_test.stan')
fit = sm.sampling(data = Data, iter =
10_000, n_jobs = 1)
print(fit)
sb.kdeplot(fit['Pd'])
plt.hist(fit['Pd'], density = True,
color = 'g', alpha = 0.5)
MAPEst = CalcMAPEst(fit['Pd'])[0]
plt.scatter([MAPEst], [0], marker =
'o', s = 200,
c = 'magenta', label = f'mode = {MAPEst:.2f}')
L_bndry, R_bndry = calcHDI(fit['Pd'],
0.05)
plt.plot([L_bndry, R_bndry], [0, 0], c
= 'r', lw = 5,
label = '95% HDI [{0:.2f}, {1:.2f}]'.format(L_bndry, R_bndry))
plt.yticks([])
plt.xlabel('Pd', fontsize = 12)
plt.title(f'Posterior distribution\nN
= {N}, k = {k}', fontsize = 16)
plt.legend()
plt.savefig('FigPd.png')
plt.show()
リスト３ ロジスティック回帰モデルのStanスクリプト（ファイル名tri_test_logistic.stan）
//
// Yasuharu
Okamoto, 2021.12
//
data {
int NB;
int kB;
int NC;
int kC;
int NBC;
int kBC;
}
parameters {
real fB;
real fC;
real fBC;
}
transformed parameters {
real<lower = 0.0, upper =
1.0> PdB;
real<lower = 0.0, upper =
1.0> PdC;
real<lower = 0.0, upper =
1.0> PdBC;
real<lower = 0.0, upper =
1.0> PcB;
real<lower = 0.0, upper =
1.0> PcC;
real<lower = 0.0, upper =
1.0> PcBC;
PdB = 1.0 / (1.0 +
exp(-fB));
PdC = 1.0 / (1.0 +
exp(-fC));
PdBC = 1.0 / (1.0 + exp(-(fB
+ fC + fBC)));
PcB = PdB + (1.0 - PdB) *
(1.0 / 3.0);
PcC = PdC + (1.0 - PdC) *
(1.0 / 3.0);
PcBC = PdBC + (1.0 - PdBC) *
(1.0 / 3.0);
}
model {
fB ~ normal(0.0, 5.0);
fC ~ normal(0.0, 5.0);
fBC ~ normal(0.0, 5.0);
kB ~ binomial(NB, PcB);
kC ~ binomial(NC, PcC);
kBC ~ binomial(NBC, PcBC);
}
リスト４ リスト３のStanスクリプトをPythonスクリプト（ファイル名TriTestLogistic.py）
import numpy as np
import scipy.stats as ss
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
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
"""
print('Calculating the
mode.')
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]
def calcHDI(sample, a = 0.05):
"""
Calculation of 100(1-a)% HDI for sample
Chen-Shao HPD Estimation algorithm
"""
print('Calculating the
HDI.')
smpl = np.sort(sample)
n = len(smpl)
ckW = max(smpl) - min(smpl)
cki = 0
RB = int(n * (1-a))
L0 = 0
U0 = n-1
while (RB + cki < n):
if
smpl[RB + cki] - smpl[cki] < ckW:
ckW = smpl[RB + cki] - smpl[cki]
L0 = cki
U0 = RB + cki
cki
+= 1
ckW = max(smpl) - min(smpl)
cki = 1
LB = n - int(n * (1-a))
L1 = 0
U1 = n-1
while (LB - cki >= 0):
if
smpl[n - cki] - smpl[LB - cki] < ckW:
ckW = smpl[n - cki] - smpl[LB - cki]
U1 = n - cki
L1 = LB - cki
cki
+= 1
return smpl[int((L1+L0)/2)],
smpl[int((U1+U0)/2)]
kB = int(input('kB = '))
NB = int(input('NB = '))
kC = int(input('kC = '))
NC = int(input('NC = '))
kBC = int(input('kBC = '))
NBC = int(input('NBC = '))
sm = pystan.StanModel(file =
'tri_test_logistic.stan')
Data = {'kB': kB, 'NB': NB, 'kC': kC,
'NC': NC, 'kBC': kBC, 'NBC': NBC}
fit = sm.sampling(data = Data, iter =
10_000, n_jobs = 1)
print(fit)
MAPfB = CalcMAPEst(fit['fB'])[0]
MAPfC = CalcMAPEst(fit['fC'])[0]
MAPfBC = CalcMAPEst(fit['fBC'])[0]
L_HDIfB, R_HDIfB = calcHDI(fit['fB'])
L_HDIfC, R_HDIfC = calcHDI(fit['fC'])
L_HDIfBC, R_HDIfBC =
calcHDI(fit['fBC'])
sb.kdeplot(fit['fB'], label = 'fB,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfB, L_HDIfB, R_HDIfB))
sb.kdeplot(fit['fC'], label = 'fC,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfC,
L_HDIfC, R_HDIfC))
sb.kdeplot(fit['fBC'], label = 'fBC,
mode = {0:.2f}, \n
95% HDI = [{1:.2f}, {2:.2f}]'.
format(MAPfBC, L_HDIfBC, R_HDIfBC))
plt.xlabel('Factor', fontsize = 14)
plt.yticks([])
plt.title('Posterior distribtion\n' +
\
('kB = {0:}, NB = {1:}, kC = {2:}, NC = {3:}\n'.format(kB, NB, kC, NC))
+\
('kBC = {0:}, NBC = {1:}'.format(kBC, NBC)), fontsize = 16)
plt.legend()
plt.tight_layout()
plt.savefig('FigFctrs.png')
plt.show()
MAPPB = CalcMAPEst(fit['PdB'])[0]
MAPPC = CalcMAPEst(fit['PdC'])[0]
MAPPBC = CalcMAPEst(fit['PdBC'])[0]
L_HDIPB, R_HDIPB = calcHDI(fit['PdB'])
L_HDIPC, R_HDIPC = calcHDI(fit['PdC'])
L_HDIPBC, R_HDIPBC =
calcHDI(fit['PdBC'])
sb.kdeplot(fit['PdB'], label = 'PdB,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPB, L_HDIPB, R_HDIPB))
sb.kdeplot(fit['PdC'], label = 'PdC,
mode = {0:.2f}, \n 95%
HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPC, L_HDIPC, R_HDIPC))
sb.kdeplot(fit['PdBC'], label = 'PdBC,
mode = {0:.2f}, \n
95% HDI = [{1:.2f}, {2:.2f}]'.
format(MAPPBC, L_HDIPBC, R_HDIPBC))
plt.xlabel('Prob. of detection',
fontsize = 14)
plt.yticks([])
plt.title('Posterior distribtion\n' +
\
('kB = {0:}, NB = {1:}, kC = {2:}, NC = {3:}\n'.format(kB, NB, kC, NC))
+\
('kBC = {0:}, NBC = {1:}'.format(kBC, NBC)), fontsize = 16)
plt.legend()
plt.tight_layout()
plt.savefig('FigPd.png')
plt.show()
リストRA1. 年齢Aを独立変数（スクリプトではAs）とするロジスティック回帰モデルのStanスクリプト（ファイル名reg_age.stan）
data {
int N;
int<lower = 0, upper =
1> Rs[N];
real As[N];
}
parameters {
real a;
real b;
}
transformed parameters {
real <lower = 0.0, upper
= 1.0> Pd[N];
real<lower = 0.0, upper =
1.0> Pc[N]; // Probability of correct response
for (i in 1:N) {
Pd[i] = inv_logit(a + b * As[i]);
// Logistic regression model
Pc[i] = Pd[i] + (1.0 - Pd[i]) * (1.0/3.0);
}
}
model {
a ~ normal(0.0, 10.0); // weakly informative priot
b ~ normal(0.0, 10.0); // weakly informative priot
Rs ~ bernoulli(Pc);
}
リストRA2．リストRA1のStanスクリプトを用いるPythonスクリプト（ファイル名TriLogRegAge.py）
import csv
from copy import deepcopy
import numpy as np
import scipy.stats as ss
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
# 日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
fin_nm = input('Input file (*.csv) =
')
with open(fin_nm, 'r') as f:
Data_in = [v for v in
csv.reader(f)]
Data_in = np.array(Data_in)
Rs = [int(v) for v in
Data_in.T[1][1:]]
As = [float(v) for v in Data_in.T[2][1:]]
print(len(Rs), len(As))
N = len(Rs)
Rs = np.array(Rs)
As = np.array(As)
Raw_As = deepcopy(As)
As_mean = As.mean()
As_sd = As.std()
As = (As - As_mean) / As_sd
# Stanスクリプトのコンパイル
sm = pystan.StanModel(file =
'reg_age.stan')
# MCMCサンプリング
fit = sm.sampling(data = {'N':N, 'Rs':
Rs, 'As': As},
pars = ['a', 'b'], iter = 10_000, n_jobs = 1)
print(fit)
#
# modeの推定
#
def est_mode(samples, a = 0.05,
n_points = 10_000):
import numpy as np
import scipy.stats as ss
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
a_mode = est_mode(fit['a'])
b_mode = est_mode(fit['b'])
# パラメータaの事後分布のグラフ
sb.kdeplot(fit['a'], label = f'mode =
{a_mode:.2f}')
plt.xlabel('a', fontsize = 14)
plt.yticks([])
plt.title('aの事後分布', fontsize = 16)
plt.legend(fontsize = 14)
plt.savefig('Fig_a.png')
plt.show()
# パラメータbの事後分布のグラフ
sb.kdeplot(fit['b'], label = f'mode =
{b_mode:.2f}')
plt.xlabel('b', fontsize = 14)
plt.yticks([])
plt.title('bの事後分布', fontsize = 16)
plt.legend(fontsize = 14)
plt.savefig('Fig_b.png')
plt.show()
#
# 値のカテゴリ化
# 0:[-inf, 29), 1:[30, 39),
2:[40, 49), 3:[50, +inf)
#
def check_cat(v):
v = int((v - 20) / 10)
if v > 3:
v =
3
return v
Ones = np.zeros(4) # 各カテゴリにおける判断１の数
Boths = np.zeros(4) # 各カテゴリにおける判断１および０の数
for r, a in zip(Rs, Raw_As):
icat = check_cat(a)
Boths[icat] += 1
if r == 1:
Ones[icat] += 1
print(Ones)
print(Boths)
# Ｘ軸上のグリッド点100個
x_coord = np.linspace(As.min(),
As.max(), 100)
def Pd(u):
return 1.0 / (1.0 +
np.exp(-u))
def Pc(x):
return Pd(x) + (1.0 - Pd(x))
/ 3
y_Pc_m = Pc(a_mode + b_mode*x_coord)
y_Pd_m = Pd(a_mode + b_mode*x_coord)
plt.plot(x_coord, y_Pc_m, label = 'Pc')
plt.plot(x_coord, y_Pd_m, label =
'Pd')
#
# 各カテゴリ区間におけるデータの正答率
#
data_x = [(v - As_mean)/As_sd for v in
[25, 35, 45, 55] ]
data_y = [Ones[i]/Boths[i] for i in
range(4)]
plt.plot(data_x, data_y, 'o', label =
'Data')
x_vs = [(v - As_mean)/As_sd for v in
[20, 30, 40, 50, 60]]
x_lbls = [f'{v}' for v in [20, 30, 40,
50, 60]]
plt.xticks(x_vs, x_lbls)
plt.xlabel('独立変数', fontsize = 14)
plt.legend(fontsize = 14)
plt.ylabel('確率', fontsize = 14)
plt.title('３点識別試験', fontsize = 16)
plt.savefig('FigRegAge.png')
plt.show()