第2回/ハンズオンの履歴(No.7)

第2回ハンズオン†

第2回ハンズオン

0. 準備†

Google Driveにアクセスして，Google Colabを新規作成してください
名前を handson2.ipynb とします

下記のコードを貼り付けて実行してください

# PandasとNumpyをインポート
import pandas as pd
import numpy as np

# 日本語化Matplotlibもインポート
import matplotlib.pyplot as plt
!pip install japanize-matplotlib
import japanize_matplotlib

#Scipy statsをインポート
from scipy import stats

↑

1. 第2回で出てきた確率分布†

↑

scipy.stats の書式†

標準正規分布： stats.norm(loc=0, scale=1)
t分布： stats.t(df=自由度)
χ二乗分布： stats.chi2(df=自由度)
F分布： stats.f(dfn=自由度f1, dfd=自由度f2)

↑

標準正規分布の概形†

Y = Z = N(0,1)

#xの範囲を適当に決める (-5～5の範囲を100等分した数列)
x=np.linspace(-5,5,100)
# 標準正規分布を定義
norm = stats.norm(loc=0, scale=1)

#グラフに描く
plt.plot(x, norm.pdf(x), label="N(0,1)")
plt.legend()
plt.grid()

↑

t分布の概形†

T = Z / sqrt(Y/f)
- 標準正規分布をカイ二乗分布(自由度で平均化)の平方根で割った分布

#xの範囲を適当に決める (-5～5の範囲を100等分した数列)
x=np.linspace(-5,5,100)

#自由度fを1～10にして確率密度関数を描画する
for f in range(10):
  t = stats.t(df=f)
  plt.plot(x, t.pdf(x), label=f"t({f})")

plt.grid()
plt.legend()

↑

カイ二乗分布の概形†

Y = Z_1^2 + Z_2^2 + ... + Z_f^2
- 標準正規分布の2乗を任意の数（自由度）だけ足し合わせた分布

#xの範囲を適当に決める (0～20の範囲を100等分した数列)
x=np.linspace(0,20,100)

#自由度fを1～10にして確率密度関数を描画する
for f in range(10):
  chi2 = stats.chi2(df=f)
  plt.plot(x, chi2.pdf(x), label=f"χ^2({f})")

plt.grid()
plt.legend()

↑

F分布†

X = (Y1/f1) / (Y2/f2)
- 2つのカイ二乗分布（自由度f1, f2）の比を取った分布

#xの範囲を適当に決める (0～5の範囲を100等分した数列)
x=np.linspace(0,5,100)

#自由度f1,f2を適当に決めて確率密度関数を描画する
for f1 in [10,20,30]:
  for f2 in [5,50]:
    F = stats.f(dfn=f1, dfd=f2)
    plt.plot(x, F.pdf(x), label=f"F({f1}, {f2})")
plt.xlim(0,3)
plt.grid()
plt.legend()

↑

2. 仮説検定の基礎†

出てきた2つの例題をPythonで試してみよう

↑

例1: このコインで行うコイントスは正当か？†

p=1/2のコインを10回投げる場合の確率は，2項分布 stats.binom(p=0.5, n=10)で計算可能

#コインを投げる回数
N=10
#2項分布を使う (コインが表裏同じ確率1/2で出る場合)
binom = stats.binom(p=0.5, n=N)
print(f"コインを{N}枚投げて i枚表が出る確率を計算します")
for i in range(N+1):
  print(f"{i} 枚が表の確率： {binom.pmf(i)}")

Nが大きい場合でも，同様に計算できる

講義であった，中心極限定理を試してみる

#ゆがんだコインの母集団を作る
all_coins = stats.bernoulli(p=0.7)

#n枚取り出してコイントス
n = 30
coins = all_coins.rvs(size=n)
print(f"【コイントス結果】 {coins} -- {n}枚中{coins.sum()}枚表")

