[Python]とにかくわかりやすく!ハンバーガー統計をやってみる③

#前回の記事
前回の記事→[こちら](https://qiita.com/hiroyuki_mrp/items/f0cb8d3b94c08588f032)

#環境
Mac OSX　mojave
Python３.7

#とりあえず超超有名な入門書
ハンバーガー統計をPythonで書いてみたいと思います。
[元はこちらの内容](https://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E5%AD%A6%E3%81%8C%E3%82%8F%E3%81%8B%E3%82%8B-%E3%83%95%E3%82%A1%E3%83%BC%E3%82%B9%E3%83%88%E3%83%96%E3%83%83%E3%82%AF-%E5%90%91%E5%BE%8C-%E5%8D%83%E6%98%A5/dp/4774131903)

分量が多くなるので
その①：平均と分散、信頼区間
その②：カイ二乗検定
その③：t検定 ←本記事はここ
その④：分散分析
のようなシリーズで行きたいと思います。

#t検定（対応なし）

ワクワクバーガーとモグモグバーガーの味について差があるのかを調べます。今回は通りすがった女子高生に**交互**に点数をつけてもらいます。平均点が5点あったが、これが意味のある差なのかどうかを考えます。

ここでは、ワクワクとモグモグの平均の差の信頼区間があるのかみてみます。

データを準備します。

<img width="292" alt="スクリーンショット 2019-04-10 4.51.00.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241505/6a30b054-5823-de0c-c7ac-2db67e71f893.png">

ライブラリをインポートします。してあれば不要です。

```python
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を引数にいれました。

```python
pointdata = pd.read_csv("wakumogupoint.csv",sep ="\t")
```
店舗ごとの売上を出してarrayにします。

```python
wakupoint = pointdata.query('stores == "waku"')["points"]
mogupoint = pointdata.query('stores == "mogu"')["points"]
wakupoint = np.array(wakupoint)
mogupoint = np.array(mogupoint)
```
計算に使いそうな統計量を出しておきます。

```python
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{差の標準誤差}$で計算を行います。

差の標準誤差ですがそれぞれの誤差は独立しているものなので、加算をします。

``` math
差の標準誤差 = \sqrt{\frac{不偏分散A}{NA}+\frac{不偏分散B}{NB}}
```

```python
dif_stder = sp.sqrt(waku_var/waku_N +mogu_var/mogu_N)
```

としたいところだが、
AとBの母分散が等しいと仮定して**推定母分散**というものを使います。
※それぞれの分散が等しいと仮定せずに行なっても問題ありません。（＝welchの検定でいい）

まずは平方和を出します。通常の不偏分散の算出と考え方は一緒です。

```math
推定母分散 = \frac{標本Aの偏差平方和＋標本Bの偏差平方和}{(NA−1)+(NB-1)}
```

```python
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
```
それでは推定母分散を元に、差の標準誤差を算出してみます。※式はどれでもいいです

``` math
差の標準誤差 = \sqrt{\frac{推定母分散}{NA}+\frac{推定母分散}{NB}}\\
差の標準誤差 = \sqrt{推定母分散\times\frac{1}{NA}+推定母分散\times\frac{1}{NB}}\\
差の標準誤差 = \sqrt{推定母分散\times\bigg(\frac{1}{NA}+\frac{1}{NB}\bigg)}

```

```python
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
```

あと信頼区間なので自由度が必要です。

```math
自由度 = （NA-1）\times(NB-1)
```

```python
df = (waku_N-1) + (mogu_N-1)
print(df)
#14
```
※14の時の自由度は2.145です

では信頼区間を計算してみます。

```math
差の信頼区間 = (平均A-平均B)\pm{t}\times{差の標準誤差}
```

```python
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値を求めます。

```math
t = 標本平均の差 \div 標本平均の差の標準誤差
```

```python
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値がでてます。

```python
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検定（対応あり）

今度は同じ人にワクワクバーガー、モグモグバーガー食べてもらって評点をつけてもらいます。
つまり独立した２つの母集団ではなく、**点数の差**自体を１つのものとして考えるのです。
さらにつまり、点数の差自体に意味があるのかどうかを考えればいいということです。

データを作ります。

<img width="295" alt="スクリーンショット 2019-04-14 15.54.49.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241505/9679c36b-c126-df3b-ec22-ac6e3a089d40.png">

データの前準備をしておきます。

```python
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}}）$

差のデータを作り、統計量を計算しておきます。

```python
dif = wakupoint-mogupoint
dif_mean = sp.mean(dif)
dif_var = sp.var(dif,ddof=1)
print(dif_var)
#17.410714285714285
``` 
サンプル数を計算します。ここでの自由度は差のN数なので以下のように計算します。２パターン。

```python
N = len(set(np.array(pointdata["num"])))　#回答者IDの重複を削除したUU数を数える
N = len(dif)　#差のarrayの数を数える
print(N)
#8
```
差の標準偏差を計算します。

```python
dif_stder = sp.sqrt(dif_var/N)
print(dif_stder)
#1.4752421108802058
```

あとはt値を計算します。

```python
t = (waku_mean-mogu_mean)/dif_stder
print(t)
#-2.965614910077132
```
ここまでの計算は一瞬でできます。差のarrayが必要です。
差の平均が０と異なるかとどうかをみます。

```python
result1 = stats.ttest_1samp(dif,0)
print("t値:"+str(result1.statistic))
print("p値:"+str(result1.pvalue))
#t値:-2.965614910077132
#p値:0.020937570206924612
```
さらにもっと簡単にできます。サンプルを２つ引数にとるだけです。

```python
result2 = stats.ttest_rel(wakupoint,mogupoint)
print("t値:"+str(result2.statistic))
print("p値:"+str(result2.pvalue))
t値:-2.965614910077132
p値:0.020937570206924612
```

#終わり

t検定をやると、統計っぽさが一気に出てきます。次回は分散分析について投稿します。