大数の法則と中心極限定理のシミュレーション

1 Law of Large Numbers 大数の法則

0.6 の確率で表、0.4 の確率で裏が出るコインを何回か投げた時、何割ぐらい表が出るだろうか。期待値は ６割であるが、実際には７割になったり、３割りになったりすることもあり、ピッタリ６割になるとは限らない。しかし、試行回数（ここではコインを投げる回数）を増やせば増やすほど、表の割合は６割に近づいていく。これが大数の法則の概要である。

これが重要なのは、実際のサンプリングに応用が利く知識だからである。例えば、犯罪被害の経験のある人の比率を知りたいとする。被害経験者の真の比率が \(p\) だとすると、母集団から無作為にサンプリングして、犯罪被害の経験を尋ねても、その比率は \(p\) とはかけ離れた値になることもある。しかし、大数の法則より、サンプルサイズ（上の例では試行回数）を増やせば増やすほど、\(p\) に近い値を得られると期待できる。この場合、コイン投げに該当するのが、母集団から無作為に一人サンプルを抽出する作業である。この時、確率 \(p\) で犯罪被害の経験者にあたるので。

これは比率だけでなく、平均値についても同じことが言える。

2 Simulation シミュレーション

大数の法則は数学的に証明されており、わざわざシミュレーションする必要はない。しかし、証明を追うよりもシミュレーションをした方がわかりやすいと思われるので、実際にやってみよう。実際に自分でコインを投げてもいいのだが、コインを何千回も投げてその結果を記録するのは手間がかかるので、コンピュータ上でそれを擬似的に実現しよう。

コイン投げは二項分布の一種なので、二項分布に従う乱数を発生させる関数 rbinom() を使う。

rbinom(発生させる乱数の数, 試行回数, 表が出る確率) # 表が出る回数を返す

例えば、５回続けて 0.4 の確率で表が出るコインを投げた場合、１回も表が出ないかもしれないし、５回表が出るかもしれない。そのような５回投げを３セット（合計 15 回）した場合の結果を rbinom() によって生成すると、以下のようになる。

rbinom(3, 5, 0.4)

## [1] 2 2 2

rbinom(3, 5, 0.4)

## [1] 2 0 3

rbinom(3, 5, 0.4)

## [1] 2 5 1

乱数によって結果を作っているので、毎回異なる結果が出るはずである（たまたま同じになる場合もある）。試行回数を 1 にして、発生させる乱数の数を 6 にすれば、6 回コインを投げた時の毎回の記録を作ることもできる。

set.seed(1986)
rbinom(6, 1, 0.4)

## [1] 0 0 1 0 0 1

この例だと最初の２回は 0 （つまり裏）で、３回めは 1 という具合に、毎回の表裏の結果を生成することができる。このようにしてコイン投げを 20 回分の結果を生成すると、以下のようになる（表が出る確率は 0.5）。

set.seed(1986)
sim20 <- rbinom(20, 1, 0.5)
sim20

##  [1] 0 0 1 0 1 1 0 1 1 1 0 0 1 1 0 0 0 1 0 0

試行回数１回めだけのデータを使うと、１回中１回表が出たので、表が出た比率は 100% 、１～２回目のデータを使うと、２回中１回表が出たので、表が出た比率は 50% 、１～３回めのデータを使うと３回中１回表が出たので、表が出た比率は 33.3%、といった計算を繰り返し行なっていくと、以下のようになる。

cumsum20 <- cumsum(sim20) # 各回までに表が出た回数 cumsum() についてはこの下の解説を参照
cumsum20

##  [1] 0 0 1 1 2 3 3 4 5 6 6 6 7 8 8 8 8 9 9 9

p20 <- cumsum20 / 1:20 # サンプルサイズ（試行回数）に対応する表が出た比率
p20

##  [1] 0.0000000 0.0000000 0.3333333 0.2500000 0.4000000 0.5000000 0.4285714
##  [8] 0.5000000 0.5555556 0.6000000 0.5454545 0.5000000 0.5384615 0.5714286
## [15] 0.5333333 0.5000000 0.4705882 0.5000000 0.4736842 0.4500000

plot(p20, type="l") # type="l" で折れ線グラフになる
abline(h = 0.5, col="red") # abline(h = yの座標) で指定した座標の高さの水平線を引く col= は色の指定

図1 0.5の確率で表が出るコインを１～２０回まで投げた時の表の比率

cumsum(ベクトル) # ベクトルの累積和を返す

累積和とは、引数のベクトルと同じ長さで、返り値の最初の要素は引数のベクトルの最初の要素、返り値の２番めの要素は引数のベクトルの１番目と２番めの要素の和、返り値の３番めの要素は引数のベクトルの１～３番目の要素の和、以下同様に最後まで計算を繰り返す。例えば以下のように。

cumsum(c(2, 7, 3, 1))

## [1]  2  9 12 13

さて、試行回数を１～２０回まで試行回数を増やしていった時の表が出た比率を示したのが、上の図１であるが、最初は 0.5 から大きくかけ離れているが、その後は 近づいていっているのがわかる。さらに試行回数を 2000 回まで伸ばして同じことをしてみよう。

