z検定とt検定をPythonで解く

前書き

大学の課題でz検定とt検定を使い間違えるという醜態を晒してしまったため、二度とそのようなことがないよう、このように記事に認めておくこととする。
未だ若輩者なので、間違っていることがあればコメントなどでお手柔らかにご指摘願いたい。
また、当記事は栗原　伸一先生の『入門統計学　ー検定から多変量解析・実験計画法までー』　の内容をメインに執筆した。引用には細心の注意を払うが、不十分な場合にはその旨をコメント欄にてご指摘願いたい。

z分布(正規分布あるいはガウス分布)

z分布(正規分布)とはフランスの数学者ド・モアブルが**「試行回数が増え、計算が難しくなった二項分布（ベルヌーイ試行）」**を数式で近似したものである。(二項分布の説明は省略)
正規分布の確率密度関数

$$\mbox{f}(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$

この式を基に、python3のnumpyを用いて関数を定義し、matplotlibにて図表化した。
(Jupyter notebookを用いたセルプログラミング環境で実行)

ライブラリのインポート

numpy=1.14.3
matplotlib=2.2.2

import numpy as np
from matplotlib import pyplot as plt

z分布（正規分布）の確率密度関数を定義する

z分布（正規分布）の確率密度関数
$$\mbox{f}(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$

def pdf_norm(x, mu, sigma):
    destiny = 1 / np.sqrt(2 * np.pi *  sigma) * np.exp(-(x - mu)**2/(2*sigma))
    return destiny

matplotlib.pyplotを用いた標準正規分布(平均０、標準偏差1)のプロット

x = np.arange(-10, 10, 0.01)
y = [pdf_norm(i, mu=0, sigma=1) for i in x]

plt.plot(x, y, label="σ=1, µ=0")
plt.legend()
plt.show()

参考

https://gcbgarden.com/2017/05/07/python-gauss/

z検定

z検定は標準化（平均を０、標準偏差を1になるように整形）した値を元に、z(標準正規)分布によって生起確率を求め、それが分布のどこに位置するのかを求め、標本平均から推定された値が母平均値に等しいかを判断するための検定。
よって分布が正規分布性を持ち、母平均か母分散が既知でなければならない。
しかし、**標本数が多く(標本数 > 30)**あるのであれば、母分散が未知の場合であっても区間推定に利用することができる。
区間推定のための信頼区間（その分布が大体母平均値に合致する区間）公式は、有意水準（おそらく母平均値ではないと考えられる確率領域）を5%とすると、
正規分布の区間推定公式（両側）
$$-1.96<\frac{\overline{x}-\mu}{\frac{\sigma}{\sqrt{n}}}<1.96$$

参考

http://www.geisya.or.jp/~mwm48961/linear_algebra/t_test2.htm

Python3によるz検定

通常、手計算などではz分布表などを用いてパーセント点(先の公式での[-1.96, 1.96])を求めるが、Pythonの理数系ライブラリScipyのstats moduleを用いることによって、より簡単に求めることができる。
これを利用し、z検定の信頼区間とz値を返す関数を作成してみる。
以下、jupyter notebookにて作成したセルプログラミング。

ライブラリのインポート

numpy=1.14.3
scipy.stats=1.1.0

import numpy as np
from scipy import stats as ss

z値と信頼区間を計算する関数を定義

# scipy.stats.norm.ppf -> パーセント点を返してくれるメソッド
def interval_nom(alpha, x, mu, sigma, n, tile="one"):
    z = (x - mu) / (sigma/np.sqrt(n))
    z_critical_value = []
    if tile == "one":
        z_critical_value.append(ss.norm.ppf(1-alpha))
    elif tile == "two":
        # 両側検定なので、棄却域を/2している
        z_critical_value.append(ss.norm.ppf(1-alpha/2))
        z_critical_value.append(-ss.norm.ppf(1-alpha/2))
    else:
        print("Tile option error. tile option is only \'{0}\' or \'{1}\' ".format("one", "two"))
        return 
    
    return z, z_critical_value
    
    

例題

高校数学の基本問題/統計/z検定, t検定

【例１】
　標本平均が——x=56，母標準偏差がσ=8，標本の大きさがn=16のとき，
帰無仮説　H0：μ=50
対立仮説　H1：μ≠50
の検定をせよ

z, z_critical_value = interval_nom(alpha=0.05,x=56, mu=50, sigma=8, n=16, tile="two")
print("z value: {:.3f}".format(z))
if len(z_critical_value) > 1:
    print("z critical velue: {1:.3f} < z < {0:.3f}".format(z_critical_value[0], z_critical_value[1]))
else:
    print("z critical velue: {0:.3f}".format(z_critical_value))

z value: 3.000
z critical velue: -1.960 < z < 1.960

解釈

z値が有意水準-1.96 < z < 1.960外にあるため、帰無仮説を棄却する。

t分布

t分布の説明については、冒頭でも紹介した栗原　伸一先生の『入門統計学　ー検定から多変量解析・実験計画法までー』、p66の内容を引用させていただくと、

母分散が未知でサイズの小さな標本から母平均などの推定・検定を行う場合に標準正規分布の代わりに用いる連続型の確率分布。自由度によって分布の形が変化するが,　たいていは正規分布よりも両裾が厚くなる。

栗原　伸一 (2011) 『入門統計学　ー検定から多変量解析・実験計画法までー』 pp. 66　オーム社

また、標本数が多くなれば多くなるほどt分布は標準正規分布に近似していくことが知られています。

t分布の確率密度関数
$$
f(t) = \frac{\Gamma(\frac{n}{2}) } {\sqrt{(n-1)\pi}\Gamma(\frac{n-1}{2})} (1+\frac{t^2}{(n-1)})^\frac{-n}{2}
$$

ライブラリのインポート

numpy=1.14.3
matplotlib=2.2.2

import numpy as np
import math
import matplotlib.pyplot as plt

t分布の確率密度関数を定義する

t分布の確率密度関数
$$
f(t) = \frac{\Gamma(\frac{n}{2}) } {\sqrt{(n-1)\pi }\Gamma(\frac{n-1}{2})} (1+\frac{t^2}{(n-1)})^\frac{-n}{2}
$$
t値
$$
t = \frac{\overline{X} - \mu}{\frac{\hat{\sigma}}{\sqrt{n-1}}}
$$

def pdf_t(x, mu, sigmma,  n):
    t = (x - mu) / (sigmma / np.sqrt(n-1))
    destiny = (math.gamma(n / 2) / (np.sqrt((n-1) * np.pi) * math.gamma((n-1) / 2))) * (1 + (t**2/(n-1)))**(-n / 2)
    return destiny

matplotlib.pyplotを用いたt分布のプロット

x = np.arange(-10, 10, 0.01)
y = [pdf_t(i, mu=0, sigmma=1, n=10) for i in x]

plt.plot(x, y, label="σ=1, µ=0, df=10")
plt.legend()
plt.show()

ライブラリのインポート

import numpy as np
import scipy.stats as ss
import math

t値と信頼区間を計算する関数を定義

# scipy.stats.t.ppf() -> t分布のパーセント点を返してくれるメソッド
def interval_t(alpha, x, mu, sigmma, n, tile="one"):
    import numpy as np
    import scipy.stats as ss
    import math
    t = (x - mu) / (sigmma/np.sqrt(n-1))
    t_critical_value = []
    if tile == "one":
        t_critical_value.append(ss.t.ppf(1-alpha, df=n-1))
    if tile == "two":
        # 両側検定なので、棄却域を/2している
        t_critical_value.append(ss.t.ppf(1-alpha/2, df=n-1))
        t_critical_value.append(-ss.t.ppf(1-alpha/2, df=n-1))
    else:
        print("Tile option error. tile option is only \'{0}\' or \'{1}\' ".format("one", "two"))
        return
    return t, t_critical_value
    

例題

統計WEB/ホーム > 統計の時間 > Step1. 基礎編 > 24. 平均値の検定 > 練習問題 (24.平均値の検定)(3)

次の表は、1つ25.5 kgの強力粉20個をサンプリングし、重量を測定した結果をまとめたものである。このデータを用いて、強力粉の重量は25.5 kgではないと言えるかどうか検定せよ。なお、有意水準は$\alpha = 0.05$とする。

項目	測定結果
サンプルサイズ	20
平均	25.29
不偏分散	2.23
標準偏差	1.49

1. 帰無仮説:H0 強力粉の重量は25.5kgである。

2. 対立仮説:H1 強力粉の重量は25.5kgではない。

t, t_critical_value = interval_t(0.05, 25.29, 25.5,1.49, 20, tile="two")
print("t value: {:.3f}".format(t))
if len(t_critical_value) > 1:
     print("t critical velue: {1:.3f} < t < {0:.3f}".format(t_critical_value[0], t_critical_value[1]))
else:
    print("t critical velue: {0:.3f}".format(t_critical_value))
        

t value: -0.614
t critical velue: -2.093 < t < 2.093

解釈

t値が信頼区間 -2.093 < t < 2.093内にあるため、
対立仮説を棄却し帰無仮説を有意水準5%で採択する。

使い分け

長々と書いてしまいましたが、

• t検定
  • 母分散未知
  • 小標本(標本数 < 30)
• z検定
  • 母分散既知
  • 大標本(標本数 > 30)
  • 大標本であれば母分散未知でも区間推定に用いることができる

ということだと解釈しました。

参考文献

• 栗原　伸一 (2011) 『入門統計学　ー検定から多変量解析・実験計画法までー』 オーム社
• 高校数学の基本問題/統計/z検定, t検定
• 統計WEB/ホーム > 統計の時間