概要とdisclaimer

統計メモ(1)の続きです．
Pythonで母分散$\sigma_0$が未知の場合の検定サンプルプログラムを書いてます．
２つの母集団の差の検定，対応のある２集団の母平均の検定も以後書いていきます．

母分散が未知の場合の一つの標本に対する母平均μの検定

ある標本を取ってきて母平均$\mu$を推定したとき，$\mu$が帰無仮説$H_0:\mu = \mu_1$に従うかを検定します．ここでは，母分散$\sigma_0^2$は未知です．

手順

手順自体は検定統計量が$u$だったときと同じです．母分散が未知のときは検定統計量$t$を計算します．（t検定）

仮説の設定
$$帰無仮説　H_0: \mu = \mu_0$$
$$対立仮説　H_1: \mu \neq \mu_0$$
有意水準$\alpha$の設定．今回は0.05に設定．
棄却域Rの設定
$$ R:|t_0| > t(\phi, \alpha)$$
ここで，$t(\phi, \alpha)$は$Pr[|t| \geq t(\phi, P)]=P$に従う__両側__100P%点
統計検定量$t$の計算．
帰無仮説$H_0$のもとで，検定統計量$t_0$は
$$t_0 = \frac{\bar{x} - \mu_0}{\sqrt{V/n}}$$
なお，$V$は不偏分散．($V=(偏差平方和)/(n-1)$)
計算した$t_0$が棄却域Rに入っていれば帰無仮説を棄却する．

検出力1-βの計算

母分散が未知のときの検出力は少々クセ者のようです．検出力$1-\beta$は対立仮説のもとで計算されますが，このとき検定統計量$t_0$は__非心t分布__に従います．「サンプルサイズの決め方(永田靖)」では手計算を想定していたのか，非常に丁寧に近似式を導出しています．（結構精度がいいようです）ただ，scipy.statsライブラリにはstats.nct()で非心t分布の厳密値が計算できるようです．よって，せっかくなので後述のサンプルプログラムでは近似値と精確値どちらも計算しています．

例えば

$n=10$のデータ例で検出力$1-\beta$を計算してみます．

$$データ：[6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3]$$

