今回は標本の平均と分散がどのような分布になるかをシミュレーションします。

今回やりたいこと

標準正規分布に独立に従うサイズnの標本について、その平均と分散がどのような分布になるか、nを動かして調べます。

例えば以下のコードを実行すると、平均μ=0で標準偏差σ=1の正規分布からサイズn=5の標本をランダムにサンプルできます。

set.seed(123)

mu <- 0
sigma <- 1

n <- 5
sample <- rnorm(n = n, mean = mu, sd = sigma)

すると、

> sample
[1] -0.56047565 -0.23017749  1.55870831  0.07050839
[5]  0.12928774

となり、この平均と分散を計算すると、

> mean(sample)
[1] 0.1935703
> var(sample)
[1] 0.6577564

となります。

今回は、このsampleを生成する過程を繰り返しsampleの平均mean(sample)と分散var(sample)の値を集計した結果がどのような分布になるかを見ていきます。またn=5以外にもn=10, 20, 50の場合にどのようになるかを見ていきます。

標本平均

まず標本平均の分布を調べます。 N=10⁵回ランダムサンプリングしてデータを生成します。

d <- data.frame(mean = c(), var = c(), n = c())
N <- 1e+5
for(n in c(5,10,20,50)){
  means <- c()
  vars <- c()
  for(i in 1:N){
    sample <- rnorm(n = n, mean = mu, sd = sigma)
    means[i] <- mean(sample)
    vars[i] <- var(sample)
  }
  d <- rbind(d, data.frame(means, vars, n = as.factor(rep(n, N))))
}

これらについて、標本平均の分布を図示すると、

ggplot(data = d, aes(x = means, fill = n)) +
  geom_histogram(position = "identity", alpha = 0.5)

f:id:unaoya:20190116090822p:plain — 標本平均

同じものを分割して表示すると

ggplot(data = d, aes(x = means, fill = n)) +
  geom_histogram() +
  facet_wrap( ~ n)

f:id:unaoya:20190116090838p:plain — 標本平均を分割して

nを増加させると幅が狭くなりますが、いずれも0を中心に左右対称な形になります。

これらの平均と標準偏差は、

> d %>%
+   dplyr::select(means, n) %>%
+   dplyr::group_by(n) %>%
+   dplyr::summarise_all(funs(mean,sd))
# A tibble: 4 x 3
  n           mean    sd
  <fct>      <dbl> <dbl>
1 5      0.0000531 0.445
2 10    -0.00197   0.316
3 20     0.0000108 0.223
4 50     0.000201  0.142

です。平均はいずれのnでもおよそ0なのに対し、標準偏差はnが増加すると減少します。

この分布は実は正規分布になることがわかります。いわゆる正規分布の再生性という性質です。このことについては次回以降詳しくみていきます。

標本分散

次に標本分散の分布です。 Rで計算されるのは(n-1)で割る方の不偏分散ですが、定数倍なので気にしないことにします。

上でサンプリングしたデータを用いて同様に図示していきます。

ggplot(data = d, aes(x=vars, fill = n)) +
  geom_histogram(position = "identity", alpha = 0.5)

f:id:unaoya:20190116091122p:plain — 標本分散

同じものを分割して表示すると

ggplot(data = d, aes(x = vars, fill = n)) +
  geom_histogram() +
  facet_wrap( ~ n)

f:id:unaoya:20190116091140p:plain — 標本分散を分割して

nを増加させると幅が狭くなり、さらにピークが右に寄っていきます。

これらの平均と標準偏差は、

> d %>%
+   dplyr::select(vars, n) %>%
+   dplyr::group_by(n) %>%
+   dplyr::summarise_all(funs(mean,sd))
# A tibble: 4 x 3
  n      mean    sd
  <fct> <dbl> <dbl>
1 5     1.00  0.711
2 10    0.997 0.468
3 20    1.00  0.325
4 50    1.000 0.202