統計検定をpythonにやってもらおう#1

初めに

データ解析とかビックデータって近未来感あって学生の僕にはかっこよく夢のある分野なんです。ですよね？笑
でも学ぶとなるとプログラミングも統計学なんかの数学知識も学ばなきゃいけないし独学だとくじけることもしばしば、、
そこで今回は、統計学の中でも基本となる仮説検定の理解を深めるために簡単な正規分布から、t検定あたりまでの基本をpythonによって行った。そう、同時に学んでやろうって

想定読者としてはpythonは文法程度、統計学も大学でかじったよくらいの理解度の初学者で、データサイエンスなどに触れてみたい学生等のちょっとしたヒントになればうれしいです、、！
前提として、pythonの使える環境を想定していますが環境構築なんかの話は僕よりわかりやすく説明してくれている記事がたくさんあるので調べてみてください。

目次
１. グラフの作り方
２. 正規分布の作成とデータ比較
３. 仮説検定の基本(SciPy)
４. t分布の作り方
５. まとめ、次回について

１. グラフの作り方

ここでは、統計学をpythonで学ぶために必須といえるグラフの作成を簡単に紹介しておきます。詳しく知りたい方は書籍やサイトを調べれば簡単に調べられます。

例として、0~5までの0.2間隔ずつの数値を2乗、3乗した結果をグラフにしてみよう。
グラフの作成にはmatplotlibを用いています。

グラフの作り方

import numpy as np
import matplotlib.pyplot as plt

# シンプルなグラフ
t = np.arange(0.,5.,0.2) # 0から5まで0.2ごとに数を生成しリストにぶち込む
plt.title('drawing example1')  # タイトル
plt.plot(t,t,'r--',label='linear')  # tの値をそのままプロット
plt.plot(t,t**2,'bs',label='square')  # tの2乗の値をプロット
plt.plot(t,t**3,'g',label='cube')  # tの3乗の値をプロット
plt.xlabel('x')  # それぞれの軸に名前を付ける
plt.ylabel('y')
plt.legend()  # 凡例
plt.show()   # 表示

出力結果

累乗が大きいほど値が大きくなることがグラフから読み取れます。もちろんグラフのサイズや色なども自由に変えられるので興味のある方は遊んでみて、、！笑

２. 正規分布の作成とデータ比較

次に正規分布の作り方をざっくりやっていきます。今回はサンプルデータとして平成28年度の17歳男性の身長データを用いた。h28学校保健統計調査よりcsv形式で保存してください。
やってみたい方は自分で好きな年度(自分の年のとか？)でやってみるのもよいかも

注意としては年度によってデータの配置が異なるので抽出箇所をもとのデータを見て調整しておきましょう。コード内だとfilenameに自分でダウンロードしたパスを入れ、[61:110[0,13]]の部分が取り出したい範囲なのでここを変えれば好きなデータで確認できます。

正規分布と統計データの比較

from scipy.stats import norm
import pandas as pd
from matplotlib import pyplot as plt

# u=170.7 ,σ=5.86
filename = "C:\python_study\python_excel\pandas_data_analysis\input_data\h28_hoken_toukei_02_men.csv"
df = pd.read_csv(filename)  # データの読み込み
df2 = df.iloc[61:110,[0,13]].astype(float)  # 抽出範囲の指定
df2 = df2.rename(columns={'Unnamed: 0':'height','Unnamed: 13':'permil'})  # 取り出した範囲の列に変数名(height,permil)を入れる→グラフ作成時に代入できるように

u = [norm.pdf(x=i,loc=170.7,scale= 5.81)*1000.0 for i in df2['height']]    # 正規分布関数値を作成（データの平均、偏差から作成し比較対象としている)

df2['norm']= u

ax = df2.plot.scatter(x = 'height',y='permil',color='black',marker = 'x')

df2.plot(x='height',y ='norm',color='black',label='N(170.7,5.81)',ax=ax)
plt.grid()
plt.xlabel('tall')
plt.ylabel('‰')
plt.legend(loc="upper left")
plt.show()

出力結果

サンプルデータはおおむね正規分布に従っていることがグラフから読み取れました。
これをバスケ部の人だけで調査したら175くらいに上がるのかな、、

というわけで、やってみるとpythonで正規分布を作るのはコードではたった一行であることがわかった。scipyってすごく便利

ちなみに、正規分布は標準正規分布をよく使うと思うので、この身長の分布を標準化するには
x' = (x-μ)/σ = (x-170.7)/5.81
の変換を行えばいいです。

３. 仮説検定の基本(SciPy)

初めにscipyについて簡単に説明しておきます。
推測統計では、SciPyというパッケージをよく使います。さっきの正規分布でも用いましたね。
SciPyの使用にあたっては、scipy.stats.xxx.oooのように使用するメソッドの指定を行い、xxxに確率分布、oooに行う処理を指定する。

