第12章 感覚と判断
信号検出理論
信号検出理論に基づくデータ分析プログラム例を、Yes/Noの強制2選択肢課題の等分散モデル、多件法評定課題の等分散を仮定しないモデルと等分散モデルの場合について用意してみた。なお、本ウェブサイトにおける式番号あるいはページ数は<岡本安晴「心理学データ分析と測定」勁草書房>におけるものである。ベイズ分析とそのソフトStanについては<岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版>で説明している。
ノイズのみの刺激に対して、「ノイズのみ(No)」と判断された回数が19回、「シグナルを含む(Yes)」と答えた回数が8回、
シグナルを含む刺激に対して、「ノイズのみ(No)」と答えた回数が3回、「シグナルを含む(Yes)」と答えた回数が24回
であったとする。このデータを図1.1のようにExcelで用意する。
図1.1 ファイル名Data2Cat.xlsx
1行目に回答の種類を表すカテゴリ名を「ノイズのみ」、「シグナルを含む」の順に設定している。「シグナルを含むかどうか」という問いに対する答えとして「Yes/No」である。第1列目に提示刺激の種類を、「ノイズのみ」を2行目に、「シグナルを含む」を3行目に設定している。プログラムは、この図1.1の形式のデータを読み込む。
このとき、次式が成り立つ(岡本、2014、pp. 179-180参照)。
これらの式に基づいて図1.1のデータ分析用のStanスクリプトを仮定
の下で、リスト1(本サイトの後半)のように用意した。
リスト1のStanスクリプトによりベイズ分析を行うPythonスクリプトをリスト2(本サイトの後半)に示す。なお、ファイルは圧縮ファイルfiles12.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト2のプログラムを実行すると入力データファイル名を聞いてくる(図1.2)。
図1.2
入力データファイル名を設定すると、計算が始まるが、Stanスクリプトのコンパイルに少し時間が掛かる。StanスクリプトのコンパイルとMCMCによるサンプリングが終わると、判断の基準値と弁別力の事後分布のグラフが表示される(図1.3)。
図1.3
図1.3のウィンドウを閉じると、ROC曲線のグラフが表示される(図1.4)。
図1.4
弁別力の事後分布における95%Ciと50%CIに対応するROC曲線の領域がそれぞれ黄色と緑の領域で表されている。青色の太い曲線で描かれたROC曲線は、弁別力の事後分布の中央値を点推定値としたときのグラフである。
図1.4のウィンドウを閉じるとプログラムの実行終了となり、端末には以下の出力が表示されている。
C(mean) = 0.549
95%CI = [0.0632, 1.06]
median = 0.544
Q1 = 0.386, Q3 = 0.712
d'(mean) = 1.8
95%CI = [1.05, 2.65]
median = 1.79
Q1 = 1.53, Q3 = 2.06
ノイズのみの刺激あるいはシグナルを含む刺激を提示されたときに、多件法による評定、例えば4件法であれば、「ノイズのみ」、「多分ノイズのみ」、「多分シグナルを含む」、「シグナルを含む」で答えるときのデータは、例えば図2.1のようになる。
図2.1 ファイル名Data4Cat.xlsx
図2.1のデータの形式は、1行目に評定回答の種類、左が「シグナルを含まない」という評定が最も強い回答であり、右に行くほど「シグナルを含む」という評定が強くなっている。第1列目には提示刺激の種類が書かれており、2行目に「ノイズのみの刺激」に対する各評定の度数が書かれ、3行目に「シグナルを含む刺激」に対する評定回答の度数が書かれている。
多件法評定課題の4件法の場合のモデルは、以下のように表される(岡本、2014、pp. 176-177参照)。
上のモデルに対するStanスクリプトを本サイト後半のリスト3に示す。ただし、評定カテゴリ数をK(≥3)として、一般化してある。リスト3のStanスクリプトを用いるPythonスクリプトを本サイトの後半にリスト4として示した。これらのファイルは圧縮ファイルfiles12.zipとしてまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト4のスクリプトを実行すると、まず入力データファイル名が聞かれる(図2.2)。
図2.2
入力データファイル名を設定すると、データの読み込みが行われて計算が始まるが、Stanスクリプトのコンパイルに多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリングが終了するとパラメータとの事後分布のグラフが表示される(図2.3)。
図2.3
図2.3のウィンドウを閉じると、の事後分布のグラフが表示される(図2.4)。
図2.4
図2.4のウィンドウを閉じると、ROC曲線のグラフが表示される(図2.5)。
図2.5
ROC曲線は、パラメータおよびの事後分布の中央値を点推定値として描いたものが太い青の曲線で描かれている。さらに、パラメータの組とをMCMCサンプリングのものから100組ランダムに選び、それぞれの組の値に対して描いたROC曲線が薄い黄色の曲線で描かれている。これらの100組の曲線の分布から、事後分布におけるROC曲線の散らばり具合(分布)を読み取ることができる。
図2.5のウィンドウを閉じるとプログラムは実行終了となり、端末画面には以下の出力が表示されている。がd’で、がsgmで表記されている。
d'
mean = 1.43
95%CI = [0.681,
2.4]
Q1 = 1.13, Med. = 1.39, Q3 = 1.68
sgm
mean = 1.14
95%CI = [0.622,
1.97]
Q1 = 0.896, Med. = 1.08, Q3 = 1.33
C1
mean = -0.341
95%CI = [-0.795,
0.124]
Q1 =
-0.499, Med. = -0.343, Q3 = -0.179
C2
mean = 0.364
95%CI = [-0.105,
0.831]
Q1 = 0.208, Med. = 0.363, Q3 = 0.521
C3
mean = 2.26
95%CI = [1.38,
3.36]
Q1 = 1.92, Med. = 2.23, Q3 = 2.57
図2.4のの事後分布、および上の端末画面に表示された、すなわちsgmの平均値、中央値、および50%CI=[Q1, Q3]から、と考えてもよさそうである。図2.1のデータを、次の等分散モデルで分析してみる。
等分散モデルは、
と置いたものである。Stanスクリプトを本サイトの後半にリスト5として示す。リスト5を用いる分析のPythonスクリプトを本サイト後半のリスト6に示す。これらのファイルは圧縮ファイルfiles12.zipにまとめたので、ダウンロード解凍すれば自由に利用できる。
リスト6のプログラムを実行すると、まず入力データファイル名の設定が求められる(図3.1)。
図3.1
入力データファイル名を設定すると、データが読み込まれ、計算が始まるが、Stanスクリプトのコンパイルに少し時間が掛かる。Stanによるサンプリングが終わると、との事後分布のグラフが表示される(図3.2)。
図3.2
図3.2のウィンドウを閉じると、ROC曲線のグラフが表示される(図3.3)。
図3.3
の95%CIおよび50%CIに対応するROC曲線の領域が、薄い黄色と緑で示されている。太い青の曲線で描かれたROC曲線は、の事後分布の中央値に対するものである。
図3.3のウィンドウを閉じると、プログラムの実行終了となり、端末画面には以下の出力が表示されている。
d'
mean = 1.42
95%CI = [0.694,
2.34]
Q1 = 1.12, Med. = 1.39, Q3 = 1.67
C1
mean = -0.327
95%CI = [-0.803,
0.139]
Q1 =
-0.475, Med. = -0.323, Q3 = -0.171
C2
mean = 0.372
95%CI =
[-0.0635, 0.843]
Q1 = 0.211, Med. = 0.37, Q3 = 0.525
C3
mean = 2.25
95%CI = [1.37,
3.28]
Q1 = 1.9, Med. = 2.21, Q3 = 2.55
岡本安晴(2019)いまさら聞けないPythonでデータ分析。丸善出版
リスト1 強制2選択肢課題(等分散モデル)のStanスクリプト:ファイル名SDT2Cat.stan
data {
int N_no;
int N_yes;
int S_no;
int S_yes;
}
parameters {
real d;
real C;
}
transformed parameters {
real<lower = 0.0, upper =
1.0> N_p;
real<lower = 0.0, upper =
1.0> S_p;
N_p = normal_cdf(C, 0.0,
1.0);
S_p = normal_cdf(C - d, 0.0,
1.0);
}
model {
d ~ normal(0.0, 1000.0);
C ~ normal(0.0, 1000.0);
N_no ~ binomial(N_no +
N_yes, N_p);
S_no ~ binomial(S_no +
S_yes, S_p);
}
リスト2 強制2選択肢課題のPythonスクリプト:ファイル名SDT2Cat.py
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
N_no = data.values[0][1]
N_yes = data.values[0][2]
S_no = data.values[1][1]
S_yes = data.values[1][2]
print('Noise no:', N_no, ' yes:', N_yes)
print('Signal no:', S_no, ' yes:', S_yes)
Data = {'N_no': N_no, 'N_yes': N_yes,
'S_no': S_no, 'S_yes': S_yes}
sm = pystan.StanModel(file =
'SDT2Cat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
C_mean = fit['C'].mean()
C_p025, C_Q1, C_med, C_Q3, C_p975 =
np.percentile(fit['C'], [2.5, 25, 50, 75, 97.5])
d_mean = fit['d'].mean()
d_p025, d_Q1, d_med, d_Q3, d_p975 =
np.percentile(fit['d'], [2.5, 25, 50, 75, 97.5])
sb.kdeplot(fit['C'], label = 'C')
sb.kdeplot(fit['d'], label =
"d'")
plt.yticks([])
plt.title("C(med) = {0:.3} d'(med) =
{1:.3}".format(C_med, d_med), fontsize = 16)
plt.legend(fontsize = 14)
plt.show()
Cs = np.arange(-5, 5, 0.01)
FAs = [1.0 - ss.norm.cdf(v, 0.0, 1.0)
for v in Cs]
Hits = [1.0 - ss.norm.cdf(v, d_med,
1.0) for v in Cs]
Hitsp025 = [1.0 - ss.norm.cdf(v,
d_p025, 1.0) for v in Cs]
HitsQ1 = [1.0 - ss.norm.cdf(v, d_Q1,
1.0) for v in Cs]
HitsQ3 = [1.0 - ss.norm.cdf(v, d_Q3,
1.0) for v in Cs]
Hitsp975 = [1.0 - ss.norm.cdf(v,
d_p975, 1.0) for v in Cs]
plt.figure(figsize = (5, 5))
plt.plot([0, 1, 1, 0, 0], [0, 0, 1, 1,
0], c = 'k')
plt.plot([0, 1], [0, 1], c = 'k')
plt.fill_between(FAs, Hitsp025,
Hitsp975, color = 'y', alpha = 0.3)
plt.fill_between(FAs, HitsQ1, HitsQ3,
color = 'g', alpha = 0.3)
plt.plot(FAs, Hits, c = 'b')
plt.plot(N_yes/(N_yes + N_no),
S_yes/(S_yes + S_no), 'ro', label = 'Data')
plt.title(("ROC curve: d'(med) =
{0:.3}".format(d_med)) +
'\n95%CI(yellow)
50%CI(green)', fontsize = 16)
plt.xlabel('FA', fontsize = 14)
plt.ylabel('Hit', fontsize = 14)
plt.legend(fontsize = 14)
plt.show()
print('\nC(mean) = {0:.3}'.format(C_mean))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(C_p025, C_p975))
print(' median =
{0:.3}'.format(C_med))
print(' Q1 = {0:.3}, Q3 = {1:.3}'.format(C_Q1,
C_Q3))
print("\nd'(mean) =
{0:.3}".format(d_mean))
print(' 95%CI = [{0:.3}, {1:.3}]'.format(d_p025,
d_p975))
print(' median =
{0:.3}'.format(d_med))
print(' Q1 = {0:.3}, Q3 = {1:.3}'.format(d_Q1,
d_Q3))
print()
リスト3 評定カテゴリ数Kが3以上の場合のStanスクリプト(ファイル名SDTKCat.stan)。
data {
int K;
int<lower = 0>
Noise[K];
int<lower = 0>
Signal[K];
vector[K] a;
}
parameters {
simplex[K] theta;
real d;
real<lower = 0.0> sgm;
}
transformed parameters {
simplex[K] p_noise;
simplex[K] p_signal;
ordered[K-1] C;
real v;
v = theta[1];
for (i in 1:(K-1)) {
C[i]
= logit(v);
v +=
theta[i+1];
}
p_noise[1] =
normal_cdf(C[1], 0.0, 1.0);
p_signal[1] =
normal_cdf(C[1], d, sgm);
for (k in 2:(K-1)) {
p_noise[k] = normal_cdf(C[k], 0.0, 1.0) - normal_cdf(C[k-1], 0.0, 1.0);
p_signal[k] = normal_cdf(C[k], d, sgm) - normal_cdf(C[k-1], d, sgm);
}
p_noise[K] = 1.0 -
normal_cdf(C[K-1], 0.0, 1.0);
p_signal[K] = 1.0 -
normal_cdf(C[K-1], d, sgm);
}
model {
d ~ normal(0.0, 1000.0);
sgm ~ exponential(0.001);
theta ~ dirichlet(a);
Noise ~
multinomial(p_noise);
Signal ~
multinomial(p_signal);
}
リスト4 リスト3のStanスクリプトによるベイズ分析のPythonスクリプト(ファイル名SDTKCat.py)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.optimize as so
import pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
X = data.values
K = len(X[0]) - 1
print('K = ', K)
Res_noise = [v for v in X[0][1:]]
Res_signal = [v for v in X[1][1:]]
print('Res_noise =\n', Res_noise)
print('Res_signal =\n', Res_signal)
print('a = \n', np.ones(K))
Data = {'K': K, 'Noise': Res_noise,
'Signal': Res_signal, 'a': np.ones(K)}
sm = pystan.StanModel(file =
'SDTKCat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
sb.kdeplot(fit['d'], color = 'b',
label = '$\mu_s$')
for k in range(K-1):
sb.kdeplot(fit['C'].T[k],
linestyle = '--', label = 'C{}'.format(k+1))
plt.yticks([])
plt.title('Posterior distributions\n'
+
'$\mu_s$ and Cs', fontsize = 16)
plt.legend(fontsize = 12)
plt.show()
sb.kdeplot(fit['sgm'])
plt.title('Posterior distribution of
$\sigma_s$', fontsize = 16)
plt.show()
d_mean = fit['d'].mean()
d_p025, d_Q1, d_med, d_Q3, d_p975 =
np.percentile(fit['d'], [2.5, 25, 50, 75, 97.5])
C_means = []
for k in range(K - 1):
C_means.append(fit['C'][:,
k].mean())
C_pest = []
for k in range(K - 1):
C_pest.append(np.percentile(fit['C'][:,
k], [2.5, 25, 50, 75, 97.5]))
sgm_mean = fit['sgm'].mean()
sgm_p025, sgm_Q1, sgm_med, sgm_Q3,
sgm_p975 = np.percentile(fit['sgm'], [2.5, 25, 50, 75, 97.5])
n = 100
Param_sgm = fit['sgm']
Param_d = fit['d']
smpl_pos =
np.random.choice(np.arange(len(Param_sgm)), size = n, replace = False)
plt.figure(figsize = (5, 5))
plt.plot([0, 1, 1, 0, 0], [0, 0, 1, 1,
0], c = 'k')
plt.plot([0, 1], [0, 1], c = 'k')
Cs = np.arange(-5, 5, 0.01)
for i in range(n):
if i % 10 == 0:
if i
== 0:
print('step-{0}/{1}'.format(i, n), end = '')
else:
print('\rstep-{0}/{1}'.format(i, n), end = '')
FAs = [1.0 - ss.norm.cdf(v,
0.0, 1.0) for v in Cs]
Hits = [1.0 - ss.norm.cdf(v,
Param_d[smpl_pos[i]], Param_sgm[smpl_pos[i]]) for v in Cs]
plt.plot(FAs, Hits, color =
'y', alpha = 0.2)
print()
FAs = [1.0 - ss.norm.cdf(v, 0.0, 1.0)
for v in Cs]
Hits = [1.0 - ss.norm.cdf(v, d_med,
sgm_med) for v in Cs]
plt.plot(FAs, Hits, color = 'b',
linewidth = 3, solid_capstyle = 'round')
CFA = [1.0 - ss.norm.cdf(C_pest[k][2],
0.0, 1.0) for k in range(K-1)]
CHit = [1.0 -
ss.norm.cdf(C_pest[k][2], d_med, sgm_med) for k in range(K-1)]
plt.plot(CFA, CHit, 'gs', label =
'Model')
Res_noise = np.array(Res_noise)
Res_signal = np.array(Res_signal)
sum_noise = Res_noise.sum()
sum_signal = Res_signal.sum()
cum_noise = 0
cum_signal = 0
for k in range(K - 1):
cum_noise +=
Res_noise[K-1-k]
cum_signal +=
Res_signal[K-1-k]
if k == 0:
plt.plot(cum_noise/sum_noise, cum_signal/sum_signal, 'ro', label =
'Data')
else:
plt.plot(cum_noise/sum_noise, cum_signal/sum_signal, 'ro')
plt.title('ROC curves\n' +
'$\mu_s$(med) = {0:.3}
$\sigma_s$(med) = {1:.3}'.format(d_med, sgm_med),
fontsize = 16)
plt.legend(fontsize = 14)
plt.show()
print("\nd'")
print(' mean =
{0:.3}'.format(d_mean))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(d_p025, d_p975))
print(' Q1 =
{0:.3}, Med. = {1:.3}, Q3 = {2:.3}'.format(d_Q1, d_med,
d_Q3))
print('\nsgm')
print(' mean =
{0:.3}'.format(sgm_mean))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(sgm_p025, sgm_p975))
print(' Q1 =
{0:.3}, Med. = {1:.3}, Q3 = {2:.3}'.format(sgm_Q1,
sgm_med, sgm_Q3))
for k in range(K-1):
print()
print('C{}'.format(k+1))
print(' mean =
{0:.3}'.format(C_means[k]))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(C_pest[k][0], C_pest[k][4]))
print(' Q1 =
{0:.3}, Med. = {1:.3}, Q3 = {2:.3}'.format(
C_pest[k][1], C_pest[k][2], C_pest[k][3]))
print()
リスト5 多件法評価課題のStanスクリプト(ファイル名SDTKCatSgmOne.stan)
data {
int K;
int<lower = 0>
Noise[K];
int<lower = 0>
Signal[K];
vector[K] a;
}
parameters {
simplex[K] theta;
real d;
}
transformed parameters {
simplex[K] p_noise;
simplex[K] p_signal;
ordered[K-1] C;
real v;
v = theta[1];
for (i in 1:(K-1)) {
C[i]
= logit(v);
v +=
theta[i+1];
}
p_noise[1] = normal_cdf(C[1],
0.0, 1.0);
p_signal[1] =
normal_cdf(C[1], d, 1.0);
for (k in 2:(K-1)) {
p_noise[k] = normal_cdf(C[k], 0.0, 1.0) - normal_cdf(C[k-1], 0.0, 1.0);
p_signal[k] = normal_cdf(C[k], d, 1.0) - normal_cdf(C[k-1], d, 1.0);
}
p_noise[K] = 1.0 -
normal_cdf(C[K-1], 0.0, 1.0);
p_signal[K] = 1.0 -
normal_cdf(C[K-1], d, 1.0);
}
model {
d ~ normal(0.0, 1000.0);
theta ~ dirichlet(a);
Noise ~
multinomial(p_noise);
Signal ~
multinomial(p_signal);
}
リスト6 多件法評価課題の分析スクリプト(ファイル名SDTKCatSgmOne.py)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.optimize as so
import pystan
import seaborn as sb
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
X = data.values
K = len(X[0]) - 1
print('K = ', K)
Res_noise = [v for v in X[0][1:]]
Res_signal = [v for v in X[1][1:]]
print('Res_noise =\n', Res_noise)
print('Res_signal =\n', Res_signal)
print('a = ', np.ones(K))
Data = {'K': K, 'Noise': Res_noise,
'Signal': Res_signal, 'a': np.ones(K)}
sm = pystan.StanModel(file =
'SDTKCat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
sb.kdeplot(fit['d'], color = 'b',
label = '$\mu_s$')
for k in range(K-1):
sb.kdeplot(fit['C'].T[k],
linestyle = '--', label = 'C{}'.format(k+1))
plt.yticks([])
plt.title('Posterior distributions\n'
+
'$\mu_s$ and Cs', fontsize = 16)
plt.legend(fontsize = 12)
plt.show()
d_mean = fit['d'].mean()
d_p025, d_Q1, d_med, d_Q3, d_p975 =
np.percentile(fit['d'], [2.5, 25, 50, 75, 97.5])
C_means = []
for k in range(K - 1):
C_means.append(fit['C'][:,
k].mean())
C_pest = []
for k in range(K - 1):
C_pest.append(np.percentile(fit['C'][:,
k], [2.5, 25, 50, 75, 97.5]))
Param_d = fit['d']
plt.figure(figsize = (5, 5))
plt.plot([0, 1, 1, 0, 0], [0, 0, 1, 1,
0], c = 'k')
plt.plot([0, 1], [0, 1], c = 'k')
Cs = np.arange(-5, 5, 0.01)
FAs = [1.0 - ss.norm.cdf(v, 0.0, 1.0)
for v in Cs]
Hits = [1.0 - ss.norm.cdf(v, d_med,
1.0) for v in Cs]
Hitsp025 = [1.0 - ss.norm.cdf(v,
d_p025, 1.0) for v in Cs]
HitsQ1 = [1.0 - ss.norm.cdf(v, d_Q1,
1.0) for v in Cs]
Hitsmed = [1.0 - ss.norm.cdf(v, d_med,
1.0) for v in Cs]
HitsQ3 = [1.0 - ss.norm.cdf(v, d_Q3,
1.0) for v in Cs]
Hitsp975 = [1.0 - ss.norm.cdf(v,
d_p975, 1.0) for v in Cs]
plt.fill_between(FAs, Hitsp025,
Hitsp975, color = 'y', alpha = 0.3)
plt.fill_between(FAs, HitsQ1, HitsQ3,
color = 'g', alpha = 0.3)
plt.plot(FAs, Hits, color = 'b',
linewidth = 3, solid_capstyle = 'round')
CFA = [1.0 - ss.norm.cdf(C_pest[k][2],
0.0, 1.0) for k in range(K-1)]
CHit = [1.0 -
ss.norm.cdf(C_pest[k][2], d_med, 1.0) for k in range(K-1)]
plt.plot(CFA, CHit, 'gs', label =
'Model')
Res_noise = np.array(Res_noise)
Res_signal = np.array(Res_signal)
sum_noise = Res_noise.sum()
sum_signal = Res_signal.sum()
cum_noise = 0
cum_signal = 0
for k in range(K - 1):
cum_noise +=
Res_noise[K-1-k]
cum_signal +=
Res_signal[K-1-k]
if k == 0:
plt.plot(cum_noise/sum_noise, cum_signal/sum_signal, 'ro', label =
'Data')
else:
plt.plot(cum_noise/sum_noise, cum_signal/sum_signal, 'ro')
plt.title('ROC curves', fontsize = 16)
plt.xlabel('FA', fontsize = 14)
plt.ylabel('Hit', fontsize = 14)
plt.legend(fontsize = 14)
plt.show()
print(fit)
print("\nd'")
print(' mean =
{0:.3}'.format(d_mean))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(d_p025, d_p975))
print(' Q1 =
{0:.3}, Med. = {1:.3}, Q3 = {2:.3}'.format(d_Q1, d_med,
d_Q3))
for k in range(K-1):
print()
print('C{}'.format(k+1))
print(' mean =
{0:.3}'.format(C_means[k]))
print(' 95%CI = [{0:.3},
{1:.3}]'.format(C_pest[k][0], C_pest[k][4]))
print(' Q1 =
{0:.3}, Med. = {1:.3}, Q3 = {2:.3}'.format(
C_pest[k][1], C_pest[k][2], C_pest[k][3]))
print()