各種のパラメータを設定
$$\Delta=0.5, \mu_0=5.0(既知)$$
統計検定量$t_0$は帰無仮説$H_0$のもとで計算される．例えば今回の場合は$t_0 = 4.195$
検出力$1-\beta$の計算
対立仮説$H_1$のもとで，
$$1-\beta = Pr[|t_0| \geq t(\phi, \alpha)] = Pr[t_0 \leq -t(\phi, \alpha)] + Pr[t_0 \geq t(\phi, \alpha)]$$
ここで，非心t分布を計算する必要があります．上で紹介した参考文献の中では，近似式が与えられてり，
$$Pr[t' \leq w] \approx Pr\Bigl[\frac{w(1 - 1/(4\phi)) - \lambda}{\sqrt{1+w^2/(2\phi)}}\Bigr]$$

これらを用いて，$1-\beta$を計算すると$1-\beta = 0.293 (精確値:0.293)$

サンプルサイズの決め方

ここでのサンプルサイズの決め方は，参考文献の近似式を使用します．ある検出力$1-\beta$を保証したい場合にサンプルサイズ$n$をいくらにすればいいかの概算値を計算できます．ここで詳しく書くと~~めんどくさい~~長くなるのでここでは導出は省略します．

帰無仮説が

H_0: \mu_0 = \mu

対立仮説が

H_1: \mu_0 \neq \mu

であるとき，$|\Delta| = (\mu-\mu_0)/\sigma > \Delta_0$を高い検出力$1-\beta$で検知したいとき，必要なサンプルサイズ$n$は

n \approx \left(
\frac{z_{\alpha/2} - z_{1-\beta}}{\Delta_0}
\right)^2

+ \frac{z_{\alpha/2}^2}{2}

Pythonでの実装

sample.py

from scipy.stats import *
import numpy as np


class Mu_1sample:

    def __init__(self, data, alpha=0.05, mu0=0.0):
        self.n = len(data)
        self.df = self.n - 1
        self.alpha = alpha
        self.mu0 = mu0
        self.data = data

    def t_test(self, mode=1):
        """
        Null Hypothesis Ho: mu = mu0

        :param mode:
        Alternative Hypothesis
        1)H1: mu != mu0 -> two-sided test
        2)H1: mu > mu0  -> one-sided test
        3)H1: mu < mu0  -> one-sided test
        :return:
        """
        V = np.var(self.data, ddof=1)
        t_0 = (np.average(self.data) - self.mu0)/np.sqrt(V/self.n)
        print("t_0: ", t_0)

        rv = t(df=self.df)

        if mode == 1:
            R = rv.isf(self.alpha/2)
            print("R: ", R)
            if abs(t_0) >= R:
                print("Ho rejected.")
            else:
                print("Ho not rejected")
        elif mode == 2:
            R = rv.isf(self.alpha)
            print("R: ", R)
            if t_0 >= R:
                print("Ho rejected")
            else:
                print("Ho not rejected")
        elif mode == 3:
            R = rv.ppf(self.alpha)
            print("R: ", R)
            if t_0 <= R:
                print("Ho rejected")
            else:
                print("Ho not rejected")

    def statistic_power(self, delta=0.5, mode=1):
        lambda_ = np.sqrt(self.n) * delta  # parameter of non-centrality

        rv1 = t(df=self.df)  # t distribution of df = n-1
        rv2 = nct(df=self.df, nc=lambda_)  # non-central t distribution of df=n-1, nc=lambda_

        if mode == 1:
            # approximation
            w = rv1.isf(self.alpha/2)
            power = norm.cdf((-1*w*(1-1/(4*self.df)) - lambda_) / np.sqrt(1 + w**2/(2*self.df))) + 1 - \
                    norm.cdf((w*(1-1/(4*self.df)) - lambda_) / np.sqrt(1 + w**2/(2*self.df)))
            print("power(approx.): ", power)

            # exact value
            power = rv2.cdf(-1*t.isf(self.alpha/2, df=self.df)) + \
                    rv2.sf(t.isf(self.alpha/2, df=self.df))
            print("power(exact): ", power)
        elif mode == 2:
            # approximation
            w = rv1.isf(self.alpha)
            power = 1 - norm.cdf((w*(1-1/(4*self.df)) - lambda_) / np.sqrt(1 + w**2/(2*self.df)))
            print("power(aprox.): ", power)

            # exact value
            power = rv2.sf(t.isf(self.alpha, df=self.df))
            print("power(exact): ", power)

        elif mode == 3:
            # approximation
            w = rv1.isf(self.alpha)
            power = norm.cdf((-1*w*(1-1/(4*self.df)) - lambda_) / np.sqrt(1 + w**2/(2*self.df)))
            print("power(approx.): ", power)

            # exact value
            power = rv2.cdf(-t.isf(self.alpha, df=self.df))
            print("power(exact): ", power)

class Sample_size_estimation:

    def __init__(self):
        pass

    def estimate(self, alpha=0.05, delta=0.5, required_power=0.90, mode=1):

        if mode == 1:
            n = ((norm.isf(alpha/2) - norm.isf(required_power)) / delta)**2 + 0.5 * norm.isf(alpha/2)**2

            n_ = np.floor(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_-1, nc=lambda_)
            power = rv.cdf(-1*t.isf(alpha/2, df=n_-1)) + rv.sf(t.isf(alpha/2, df=n_-1))
            if power >= required_power:
                return n_

            n_ = np.ceil(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_ - 1, nc=lambda_)
            power = rv.cdf(-1 * t.isf(alpha/2, df=n_-1)) + rv.sf(t.isf(alpha/2, df=n_-1))
            if power >= required_power:
                return n_

        elif mode == 2:
            n = ((norm.isf(alpha) - norm.isf(required_power)) / delta)**2 + 0.5 * norm.isf(alpha)**2

            n_ = np.floor(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_ - 1, nc=lambda_)
            power = rv.sf(t.isf(alpha, df=n_ - 1))
            if power >= required_power:
                return n_

            n_ = np.ceil(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_ - 1, nc=lambda_)
            power = rv.sf(t.isf(alpha, df=n_ - 1))
            if power >= required_power:
                return n_

            return np.floor(n) + 2

        elif mode == 3:
            n = ((norm.isf(alpha) - norm.isf(required_power)) / delta) ** 2 + 0.5 * norm.isf(alpha) ** 2

            n_ = np.floor(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_ - 1, nc=lambda_)
            power = rv.cdf(-t.isf(alpha, df=n_ - 1))
            if power >= required_power:
                return n_

            n_ = np.ceil(n)
            lambda_ = np.sqrt(n_) * delta
            rv = nct(df=n_ - 1, nc=lambda_)
            power = rv.cdf(-t.isf(alpha, df=n_ - 1))
            if power >= required_power:
                return n_

            return np.floor(n) + 2


if __name__ == '__main__':
    # data = [6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3]
    # data = [10.8, 11.2, 9.7, 9.9, 12.0, 9.6, 10.5, 10.7, 10.1]
    data = [21, 19, 16, 19, 22, 18, 20, 21]
    test = Mu_1sample(data, alpha=0.05, mu0=20.0)
    test.t_test(mode=3)
    test.statistic_power(delta=-1.0, mode=3)

    n = Sample_size_estimation().estimate(delta=-1.5, required_power=0.95, mode=3)
    print(n)