拡張２肢強制選択評定法データの分析
PyStan3対応スクリプト
２AFC課題については、このウェブサイトで説明している。
どちらの条件の方が記憶が強いかという基本的な科学的記憶研究においては２AFCが用いられるべきであるとBradyら（2023）は言っている（p. 422）。しかし、２AFC課題では等分散モデルが仮定されていると説明されている（p. 442）が、ROC曲線は不等分散モデルがよく合うと言っている（p. 442）。拡張2AFC評定法（https://osf.io/ts2eq/）を用いると不等分散モデルによる分析が可能であるので、Bradyら（2023）の主張に従えば、記憶研究においては拡張２肢強制選択評価法が用いられるべきであると言える。不確かな仮定（等分散仮定）に基づいた推定値
は、等分散仮定が成り立たない度合いに応じての推定値も誤ったものになる（下図１参照）。等分散仮定が不確かなときは、等分散仮定を必要としない拡張２肢強制選択評価法を用いれば、間違った等分散仮説の影響を防ぐことができる。
２AFC課題では、シグナル刺激感覚の平均値と分散を個別に推定することはできない。このことをグラフ（図１）で見てみる。
いま、ノイズ刺激感覚
とシグナル刺激感覚
が次の正規分布に従うとする。
このとき
となるので、正答確率は
で与えられる。
正答確率
ののグラフは、図１のようである。
図１
図１のグラフはリスト1のスクリプトで描いている。
シグナル感覚の平均値が大きくなると、正答確率も上がるが、標準偏差
が大きくなると正答確率は下がる。
同じ等高線上のとの値の組は同じ正答確率Pcを与えるので、正答確率Pcからとの値の組を決めることはできない。との値の組は、２AFC課題を以下のように拡張すると求めることができる。
２AFCでは、提示刺激対は＜ノイズ、シグナル＞と＜シグナル、ノイズ＞の2種類である。拡張２肢強制選択評定課題では、さらに＜ノイズ、ノイズ＞と＜シグナル、シグナル＞の２種類を加えて、４種類の刺激対に対して評定法を用いるので、シグナル刺激感覚の平均値と分散の両方の推定を行うことができる。
２AFC課題では、各試行において２つの刺激（刺激対）が与えられる。各刺激の提示位置（空間的、時間的）を＜S1, S2＞で表す。第1提示位置がS1、第２提示位置がS2である。＜n, s＞（あるいは、n-s）は、第１位置にノイズ刺激、第２位置にシグナル刺激が提示されたことを表す。
第１提示位置の刺激感覚を
、第２刺激提示位置の刺激感覚を
で表す。
第１刺激位置に対する第２刺激位置のバイアスをbで表す。このバイアスは恒常誤差のように感覚のものと反応のものの両者を含むものである。
このとき、以下のようになる。
のとき、
のとき、
のとき、
のとき、
いま、判断がK個のカテゴリによる評定で回答されるものとし、判断のカテゴリ境界を以下のようにおく。
このとき、
であれば、提示刺激対<
>に対する判断
は
が与えられるとする。すなわち、
である。
例えば、判断カテゴリー数が
個のときは、評定判断は以下のようなものが考えられる。
： シグナル刺激は確かに第１位置である。
： シグナル刺激は第１位置である。
： シグナル刺激はたぶん第１刺激である。
： わからない／どちらも同じと思われる。
： シグナル刺激はたぶん第２位置である。
： シグナル刺激は第２位置である。
： シグナル刺激は確かに第２位置である。
平均
、分散
の正規分布の累積分布を
で表し、以下のように便宜的に設定する。
このとき、差
の分布は以下のように与えられる。
刺激対 
が提示されとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
刺激対
が提示されたとき、
したがって、
いま、カテゴリ境界は、原点を中心として対称に設定されているとする。
したがって、Kが偶数であれば、
奇数であれば、
である。
上のモデルに基づいてベイズ分析するためのStanスクリプトを、Kが偶数の場合をリスト２に、Kが奇数の場合をリスト3に示す。これらのスクリプトファイルを用いて拡張２AFC評定法によるデータをベイズ分析するスクリプトをリスト4とリスト５に示す。これらのファイルは圧縮ファイルscriptfiles.zipにまとめた。ダウンロード展開して自由に利用されたい。
なお、PyStan ２用のスクリプトはこのウェブサイトに用意している。
スクリプトの実行は、リスト４を行ってからリスト５を行う。リスト４では、MCMCサンプリングが実行され、得られたサンプルの分析をリスト５で行う。
ただし、Windowsの場合は、リスト４の実行は、WSL（Windows SubSystem for Linux）で行う。WSL上のUbuntuにインストールされたAnacondaにおいてPython 3.11の仮想環境を用意し、PyStan 3.7を利用した。PyStan 3はコマンド
pip install
pystan
でインストールできる。
古いPyStan 3が既にインストールされているときは、一旦アンインストールしてから、再度インストールすると新しいPyStan 3がインストールされる。
PyStan 3の再インストールは、コマンド
pip uninstall
pystan
の実行後、
pip install
pystan
を実行する。
Python 3.11仮想環境も、新しく作成した方が良い。
PyStan 3をプログラム内で使用するときは
import stan
である。
リスト４の実行後、Windowsの場合は、WSLを離れて、Windows上でリスト５を実行する。WSLでは、画像表示ができない。
Ubuntuの場合は、Anaconda上のPython3.11仮想環境で続けてリスト５を実行すればよい。
分析例を、まず、カテゴリ数が６個の偶数の場合について示し、続いてカテゴリ数が７個の奇数の場合を示す。
表1のようなデータが与えられたとする。
|
表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 |
まず、リスト４のファイルsampling.pyを実行する。コマンドで実行する場合は、以下のようである。
(py311) …/scriptfiles$ python sampling.py
プロンプト “n_n =” が表示されたら、カテゴリ１から６までの提示刺激対条件<noise, noise>に対する判断回答数を入力する。表１のデータであれば、次のように数値を半角空白文字で区切って並べて最後にEnterキーを押して入力する。
n_n = 2 5 15
21 6 1
Enterキーを押すと、次の入力プロンプトが表示されるので、同様に入力していく。
n_n = 2 5 15
21 6 1
n_s =
表１のデータをすべて入力すると、端末画面は以下のようになっている。
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サンプリングが始まる。
サンプリングが終了すると、端末画面は以下のようになっている。
・
・
・
ss_cum_p[7] 40000.0 NaN
C[0]
30166.0 1.0
C[1] 28251.0 1.0
C[2]
28233.0 1.0
C[3]
28233.0 1.0
C[4]
28251.0 1.0
C[5]
30166.0 1.0
vsum
30166.0 1.0
sgm_nn
40000.0 NaN
sgm_ns
27326.0 1.0
sgm_ss
27326.0 1.0
i_data.pkl was saved.
d_frame.pkl was saved.
K.pkl was saved.
End.
(py311) …/scriptfiles$
MCMCサンプリングの実行結果などが、ファイルi_data.pkl、d_frame.pkl、K.pklに保存されている。
上の出力において、NaNとあるのは、定数として設定された変数（
など）はRhatの計算において分母が０になることによる。
続いてリスト５のスクリプトを実行するが、Windowsの場合はWSLを抜け出て、Windows上で実行する。Ubuntuの場合は、そのままでよい。
(py311) …3\scriptfiles> python
analysis.py
上のコマンドを実行すると、リスト４の実行で作成されたファイルが読み込まれて、分析が行われる。
まず、の事後分布とMAP推定値が表示される（図２）。
図２
MAP推定値が良いことは、岡本（2022；https://doi.org/10.14947/psychono.41.1）に報告されている。
図２のフォームを閉じると、の事後分布とMAP推定値が表示される（図３）。
図３
図３のフォームを閉じると、カテゴリ境界値C1, C2, C4, C5の事後分布とMAP推定値のグラフが表示される（図４）。
図４
カテゴリ数が偶数であるので、C3は０に固定されているが、点（０，０）の赤い半円で表されている。
図４のフォームを閉じると、パラメータｂの事後分布とMAP推定値のグラフが表示される（図５）。
図５
図５のフォームを閉じると、図６のグラフが表示される。
図６
各提示刺激対に対して、評定カテゴリ判断の累積確率のデータ値（縦軸）とモデルによる予測値（横軸）のグラフが表示されている。破線の対角線は、データ値とモデル予測値が完全に一致した場合のグラフである。データ値と予測値の関係は、この対角線にほぼ沿っているので、モデルのデータに対する当てはまりは良いといえる。
図６のフォームを閉じると、プログラムの終了である。
グラフの図は、スクリプトファイルと同じフォルダにファイルとして保存されている（図７）。
図７
リスト4のスクリプトsampling.pyを次のように実行する。ただし、Windowsの場合は、WSL上のUbuntuにおいて実行する。
(py311) …/scriptfiles$ python
sampling.py
まず、提示刺激対＜ノイズ、ノイズ＞に対するデータの入力が求められる。
表２の場合のデータを以下のように入力する。
|
表２． カテゴリ数７個の場合のデータ |
||||
|
評定カテゴリ |
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
＜シグナル、シグナル＞に対する入力を終えると、ＭＣＭＣサンプリングが始まる。
サンプリングの終了後、端末画面は以下のようになる。
ss_cum_p[7] 40000.0 NaN
C[0]
29896.0 1.0
C[1]
27687.0 1.0
C[2]
28573.0 1.0
C[3]
28573.0 1.0
C[4]
27687.0 1.0
C[5]
29896.0 1.0
vsum
29896.0 1.0
sgm_nn 40000.0 NaN
sgm_ns
27567.0 1.0
sgm_ss
27567.0 1.0
i_data.pkl was saved.
d_frame.pkl was saved.
K.pkl was saved.
End.
(py311) …/scriptfiles$
MCMCサンプリングの実行結果などが、ファイルi_data.pkl、d_frame.pkl、K.pklに保存されている。
上の出力において、NaNとあるのは、定数に設定された変数（など）はRhatの計算において分母が０になることによる。
続いてリスト５のスクリプトを実行するが、Windowsの場合はWSLを抜け出て、Windows上で実行する。Ubuntuの場合は、そのままでよい。
(py311) …3\scriptfiles> python
analysis.py
上のコマンドを実行すると、リスト４の実行で作成されたファイルが読み込まれて、分析が行われる。
まず、の事後分布とMAP推定値が表示される（図８）。
図８
MAP推定値が良いという報告がある（岡本、2022；https://doi.org/10.14947/psychono.41.1）
図８のフォームを閉じると、の事後分布とMAP推定値が表示される（図９）。
図９
図９を閉じると、評定境界値C1, C2, C3, C4, C5, C6の事後分布とMAP推定値が表示される（図10）。
図10
図10のフォームを閉じると、パラメータｂの事後分布とＭＡＰ推定値が表示される（図11）。
図11
図11のフォームを閉じると、図12のグラフが表示される。
図12
各提示刺激対に対して、評定カテゴリ判断の累積確率のデータ値（縦軸）とモデルによる予測値（横軸）のグラフが表示されている。破線の対角線は、データ値とモデル予測値が完全に一致した場合のグラフである。データ値と予測値の関係は、この対角線にほぼ沿っているので、モデルのデータに対する当てはまりは良いといえる。
図Ex.10のフォームを閉じると、プログラムの終了である。
グラフの図は、スクリプトファイルと同じフォルダにファイルとして保存されている（図13）。
図13
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;
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);
}
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);
}
リスト４．拡張２肢強制選択評定法データ分析用Pythonスクリプト：MCMCサンプリング（sampling.py）
import numpy as np
import stan
import arviz as az
import pandas as pd
import pickle
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
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)
Data = {'K':K, 'nn_cond':n_n,
'ns_cond':n_s, 'sn_cond':s_n, 'ss_cond':s_s}
with open('SDT_2AFC_KCat_odd.stan' if
K % 2 == 1 else
'SDT_2AFC_KCat_even.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
fit = sm.sample(num_samples=10000)
print('fit\n', fit)
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
with open('i_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print('i_data.pkl was saved.')
df = fit.to_frame()
with open('d_frame.pkl', 'wb') as f:
pickle.dump(df, f)
print('d_frame.pkl was saved.')
with open('K.pkl', 'wb') as f:
pickle.dump((K, n_n, n_s,
s_n, s_s), f)
print('K.pkl was saved.')
print('\nEnd.\n')
リスト５．拡張２肢強制選択評定法データ分析用Pythonスクリプト：MCMCサンプルの分析（analysis.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pandas as pd
import pickle
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
def CalcMAPEst(samples, a = 0.05,
n_points = 2000):
"""
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]
with open('i_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
with open('d_frame.pkl', 'rb') as f:
fit = pickle.load(f)
with open('K.pkl', 'rb') as f:
(K, n_n, n_s, s_n, s_s) =
pickle.load(f)
smry = az.summary(inference_data, var_names
= ['mu_s', 'sgm_s'])
print('Summary =\n',smry)
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()
print('\n I am calculating.')
s_title = ''
if K % 2 == 1:
for k in range(K-1):
Csmpl
= fit[f'C.{k+1}']
sb.kdeplot(Csmpl)
c_map = CalcMAPEst(Csmpl)[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:
Csmpl = fit[f'C.{k+1}']
sb.kdeplot(Csmpl)
c_map = CalcMAPEst(Csmpl)[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.empty(K-1)
for k in range(K-1):
est_pcum_n_n[k] =
np.median(fit[f'nn_cum_p.{k+2}'])
est_pcum_n_s = np.empty(K-1)
for k in range(K-1):
est_pcum_n_s[k] =
np.median(fit[f'ns_cum_p.{k+2}'])
est_pcum_s_n = np.empty(K-1)
for k in range(K-1):
est_pcum_s_n[k] =
np.median(fit[f'sn_cum_p.{k+2}'])
est_pcum_s_s = np.empty(K-1)
for k in range(K-1):
est_pcum_s_s[k] =
np.median(fit[f'ss_cum_p.{k+2}'])
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()
print('End.')
データ分析用のPython解説書＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を用意した。