#概要（というかdisclaimer)
-
Pythonで母分散$\sigma_0 ^2$が既知の場合の検定のサンプルプログラムを実装しています．scipy.statsライブラリを使っていますが，$N(0, 1^2)$の確率密度関数や$\chi^2$の値を求めるための利用にとどめ，ttest()などのデータを渡せば検定が終わってしまう関数は使いません．
-
自分向けのメモなので，細かい説明はかなり省略してます．申し訳ない．
-
対立仮説$H_1$は
3種類すべてについて書くとめんどくさいので$H_1:\mu \neq \mu_0$のケースのみを考えます． -
参考文献は「サンプルサイズの決め方(著：永田靖) 朝倉書店」，および「keisukeのブログ」さんです．
#母分散が既知の場合の1つの標本に対する母平均の検定
ある標本を取ってきて母平均$\mu$を推定したとき，$\mu$が帰無仮説$H_0:\mu = \mu_0$に従うかどうか検定します．ここでは，母分散$\sigma_0^2$が既知であると考え話を進めます．
実際の場面で母平均がわからないのに母分散が既知などという状態は普通ありえないので，この検定法は「実務的には」役立ちません．しかし，統計的検定の考え方を理解するのに非常にいい練習になりました．
#手順
-
仮説の設定
$H_0: \mu = \mu_0$ （なお，$\mu_0$は指定の値. 例えば，$\mu$が新製品バッテリーの持続時間で，$\mu_0$が旧製品バッテリーのそれ.$H_0$が棄却できれば，新しいバッテリーは持続時間は旧製品のそれから違うものとなっている．）
$H_1: \mu \neq \mu_0$ -
有意水準$\alpha$の設定．ここでは$\alpha = 0.05$と設定
-
棄却域Rの設定
$|u_0| \geq z_{\alpha/2}$
$z_P$は標準正規分布$N(0, 1^2)$の上側100P%点で，$Pr(u \geq z_P)=P$が成り立つ．
実際に計算すると，$z_{\alpha/2}$は$z_{\alpha/2}=1.960$ -
データ$x_1, x_2, ..., x_n$に対して，帰無仮説$H_0:\mu=\mu_0$のもとで検定統計量$u_0$の値を計算する．
$u_0 = \frac{\bar{x} - \mu_0}{\sqrt{\sigma_0^2/n}}$
例えば，以下のデータ($n=10, \mu_0=5.0, \sigma_0^2=1.0^2$)を取得したと仮定する．
データ：6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3
このとき，$u_0 = \frac{\bar{x} - \mu_0}{\sqrt{\sigma_0^2/n}} = 4.364$ -
計算した$u_0$が棄却域Rに入っていれば$H_0$は棄却され，対立仮説$H_1$を採択する．
いま，$u_0 = 4.364 > u_{\alpha/2} = 1.960$なので，帰無仮説は棄却される．
#検出力1-βの計算
t検定や$\chi^2$検定における第一種の過誤($\alpha$)について述べている入門書は多くありますが，検出力(power)についても述べているものは少ないようです．なので自分も「サンプルサイズの決め方」を読んでなんとなーくの理解しかありませんが，要は検出力$1-\beta$は$H_0$の状態から$\Delta_0$だけの差がある場合に「$H_1$が本当は真であるのに，$H_0$を棄却しない（第2種の過誤）」ミスが起こらないように検出できる確率であると解釈できるものです．
#例えば
上での$n=10$のデータ例で検出力$1-\beta$を計算してみます．
- $\Delta = \frac{\mu - \mu_0}{\sigma_0} = 0.5$，$\mu_0 = 5.0, \sigma_0^2=1.0^2$と設定する
- 検定統計量は帰無仮説$H_0:\mu = \mu_0$のもとで計算される．先に計算したように，$u_0 = 4.364$である．
- 検出力$1-\beta$を計算する．これは第二種の過誤を起こさない確率である．検出力は対立仮説$H_1:\mu \neq \mu_0$のもとで計算される．
$1 - \beta = Pr(|u_0| \geq z_{\alpha/2}) = Pr(\frac{\bar{x} - \mu_0}{\sqrt{\sigma_0^2/n}} \leq -z_{\alpha/2}) + Pr(\frac{\bar{x} - \mu_0}{\sqrt{\sigma_0^2/n}} \geq z_{\alpha/2})$
ここで，
$1 - \beta = Pr(\frac{\bar{x} - \mu}{\sqrt{\sigma_0^2/n}} + \frac{\mu - \mu_0}{\sqrt{\sigma_0^2/n}} \leq -z_{\alpha/2}) + Pr(\frac{\bar{x} - \mu}{\sqrt{\sigma_0^2/n}} + \frac{\mu - \mu_0}{\sqrt{\sigma_0^2/n}} \geq z_{\alpha/2})$
と変形できるので，$\Delta = (\mu - \mu_0)/\sigma_0$を考慮して，
$1 - \beta = Pr(u \leq -z_{\alpha/2} - \sqrt{n}\Delta) + Pr(z_{\alpha/2} - \sqrt{n}\Delta)$
と最終的に計算され，各種の値を代入すると$1-\beta = 0.353$とわかる．
#サンプルサイズの決め方
以上の議論を踏まえると，サンプルサイズを適切に設定することで検出力$1-\beta$をコントロールすることが出来ることがわかります．実際に，
$1 - \beta = Pr(u \leq -z_{\alpha/2} - \sqrt{n}\Delta) + Pr(z_{\alpha/2} - \sqrt{n}\Delta)$
を考えると，$(右辺第1項目)\approx 0$であるので，
$1 - \beta \approx Pr(z_{\alpha/2} - \sqrt{n}\Delta)$となります．
すなわち，標準正規分布の確率密度関数を$z_{\alpha/2} - \sqrt{n}\Delta$から∞まで積分すれば$1-\beta$に等しいということです．よって，上側100P%点の定義より，
$z_{1-\beta} \approx z_{\alpha/2} - \sqrt{n}\Delta$
検出したい最小の$\Delta$という意味で，$\Delta = \Delta_0$とおいて，この近似式をｎについて解くと，
$n \approx \Bigl(\frac{z_{\alpha/2} - z_{1-\beta}}{\Delta_0}\Bigr)^2$
を得ることができます．
#Pythonでの実装
from scipy.stats import *
import numpy as np
class Sigma0_Known:
def __init__(self, data, alpha=0.05, sigma0=1.0, mu0=0.0):
self.n = len(data)
self.df = self.n - 1
self.alpha = alpha
self.sigma0 = sigma0
self.mu0 = mu0
self.data = data
def statistic_test(self, mode=1):
"""
Null Hypothesis Ho: mu = mu0
:param mode:
Alternative Hypothesis
1)H1: mu != mu0
2)H1: mu > mu0
3)H1: mu < mu0
:return:
"""
u_0 = (np.average(self.data) - self.mu0)/np.sqrt((self.sigma0**2/self.n))
print("u_0: ", u_0)
# mode = 1: two-sided test
# mode = 2, 3: one-sided test
if mode == 1:
u = norm.isf(self.alpha/2, loc=0, scale=1)
print("u: ", u)
if abs(u_0) >= u:
print("Ho rejected. mu not equals mu_0")
else:
print("Ho not rejected")
elif mode == 2:
u = norm.isf(self.alpha, loc=0, scale=1)
print("u: ", u)
if u <= u_0:
print("Ho rejected. mu not equals mu_0")
else:
print("Ho not rejected")
elif mode == 3:
u = norm.ppf(self.alpha, loc=0, scale=1)
print("u: ", u)
if u_0 <= u:
print("Ho rejected. mu not equals mu_0")
else:
print("Ho not rejected")
def statistic_power(self, delta=0.5, mode=1):
"""
power = 1 - beta
:param delta:
:return:
"""
if mode == 1:
power = norm.cdf(-norm.isf(self.alpha/2)-np.sqrt(self.n)*delta) + \
norm.sf(norm.isf(self.alpha/2)-np.sqrt(self.n)*delta)
print("power: ", power)
elif mode == 2:
power = norm.sf(norm.isf(self.alpha)-np.sqrt(self.n)*delta)
print("power: ", power)
elif mode == 3:
power = norm.cdf(-norm.isf(self.alpha)-np.sqrt(self.n)*delta)
print("power: ", power)
class Sample_size_estimation:
def __init__(self):
pass
def estimate(self, alpha=0.05, delta0=0.5, required_power=0.90, mode=1):
if mode == 1:
n = ((norm.isf(alpha/2) - norm.isf(required_power))/delta0)**2
n_ = np.floor(n)
if norm.cdf(-norm.isf(alpha/2)-np.sqrt(n_)*delta0)+norm.sf(norm.isf(alpha/2)-np.sqrt(n_)*delta0) >= required_power:
return n_
n_ = np.ceil(n)
if norm.cdf(-norm.isf(alpha/2)-np.sqrt(n_)*delta0)+norm.sf(norm.isf(alpha/2)-np.sqrt(n_)*delta0) >= required_power:
return n_
elif mode == 2:
n =((norm.isf(alpha) - norm.isf(required_power))/delta0)**2
n_ = np.floor(n)
if norm.sf(norm.isf(alpha)-np.sqrt(n_)*delta0) >= required_power:
return n_
n_ = np.ceil(n)
if norm.sf(-norm.isf(alpha)-np.sqrt(n_)*delta0) >= required_power:
return n_
elif mode == 3:
n = ((norm.isf(alpha) - norm.isf(required_power)) / delta0) ** 2
n_ = np.floor(n)
if norm.cdf(-norm.isf(alpha)-np.sqrt(n_)*delta0) >= required_power:
return n_
n_ = np.ceil(n)
if norm.cdf(-norm.isf(alpha)-np.sqrt(n_)*delta0) >= required_power:
return n_
if __name__ == '__main__':
data = [21, 19, 16, 19, 22, 18, 20, 21]
test = Sigma0_Known(data, sigma0=2.0, mu0=20.0)
test.statistic_test(mode=3)
test.statistic_power(delta=-1.0, mode=3)
n = Sample_size_estimation().estimate(alpha=0.05, delta0=-1.5, required_power=0.95, mode=3)
print(n)
一応検算したのでミス等は無いと思います．