演習問題

*標本平均を \(\bar{X}\) 、標準偏差の不偏推定量を \(\bar{S}\)、データサイズを \(N \ge 30\) とするとき、未知の母集団平均 \(\mu\) の95%信頼区間、99%信頼区間を計算せよ (上記の計算を確認せよ)

*上記にて、\(\bar{X} = 155\) 、標準偏差の不偏推定量を \(\bar{S} = 5\)、データサイズを \(N = 100\) とするとき、未知の母集団平均 \(\mu\) の95%信頼区間、99%信頼区間を計算せよ (上記の計算を確認せよ)

統計の概念図 (再掲)

需要な点なので、再掲する。

母集団からランダムにサンプリングされるデータを確率変数とみなして、計測されるデータから未知の母集団の平均や分散を推定する。これが統計の枠組みである。
具体的に計測されたデータ (e.g., 150cm, 165cm, … などの身長データ) の平均を計算すると、期待値としては真の平均に一致することを示した。ただし、あくまで期待値として一致するのみであり、分散 \(\frac{\sigma^2}{N}\) と伴いばらつく (\(\sigma^2\)も未知)。

差の信頼区間

引き続き簡略化のため、\(N \ge 30\) においては中心極限定理はおおよそ成り立ち、標準偏差(or 分散)の不偏推定量は、おおよそ母集団の標準偏差 (or 分散)と近いものと仮定する。

中心極限定理は、\(N \to \infty\) にて、標準化したデータ平均がガウス分布になることを保証してくれるものである。
つまり、データ \(x_1, ..., x_N\) にて (すべて独立かつ同じ分布から生成、平均 \(\mu\)、分散 \(\sigma^2\))、 \[\begin{equation} z = \frac{\frac{1}{N}\sum_i x_i - \mu}{\frac{\sigma}{\sqrt{N}}} \end{equation}\] という確率変数 \(z\) は、\(N\to\infty\)にて平均 0、分散 1のガウス分布へと収束していくことを意味する。

しばしば統計で扱う対象は、AとBとでどの程度差があるかという問題である。例えば、男性と女性とで身長が異なるか、という問題を考えてみよう。女性の身長データを \(x_{11}, ..., x_{1N_1}\)、男性の身長データを \(x_{21}, ..., x_{2N_2}\) とする。考えるべき問題は、女性の身長データの平均 \(\bar{x}_1 = \frac{1}{N_1}\sum_{i=1}^{N_1}x_{1i}\)、 男性の身長データの平均 \(\bar{x}_2 = \frac{1}{N_2}\sum_{i=1}^{N_2}x_{2i}\) の差の信頼区間であろう。これが0と大きく異なるようであれば差があると言えるであろう。

女性の身長データは未知の母集団分布に従う確率変数であり、平均 \(\mu_1\)、分散 \(\sigma_1^2\) をもつものとする。男性の身長データは未知の母集団分布に従う確率変数であり、平均 \(\mu_2\)、分散 \(\sigma_2^2\) をもつものとする。すると、変数変換と平均・分散の関係から、 \[\begin{equation} E[\bar{x}_1] = \mu_1 \end{equation}\] \[\begin{equation} \mathrm{Var}[\bar{x}_1] = \frac{\sigma_1^2}{N_1} \end{equation}\] \[\begin{equation} E[\bar{x}_2] = \mu_2 \end{equation}\] \[\begin{equation} \mathrm{Var}[\bar{x}_2] = \frac{\sigma_2^2}{N_2} \end{equation}\] が成り立つ。

さらに、\(x_{1i}\)と\(x_{2j}\)は独立であると想定すると(ランダムサンプリングを想定)、 \[\begin{equation} E[\bar{x}_1 \pm \bar{x}_2] = \mu_1 \pm \mu_2 \end{equation}\] \[\begin{align} \mathrm{Var}[\bar{x}_1 \pm \bar{x}_2] &= E\left[ \left(\sum_i\frac{1}{N_1}x_{1i} \pm \sum_i\frac{1}{N_2}x_{2i} - (\mu_1 \pm \mu_2) \right)^2 \right] \nonumber\\ &= E\left[ \left(\sum_i\frac{1}{N_1}x_{1i} - \mu_1\right)^2\right] + E\left[ \left(\sum_i\frac{1}{N_2}x_{2i} - \mu_2\right)^2\right] \nonumber\\ &= \mathrm{Var}(\bar{x}_1) + \mathrm{Var}(\bar{x}_2)\nonumber\\ &= \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2} \end{align}\] を得る。つまり、平均の和を考えようが差を考えようが、分散は何かを足すと大きくなる。

