有意水準と期待効用
有意水準を、検定力との関係で決めるときは、統計的結論の効用は考慮されていない。結論の効用を考えた有意水準の検討を、期待効用の観点から試みる。
期待効用算出のために、以下の設定を行う。
第一種の誤りの確率、すなわち有意水準
帰無仮説が正しいときに、効果有りと判断することの効用
帰無仮説が正しいときに、効果無しと判断することの効用
: 対立仮説が正しい確率
第二種の誤りの確率、すなわち
: 検定力
対立仮説が正しいときに、効果有りと判断することの効用
対立仮説が正しいときに、効果無しと判断することの効用
簡単のために、検定統計量は正規分布に従うものとする（図１）。その後、F分布、ｔ分布、カイ二乗分布の場合についても扱う。
図１
このとき、期待効用は有意水準
の関数として、次式で表すことができる。
パラメータの幾つかのパターンについて、以下のように計算してみた。スクリプトは、リスト１に示すものである。スクリプトファイルなどは、まとめて圧縮ファイルoptalphafiles.zipとした。期待効用の合理的な観点からは、最適な有意水準としてパターン１では30.9％、パターン２では1.4％が得られた。
図２に、使用した入力データファイルを示す。
図２
入力データファイルは、CSV形式のファイルで、表１に示す形式で用意されている。
|
|
mu |
sigma |
Prob |
Yes |
No |
|
帰無仮説 |
帰無仮説の 正規分布の 平均値 |
帰無仮説の 正規分布の 標準偏差 |
帰無仮説が 正しい確率 |
帰無仮説が 正しいときに 有意になった効用 |
帰無仮説が 正しいときに 有意でない効用 |
|
対立仮設 |
対立仮設の 正規分布の 平均値 |
対立仮設の 正規分布の 標準偏差 |
対立仮設が 正しい確率 |
対立仮設が 正しいときに 有意になった効用 |
対立仮設が 正しいときに 有意でない効用 |
２行目が帰無仮説についてのパラメータであり、３行目が対立仮設についてのパラメータである。
2列目に仮説に対応する正規分布の平均値を、３列目に標準偏差を書く。４列目にそれぞれの仮説が正しい（成り立っている）確率を書く。
５列目に、有意な結果が得られたときの効用、６列目に有意でない結果が得られたときの効用を設定する。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
次の例では、図２のファイルsample.csvが設定されている。データファイルとスクリプトファイルは同じフォルダ内に置かれているものとする。
=
RESTART: C:\*****\optalphafiles\utility_curve.py
Input
data file (*.csv) = sample.csv
データファイル名を設定してEnterキーを押すと計算が始まり、計算結果が表示される（図３）。
図３
左側のグラフに期待効用値のグラフが最大値を与える値とともに表示されている。期待効用を最大にする有意水準は、
であると示されている。
この有意水準の値は、通常用いられているものより高い値である。次のパターン２では、帰無仮説が正しいにもかかわらず有意差が認められたときの効用を高い負の効用として設定してみる。
帰無仮説が正しいときに有意な差が出たときの効用を、大きい負の値
に設定した（図４）。
図４
リスト１のスクリプトを実行して、図４のファイルを入力データファイルとすると、図５のグラフを得た。
図５
期待効用が最大となる有意水準が、
であることが示されている。
F分布の場合のグラフは、図６のようである。
図６
帰無仮説は非心度０のF分布で表されるが、対立仮設は非心度が正の値である非心F分布で表される。
このとき、期待効用は有意水準の関数として、次式で表すことができる（再掲）。
パラメータの幾つかのパターンについて、以下のように計算してみた。スクリプトは、リスト２に示すものである。スクリプトファイルなどは、まとめて圧縮ファイルoptalphafiles.zipとした。期待効用の合理的な観点からは、最適な有意水準としてパターン３では27.2％、パターン４では2.3％が得られた。
図７に、使用した入力データファイルを示す。
図７
入力データファイルは、CSV形式のファイルで、表２に示す形式で用意されている。
|
表２ |
||||
|
|
非心度 |
Prob |
Yes |
No |
|
帰無仮説 |
このセルの値は用いない。空セルでよい |
帰無仮説が 正しい確率 |
帰無仮説が 正しいときに 有意になった効用 |
帰無仮説が 正しいときに 有意でない効用 |
|
対立仮設 |
対立仮設の 非心度 |
対立仮設が 正しい確率 |
対立仮設が 正しいときに 有意になった効用 |
対立仮設が 正しいときに 有意でない効用 |
|
分子の自由度 |
df1 |
|
|
|
|
分母の自由度 |
df2 |
|
|
|
２行目が３列目以降に帰無仮説についてのパラメータが設定されている。
３行目が対立仮設についてのパラメータであり、２列目に非心度が設定されている。
2行目、３行目とも、３列目にそれぞれの仮説が正しい確率、４列目に有意な結果が得られたときの効用、５列目に有意な結果が得られなかったときの効用が設定されている。
4行目、５行目の２列目に自由度df1とdf2が設定されている。
る。
リスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。
次の例では、図７のファイルFsample.csvが設定されている。データファイルとスクリプトファイルは同じフォルダ内に置かれているものとする。
=
RESTART: C:\*****\optalphafiles\utility_curve.py
Input
data file (*.csv) = Fsample.csv
データファイル名を設定してEnterキーを押すと計算が始まり、計算結果が表示される（図８）。
図８
左側のグラフに期待効用値のグラフが最大値を与える値とともに表示されている。期待効用を最大にする有意水準は、
であると示されている。
この有意水準の値は、通常用いられているものより高い値である。次のパターン４では、帰無仮説が正しいにもかかわらず有意差が認められたときの効用を高い負の効用として設定してみた。
帰無仮説が正しいときに有意な結果が得られたときの効用をとした設定を用意した（図９）。
図９
図９のファイルを入力データファイルとしてリスト２のスクリプトを実行すると、図10の結果を得た。
図10
最適な有意水準として
を得た。
図T1
帰無仮説は非心度０のｔ分布で表されるが、対立仮設は非心度が正の値である非心ｔ分布で表される。
このとき、期待効用は有意水準の関数として、次式で表すことができる（再掲）。
パラメータの幾つかのパターンについて、以下のように計算してみた。スクリプトは、リスト３に示すものである。スクリプトファイルなどは、まとめて圧縮ファイルoptalphafiles.zipとした。期待効用の合理的な観点からは、最適な有意水準としてパターンT1では31.8％、パターンT2では1.1％が得られた。
図T２
入力データファイルは、CSV形式のファイルで、表T1に示す形式で用意されている。
|
表T1 |
||||
|
|
パラメータ |
Prob |
Yes |
No |
|
帰無仮説 |
このセルの値は用いない。空セルでよい |
帰無仮説が 正しい確率 |
帰無仮説が 正しいときに 有意になった効用 |
帰無仮説が 正しいときに 有意でない効用 |
|
対立仮設(nc) |
対立仮設の 非心度 |
対立仮設が 正しい確率 |
対立仮設が 正しいときに 有意になった効用 |
対立仮設が 正しいときに 有意でない効用 |
|
自由度 |
df |
|
|
|
２行目が３列目以降に帰無仮説についてのパラメータが設定されている。
３行目が対立仮設についてのパラメータであり、２列目に非心度が設定されている。
2行目、３行目とも、３列目にそれぞれの仮説が正しい確率、４列目に有意な結果が得られたときの効用、５列目に有意な結果が得られなかったときの効用が設定されている。
4行目の２列目に自由度dfが設定されている。
リスト３のスクリプトを実行すると、入力データファイル名の設定が求められる。
次の例では、図T2のファイルt_sample.csvが設定されている。データファイルとスクリプトファイルは同じフォルダ内に置かれているものとする。
= RESTART:
C:\*****\optalphafiles\t_utility_curve.py
Input data file (*.csv) = t_sample.csv
データファイル名を設定してEnterキーを押すと計算が始まり、計算結果が表示される（図T3）。
図T3
左側のグラフに期待効用値のグラフが最大値を与える値とともに表示されている。期待効用を最大にする有意水準は、
であると示されている。
この有意水準の値は、通常用いられているものより高い値である。次のパターンT2では、帰無仮説が正しいにもかかわらず有意差が認められたときの効用を高い負の効用として設定してみた。
帰無仮説が正しいときに有意な結果が得られたときの効用をとした設定を用意した（図T4）。
図T4
図T4のファイルを入力データファイルとしてリスト３のスクリプトを実行すると、図T5の結果を得た。
図T5
最適な有意水準として
を得た。
カイ二乗分布の場合のグラフは、図C1のようである。
図C1
帰無仮説は非心度０のカイ二乗分布（上のグラフ）で表されるが、対立仮設は非心度が正の値である非心カイ二乗分布（下のグラフ）で表される。
このとき、期待効用は有意水準の関数として、次式で表すことができる（再掲）。
パラメータの幾つかのパターンについて、以下のように計算してみた。スクリプトは、リスト４に示すものである。スクリプトファイルなどは、まとめて圧縮ファイルoptalphafiles.zipとした。期待効用の合理的な観点からは、最適な有意水準としてパターンC1では23.4％、パターンC2では2.5％が得られた。
図C2に、使用した入力データファイルを示す。
図C２
入力データファイルは、CSV形式のファイルで、表C1に示す形式で用意されている。
|
表C1 |
||||
|
|
パラメータ |
Prob |
Yes |
No |
|
帰無仮説 |
このセルの値は用いない。空セルでよい |
帰無仮説が 正しい確率 |
帰無仮説が 正しいときに 有意になった効用 |
帰無仮説が 正しいときに 有意でない効用 |
|
対立仮設(nc) |
対立仮設の 非心度 |
対立仮設が 正しい確率 |
対立仮設が 正しいときに 有意になった効用 |
対立仮設が 正しいときに 有意でない効用 |
|
自由度 |
df |
|
|
|
２行目が３列目以降に帰無仮説についてのパラメータが設定されている。
３行目が対立仮設についてのパラメータであり、２列目に非心度が設定されている。
2行目、３行目とも、３列目にそれぞれの仮説が正しい確率、４列目に有意な結果が得られたときの効用、５列目に有意な結果が得られなかったときの効用が設定されている。
4行目の２列目に自由度dfが設定されている。
リスト４のスクリプトを実行すると、入力データファイル名の設定が求められる。
次の例では、図C2のファイルchi2_sample.csvが設定されている。データファイルとスクリプトファイルは同じフォルダ内に置かれているものとする。
= RESTART:
C:\*****\optalphafiles\chi_utility_curve.py
Input data file (*.csv) =
chi2_sample.csv
データファイル名を設定してEnterキーを押すと計算が始まり、計算結果が表示される（図C3）。
図C3
左側のグラフに期待効用値のグラフが最大値を与える値とともに表示されている。期待効用を最大にする有意水準は、
であると示されている。
この有意水準の値は、通常用いられているものより高い値である。次のパターンC2では、帰無仮説が正しいにもかかわらず有意差が認められたときの効用を高い負の効用として設定してみた。
帰無仮説が正しいときに有意な結果が得られたときの効用をとした設定を用意した（図C4）。
図C4
図C4のファイルを入力データファイルとしてリスト４のスクリプトを実行すると、図C5の結果を得た。
図C5
最適な有意水準として
を得た。
参考文献
岡本安晴「いまさら聞けないPythonでデータ分析」、2019、丸善出版
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
class ExpectedUtil():
"""
期待効用値算出クラス
"""
def __init__(self, p_H0,
p_H1, U0sig, U0nsig, U1sig, U1nsig,
mu0, sgm0, mu1, sgm1):
self.p_H0 = p_H0
self.p_H1 = p_H1
self.U0sig = U0sig
self.U0nsig = U0nsig
self.U1sig = U1sig
self.U1nsig = U1nsig
self.mu0 = mu0
self.sgm0 = sgm0
self.mu1 = mu1
self.sgm1 = sgm1
def calc_criterion(self, a):
return
ss.norm.ppf(1.0-a, self.mu0, self.sgm0)
def expectedU(self, p):
a =
p
c =
self.calc_criterion(a)
b =
ss.norm.cdf(c, self.mu1, self.sgm1)
exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \
self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)
return exp_util
def calc_beta(self, a):
c =
self.calc_criterion(a)
return ss.norm.cdf(c, self.mu1, self.sgm1)
pd_df = pd.read_csv(input('Input data
file (*.csv) = '))
print(pd_df)
vals = pd_df.values
print(vals)
mu0 = vals[0][1]
sgm0 = vals[0][2]
mu1 = vals[1][1]
sgm1 = vals[1][2]
p_H0 = vals[0][3] / (vals[0][3] +
vals[1][3])
p_H1 = vals[1][3] / (vals[0][3] +
vals[1][3])
U0sig = vals[0][4]
U0nsig = vals[0][5]
U1sig = vals[1][4]
U1nsig = vals[1][5]
plt.figure(figsize=(12,5))
plt.subplot(1,3,(1,2)) # 左２区画に期待効用曲線を描画
e_u = ExpectedUtil(p_H0, p_H1, U0sig,
U0nsig, U1sig, U1nsig,
mu0, sgm0, mu1, sgm1)
ps = np.linspace(0.0001, 0.9999,
10000)
expected_utils = e_u.expectedU(ps)
ck_idx =
np.argsort(-expected_utils)
# 最大値を与えるインデックス
opt_a = ps[ck_idx[0]]
# 最大期待効用を与える有意水準
print('opt_a =\n', opt_a)
plt.plot(ps, expected_utils,
label=f'Exp.U')
plt.plot([opt_a],
[expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,
label=fr'opt-$\alpha$={opt_a:.3f}')
plt.xlabel(r'$\alpha$', fontsize= 18)
plt.ylabel('Expect.Utility',
fontsize=18)
plt.legend()
plt.title(fr'$\mu0$={mu0},
$\sigma0$={sgm0}, p0={p_H0:.3f},' +
f'U0sig={U0sig}, U0nsig={U0nsig}' + '\n' +
fr'$\mu1$={mu1}, $\sigma1$={sgm1}, p1={p_H1:.3f},' +
f'U1sig={U1sig}, U1nsig={U1nsig}', fontsize=16)
plt.tight_layout()
plt.subplot(1, 3, 3)
# αに対するβと検出力のグラフ
alphas = np.linspace(0.0001, 0.9999,
1000)
betas = e_u.calc_beta(alphas)
pwrs = 1.0 - betas
plt.plot(alphas, betas,
label=r'$\beta$')
plt.plot(alphas, pwrs, label='power')
plt.plot(opt_a, e_u.calc_beta(opt_a),
'ro', markersize=12,
label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}')
plt.xlabel(r'$\alpha$', fontsize=18)
plt.ylabel(r'$\beta$, power',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.special as sspec
class ExpectedUtil():
"""
期待効用値算出クラス、F分布対応
"""
def __init__(self, p_H0,
p_H1, U0sig, U0nsig, U1sig, U1nsig,
df1, df2, nc):
self.p_H0
= p_H0
self.p_H1 = p_H1
self.U0sig = U0sig
self.U0nsig = U0nsig
self.U1sig = U1sig
self.U1nsig = U1nsig
self.df1 = df1
self.df2 = df2
self.nc = nc
def calc_criterion(self, a):
c =
sspec.fdtri(self.df1, self.df2, 1.0 - a)
return c
def expectedU(self, p):
a =
p
c =
self.calc_criterion(a)
b =
sspec.ncfdtr(self.df1, self.df2, self.nc, c)
exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \
self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)
return exp_util
def calc_beta(self, a):
c =
sspec.fdtri(self.df1, self.df2, 1.0 - a)
return
sspec.ncfdtr(self.df1, self.df2, self.nc, c)
pd_df = pd.read_csv(input('Input data
file (*.csv) = '))
print(pd_df)
vals = pd_df.values
print(vals)
nc = vals[1][1]
df1 = vals[2][1]
df2 = vals[3][1]
p_H0 = vals[0][2] / (vals[0][2] +
vals[1][2])
p_H1 = vals[1][2] / (vals[0][2] +
vals[1][2])
U0sig = vals[0][3]
U0nsig = vals[0][4]
U1sig = vals[1][3]
U1nsig = vals[1][4]
print('nc =', nc)
print('df1 =', df1)
print('df2 =', df2)
print('p_H0 =', p_H0)
print('p_H1 =', p_H1)
print('U0sig =', U0sig)
print('U0nsig =', U0nsig)
print('U1sig =', U1sig)
print('U1nsig =', U1nsig)
plt.figure(figsize=(12,5))
plt.subplot(1,3,(1,2))
# 左２区画に期待効用曲線を描画
e_u = ExpectedUtil(p_H0, p_H1, U0sig,
U0nsig, U1sig, U1nsig,
df1, df2, nc)
ps = np.linspace(1.0e-9, 1.0 - 1.0e-9,
10000)
expected_utils = e_u.expectedU(ps)
ck_idx =
np.argsort(-expected_utils)
# 最大値を与えるインデックス
opt_a = ps[ck_idx[0]]
# 最大期待効用を与える有意水準
print('opt_a =\n', opt_a)
plt.plot(ps, expected_utils,
label=f'Exp.U')
plt.plot([opt_a],
[expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,
label=fr'opt-$\alpha$={opt_a:.3f}')
plt.xlabel(r'$\alpha$', fontsize= 18)
plt.ylabel('Expect.Utility',
fontsize=18)
plt.legend()
plt.title(fr'nc={nc}, df1={df1},
df2={df2}' + '\n' +
fr'U0sig={U0sig}, U0nsig={U0nsig}, p0={p_H0:.3f}' + '\n' +
fr'U1sig={U1sig}, U1nsig={U1nsig}, p1={p_H1:.3f}', fontsize=16)
plt.tight_layout()
plt.subplot(1, 3, 3)
# αに対するβと検出力のグラフ
alphas = np.linspace(1.0e-9, 1.0 -
1.0e-9, 10_000)
betas = e_u.calc_beta(alphas)
pwrs = 1.0 - betas
plt.plot(alphas, betas,
label=r'$\beta$')
plt.plot(alphas, pwrs, label='power')
plt.plot(opt_a, e_u.calc_beta(opt_a),
'ro', markersize=12,
label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}')
plt.xlabel(r'$\alpha$', fontsize=18)
plt.ylabel(r'$\beta$, power',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.special as sspec
class ExpectedUtil():
"""
期待効用値算出クラス、t分布対応
"""
def __init__(self, p_H0, p_H1,
U0sig, U0nsig, U1sig, U1nsig,
df, nc):
self.p_H0 = p_H0
self.p_H1 = p_H1
self.U0sig = U0sig
self.U0nsig = U0nsig
self.U1sig = U1sig
self.U1nsig = U1nsig
self.df = df
self.nc = nc
def calc_criterion(self, a):
c =
sspec.stdtrit(self.df, 1.0 - a)
return c
def expectedU(self, p):
a =
p
c =
self.calc_criterion(a)
b =
sspec.nctdtr(self.df, self.nc, c)
exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \
self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)
return exp_util
def calc_beta(self, a):
c =
sspec.stdtrit(self.df, 1.0 - a)
return sspec.nctdtr(self.df, self.nc, c)
pd_df = pd.read_csv(input('Input data
file (*.csv) = '))
print(pd_df)
vals = pd_df.values
print(vals)
nc = vals[1][1]
df = vals[2][1]
p_H0 = vals[0][2] / (vals[0][2] +
vals[1][2])
p_H1 = vals[1][2] / (vals[0][2] +
vals[1][2])
U0sig = vals[0][3]
U0nsig = vals[0][4]
U1sig = vals[1][3]
U1nsig = vals[1][4]
print('nc =', nc)
print('df =', df)
print('p_H0 =', p_H0)
print('p_H1 =', p_H1)
print('U0sig =', U0sig)
print('U0nsig =', U0nsig)
print('U1sig =', U1sig)
print('U1nsig =', U1nsig)
import sys
#sys.exitc()
plt.figure(figsize=(12,5))
plt.subplot(1,3,(1,2))
# 左２区画に期待効用曲線を描画
e_u = ExpectedUtil(p_H0, p_H1, U0sig, U0nsig,
U1sig, U1nsig,
df, nc)
ps = np.linspace(1.0e-9, 1.0 - 1.0e-9,
10000)
expected_utils = e_u.expectedU(ps)
ck_idx =
np.argsort(-expected_utils)
# 最大値を与えるインデックス
opt_a = ps[ck_idx[0]]
# 最大期待効用を与える有意水準
print('opt_a =\n', opt_a)
plt.plot(ps, expected_utils,
label=f'Exp.U')
plt.plot([opt_a],
[expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,
label=fr'opt-$\alpha$={opt_a:.3f}')
plt.xlabel(r'$\alpha$', fontsize= 18)
plt.ylabel('Expect.Utility',
fontsize=18)
plt.legend()
plt.title(fr'nc={nc}, df={df}' + '\n'
+
fr'U0sig={U0sig}, U0nsig={U0nsig}, p0={p_H0:.3f}' + '\n' +
fr'U1sig={U1sig}, U1nsig={U1nsig}, p1={p_H1:.3f}', fontsize=16)
plt.tight_layout()
plt.subplot(1, 3, 3)
# αに対するβと検出力のグラフ
alphas = np.linspace(1.0e-9, 1.0 -
1.0e-9, 10_000)
betas = e_u.calc_beta(alphas)
pwrs = 1.0 - betas
plt.plot(alphas, betas,
label=r'$\beta$')
plt.plot(alphas, pwrs, label='power')
plt.plot(opt_a, e_u.calc_beta(opt_a),
'ro', markersize=12,
label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}')
plt.xlabel(r'$\alpha$', fontsize=18)
plt.ylabel(r'$\beta$, power',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
import scipy.special as sspec
class ExpectedUtil():
"""
期待効用値算出クラス、カイ二乗分布対応
"""
def __init__(self, p_H0,
p_H1, U0sig, U0nsig, U1sig, U1nsig,
df, nc):
self.p_H0 = p_H0
self.p_H1 = p_H1
self.U0sig = U0sig
self.U0nsig = U0nsig
self.U1sig = U1sig
self.U1nsig = U1nsig
self.df = df
self.nc = nc
def calc_criterion(self, a):
c =
ss.chi2.ppf(1.0 - a, self.df) # sspec.stdtrit(self.df, 1.0 - a)
return c
def expectedU(self, p):
a =
p
c =
self.calc_criterion(a)
b =
ss.ncx2.cdf(c, self.df, self.nc) # sspec.nctdtr(self.df, self.nc, c)
exp_util = self.p_H0 * (a * self.U0sig + (1 - a) * self.U0nsig) + \
self.p_H1 * ((1 - b) * self.U1sig + b * self.U1nsig)
return exp_util
def calc_beta(self, a):
c =
ss.chi2.ppf(1.0 - a, self.df) #
sspec.stdtrit(self.df, 1.0 - a)
return ss.ncx2.cdf(c, self.df, self.nc) # sspec.nctdtr(self.df, self.nc, c)
pd_df = pd.read_csv(input('Input data
file (*.csv) = '))
print(pd_df)
vals = pd_df.values
print(vals)
nc = vals[1][1]
df = vals[2][1]
p_H0 = vals[0][2] / (vals[0][2] +
vals[1][2])
p_H1 = vals[1][2] / (vals[0][2] +
vals[1][2])
U0sig = vals[0][3]
U0nsig = vals[0][4]
U1sig = vals[1][3]
U1nsig = vals[1][4]
print('nc =', nc)
print('df =', df)
print('p_H0 =', p_H0)
print('p_H1 =', p_H1)
print('U0sig =', U0sig)
print('U0nsig =', U0nsig)
print('U1sig =', U1sig)
print('U1nsig =', U1nsig)
import sys
#sys.exitc()
plt.figure(figsize=(12,5))
plt.subplot(1,3,(1,2))
# 左２区画に期待効用曲線を描画
e_u = ExpectedUtil(p_H0, p_H1, U0sig,
U0nsig, U1sig, U1nsig,
df, nc)
ps = np.linspace(1.0e-9, 1.0 - 1.0e-9,
10000)
expected_utils = e_u.expectedU(ps)
ck_idx =
np.argsort(-expected_utils)
# 最大値を与えるインデックス
opt_a = ps[ck_idx[0]]
# 最大期待効用を与える有意水準
print('opt_a =\n', opt_a)
plt.plot(ps, expected_utils,
label=f'Exp.U')
plt.plot([opt_a],
[expected_utils[ck_idx[0]]], 'ro', markersize=16, lw=0,
label=fr'opt-$\alpha$={opt_a:.3f}')
plt.xlabel(r'$\alpha$', fontsize= 18)
plt.ylabel('Expect.Utility',
fontsize=18)
plt.legend()
plt.title(fr'nc={nc}, df={df}' + '\n'
+
fr'U0sig={U0sig}, U0nsig={U0nsig}, p0={p_H0:.3f}' + '\n' +
fr'U1sig={U1sig}, U1nsig={U1nsig}, p1={p_H1:.3f}', fontsize=16)
plt.tight_layout()
plt.subplot(1, 3, 3)
# αに対するβと検出力のグラフ
alphas = np.linspace(1.0e-9, 1.0 -
1.0e-9, 10_000)
betas = e_u.calc_beta(alphas)
pwrs = 1.0 - betas
plt.plot(alphas, betas,
label=r'$\beta$')
plt.plot(alphas, pwrs, label='power')
plt.plot(opt_a, e_u.calc_beta(opt_a),
'ro', markersize=12,
label=fr'$\beta$={e_u.calc_beta(opt_a):.3f}')
plt.xlabel(r'$\alpha$', fontsize=18)
plt.ylabel(r'$\beta$, power',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()