推測統計ではxxxに対して主に下記などの確率分布の指定を行う。

norm: 正規分布
t:t 分布
chi2: χ2分布

同様に推測統計ではoooに対して、主に下記の処理の指定を行う。

cdf: 統計量から有意水準・値の計算
ppf: 有意水準・値から統計量の計算
pdf: 統計量から確率密度関数の計算

特にcdfとppfは区間推定や検定で統計数値表の代わりに用いることができる。

というわけで使い方の確認ついでに、正規分布のpdfとcdfを作ってみましょう。

pdfとcdf

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

x = np.arange(-3,3.01,0.01)
f_pdf = stats.norm.pdf(x)
F_cdf = stats.norm.cdf(x)

plt.plot(x,f_pdf,label="pdf")
plt.plot(x,F_cdf,label="cdf")

出力結果

確率密度関数は、正規分布の積分なので1に近づいていくことがわかる。

次に、確率密度関数からどんな結果が得られるか見てみよう。

確率密度関数

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

alpha = np.array([0.1, 0.05, 0.025, 0.01, 0.005]) # 複数の確率で行ってみる
x = np.arange(-3.,3.01,0.01)
F_cdf = stats.norm.cdf(x)

plt.plot(x,F_cdf,label="cdf")
for i in range(alpha.shape[0]):
    print("alpha: {}, statistic_upper: {:.2f}".format(alpha[i], stats.norm.ppf(1.-alpha[i])))  # 結果を表示
    x_upper = np.repeat(stats.norm.ppf(1.-alpha[i]), 100)
    y = np.linspace(0,1.-alpha[i],100)
    plt.plot(x_upper,y,label="alpha: {}".format(alpha[i]))

plt.legend(loc=2)

出力結果

alpha: 0.1, statistic_upper: 1.28
alpha: 0.05, statistic_upper: 1.64
alpha: 0.025, statistic_upper: 1.96
alpha: 0.01, statistic_upper: 2.33
alpha: 0.005, statistic_upper: 2.58

結果を見ると、例えば
alpha(確率)：0.1,statistic_upper(上限確率値):1.28　であれば、-3から3までの0.01刻みのランダムな数値の分布において1.28以上の数値が出る確率が10％であることがわかる。

後に区間推定なんかを触れていくのでめっちゃざっくり活用法をみておくと、信頼区間95％から外れる点がalpha:0.025(両側5%)の点で1.96を実測値が超えてれば帰無仮説を棄却”！みたいな感じですね(次回以降詳しくやっていきます!)

４. t分布の作り方

やっと統計学！！という感じがしてきましたが、もともと統計をしっかりしてる人からみたらまだ入り口ですけどって感じでしょうか、、？笑

というわけで、t分布についてですが
そもそもどんな時に使うかというと、母集団が大きすぎて母平均がわからない！みたいなときです(ざっくりすぎでしょうか笑)
こういう時は、何個か標本を取ってきて標本平均を代わりに使えば何とか推定できないかというわけです。
しかし、その場合取ってくるサンプルサイズによって偏りが生じます。自由度ってやつです。詳しくやってると日が暮れそうなので省略します、、詳しくはこちらのサイトさんがわかりやすかったのでぜひ (t分布と自由度)

それでは、自由度が2,20,100のときで分布がどうなるか見ていきましょう

t分布

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st

dfs = [2,20,100]
mk =['.','x','o']
x = np.arange(-6,6,0.1)

plt.figure(figsize=(10,5))
for i,df in enumerate(dfs):
  y = st.t.pdf(x,df)
  plt.plot(x,y,color = 'black',label = 'df='+str(df),marker=mk[i])
  plt.xlabel('x')
  plt.ylabel('y')
  plt.legend()
  plt.grid(color ='black')
plt.show()

出力結果

結果を見ると、自由度が2と20では結構変わりますが20と100では差がかなり少なくなっていき正規分布に近づくことがわかります。
つまり、データが複雑な要素の絡みがないシンプルなものであれば、標本を100くらい用意しておけば十分有意な結果が得られるというわけです。なんか思ったより少なくていいという感じ笑

ちなみにt分布の作成はst.t.pdf(x,df)によって作れるのでこれまた簡単に作れます。

５. まとめと次回について

最後まで読んでいただきありがとうございました！
統計学の基礎の基礎をpythonを使って簡単に解説しましたがどうでしたでしょうか？
かなり端折ったり雑な説明になってしまっている個所もあるとは思いますが、まずは触れてみるのがだいじだと思っているので、ぜひ自分でコードをカスタムして使ってみてください！！
次回は、今回の知識を使って区間推定について書こうと考えています。次回からはデータを分析してるって思える内容になると思います多分、、( ´∀｀ )

なお使うデータなどをホームページかなにかにまとめておきたいですが，作り方がわかんないので勉強します泣

参考：データサイエンスのための統計学入門