2013年5月26日日曜日

中心極限定理をRでシミュレーション

中心極限定理とは、大まかに言えば、
母集団の確率変数\(X\)がどんな確率分布であっても、その平均と分散が\(\mu\)と\(\sigma^2\)であれば、その標本(サンプル数\(n\))の確率変数\(X\)の和や平均は、それぞれ正規分布 \(N(n\mu, n\sigma^2)\)、 \(N(\mu, \sigma^2/n)\)に従うということである。

それが本当かどうか確かめたかったので、Rのプログラムを書いて検証してみた。
ここでは、「指数分布に従う確率変数の平均」で検証した。

以下にコードを置いておく。
実際にCONSTANTSの部分をいろいろ変えて
  • 分布が正規分布に従い、その平均と分散が中心極限定理の予言通りか?
  • 標本数が大きくなれば、分散は実際に小さくなっていくのか?
などを確認してみるとよい。


#中心極限定理が成り立つかをシミュレーションするRスクリプト
#ここでは、指数分布の確率分布に従う確率変数Xの平均値について検証。

## CONSTANTS
SAMPLE.NUM <- 10 #1回の試行での標本数。
TEST.NUM <- 50000 #試行回数(この数の試行を繰り返して分布が正規分布になるかをみる。)
EXP.LAMBDA <- 3 # 指数関数パラメータλ
PLOT.MAX.X <- 1 #こ

#各行が指数分布に従う確率変数の一つの標本を表わす行列を生成。
data.matrix <- matrix(rexp(TEST.NUM, rate = EXP.LAMBDA), ncol = SAMPLE.NUM)

#各試行(= 各行)の平均をとる。
means.of.each.test <- apply(data.matrix,MARGIN=1,mean)

#各試行の平均の、相対度数ヒストグラムを描画
hist(means.of.each.test, breaks=seq(0, PLOT.MAX.X, by=0.05), xlim=c(0, PLOT.MAX.X), probability = TRUE)

#中心極限定理が予言する平均と分散を求める。
#ここで母集団分布がλ=EXP.LAMBDAの指数分布であり、
#その平均と分散(母平均&母分散)がそれぞれ、
#母平均=1/λ、母分散=1/λ^2なので・・・
theorem.mean <- 1/EXP.LAMBDA
print(theorem.mean)
theorem.var <- 1/EXP.LAMBDA^2/SAMPLE.NUM
print(theorem.var)

#中心極限定理が予言する正規分布を描画する。
#ヒストグラムと正規分布曲線が一致すれば、中心極限定理は正しい!
x <- seq(0, PLOT.MAX.X, by=0.01)
lines(x, dnorm(x, theorem.mean, sqrt(theorem.var)), col="red")

0 件のコメント:

コメントを投稿