#author("2023-10-09T21:09:00+09:00","default:cmdsadmin","cmdsadmin")
[[第2回]]
#author("2024-05-07T10:36:42+09:00","default:cmdsadmin","cmdsadmin")

* 第2回ハンズオン [#r15dffde]

#contents

** 0. 準備 [#fd9a7c5f]
- [[Google Drive:https://drive.google.com/]]にアクセスして，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回で出てきた確率分布 [#ndd02cea]
*** scipy.stats の書式 [#u9a47ff5]
- 標準正規分布： stats.norm(loc=0, scale=1)
- t分布： stats.t(df=自由度)
- χ二乗分布： stats.chi2(df=自由度)
- F分布： stats.f(dfn=自由度f1, dfd=自由度f2)

*** 標準正規分布の概形 [#pf2c413c]
- 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分布の概形 [#p084ccfe]
- 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()


*** カイ二乗分布の概形 [#z4767d36]
- 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分布 [#z4bd240f]
- 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. 仮説検定の基礎 [#m875ff71]
講義に出てきた2つの例題をPythonで試してみよう

*** 例1: このコインで行うコイントスは正当か？ [#s47f7961]

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(0,1)(X>=Z) + N(0,1)(X<=-Z) (両側) を計算すればよい
 absZ = np.abs(Z) #絶対値をとる
 p_value= norm.sf(absZ) + norm.cdf(-absZ) 
 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()


上記は標準正規分布を使った近似検定だが，scipy.statsの''2項検定''を使えば，正確な2項分布を使って検定できる
- stats.binomtest(k=成功数, n=総試行数, p=想定する母平均)
 stats.binomtest(k=21, n=30, p=0.5)

*** 例2：日本人男性の平均身長は160cmか？ [#y798484d]

母集団が正規分布であることを仮定し，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 = (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) + t(n-1)(X<=-T) を計算すればよい
 absT = np.abs(T) #絶対値を取る
 p_value= tdist.sf(absT) + tdist.cdf(-absT)
 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()

上記の検定は，scipy.statsの''1変数のt検定''で1行で書ける
- stats.ttest_1samp(標本, popmean=母平均の値)
 stats.ttest_1samp(boys,popmean=160)

** 3. 各種検定手法 [#i4e8e575]

*** Q-Qプロット [#hc6b94b2]
サンプルした男性の身長 boys が正規分布に沿っているか？
 stats.probplot(boys, dist="norm", plot=plt)

*** 正規性判定 [#w855c85a]
Shapiro-Wilk検定
 stats.shapiro(boys)

Kolmogorov-Smirnov検定
 stats.kstest(boys, "norm")




*** ブートストラップ検定 [#wbac2364]

【例題】ある工場で製造される電球は，平均寿命が3000時間であることを公表している．

いま，完成品からランダムに50個の電球を取り出し，寿命を調べたところ，以下の通りであった．

,1167,7295,1154, 537,1498, 327,1372,1809,2970, 844
,5604,2075,9119,3212,1634, 309,  28,2852, 804, 273
,5075, 373, 457,2701,4451, 836, 474, 391,1902,1926
,2210,1545, 445, 408, 299, 393,1964,1554, 595,5680
,3777,7511, 783,9516, 817,4359, 207,6039, 609,3069

データより50個の平均値は2305時間であり，一見3000時間より短いように見える．
工場内のすべての電球の寿命が3000時間より短いかどうかを仮説検定しなさい．

- H0: μ = 3000
- H1: μ < 3000
- α=0.05

 #50個の電球の故障時間
 data = [1167,7295,1154, 537,1498, 327,1372,1809,2970, 844,
         5604,2075,9119,3212,1634, 309,  28,2852, 804, 273,
         5075, 373, 457,2701,4451, 836, 474, 391,1902,1926,
         2210,1545, 445, 408, 299, 393,1964,1554, 595,5680,
         3777,7511, 783,9516, 817,4359, 207,6039, 609,3069]
 
 #配列に入れる
 samples = np.array(data)
 print(f"【初期標本】 {samples}")
 N=len(samples)
 print(f"【標本数】 {N}")
 print(f"【初期標本の平均】 {samples.mean()}")
 
 #Q-Qプロットを描いてみる
 stats.probplot(samples, dist="norm", plot=plt)

 #ブートストラップ標本を作る
 #離散一様分布 (0～N-1の一様な分布)
 uni = stats.randint(0,N)
 #復元抽出する元標本のインデクスを作る
 idx = uni.rvs(size=N)
 #わかりやすいようにソート
 idx = np.sort(idx)
 print(f"【インデクス】{idx}")
 #待ち時間の標本からリサンプリングする
 resamples = samples[idx]
 print(f"【初期標本】{samples}")
 print(f"【ブートストラップ標本】{resamples}")


 #ブートストラップ検定を行う
 
 #反復回数B
 B=5000
 #帰無仮説H0の母平均μ0
 mu0 = 3000
 
 #有意水準
 alpha = 0.05
 
 #ブートストラップ標本をいれる辞書
 b_samples = {}
 
 #初期標本XからT(X)を計算する
 T = (samples.mean() - mu0) / np.sqrt(samples.var(ddof=0)/N)
 print(f"【初期標本の検定統計量T】 {T}") 
 
 #離散一様分布 (0～N-1の一様な分布)
 uni = stats.randint(0,N)
 
 #初期標本からB個のブートストラップ標本をつくる
 for i in range(B):
   #復元抽出する標本のインデクスをランダムに決める
   idx = np.sort(uni.rvs(size=N))
   #初期標本からブートストラップ標本を作る
   b_samples[i] = samples[idx]
 
 #作成したB個のブートストラップ標本をデータフレームに入れる
 df_bs = pd.DataFrame(b_samples)
 #print("【ブートストラップ標本】\n", df_bs)
 
 #帰無仮説H0下での検定統計量を推定するため，値をシフトする
 df_bs_z = df_bs - samples.mean() + mu0
 #print("【H0下での標本にシフト】\n", df_bs_z) 
 
 #シフトしたブートストラップ標本から，検定統計量を求める
 T_z = (df_bs_z.mean() - mu0) / np.sqrt(df_bs_z.var(ddof=0)/N)
 #検定統計量の分布を描く
 ax = T_z.plot.hist(bins=50, title="【H0を仮定した場合の検定統計量の分布】")
 #初期標本の検定統計量の位置
 ax.axvline(x=T, c="green")
 
 #p値は，T_zの中で，Tより小さい部分の数が全体に占める割合
 p_value = T_z[T_z<=T].count() / T_z.count()
 
 print("【p値】", p_value)
 if (p_value < alpha):
   print(f"帰無仮説H0: μ0 = {mu0} を棄却します．")
 else:
   print(f"帰無仮説H0: μ0 = {mu0} を保留します．")

*** 2群比較のための仮説検定 (対応がない場合) [#t414d746]

【例題】ある企業において20代の従業員10名と，50代の従業員10名を無作為抽出し，一般教養クイズをしてもらったところ，以下の点数を得た．

|No.|#1|#2|#3|#4|#5|#6|#7|#8|#9|#10|平均|h
|20代|70|64|83|67|75|76|71|76|52|69|70.3|
|50代|67|80|86|71|79|67|71|84|78|83|76.6|


この時，この企業における20代と50代の一般教養に差があるかどうか，有意水準α=0.05で検定しなさい．なお，母集団は正規分布を仮定してよい．

■仮説の設定
- H0: μ1 ＝ μ2
- H1: μ1 ≠ μ2
- alpha: 0.05

 #有意水準
 alpha = 0.05
 
 #20代10名の得点
 a = np.array([70, 64, 83, 67, 75, 76, 71, 76, 52, 69])
 #50代10名の得点
 b = np.array([67, 80, 86, 71, 79, 67, 71, 84, 78, 83])
 
 #検定統計量
 T = (a.mean() - b.mean()) / np.sqrt(a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b))
 
 print(f"【検定統計量 T】{T}")
 
 #t分布の自由度
 d = ((a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b)))**2 / \
     (((a.var(ddof=1)/len(a))**2 / (len(a)-1)) + 
     ((b.var(ddof=1)/len(b))**2 / (len(b)-1)) ) 
 
 print(f"【t分布の自由度 d】 {d}")
 
 #自由度dのt分布
 tdist = stats.t(df=d)
 
 #p値を求める．Tの絶対値を取り，t(d)(X<=-T) + t(d)(X>=T)を求める
 absT = np.abs(T)
 p_value = tdist.cdf(-absT) + tdist.sf(absT)
 print(f"【p値】{p_value}")
 
 if (p_value < alpha):
   print("H0を棄却し，両グループに差を認めます")
 else:
   print("H0を保留します．両グループに有意な差が認められません")
 
上記のt検定は，以下の1行で実行できる
- stats.ttest_ind(群1,群2, equal_var=等分散を仮定するか否か)

 stats.ttest_ind(a,b, equal_var=False)

*** 2群比較のための仮説検定 (対応がある場合) [#x0aaf1a2]

数値的に全く同じデータでも，''2群間に対応がある''場合，異なる検定方法になる
- 1群の2つのデータの差を検定することになる

【例題】ある大企業において20代の従業員に一般教養の研修を受講してもらったところ，受講前と受講後のテストの成績が以下のようになった．

|No.|#1|#2|#3|#4|#5|#6|#7|#8|#9|#10|平均|h
|受講前|70|64|83|67|75|76|71|76|52|69|70.3|
|受講後|67|80|86|71|79|67|71|84|78|83|76.6|

この時，この研修の効果があったかどうか，有意水準α=0.05で検定しなさい．なお，母集団は正規分布を仮定してよい．

■仮説の設定
- H0: μ1 ＝ μ2
- H1: μ1 ＜ μ2 (受講後の方が点数が高い)
- alpha: 0.05　片側検定
- 2群に対応がある

 stats.ttest_rel(a,b, alternative="less")

あるいは，受講後と受講前の差分を取り，それがプラスかどうかを1群の検定を行う
|No.|#1|#2|#3|#4|#5|#6|#7|#8|#9|#10|h
|受講後-受講前|-3|+16|+3|+4|+4|+11|0|+8|+26|+14|

■仮説の設定
- H0: μ2 - μ1 ＝ 0
- H1: μ2 - μ1 ＞ 0 (受講後の点数増が0より大きい)
- alpha: 0.05　片側検定
- 1群


 stats.ttest_1samp(b-a, popmean=0, alternative="greater")

** さらなる勉強のために [#v65b80eb]
- [[【Pythonで行う】仮説検定:https://corvus-window.com/python_hypothesis-test/]]
- [[統計Web:https://bellcurve.jp/statistics/course/]]


** 第2回演習へ [#o4f8d75d]
- [[第2回/演習]]
トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS