中心極限定理を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 #１回の試行での標本数。
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")

R言語のcut関数の使い方

R初心者として、cut関数がいまいち分かりにくかったので、ここで少しまとめておく。

■cut関数は何をする関数？

一言でいうと「数値データを、指定した分割基準でカテゴリに変換する関数」だ。

もう少し詳しくいうと、例えば英語の試験を行い各人の試験結果を

> x <- c(90, 55, 79, 80, 100)

とする。80点未満を「不合格」、80点以上を「合格」と分けるとすると、xの変数の内容を「合格、不合格、不合格, 合格、合格」と変換したfactor型のオブジェクトを返すのがcut関数だ。

■実際にcut関数を動かしてみる。

上記の例を実際にR上で行ったのが以下。

> x <- c(90, 55, 79, 80, 100)
> cut(x,breaks=c(0,80,100), labels=c("不合格","合格"), right = FALSE, include.lowest = TRUE)

[1] 合格   不合格 不合格 合格   合格  
Levels: 不合格 合格

点数が先ほどの基準に従ってカテゴリ（合格、不合格）に変換されているのが分かる。

■right = FALSEのオプションは？
先のcut使用例で「right=FALSE」のオプションをしている。この説明をしなければならない。
cutのデフォルト動作は、例えばbreaks=c(0,80,100)と指定した場合、

• (0,80]の区間を不合格
• (80,100]の区間を合格

とカテゴリ分けする。
ここで丸括弧は開区間、角括弧は閉区間を示す。でも今回は
80点未満を不合格、80点以上を合格としたいから、

• [0,80)の区間を不合格
• [80,100)の区間を合格 （注：ここで100がカテゴリ分けに入らないのは後述）

であるべき。このように区間を分けるためにright = FALSE とする。

■include.lowest = TRUEのオプションは？
先のright=FALSEだけでは「100」がカテゴリ分けに含まれない。そこでカテゴリ分けの末端の開区間を閉区間にし、今回の例では100をカテゴリに含むようにするのが、このオプション。

「right = FALSE」「include.lowest = TRUE」を指定することで、今回我々が望む

• [0,80)の区間を不合格
• [80,100]の区間を合格

というカテゴリ分けが実現できる。