確かめてみると、

Ite = 1000
mu1 = matrix(0, Ite, 1)
mu2 = matrix(0, Ite, 1)
mud = matrix(0, Ite, 1)

for(i in 1:Ite){
mu1[i] = mean(rnorm(1000, 0, 1))
mu2[i] = mean(rnorm(1000, 0.5, 1))
mud[i] = mu1[i] - mu2[i]
}

hist(mu1, prob = T, ylim = c(0, 12), xlim = c(-1.0, 0.75), col="#ff000020")
hist(mu2, prob = T, ylim = c(0, 12), xlim = c(-1.0, 0.75), col="#00ff0020", add = T)
hist(mud, prob = T, ylim = c(0, 12), xlim = c(-1.0, 0.75), col="#00000020", add = T)

確かに黒いヒストグラム (赤 - 緑) は分散が大きいことがわかる。そして、

print(c(var(mu1), var(mu2), var(mud), 1/1000 + 1/1000))
## [1] 0.0010082505 0.0009636841 0.0020066476 0.0020000000

理論解も正しいことがわかる。以上より、\(\bar{x}_1 \pm \bar{x}_2\)は、 \[\begin{equation} z = \frac{\bar{x}_1 \pm \bar{x}_2 - (\mu_1 \pm \mu_2)}{\sqrt{ \frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2} }} \end{equation}\] とすると、中心極限定理より \(z\) は平均 0、分散 1 のガウス分布に従う。確かめてみると、

Ite = 1000
mu1 = matrix(0, Ite, 1)
mu2 = matrix(0, Ite, 1)
mud = matrix(0, Ite, 1)
z = matrix(0, Ite, 1)


for(i in 1:Ite){
  dat1 = rnorm(1000, 0, 1)
  dat2 = rnorm(1000, 0.5, 1)
  mu1[i] = mean(dat1)
  mu2[i] = mean(dat2)
  mud[i] = mu1[i] - mu2[i]

  z[i] = (mud[i] - (0 - 0.5))/sqrt((var(dat1)/1000 + var(dat2)/1000))

}

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

となり、たしかに中心極限定理が成立してそうである。

下記の表を利用すると、

信頼性 95% 99%
\(\alpha\) 1.96 2.58

95%信頼区間は、 \[\begin{equation} -1.96 \le z \le 1.96 \end{equation}\] より、 \[\begin{equation} \bar{x}_1 - \bar{x}_2 - 1.96\sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}} \le \mu_1 - \mu_2 \le \bar{x}_1 - \bar{x}_2 + 1.96\sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}} \end{equation}\] となる。同様に、99%信頼区間は、 \[\begin{equation} \bar{x}_1 - \bar{x}_2 - 2.58\sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}} \le \mu_1 - \mu_2 \le \bar{x}_1 - \bar{x}_2 + 2.58\sqrt{\frac{\sigma_1^2}{N_1} + \frac{\sigma_2^2}{N_2}} \end{equation}\] となる。

N < 30の場合

これまで、\(N \ge 30\) を想定して、分散の不偏推定量が母集団の分散とおおかた一致すること、そして中心極限定理がおおよそ成り立つことを想定した。それでは \(N < 30\)ではどのようになるのか検証してみよう。 これまでの議論から、 \[\begin{equation} z = \frac{\bar{x} - \mu}{\sqrt{ \frac{\sigma^2}{N}}} \end{equation}\] は平均0、分散1になる。\(N\)を色々変えて検証してみよう。また、

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)の例ではどうであろうか？もう少し拡大してみてみよう。

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)では適合度が高いように見える。加えて、スチューデントのt分布は自由度が十分大きいときはガウス分布と重なるため、常にスチューデントのt分布を利用したほうが妥当性が高いこともわかる。下記では、スチューデントのt分布は赤線、ガウス分布は黒線で表している。

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\) と対応して提示されることが多い。

演習問題