のまとめやコードを実行してみた際のメモです．統計的仮説検定
統計的仮説検定（検定）：標本を使って母集団に関する統計的な"判断"を下す方法．判断の基準として，確率的な表現を使う．
統計的仮説検定には様々な種類があり，判断する対象も手法によって様々です．
1変量データのt検定
1変量データのt検定の特徴
p値の計算方法
両側検定を行うことを前提としてp値を計算します．
標本から計算されたt値を$t_{標本}$と呼ぶことにします．t分布の累積分布関数を使うと，母平均を仮定した時にt値が$t_{標本}$を下回る確率を計算することができます．この確率を$\alpha$とするとp値は以下のように計算されます．
t検定の実装例
統計的仮説検定には様々な種類があり，判断する対象も手法によって様々です．
1変量データのt検定
1変量データのt検定の特徴
- 対象：平均値
- 判断すること：平均値が"ある値"と異なると言えるかどうか
ここでの"ある値"は，各々が設定する値です．
t検定の直感的な考え方
t検定では，以下の3つの条件が満たされているときに有意差があると判断します．
- サンプルサイズが大きい
- データのばらつき（分散）が小さい
- 平均値の差が大きい
t値
$t値 = \cfrac{標本平均 - 比較対象値}{標準偏差 \div \sqrt{サンプルサイズ}} = \cfrac{標本平均 - 比較対象値}{標準誤差}$
標本平均が比較対象と比べてとても小さかった場合は，t値が小さな値になります．t値はその絶対値に意味があると考えることができます．
統計的仮説検定の枠組み：帰無仮説・対立仮説
統計的仮説検定では，ある仮説を立てて，その仮説を棄却するか，しないかという判断を下すことで，データから客観的な判断を試みます．
帰無仮説：棄却される対象となる最初の仮説
対立仮説：帰無仮説と対立する仮説
帰無仮説が棄却された（帰無仮説が間違っていると判断された）ならば，有意性があると判断することができます．
p値
p値：標本と帰無仮説との矛盾となる指標
p値が小さいほど帰無仮説と標本が矛盾していると考えます．p値は確率として表現され，信頼区間と同様に，全く同じ条件で何度も標本抽出〜t値の計算を繰り返すことで確率を解釈します．
有意水準
p値が有意水準を下回った時に帰無仮説を棄却します．
有意水準：帰無仮説を棄却する基準となる値．5%が使われることが多く，危険率とも呼ばれる．
t検定とt分布の関係
t値の定義は以下のようになります．
$t値 = \cfrac{標本平均 - 母平均}{標準誤差} = \cfrac{{\hat{\mu}} - \mu}{{\hat{\mu}} / {\sqrt{N}}}$
ここに，$\hat{\mu}$：標本平均，$\mu$：母平均，$\hat{\sigma}$：標本から計算された標準偏差（不偏分散の平方根），$N$：サンプルサイズです．
p値の計算方法
両側検定を行うことを前提としてp値を計算します．
標本から計算されたt値を$t_{標本}$と呼ぶことにします．t分布の累積分布関数を使うと，母平均を仮定した時にt値が$t_{標本}$を下回る確率を計算することができます．この確率を$\alpha$とするとp値は以下のように計算されます．
p値 = $(1 - \alpha) \times 2$
ここで，最後に2倍しているのでは両側検定だからです．t検定の実装例
>>> import numpy as np
>>> import pandas as pd
>>> import scipy as sp
>>> from scipy import stats
>>>
>>> from matplotlib import pyplot as plt
>>> import seaborn as sns
>>> sns.set()
分析対象となるデータ（サンプルサイズ：20）をシリーズ形式で読み込みます．このデータはジャンクフードの重量を測定したデータです．
>>> junk_food = pd.read_csv("3-8-1-junk-food-weight.csv")["weight"]
>>> junk_food.head()
0 58.529820
1 52.353039
2 74.446169
3 52.983263
4 55.876879
Name: weight, dtype: float64
このデータを対象として，1標本のt検定を実行します．検定は以下の要領で行います．
帰無仮説：ジャンクフードの平均重量は50gである
対立仮説：ジャンクフードの平均重量は50gと異なる
有意水準は5%とし，p値が0.05を下回れば，帰無仮説は棄却され，ジャンクフードの重量は有意に50gと異なると主張ができることになります．
t値の計算式は
$t値 = \cfrac{標本平均 - 比較対象値}{標準誤差}$
です．最初に標本平均を求めます（以下の例ではScipyのmean関数を用いていますが，scipyのmean関数はSciPy 2.0.0で削除されるそうなので，代わりにnumpy.meanを使用しても同じ結果が得られます）．
>>> mu = sp.mean(junk_food) # np.mean(junk_food)
>>> mu
55.38496619666667
次に自由度を求めます．サンプルサイズから1を引くだけです．
>>> df = len(junk_food) - 1
>>> df
19
続いて標準誤差を求めます．標準誤差は
標準偏差 / サンプルサイズの平方根
で求められるので，以下のようになります（以下の例ではScipyのstd, sqrt関数を用いていますが，これらはSciPy 2.0.0で削除されるそうなので，代わりにNumpyのstd, sqrt関数を使用しても同じ結果が得られます）．
>>> sigma = sp.std(junk_food, ddof = 1) # sigma = np.std(junk_food, ddof = 1)
>>> se = sigma / sp.sqrt(len(junk_food)) # se = sigma /
>>> se
1.957927680575589
最後にt値を計算します．
>>> t_value = (mu - 50) / se
>>> t_value
2.750339683171343
p値の計算
母集団に正規分布を仮定すると，t値はt分布に従うと考えられるので，t分布の累積分布関数が使えます．
p値の計算方法は，t分布の累積分布関数を使うと，母平均を50と仮定した時に，t値が$t_{標本}$を下回る確率を計算することができ，この確率をaとします．1 - a を求めると，母集団を50と仮定した時に，t値が$t_{標本}$を上回る確率を計算できます．1-aが小さければ，t値が$t_{標本}$を上回る確率が高い（$t_{標本}$は十分に大きい）ことになり有意差が得られそうです．
両側検定を実施した場合は，p値は$(1 - a) \times 2$として計算されます．
>>> alpha = stats.t.cdf(t_value, df = df)
>>> (1 - alpha) * 2
0.012725590012524046
p値が有意水準5%（0.05）を下回っているので，有意差ありとみなせます（この場合は，ジャンクフードの重量が平均重量50gと有意に異なっていると判断することができます．
stats.ttest_1samp を使うと，以下のように簡単に1標本のt検定を行うことができます．
>>> stats.ttest_1samp(junk_food, 50)
Ttest_1sampResult(statistic=2.750339683171343, pvalue=0.012725590012524182)
この例でstatisticsはt値で，pvalueはp値です．
p値は，帰無仮説が正しいと仮定して，何度も標本抽出〜t値計算を繰り返した時に，$t_{標本}$と同じかそれより大きなt値が得られる割合と解釈できます．両側検定の場合は，この割合を2倍したものがp値となります．この割合が小さいということは$t_{標本}$を超えることがめったにないということになり，$t_{標本}$が十分に大きいと考えられ，有意さが得られると判断することができます．
以下に，p値をシミュレーションで求めてみます．まずは，標本の情報（サンプルサイズと標準偏差）を変数に格納します．
>>> size = len(junk_food)
>>> sigma = sp.std(junk_food, ddof = 1) # sigma = sp.std(junk_food, ddof = 1)
シミュレーション（50000回t値を計算）を行なってみます．t_value_arrayは，50000個のt値を格納する変数です．stats.normに指定する平均値を50とすることで，帰無仮説が正しいと仮定して，50000回，標本抽出〜t値計算を繰り返します．そして，50000個のt値のうち，$t_{標本}$を上回った割合を求めて，これに2をかけるとp値となります．
>>> t_value_array = np.zeros(50000)
>>>
>>> np.random.seed(1)
>>> norm_dist = stats.norm(loc = 50, scale = sigma)
>>> for i in range(0, 50000):
... sample = norm_dist.rvs(size = size)
... sample_mean = sp.mean(sample) # sample_mean = np.mean(sample)
... sample_std = sp.std(sample, ddof = 1) # sample_std = np.std(sample, ddof = 1)
... sample_se = sample_std / sp.sqrt(size) # sample_se = sample_std / np.sqrt(size)
... t_value_array[i] = (sample_mean - 50) / sample_se
...
>>>
>>> (sum(t_value_array > t_value) / 50000) * 2
0.01324
0 件のコメント :
コメントを投稿