第2回

第2回ハンズオン†

↑

0. 準備†

↑

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

↑

scipy.stats の書式†

↑

標準正規分布の概形†

#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分布の概形†

#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()
↑

カイ二乗分布の概形†

#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の範囲を適当に決める (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()

上記は標準正規分布を使った近似検定だが，scipy.statsの2項検定を使えば，正確な2項分布を使って検定できる

↑

例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 = (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行で書ける

↑

3. 各種検定手法†

↑

正規性判定†

Shapiro-Wilk検定

stats.shapiro(boys)

Kolmogorov-Smirnov検定

stats.kstest(boys, "norm")
↑

Q-Qプロット†

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

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

ブートストラップ検定†

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

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

1167729511545371498327137218092970844
56042075911932121634309282852804273
50753734572701445183647439119021926
22101545445408299393196415545955680
377775117839516817435920760396093069

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

#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群比較のための仮説検定 (対応がない場合)†

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

No.#1#2#3#4#5#6#7#8#9#10平均
20代7064836775767176526970.3
50代6780867179677184788376.6

この時，この企業における20代と50代の一般教養に差があるかどうか，有意水準α=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(a,b, equal_var=False)
↑

2群比較のための仮説検定 (対応がある場合)†

数値的に全く同じデータでも，2群間に対応がある場合，異なる検定方法になる

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

No.#1#2#3#4#5#6#7#8#9#10平均
受講前7064836775767176526970.3
受講後6780867179677184788376.6

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

■仮説の設定

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

あるいは，受講後と受講前の差分を取り，それがプラスかどうかを1群の検定を行う

No.#1#2#3#4#5#6#7#8#9#10
受講後-受講前-3+16+3+4+4+110+8+26+14

■仮説の設定

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

さらなる勉強のために†

トップ   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS