先週のおさらい + 演習問題の確認

男性の身長に関する母集団平均を\(\mu_1\) (データ数\(N_1\)、母集団分散\(\sigma_1^2\))、女性の身長に関する母集団平均を\(\mu_2\)とすると(データ数\(N_2\)、母集団分散\(\sigma_2^2\))、知りたい値は\(\mu_1 - \mu_2\)。\(\mu_1 - \mu_2\)に関する信頼区間を計算するために、 \[\begin{equation} \frac{1}{N_1}\sum_{i = 1}^{N_1}x_{1i} - \frac{1}{N_2}\sum_{j = 1}^{N_2}x_{2j} \end{equation}\] を標準化する(流れは先週の講義資料を参照)。ただし、\(x_{1i}\)は男性の身長に関する母集団分布に従う確率変数\(X_{1i}\)の実現値、\(x_{2i}\)は女性の身長に関する母集団分布に従う確率変数\(X_{2i}\)の実現値である。





標準化に必要な値を計算する。 \[\begin{equation} E\left[ \frac{1}{N_1}\sum_{i = 1}^{N_1}x_{1i} - \frac{1}{N_2}\sum_{j = 1}^{N_2}x_{2j} \right] = \mu_1 - \mu_2 \end{equation}\] \[\begin{equation} \mathrm{Var}\left[ \frac{1}{N_1}\sum_{i = 1}^{N_1}x_{1i} - \frac{1}{N_2}\sum_{j = 1}^{N_2}x_{2j} \right] = \sigma_1^2/N_1 + \sigma_2^2/N_2 \end{equation}\] 標準化されたデータ平均の差分は、 \[\begin{equation} z = \frac{\frac{1}{N_1}\sum_{i = 1}^{N_1}x_{1i} - \frac{1}{N_2}\sum_{j = 1}^{N_2}x_{2j} - (\mu_1 - \mu_2)}{\sqrt{\sigma_1^2/N_1 + \sigma_2^2/N_2}} \end{equation}\] である。これはデータ数が十分多ければ平均0、分散1のガウス分布に従うため、95%の信頼度で \[\begin{equation} -1.96 \le z \le 1.96 \end{equation}\] 99%の信頼度で \[\begin{equation} -2.58 \le z \le 2.58 \end{equation}\] を満たす。あとはそれぞれ値を代入して、\(\mu_1 - \mu_2\)について解けばOK。





統計的検定を行う場合、第一に帰無仮説\(\mathrm{H}_0: \mu = 1600\)、対立仮説\(\mathrm{H}_1: \mu \neq 1600\)をたてる。そして、帰無仮説がどれほど妥当かを検証する。ただし、\(\mu\)は、電球の使用可能時間の母集団平均。





第一に、電球の使用可能時間のデータ平均 \[\begin{equation} \bar{x} = \frac{1}{N}\sum_{i = 1}^Nx_i \end{equation}\] を標準化すると、 \[\begin{equation} z = \frac{\bar{x} - \mu}{\sigma / \sqrt{N}} \end{equation}\] を得る。値を代入していくと、 \[\begin{equation} z = 1598 - \mu \end{equation}\] を得る。データ数が十分多ければ平均0、分散1のガウス分布に従うため、95%の信頼度で \[\begin{equation} -1.96 \le z \le 1.96 \end{equation}\] 99%の信頼度で \[\begin{equation} -2.58 \le z \le 2.58 \end{equation}\] を満たす(ここまでは信頼区間の計算方法と同じ)。





そして、帰無仮説\(\mathrm{H}_0: \mu = 1600\)がどれほど妥当かを確認する。具体的には、 \[\begin{equation} z = 1598 - \mu \end{equation}\] に\(\mu = 1600\)を代入する。すると、 \[\begin{equation} z = -2 \end{equation}\] であり、95%信頼区間には含まれない。このとき、5%有意水準にて帰無仮説は棄却されるとみなす(つまり、\(\mu = 1600\)という表示は、5%有意水準[95%信頼区間]を基準にすると妥当ではない)。一方、\(z = -2\)は99%信頼区間には含まれる。このとき、1%有意水準にて帰無仮説は棄却されない(つまり、\(\mu = 1600\)という表示は、1%有意水準[99%信頼区間]を基準にすると棄却できない[\(\mu = 1600\)]が正しいという意味ではない、cf. \(\mu = 1600.1\))。一般的には、このような場合でも5%有意水準(95%信頼区間)を基準にすることが多いように思う。





  新しい用語として、有意水準という言葉が出てきた。これは、\(100(1 - \alpha)\)%信頼区間を元に帰無仮説を検証する場合、有意水準 \(\alpha\) にて帰無仮説が棄却できるかどうかを判断するものである。つまり、帰無仮説を棄却する基準を意味する。これまでに、95%信頼区間よりも、99%信頼区間の方が保守的であり範囲が広いことを示してきた。つまり、より広い99%信頼区間にさえ \(z\) が含まれないようであれば、設定した帰無仮説はより“強く”棄却されるべきであろう。狭い95%信頼区間に \(z\) が含まれないならば、設定した帰無仮説の棄却はより保守的であるべきだと考えられる。つまり、語弊を恐れなければ、有意水準は帰無仮説を棄却する“強さ”のようなものを表しているといえる。






  近年では、(分野にもよるが) 5%有意水準、1%有意水準の議論では不十分であることが多い。例えば、帰無仮説に従い計算した\(z\)が、\(z = -2.581\)のときと\(z = -10\)ではどちらも99%信頼区間に含まず、帰無仮説が1%有意水準で棄却される。しかし、99%信頼区間からの外れ具合は大分異なる。つまり、\(z = -10\)の方が、より強く帰無仮説が棄却されるべきであろう。この違いはどのように定量化すべきであろうか？





p値(研究において頻出)

もう一度、95%信頼区間について確認する。95%信頼区間は、平均0、分散1のガウス分布が描く、以下の灰色の領域である。標準化したデータ平均\(z\)についての式で書くならば、 \[\begin{equation} P(-1.96 \le z \le 1.96) = 0.95 \end{equation}\] である。言い換えるならば、\(P(z \ge 1.96) = (1 - 0.95)/0.25\)である。もしくは、\(P(|z| \ge 1.96) = 0.05\)。

curve(dnorm, -5, 5)

a = -1.96
b = 1.96
i = 200 #分割数

xx = seq(a, b, length=i) 
yy = dnorm(xx) 

xx<- c(a, xx, b, a)
yy<- c(0, yy, 0, 0)

polygon(xx, yy, col="gray")

別の場合についても見てみよう。以下は、 \[\begin{equation} P(-3 \le z \le 3) = 0.9973 \end{equation}\] を示している。言い換えるならば、\(P(z \ge 3.00) = (1 - 0.9973)/0.25\)である。もしくは、\(P(|z| \ge 3.00) = 0.0027\)。である。

curve(dnorm, -5, 5)

a = -3.00
b = 3.00
i = 200 #分割数

xx = seq(a, b, length=i) 
yy = dnorm(xx) 

xx<- c(a, xx, b, a)
yy<- c(0, yy, 0, 0)

polygon(xx, yy, col="gray")

a = -3.00
b = 3.00


f <- function(x) exp(-0.5*x^2)/sqrt(2*pi) # 積分する関数を定義
integrate(f, a, b) # 積分実行
## 0.9973002 with absolute error < 9.3e-07

いま、帰無仮説にもとづき計算した\(z\)は、\(z = y\)を満たすとする。このとき、 \[\begin{equation} P(|z| > |y|) \end{equation}\] を(両側検定の)p値と呼ぶ。例えば\(y = 1.96\)ならばp値は0.05(95%信頼区間、5%有意水準)、\(y = 2.58\)ならばp値は0.01(99%信頼区間、1%有意水準)、\(y = 3.00\)ならばp値は0.0027となる。つまり、p値が0.05より小さければ、5%有意水準にて帰無仮説は棄却されると判断される。加えて、\(p = 0.0499\)よりも\(p = 0.0027\)の方が強く帰無仮説が棄却されると判断できる。したがって、p値をしばしば統計的検定の基準値として利用する。





(補足) 上記にて両側検定という言葉があるが、これは対立仮説として\(\mathrm{H_1}: \mu_0 \neq \mu_1\)を想定することを意味する。つまり、比較したい\(\mu_0\)と\(\mu_1\)に対して、どちらが大きいか、どちらが小さいかは想定しない。対立仮説として\(\mathrm{H_1}: \mu_0 > \mu_1\)を想定する統計的検定を片側検定と呼ぶ。場合にもよるが、片側検定の方が帰無仮説が棄却されやすいため、基本はおすすめしない。よほど片側検定が妥当な場面ならば利用することもあるかもしれない。





N < 30の場合(統計をツールとして利用するならば、超重要)

  これまで、\(N \ge 30\) を想定して、分散の不偏推定量が母集団の分散とおおかた一致すること、そして中心極限定理がおおよそ成り立つことを想定した。それでは \(N < 30\)ではどのようになるのか検証してみよう。 これまでの議論から、 \[\begin{equation} z = \frac{\bar{x} - \mu}{\sqrt{ \frac{\sigma^2}{N}}} \end{equation}\] は平均0、分散1になる。\(N\)を色々変えて検証してみよう。下の図は、赤(\(N=3\))、緑(\(N=5\))、青(\(N=10\))、黒(\(N=30\))の標本分布(データ平均を何度も計算してヒストグラムを描いたもの、直線は平均0、分散1のガウス分布)。 





Ite = 1000
z1 = matrix(0, Ite, 1)
z2 = matrix(0, Ite, 1)
z3 = matrix(0, Ite, 1)
z4 = matrix(0, Ite, 1)

for(i in 1:Ite){
  dat = rnorm(30, 0, 1)
  z1[i] = (mean(dat[1:3]) - 0)/(sqrt(var(dat[1:3])/1))
  z2[i] = (mean(dat[1:5]) - 0)/(sqrt(var(dat[1:5])/5))
  z3[i] = (mean(dat[1:10]) - 0)/(sqrt(var(dat[1:10])/10))
  z4[i] = (mean(dat[1:30]) - 0)/(sqrt(var(dat[1:30])/30))
}

par(mfrow=c(1,4)) 
hist(z1, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 20, col="#ff000020")
curve(dnorm, -6, 6, add = T)

hist(z2, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 20, col="#00ff0020")
curve(dnorm, -6, 6, add = T)

hist(z3, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 20, col="#0000ff20")
curve(dnorm, -6, 6, add = T)

hist(z4, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 20, col="#00000020")
curve(dnorm, -6, 6, add = T)

黒のように\(N=30\)のときは、確かに中心極限定理が成り立っているように見える。ただし、赤(N=3)、緑(N=5)、青(N=10)の例ではどうであろうか？特に赤(N=3)、緑(N=5)をもう少し拡大してみてみよう。 





Ite = 10000
z1 = matrix(0, Ite, 1)
z2 = matrix(0, Ite, 1)
z3 = matrix(0, Ite, 1)
z4 = matrix(0, Ite, 1)

for(i in 1:Ite){
  dat = rnorm(30, 0, 1)
  z1[i] = (mean(dat[1:3]) - 0)/(sqrt(var(dat[1:3])/1))
  z2[i] = (mean(dat[1:5]) - 0)/(sqrt(var(dat[1:5])/5))
  z3[i] = (mean(dat[1:10]) - 0)/(sqrt(var(dat[1:10])/10))
  z4[i] = (mean(dat[1:30]) - 0)/(sqrt(var(dat[1:30])/30))
}

par(mfrow=c(1,2)) 
hist(z1, prob = T, ylim = c(0, 0.6), xlim = c(-6, 6), 400, col="#ff000020")
curve(dnorm, -6, 6, add = T)

hist(z2, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 100, col="#00ff0020")
curve(dnorm, -6, 6, add = T)

ややガウス分布からずれていることがわかる。特に、ガウス分布ではほとんど出現しない-3より小さい範囲、3より大きい範囲にて複数の値が観測される。95%信頼区間、99%信頼区間の計算に用いてきた表はガウス分布の想定の元に妥当なものであり、\(N\)が小さい範囲では妥当性が低いことがわかる。ちなみにこのズレは、\(N\)が小さいときには分散の不偏推定量が母集団の分散と大きく異なる確率が高いことに起因する。






  \(N < 30\)のとき、より妥当な確率分布は自由度 \(N-1\) のスチューデントのt分布であることが知られている。重ね書きしてみると、

Ite = 10000
z1 = matrix(0, Ite, 1)
z2 = matrix(0, Ite, 1)
z3 = matrix(0, Ite, 1)
z4 = matrix(0, Ite, 1)

for(i in 1:Ite){
  dat = rnorm(30, 0, 1)
  z1[i] = (mean(dat[1:3]) - 0)/(sqrt(var(dat[1:3])/1))
  z2[i] = (mean(dat[1:5]) - 0)/(sqrt(var(dat[1:5])/5))
  z3[i] = (mean(dat[1:10]) - 0)/(sqrt(var(dat[1:10])/10))
  z4[i] = (mean(dat[1:30]) - 0)/(sqrt(var(dat[1:30])/30))
}

par(mfrow=c(1,2)) 
hist(z1, prob = T, ylim = c(0, 0.6), xlim = c(-6, 6), 400, col="#ff000020")
curve(dt(x, 2), -6, 6, add = T)

hist(z2, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 100, col="#00ff0020")
curve(dt(x, 4), -6, 6, add = T)

となり、赤(\(N=3\))ではさすがにまだずれがあるものの、緑(\(N=5\))では適合度が高いように見える。\(N = 10\)、\(N = 30\)のときも見てみると、適合度が高いように見える。

Ite = 10000
z1 = matrix(0, Ite, 1)
z2 = matrix(0, Ite, 1)
z3 = matrix(0, Ite, 1)
z4 = matrix(0, Ite, 1)

for(i in 1:Ite){
  dat = rnorm(30, 0, 1)
  z1[i] = (mean(dat[1:3]) - 0)/(sqrt(var(dat[1:3])/1))
  z2[i] = (mean(dat[1:5]) - 0)/(sqrt(var(dat[1:5])/5))
  z3[i] = (mean(dat[1:10]) - 0)/(sqrt(var(dat[1:10])/10))
  z4[i] = (mean(dat[1:30]) - 0)/(sqrt(var(dat[1:30])/30))
}

par(mfrow=c(1,2)) 
hist(z3, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 100, col="#0000ff20")
curve(dt(x, 9), -6, 6, add = T)

hist(z4, prob = T, ylim = c(0, 0.45), xlim = c(-6, 6), 100, col="#00000020")
curve(dt(x, 29), -6, 6, add = T)

加えて、スチューデントのt分布は自由度が十分大きいときはガウス分布と重なるため、常にスチューデントのt分布を利用したほうが妥当性が高いこともわかる。下記では、スチューデントのt分布は赤線、ガウス分布は黒線で表している。左上が\(N = 3\), 右上が\(N = 5\), 左下が\(N = 10\)、右下が\(N = 30\)。 





par(mfrow=c(2,2)) 
curve(dnorm, -6, 6)
curve(dt(x, 2), -6, 6, add = T, col = "red")

curve(dnorm, -6, 6)
curve(dt(x, 4), -6, 6, add = T, col = "red")

curve(dnorm, -6, 6)
curve(dt(x, 9), -6, 6, add = T, col = "red")

curve(dnorm, -6, 6)
curve(dt(x, 29), -6, 6, add = T, col = "red")

N=30 (右下) でほぼ一致することがわかる。また、自由度\(\nu\)のスチューデントのt分布は、 \[\begin{equation} p(x | \nu) = \frac{\Gamma(\frac{\nu + 1}{2})}{\sqrt{\nu\pi} \Gamma(\frac{\nu}{2}} \left(1 + \frac{x^2}{\nu}\right)^{- \frac{\nu + 1}{2}} \end{equation}\] と定義される。ただし、\(\Gamma(x) = \int_0^{\infty} t^{x-1}e^{-t}dt\) はガンマ関数である。

N < 30の場合の信頼区間

  上記の議論をまとめると、データ数\(N\)、データ平均を\(\bar{x}\)、分散の不偏推定量を\(\bar{\sigma}^2\)、未知の母集団の平均を\(\mu\)とすると、 \[\begin{equation} z = \frac{\bar{x} - \mu}{\frac{\bar{\sigma}}{\sqrt{N}}} \end{equation}\] は自由度 \(N-1\) のスチューデントのt分布に従う。すると、信頼区間を計算するための表を改めて作成する必要がでてくる。実際に計算してみると、

Nmax = 50

conf95 = matrix(0, 1, Nmax)
conf99 = matrix(0, 1, Nmax)

for(N in 2:Nmax){
  nu = N - 1
  f <- function(x) gamma((nu+1)*0.5)/(sqrt(nu*pi)*gamma((nu)*0.5))*(1+x^2/nu)^(-(nu+1)*0.5) # 積分する関数を定義
  
  dseq = seq(0, 7, by = 0.001)
  f95 = matrix(0, 1, length(seq))
  f99 = matrix(0, 1, length(seq))
  
  for(i in 1:length(dseq)){
    fv = integrate(f, -dseq[i], dseq[i]) # 積分実行
    f95[i] = abs(0.95 - fv$value)
    f99[i] = abs(0.99 - fv$value)
  }
  
  conf95[N] = dseq[which.min(f95)]
  conf99[N] = dseq[which.min(f99)]
}

print(conf95)
##      [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## [1,]    0    7 4.303 3.182 2.776 2.571 2.447 2.365 2.306 2.262 2.228 2.201
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 2.179  2.16 2.145 2.131  2.12  2.11 2.101 2.093 2.086  2.08 2.074 2.069
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## [1,] 2.064  2.06 2.056 2.052 2.048 2.045 2.042  2.04 2.037 2.035 2.032  2.03
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 2.028 2.026 2.024 2.023 2.021  2.02 2.018 2.017 2.015 2.014 2.013 2.012
##      [,49] [,50]
## [1,] 2.011  2.01
print(conf99)
##      [,1] [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## [1,]    0    7    7 5.841 4.604 4.032 3.707 3.499 3.355  3.25 3.169 3.106 3.055
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## [1,] 3.012 2.977 2.947 2.921 2.898 2.878 2.861 2.845 2.831 2.819 2.807 2.797
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## [1,] 2.787 2.779 2.771 2.763 2.756  2.75 2.744 2.738 2.733 2.728 2.724 2.719
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
## [1,] 2.715 2.712 2.708 2.704 2.701 2.698 2.695 2.692  2.69 2.687 2.685 2.682
##      [,50]
## [1,]  2.68
x = data.frame(N = 1:Nmax, conf_95 = t(conf95), conf_99 = t(conf99))

x
##     N conf_95 conf_99
## 1   1   0.000   0.000
## 2   2   7.000   7.000
## 3   3   4.303   7.000
## 4   4   3.182   5.841
## 5   5   2.776   4.604
## 6   6   2.571   4.032
## 7   7   2.447   3.707
## 8   8   2.365   3.499
## 9   9   2.306   3.355
## 10 10   2.262   3.250
## 11 11   2.228   3.169
## 12 12   2.201   3.106
## 13 13   2.179   3.055
## 14 14   2.160   3.012
## 15 15   2.145   2.977
## 16 16   2.131   2.947
## 17 17   2.120   2.921
## 18 18   2.110   2.898
## 19 19   2.101   2.878
## 20 20   2.093   2.861
## 21 21   2.086   2.845
## 22 22   2.080   2.831
## 23 23   2.074   2.819
## 24 24   2.069   2.807
## 25 25   2.064   2.797
## 26 26   2.060   2.787
## 27 27   2.056   2.779
## 28 28   2.052   2.771
## 29 29   2.048   2.763
## 30 30   2.045   2.756
## 31 31   2.042   2.750
## 32 32   2.040   2.744
## 33 33   2.037   2.738
## 34 34   2.035   2.733
## 35 35   2.032   2.728
## 36 36   2.030   2.724
## 37 37   2.028   2.719
## 38 38   2.026   2.715
## 39 39   2.024   2.712
## 40 40   2.023   2.708
## 41 41   2.021   2.704
## 42 42   2.020   2.701
## 43 43   2.018   2.698
## 44 44   2.017   2.695
## 45 45   2.015   2.692
## 46 46   2.014   2.690
## 47 47   2.013   2.687
## 48 48   2.012   2.685
## 49 49   2.011   2.682
## 50 50   2.010   2.680

のようになり、データ数 \(N\) (自由度 \(N-1\)) に依存して、95%信頼区間、99%信頼区間に利用する値が異なる。例えば、\(N=10\)のとき、

信頼性(N=10) 95% 99%
\(\alpha\) 2.26 3.25

となる。そして、95%信頼区間は、 \[\begin{equation} \bar{x} - 2.26\frac{\bar{\sigma}}{\sqrt{N}} \le \mu \le \bar{x} + 2.26\frac{\bar{\sigma}}{\sqrt{N}} \end{equation}\] となる。

  通常、スチューデントのt分布の信頼区間に関する表は、自由度 \(\nu = N-1\) と対応して提示されることが多い。