set.seed(1986)
sim2000 <- rbinom(2000, 1, 0.5)
cumsum2000 <- cumsum(sim2000)
p2000 <- cumsum2000 / 1:2000 
plot(p2000, type="l") 
abline(h = 0.5, col="red") 

図2: 試行回数を2000回まで

このように試行回数を増やせばかなり 0.5 に近づいていくことがわかる。ただし、このような接近は単調に生じるのではなく、1000回目あたりで、それ以前よりも 0.5 から離れてしまっている。また、この例では 300回ぐらいまでで、ほとんど 0.5 になっているが、これは偶然によってもっと早く生じるかもしれないし、もっと遅く生じるかもしれない。もう一回同じシミュレーションをやってみよう。

set.seed(1905) # 前回とは同じ数字にしないと同じ結果になる
sim2000 <- rbinom(2000, 1, 0.5)
cumsum2000 <- cumsum(sim2000)
p2000 <- cumsum2000 / 1:2000 
plot(p2000, type="l") 
abline(h = 0.5, col="red") 

図3: 試行回数を2000回までのシミュレーションをもう一度

少し違ったグラフのかたちになるのがわかるだろう。

3 Central Limit Theorem: 中心極限定理

中心極限定理も統計学の基本定理で、平均値の区間推定をはじめとして、さまざまな統計技法の基礎である。無作為抽出したサンプルのある変数\(x\)の平均値（標本平均）は、サンプリングの際の偶然によって、母集団の平均値（母平均）よりも大きくなったり、小さくなったりする。つまり確率的に分布する。しかし、サンプルサイズを大きくするに従って、標本平均の分布は、平均＝母平均、分散\(=N/\sigma^2\)（\(N\)はサンプルサイズ、\(\sigma^2\)は\(x\)の母分散）の正規分布に近似する。この知識によって、さまざまな推測が可能になっているのである。

中心極限定理も数学的に証明されており、シミュレーションの必要はない。しかし、その意味を理解するためにシミュレーションしてみよう。

# 架空の10万人分の収入（単位：百万円）のデータを作り、これを母集団とみなす
set.seed(1986)
population <- round(3 * rexp(100000), 1) # rexp(乱数の数) で指数分布に従う乱数を生成する。これが母集団の収入
population[1:30]

##  [1] 2.5 3.1 0.6 1.0 5.2 1.4 1.1 0.5 1.8 0.5 4.3 0.2 1.8 5.0 1.7 5.7 0.6
## [18] 1.1 2.6 4.4 0.4 0.4 2.4 0.4 1.6 0.1 2.2 2.6 0.6 2.4

pop.mean <- mean(population) # 母平均
pop.mean

## [1] 2.983785

pop.var <- var(population) # 母分散

sample1 <- sample(population, 20) # 上の母集団から 20 人無作為抽出
sample1

##  [1] 2.8 3.7 2.6 1.7 0.2 0.6 0.9 1.2 1.9 1.5 0.0 1.1 0.3 0.0 4.9 4.2 3.2
## [18] 2.7 1.8 0.4

mean1 <- mean(sample1) # 標本平均
mean1

## [1] 1.785

var(sample1)

## [1] 2.097132

20人無作為抽出して標本平均を計算すると、母平均とピッタリとは一致しない。もう一度サンプリングして標本平均と標本分散を計算すると以下のようになる。

sample2 <- sample(population, 20) # 上の母集団から 20 人無作為抽出
sample2

##  [1] 1.8 5.2 0.1 4.0 3.7 7.9 3.0 1.8 1.8 2.6 3.9 1.8 2.1 2.5 1.9 4.0 0.7
## [18] 1.6 2.8 1.6

mean2 <- mean(sample2) # 標本平均
mean2

## [1] 2.74

var(sample2)

## [1] 2.981474

やはり、サンプリングの際の偶然によって、違った平均と分散が得られる。中心極限定理によれば、このようなサンプリングを何回も繰り返すと標本平均の平均は母平均（上の例では 2.98）に近づいていくはずであり、標本平均の分散は \(N/\sigma^2\)（上の例では 8.95 \(/20=\) 0.45にちかづいていくはずである。やってみよう。

sample.mean <- rep(NA, 2000) # 2000回サンプリングを繰り返した時の標本平均を記録するためのベクトル
for(i in 1:2000){
  sample3 <- sample(population, 20)
  sample.mean[i] <- mean(sample3)
}
sample.mean[1:10] # 20回分の標本平均だけ表示

##  [1] 2.955 3.125 4.925 3.710 3.490 4.000 2.530 2.575 2.920 2.710

mean(sample.mean)

## [1] 2.983808

var(sample.mean)

## [1] 0.4555384

pop.var/20

## [1] 0.4476377

hist(sample.mean, breaks=(4:25)/4)

標本平均の平均と分散は概ね中心極限定理の予測通りの値になる。正規分布になっているかどうかは、ヒストグラムは正規分布らしいかたちになっている。詳細は省略するが、概ね正規分布である。下は 分位点－分位点プロット(Quantile Quantile Plot)。

qqnorm(scale(sample.mean))
abline(a=0, b=1, col="red")

4 Practice: 練習問題

上の2種類のシミュレーションを自分でも実行し、概ね同じ結果になるかどうか確認せよ。