第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()

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

ブートストラップ検定†

ある工場で製造される電球は，平均寿命が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} を保留します．")
トップ   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS