#前回の記事
前回の記事→こちら
#環境
Mac OSX mojave
Python３.7
#とりあえず超超有名な入門書
ハンバーガー統計をPythonで書いてみたいと思います。
元はこちらの内容
分量が多くなるので
その①：平均と分散、信頼区間
その②：カイ二乗検定
その③：t検定 ←本記事はここ
その④：分散分析
のようなシリーズで行きたいと思います。
#t検定（対応なし）
ワクワクバーガーとモグモグバーガーの味について差があるのかを調べます。今回は通りすがった女子高生に交互に点数をつけてもらいます。平均点が5点あったが、これが意味のある差なのかどうかを考えます。
ここでは、ワクワクとモグモグの平均の差の信頼区間があるのかみてみます。
データを準備します。
ライブラリをインポートします。してあれば不要です。
import pandas as pd
import numpy as np
import scipy as sp
from scipy import stats
import matplotlib.pyplot as plt
%matplotlib inline
%precision 3
import seaborn as sns
sns.set()
データを読み込みます。タブ区切りになっていたので、sepを引数にいれました。
pointdata = pd.read_csv("wakumogupoint.csv",sep ="\t")
店舗ごとの売上を出してarrayにします。
wakupoint = pointdata.query('stores == "waku"')["points"]
mogupoint = pointdata.query('stores == "mogu"')["points"]
wakupoint = np.array(wakupoint)
mogupoint = np.array(mogupoint)
計算に使いそうな統計量を出しておきます。
waku_N = len(wakupoint)
mogu_N = len(mogupoint)
waku_mean = sp.mean(wakupoint)
mogu_mean = sp.mean(mogupoint)
waku_var = sp.var(wakupoint,ddof=1)
mogu_var = sp.var(mogupoint,ddof=1)
さてt検定を行っていくにあたって、考えるのが信頼区間です。
信頼区間は
$標本平均\pm{t}\times{標準誤差}$
というの計算式でしたが、今回は差がテーマなので
$(標本平均A-標本平均B)\pm{t}\times{差の標準誤差}$で計算を行います。
差の標準誤差ですがそれぞれの誤差は独立しているものなので、加算をします。
差の標準誤差 = \sqrt{\frac{不偏分散A}{NA}+\frac{不偏分散B}{NB}}
dif_stder = sp.sqrt(waku_var/waku_N +mogu_var/mogu_N)
としたいところだが、
AとBの母分散が等しいと仮定して推定母分散というものを使います。
※それぞれの分散が等しいと仮定せずに行なっても問題ありません。（＝welchの検定でいい）
まずは平方和を出します。通常の不偏分散の算出と考え方は一緒です。
推定母分散 = \frac{標本Aの偏差平方和＋標本Bの偏差平方和}{(NA−1)+(NB-1)}
sumsqu_waku = sp.sum((wakupoint-waku_mean)**2)
sumsqu_mogu = sp.sum((mogupoint-mogu_mean)**2)
espopu_var = (sumsqu_waku+sumsqu_mogu)/((waku_N-1)+(mogu_N-1))
print(espopu_var)
#60.267857142857146
それでは推定母分散を元に、差の標準誤差を算出してみます。※式はどれでもいいです
差の標準誤差 = \sqrt{\frac{推定母分散}{NA}+\frac{推定母分散}{NB}}\\
差の標準誤差 = \sqrt{推定母分散\times\frac{1}{NA}+推定母分散\times\frac{1}{NB}}\\
差の標準誤差 = \sqrt{推定母分散\times\bigg(\frac{1}{NA}+\frac{1}{NB}\bigg)}
dif_stder = sp.sqrt(espopu_var/(waku_N)+espopu_var/(mogu_N))
dif_stder = sp.sqrt(espopu_var*(1/waku_N)+espopu_var*(1/mogu_N))
dif_stder = sp.sqrt(espopu_var*(1/waku_N+1/mogu_N))
print(dif_stder)
#3.8816187713007424
あと信頼区間なので自由度が必要です。
自由度 = （NA-1）\times(NB-1)
df = (waku_N-1) + (mogu_N-1)
print(df)
#14
※14の時の自由度は2.145です
では信頼区間を計算してみます。
差の信頼区間 = (平均A-平均B)\pm{t}\times{差の標準誤差}
min_kukan = waku_mean-mogu_mean-2.145*dif_stder
max_kukan = waku_mean-mogu_mean+2.145*dif_stder
print(str(min_kukan)+"~"+str(max_kukan))
#-13.326072264440093~3.326072264440093
つまり今回(標本)の平均点の差は5点だったが、母集団の平均の差は-13.326〜3.326であり、その範囲には0が含まれる（差がない）ため、有意差は認められない。となりました。
とここまでのことは信頼区間からの解釈でしかないので、最後にt検定を行います。
まずt値を求めます。
t = 標本平均の差 \div 標本平均の差の標準誤差
t = (waku_mean-mogu_mean)/dif_stder
print(t)
#-1.2881223774390613
先ほどの通り,自由度が14の時のtは２.１４５でした。
つまりt値が、２.１４５より大きい、−２.１４５より小さいというのは（正規分布のためプラスマイナス）、5%未満の確率でしかおこらないことであると解釈でき、そこ（棄却域）に入るときは帰無仮説を棄却する。
上記の通り今回は−1.288なので帰無仮説を採用する。
ちなみに
自由度が14の時の有意水準1%の時は、tは２.977です。
同様にt値が、２.977より大きい、−２.977より小さいというのは、1%未満の確率でしかおこらないことであると解釈でき、そこ（棄却域）に入るときは帰無仮説を棄却する。
ここまでのt検定も一発でできます。
上は等分散を仮定したもの、下がwelchの検定になります。
t値とP値がでてます。
result1 = stats.ttest_ind(wakupoint, mogupoint,equal_var = True)
print("t値:"+str(result1.statistic))
print("p値:"+str(result1.pvalue))
#t値:-1.299867367239363
#p値:0.21464181738984914
result2 = stats.ttest_ind(wakupoint, mogupoint,equal_var = False)
print("t値:"+str(result2.statistic))
print("p値:"+str(result2.pvalue))
#t値:-1.299867367239363
#p値:0.21482211719824498
#t検定（対応あり）
今度は同じ人にワクワクバーガー、モグモグバーガー食べてもらって評点をつけてもらいます。
つまり独立した２つの母集団ではなく、点数の差自体を１つのものとして考えるのです。
さらにつまり、点数の差自体に意味があるのかどうかを考えればいいということです。
データを作ります。
データの前準備をしておきます。
pointdata = pd.read_csv("wakumogupoint2.csv",sep ="\t")
wakupoint = pointdata.query('stores == "waku"')["points"]
mogupoint = pointdata.query('stores == "mogu"')["points"]
wakupoint = np.array(wakupoint )
mogupoint = np.array(mogupoint )
waku_N = len(wakupoint)
mogu_N = len(mogupoint)
waku_mean = sp.mean(wakupoint)
mogu_mean = sp.mean(mogupoint)
waku_var = sp.var(wakupoint,ddof=1)
mogu_var = sp.var(mogupoint,ddof=1)
ここからが「対応なし」と違うところです。
対応なし：
$t = 標本平均の差\div標本平均の差の標準誤差$
対応あり：
$t = 標本平均の差（差の平均）\div差の標準誤差（\sqrt{不偏分散\div{N}}）$
差のデータを作り、統計量を計算しておきます。
dif = wakupoint-mogupoint
dif_mean = sp.mean(dif)
dif_var = sp.var(dif,ddof=1)
print(dif_var)
#17.410714285714285
サンプル数を計算します。ここでの自由度は差のN数なので以下のように計算します。２パターン。
N = len(set(np.array(pointdata["num"]))) #回答者IDの重複を削除したUU数を数える
N = len(dif) #差のarrayの数を数える
print(N)
#8
差の標準偏差を計算します。
dif_stder = sp.sqrt(dif_var/N)
print(dif_stder)
#1.4752421108802058
あとはt値を計算します。
t = (waku_mean-mogu_mean)/dif_stder
print(t)
#-2.965614910077132
ここまでの計算は一瞬でできます。差のarrayが必要です。
差の平均が０と異なるかとどうかをみます。
result1 = stats.ttest_1samp(dif,0)
print("t値:"+str(result1.statistic))
print("p値:"+str(result1.pvalue))
#t値:-2.965614910077132
#p値:0.020937570206924612
さらにもっと簡単にできます。サンプルを２つ引数にとるだけです。
result2 = stats.ttest_rel(wakupoint,mogupoint)
print("t値:"+str(result2.statistic))
print("p値:"+str(result2.pvalue))
t値:-2.965614910077132
p値:0.020937570206924612
#終わり
t検定をやると、統計っぽさが一気に出てきます。次回は分散分析について投稿します。