２選択肢強制法3対版
Two Alternative Forced Choice (2AFC)
Rating Method with Tri-Pair
シグナル刺激感覚の分散（標準偏差）の推定
Bradyら（2023）は、記憶のパフォーマンスの研究においては２AFCを用いることを推奨しているが、さらに、不等分散モデルがROC曲線によく合うと説明している（p. 442）。Decarlo(2010)は、不等分散モデルは、広く用いられていると説明している。
標準的な２AFC課題では等分散モデルが仮定されているが、2AFCを拡張して不等分散モデルが適用できるように工夫された実験法と分析モデルが提案されている（Okamoto, 2023）。不等分散が予想される場合は、この拡張2肢強制選択評定法を用いることができる。また、不等分散とは等分散を仮定しないということであり、データが等分散のものであっても不等分散モデルを適用することができる。この場合は、シグナル刺激の分散がノイズ刺激と同じ１であるという分析結果が予想される。
Okamoto(2023)では、提示刺激対は、標準2AFCの＜n,s＞と＜s,n＞に加えて、シグナル刺激の分散を推定するために＜n,n＞と＜s,s＞の２対が加えられている。しかし、＜s,s＞を除いた３対でもパラメータの推定が可能であるので、ここでは、＜n,n＞、＜n,s＞、＜s,n＞の３対を用いた場合について説明する。
まず、モデルについて説明を行い、その後、分析例を示す。ファイルは、tripfiles.zipにまとめた。
モデル
ノイズ刺激
の感覚を
で、シグナル刺激
の感覚を
で表し、それぞれ次の正規分布に従うとする。
ここで、感覚次元の原点と単位をノイズ刺激の感覚の平均値と標準偏差（分散の平方根）が次式を満たすように設定する（Wickens, 2002）。
検出力は、シグナル刺激感覚の平均値
とノイズ刺激感覚の平均値
の差で表される（Green & Swets, 1988）。
である。
シグナル刺激感覚の分散がノイズ刺激感覚の分散と等しいとき、すなわち、
のときは、検出力は
で表される。すなわち、
である。
と
の値は、２AFC課題を以下のように拡張すると求めることができる。
２AFCでは、提示刺激対は＜ノイズ、シグナル＞と＜シグナル、ノイズ＞の2種類である。これに、＜ノイズ、ノイズ＞を加えて、３種類の刺激対に対して評定法を用いると、シグナル刺激感覚の平均値と分散の推定を行うことができる。
２AFC課題では、各試行において２つの刺激（刺激対）が与えられる。各刺激の提示位置（空間的、時間的）を＜S1, S2＞で表す。第1提示位置がS1、第２提示位置がS2である。＜n, s＞（あるいは、n-s）は、第１位置にノイズ刺激、第２位置にシグナル刺激が提示されたことを表す。
第１提示位置の刺激感覚を
、第２刺激提示位置の刺激感覚を
で表す。
第１刺激位置に対する第２刺激位置のバイアスをbで表す。このバイアスは恒常誤差のように感覚のものと反応のものの両者を含むものである。
このとき、以下のようになる。
のとき、
のとき、
のとき、
のとき、
いま、判断がK個のカテゴリによる評定で回答されるものとし、判断のカテゴリ境界を以下のようにおく。
このとき、
であれば、提示刺激対<
>に対する判断
は
が与えられるとする。すなわち、
である。
例えば、判断カテゴリー数が
個のときは、評定判断は以下のようなものが考えられる。
： シグナル刺激は確かに第１位置である。
： シグナル刺激は第１位置である。
： シグナル刺激はたぶん第１刺激である。
： わからない／どちらも同じと思われる。
： シグナル刺激はたぶん第２位置である。
： シグナル刺激は第２位置である。
： シグナル刺激は確かに第２位置である。
平均
、分散
の正規分布の累積分布を
で表し、以下のように便宜的に設定する。
このとき、差
の分布は以下のように与えられる。
刺激対 
が提示されとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
いま、カテゴリ境界が、原点を中心として対称に設定されているとする。バイアスパラメータ
の使用により、カテゴリ境界は原点を中心として設定されていると考えることができる。
したがって、Kが偶数であれば、
奇数であれば、
上記のモデルのStanスクリプトを、
が偶数の場合をリスト１に、が奇数の場合をリスト２に示す。スクリプトファイルは、tripfiles.zipにまとめた。
これらのスクリプトを用いた分析例を次に示す。
リスト１のStanスクリプトファイルと、リスト１を使用するPythonスクリプトファイル（リスト3）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリを移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyのインストールについては、このウェブサイト、あるいはこのウェブサイトなどを参照されたい。
(stan) ****/evencatfiles$ python AnalTriP2AFCReven.py
上のコマンドを実行すると、入力データファイル名を聞いてくる。
(stan) ****/evencatfiles$
python AnalTriP2AFCReven.py
Input data file (*.xlsx) =
tripeven100.xlsx
上の例では、ファイル名tripeven100.xlsxが設定されている。
図A1に、ファイルtripeven100.xlsxをExcelで開いたものを示す。
図A1 Excelファイルtripeven100.xlsx
カテゴリ数がK=4の偶数の場合のデータ例である。
１行目に、変数名が置かれている。１列目は、判断カテゴリを表し、１から４までの判断カテゴリを表す整数値が昇順に置かれている。２列目は、刺激対＜ノイズ刺激、ノイズ刺激＞に対する頻度である。３列目に＜ノイズ刺激、シグナル刺激＞に対する頻度、４列目に＜シグナル刺激、ノイズ刺激＞に対する頻度が設定されている。判断カテゴリ数がＫのときは、図Ａ１は１＋Ｋ行のフォーマットになる。
入力データファイル名を設定すると、Stanスクリプトのビルドが始まる。MCMCサンプリングが終了すると、トレース図が表示される（図A2）。
図A2
右側に各サンプリングのトレースが表示され、左側にサンプリングの事後分布が表示されている。
図A2のWindowを閉じると、図A3のグラフが表示される。
図A3
パラメータすなわち
の事後分布である。タイトル欄に、第1四分位数、中央値（第２四分位数）、第３四分位数が表示されている。
図A3のWindowを閉じると、図A4のグラフが表示される。
図A4
シグナル刺激感覚の標準偏差
の事後分布のグラフである。タイトル欄に、第１四分位数、中央値（第２四分位数）、第３四分位数が表示されている。
図A4のWindowを閉じると、図A5のグラフが表示される。
図A5
判断のカテゴリ境界値の事後分布である。タイトル欄には、それぞれの中央値が表示されている。カテゴリ数が４の場合、カテゴリ境界はC1、C2、C3の３点であるが、制約（６）により、
なので定数である。C2が定数０であることを、図A5では赤色の半円で示している。
図A5のWindowを閉じると、図A6のグラフが表示される。
図A6
バイアスパラメータの事後分布である。第１四分位数
と第３四分位数
の間に０があるので、バイアス
であると考えてよいと思われる。
図A6のWindowを閉じると、図A7のグラフが表示される。
図A7
提示刺激対別の判断カテゴリーの比率データとモデルによる予測確率のグラフである。予測確率は、MCMCサンプリングによる事後分布の中央値である。データの方はランダムな変動が見られるが、モデルによる予測値は、事後分布の中央値を採用しているのでサンプリングに対して安定していると考えられる。図A7のグラフから、データに対してモデルは適合していると考えられる。
図A７のWindowを閉じると、図A8のグラフが表示される。
図A８
刺激対ごとのカテゴリ累積比率を、モデル予測値に対してデータの値のグラフを描いたものである。破線は、モデル予測値がデータと一致する場合のグラフである。モデル予測値は、事後分布の中央値をパラメータ値として求められている。モデル予測値のグラフは破線によく合っているので、モデルはデータに適合していると考えられる。
図A８のWindowを閉じると、スクリプトの実行終了である。端末にMCMCサンプルの統計量などが表示されているが、「NAN」などの値はパラメータに定数が設定されているものがあることによる。
リスト２のStanスクリプトファイルと、リスト２を使用するPythonスクリプトファイル（リスト４）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリを移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyのインストールについては、このウェブサイト、あるいはこのウェブサイトなどを参照されたい。
(stan) ****/oddcatfiles$ python AnalTriP2AFCRodd.py
上のコマンドを実行すると、入力データファイル名を聞いてくる。
(stan) ****/oddcatfiles$
python AnalTriP2AFCRodd.py
Input data file (*.xlsx) =
tripodd100.xlsx
上の例では、ファイル名tripodd100.xlsxが設定されている。
図B1に、ファイルtripodd100.xlsxをExcelで開いたものを示す。
図B1 Excelファイルtripodd100.xlsx
カテゴリ数がK=5の奇数の場合のデータ例である。
１行目に、変数名が置かれている。１列目は、判断カテゴリを表し、１から５までの判断カテゴリを表す整数値が昇順に置かれている。２列目は、刺激対＜ノイズ刺激、ノイズ刺激＞に対する頻度である。３列目に＜ノイズ刺激、シグナル刺激＞に対する頻度、４列目に＜シグナル刺激、ノイズ刺激＞に対する頻度が設定されている。判断カテゴリ数がＫのときは、図B１は１＋Ｋ行のフォーマットになる。
入力データファイル名を設定すると、Stanスクリプトのビルドが始まる。MCMCサンプリングが終了すると、トレース図が表示される（図B2）。
図B2
右側に各サンプリングのトレースが表示され、左側にサンプリングの事後分布が表示されている。
図B2のWindowを閉じると、図B3のグラフが表示される。
図B3
パラメータすなわちの事後分布である。タイトル欄に、第1四分位数、中央値（第２四分位数）、第３四分位数が表示されている。
図B3のWindowを閉じると、図B4のグラフが表示される。
図B4
シグナル刺激感覚の標準偏差の事後分布のグラフである。タイトル欄に、第１四分位数、中央値（第２四分位数）、第３四分位数が表示されている。
図B4のWindowを閉じると、図B5のグラフが表示される。
図B5
判断のカテゴリ境界値の事後分布である。タイトル欄には、それぞれの中央値が表示されている。
図B5のWindowを閉じると、図B6のグラフが表示される。
図B6
バイアスパラメータの事後分布である。第１四分位数
と第３四分位数
の間の区間の右端に０がある。負のバイアスが示唆されている。
図B6のWindowを閉じると、図B7のグラフが表示される。
図B7
提示刺激対別の判断カテゴリーの比率データとモデルによる予測確率のグラフである。予測確率は、MCMCサンプリングによる事後分布の中央値である。データの方はランダムな変動が見られるが、モデルによる予測値は、事後分布の中央値を採用しているのでサンプリングに対して安定していると考えられる。図B7のグラフから、データに対してモデルは適合していると考えられる。
図B７のWindowを閉じると、図B8のグラフが表示される。
図B８
刺激対ごとのカテゴリ累積比率を、モデル予測値に対してデータの値のグラフを描いたものである。破線は、モデル予測値がデータと一致する場合のグラフである。モデル予測値は、事後分布の中央値をパラメータ値として求められている。モデル予測値のグラフは破線によく合っているので、モデルはデータに適合していると考えられる。
図B８のWindowを閉じると、スクリプトの実行終了である。端末にMCMCサンプルの統計量などが表示されているが、「NAN」などの値はパラメータに定数が設定されているものがあることによる。
Brady, T. F., Robinson, M. M., Williams, J. R., & Wixted, J. T. (2023). Measuring memory is harder than you think: How to avoid problematic measurement practices in memory research. Psychonomic Bulletin & Review, 30, 421-449.
DeCarlo, L. T. (2010). On the statistical and theoretical basis of signal detection theory and extensions: Unequal variance, random coefficients, and mixture models. Journal of Mathematical Psychology, 54, 304-313.
Green, D. M. & Swets, J. A. (1988). Signal detection theory and psychophysics. Peninsula Publishing.
Hautus, M. J., Macmillan, N. A., & Creelman, C. D. (2022) Detection theory: A user’s guide. Rootledge.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
Okamoto, Y. (2023). Extended 2AFC Rating Task. https://doi.org/10.17605/OSF.IO/TS2EQ
岡本安晴（2025）感覚・知覚測定法．和氣典二・重野純・村上郁也（編）感覚・知覚心理学ハンドブック 第三版（第2章）、誠信書房
Wickens, T. D. (2002). Elementary Signal Detection Theory. Oxford University Press.
Wixted, J. T. (2019). The forgotten history of signal detection theory. Journal of Experimental Psychology: Learning, Memory, & Cognition, 46, 201-233.
附記
リスト１ カテゴリ数が偶数の場合（trip2AFC_KCat_even.stan）
//
//
Yasuharu Okamoto, 2026
//
functions {
real my_phi(real
x) {
if
(x <= -37.5) {
return 0.0;
}
else
if (x >= 8.25) {
return 1.0;
}
else
{
return std_normal_cdf(x);
}
}
}
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_cond;
}
transformed data {
int hK;
vector[K%/%2] alpha;
hK
= K %/% 2;
for (i
in 1:hK) {
alpha[i] = 1.0;
}
}
parameters {
real<lower=-10,
upper=10> mu_s;
real<lower = 0.01,
upper=10.0> sgm_s;
real b;
simplex[hK]
preC;
}
transformed parameters {
simplex[K] nn_theta;
simplex[K] ns_theta;
simplex[K] sn_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
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));
nn_cum_p[1]
= 0.0;
ns_cum_p[1]
= 0.0;
sn_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;
for (i
in 2:K) {
nn_cum_p[i] = my_phi((C[i-1]
- b) / sgm_nn);
ns_cum_p[i] = my_phi((C[i-1]
- (mu_s + b)) / sgm_ns);
sn_cum_p[i] = my_phi((C[i-1]
- (-mu_s + b)) / sgm_ns);
}
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];
}
}
model {
mu_s
~ uniform(-10, 10);
sgm_s
~ uniform(0.01, 10.0);
preC
~ dirichlet(alpha);
b ~ normal(0.0, 1000.0);
nn_cond
~ multinomial(nn_theta);
ns_cond
~ multinomial(ns_theta);
sn_cond
~ multinomial(sn_theta);
}
リスト２ カテゴリ数が奇数の場合（trip2AFC_KCat_odd.stan）
//
//
Yasuharu Okamoto, 2026
//
functions {
real my_phi(real
x) {
if
(x <= -37.5) {
return 0.0;
}
else
if (x >= 8.25) {
return 1.0;
}
else
{
return std_normal_cdf(x);
}
}
}
data {
int K;
array[K] int nn_cond;
array[K] int ns_cond;
array[K] int sn_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<lower=-10,
upper=10> mu_s;
real<lower = 0.01,
upper=10.0> sgm_s;
real b;
simplex[hK+1] preC;
}
transformed parameters {
simplex[K] nn_theta;
simplex[K] ns_theta;
simplex[K] sn_theta;
vector[K+1] nn_cum_p;
vector[K+1] ns_cum_p;
vector[K+1] sn_cum_p;
vector[K-1] C;
real vsum;
real sgm_nn;
real sgm_ns;
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));
nn_cum_p[1]
= 0.0;
ns_cum_p[1]
= 0.0;
sn_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;
for (i
in 2:K) {
nn_cum_p[i] = my_phi((C[i-1]
- b) / sgm_nn);
ns_cum_p[i] = my_phi((C[i-1]
- (mu_s + b)) / sgm_ns);
sn_cum_p[i] = my_phi((C[i-1]
- (-mu_s + b)) / sgm_ns);
}
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];
}
}
model {
mu_s
~ uniform(-10, 10);
sgm_s
~ uniform(0.01, 10.0);
preC
~ dirichlet(alpha);
b ~ normal(0.0, 1000.0);
nn_cond
~ multinomial(nn_theta);
ns_cond
~ multinomial(ns_theta);
sn_cond
~ multinomial(sn_theta);
}
リスト３ リスト１のStanスクリプトファイルを利用するPythonスクリプト（AnalTriP2AFCReven.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
import arviz
as az
inflnm = input('Input data file (*.xlsx) =
')
outflnm = 'Results.txt' # input('Output file (*.txt) = ')
fout = open(outflnm,
'w')
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]
if K % 2 == 1:
print('K should be even.')
import sys
sys.exit()
print('n_n
=', n_n)
print('n_s
=', n_s)
print('s_n
=', s_n)
model = CmdStanModel(stan_file = 'trip2AFC_KCat_even.stan')
Data = {'K':K, 'nn_cond':n_n, 'ns_cond':n_s, 'sn_cond':s_n}
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
inf_data = az.from_cmdstanpy(fit)
# Arviz(az) InferenceData
# トレースプロット
az.plot_trace(inf_data, var_names=['mu_s','sgm_s','C'])
plt.tight_layout()
plt.show()
fit = fit.draws_pd() # Pandas DataFrame
q1, mu_s_med,
q3 = np.percentile(fit['mu_s'],
[25, 50, 75])
sb.kdeplot(fit['mu_s'])
plt.xlabel(r'$\mu_s$',
fontsize=16)
plt.title(r'Posterior
Distribution of $\mu_s$' +
f'\nq1={q1:.3f},
Med.={mu_s_med:.3f},
q3={q3:.3f}',
fontsize=16)
plt.show()
q1, sgm_s_med,
q3 = np.percentile(fit['sgm_s'],
[25, 50, 75])
sb.kdeplot(fit['sgm_s'])
plt.title(r'Posterior
Distribution of $\sigma_s$' +
f'\nq1={q1:.3f},
Med.={sgm_s_med:.3f},
q3={q3:.3f}',
fontsize=16)
plt.xlabel(r'$\sigma_s$',
fontsize=16)
plt.show()
C = []
for k in range(K-1):
C.append(fit[f'C[{k+1}]'])
C = np.array(C).T
s_title = ''
for k in range(K-1):
if k != (K//2) - 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:
plt.plot([0], [0], marker = 'o', markersize
= 15, c = 'r')
s_title += f'C{k+1}=0, '
plt.title('Med. estimates\n' + s_title, fontsize=16)
plt.show()
q1, b_med,
q3 = np.percentile(fit['b'], [25, 50, 75])
sb.kdeplot(fit['b'])
plt.xlabel('b', fontsize=16)
plt.title('Posterior Distribution of b' +
f'\nq1={q1:.3f},
Med.={b_med:.3f},
q3={q3:.3f}',
fontsize=16)
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]
fit_nn_cum_p = []
fit_ns_cum_p = []
fit_sn_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_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
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]
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([0,1], [0,1], c = 'k', ls = '--',
label = 'Obs.=Est.')
plt.xlabel('Est.Cum.P',
fontsize=16)
plt.ylabel('Obs.Cum.P',
fontsize=16)
plt.legend()
plt.title('Cumulative Proportions', fontsize=18)
plt.show()
fout.close()
print('\n', outflnm,
'was saved.\n')
リスト４ リスト２のStanスクリプトファイルを利用するPythonスクリプト（AnalTriP2AFCRodd.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
import arviz
as az
inflnm = input('Input data file (*.xlsx) =
')
outflnm = 'Results.txt' # input('Output file (*.txt) = ')
fout = open(outflnm,
'w')
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]
if K % 2 == 0:
print('K should be odd.')
import sys
sys.exit()
print('n_n
=', n_n)
print('n_s
=', n_s)
print('s_n
=', s_n)
model = CmdStanModel(stan_file = 'trip2AFC_KCat_odd.stan')
Data = {'K':K, 'nn_cond':n_n, 'ns_cond':n_s, 'sn_cond':s_n}
fit = model.sample(data=Data)
print(fit.summary())
fout.write(f'\n{fit.summary()}\n')
inf_data = az.from_cmdstanpy(fit)
# Arviz(az) InferenceData
# トレースプロット
az.plot_trace(inf_data, var_names=['mu_s','sgm_s','C'])
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame
q1, mu_s_med,
q3 = np.percentile(fit['mu_s'],
[25, 50, 75])
sb.kdeplot(fit['mu_s'])
plt.xlabel(r'$\mu_s$')
plt.title(r'Posterior
Distribution of $\mu_s$' +
f'\nq1={q1:.3f},
Med.={mu_s_med:.3f},
q3={q3:.3f}',
fontsize=16)
plt.show()
q1, sgm_s_med,
q3 = np.percentile(fit['sgm_s'],
[25, 50, 75])
sb.kdeplot(fit['sgm_s'])
plt.title(r'Posterior
Distribution of $\sigma_s$' +
f'\nq1={q1:.3f},
Med.={sgm_s_med:.3f},
q3={q3:.3f}',
fontsize=16)
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 = ''
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 += ', '
plt.title('Med. estimates\n' + s_title)
plt.show()
q1, b_med,
q3 = np.percentile(fit['b'], [25, 50, 75])
sb.kdeplot(fit['b'])
plt.title('Posterior Distribution of b' +
f'\nq1={q1:.3f},
Med.={b_med:.3f},
q3={q3:.3f}',
fontsize=16)
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]
fit_nn_cum_p = []
fit_ns_cum_p = []
fit_sn_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_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
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]
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([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')