#コインが公正である(p0=0.5)とした場合の検定量を計算する
X_bar = coins.mean()
p0 = 0.5
Z = (X_bar - p0) / np.sqrt(p0*(1-p0)/n)
print(f"【検定量】 {Z}")

#有意水準α=0.05として検定量Zを検定してみる
alpha = 0.05
print(f"【alphaの値】{alpha}")

#標準正規分布
norm = stats.norm(loc=0, scale=1)

#N(0,1)の上側下側 alpha/2 点を調べる
pos_a_upper = norm.isf(alpha/2)
pos_a_lower = norm.ppf(alpha/2)
print(f"【alpha/2の位置】{pos_a_lower}, {pos_a_upper}")

#検定
if ((Z > pos_a_upper) or (Z < pos_a_lower)):
  print(f"【検定】Z={Z} がalphaの位置を超えるため，p={p0} を棄却します")
else:
  print(f"【検定】Z={Z} がalphaの位置を超えないため，p={p0} を保留します")

#p値の計算 N(1,0)(X>=Z) * 2 (両側) を計算すればよい
p_value= norm.sf(Z) * 2
print(f"【p値】{p_value}")
 
#正規分布をプロット
x=np.linspace(-5,5,100)
plt.plot(x, norm.pdf(x))
plt.ylim(0,0.5)

#alphaの位置をプロット（赤線）
plt.axvline(x=pos_a_upper, c="red")
plt.axvline(x=pos_a_lower, c="red")

#Zの位置をプロット（緑線）
plt.axvline(x=Z, c="green")
plt.show()

Pythonの2項検定 stats.binomtest(k=成功した数, n=試行回数, p=成功確率)を使えば，2項分布のまま1行で検定可能．

stats.binomtest(k=21, n=30, p=0.5)

↑

例2：日本人男性の平均身長は160cmか？†

母集団が正規分布であることを仮定し，Studentのt分布を使って検定する（t-検定）

#母集団．平均を160cmから少しずらしておく
all_boys = stats.norm(loc=163, scale=5)

#n人無作為抽出
n=20
boys = all_boys.rvs(size=n)
print(f"【{n}人の身長サンプル】\n{boys}")

#帰無仮説の平均
mu0 = 160
#t検定量
T = (boys.mean()-mu0)/np.sqrt(boys.var(ddof=1)/n)
print(f"【検定量】: T={T}")

#有意水準
alpha = 0.05
print(f"【alphaの値】{alpha}")

#自由度n-1のT分布
tdist = stats.t(df=n-1)

#t(n-1)の上側下側 alpha/2 点を調べる
pos_a_upper = tdist.isf(alpha/2)
pos_a_lower = tdist.ppf(alpha/2)
print(f"【alpha/2の位置】{pos_a_lower}, {pos_a_upper}")

#検定
if ((T > pos_a_upper) or (T < pos_a_lower)):
  print(f"【検定】T={T} がalphaの位置を超えるため，mu0={mu0} を棄却します")
else:
  print(f"【検定】T={T} がalphaの位置を超えないため，mu0={mu0} を保留します") 

#p値の計算 t(n-1)(X>=T) * 2 (両側) を計算すればよい
p_value= tdist.sf(T) * 2
print(f"【p値】{p_value}")

#プロット
x = np.linspace(-5,5,100)
plt.plot(x,tdist.pdf(x))
plt.axvline(x=pos_a_upper ,c="red")
plt.axvline(x=pos_a_lower ,c="red")
plt.axvline(x=T,c="green")
plt.grid()

上記の検定は，1変数のt検定 stats.ttest_1samp(標本, 平均値) で1行で書ける

stats.ttest_1samp(boys,160)

↑

3. 様々な仮説検定†

↑

正規性判定†

Shapiro-Wilk検定

stats.shapiro(boys)

Kolmogorov-Smirnov検定

stats.kstest(boys, "norm")

↑

Q-Qプロット†

サンプルした男性の身長 boys が正規分布に沿っているか？

stats.probplot(boys, dist="norm", plot=plt)

↑

第2回/ハンズオン の履歴(No.7)