情報科学と人工知能のノート

初等的な知識から最新論文の解説まで色々集めていきます.備忘録兼用.

R 備忘録 2013/05/10-3

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

平均 \mu、分散 \sigma の分布に従う独立した  n 個の確率変数  X_1, \cdots, X_n について  \bar{X} = \frac{1}{n} \sum_{i=1}^n X_i と定義する。
 n が大きいとき、 \bar{X} は平均 \mu、分散 \frac{\sigma^2}{n}正規分布に従う。

とのことです。つまり  \bar{X} をたくさん作ってヒストグラムを作ったら正規分布の密度関数の曲線が現れるわけです。
0-1 一様分布で n=10 でシミュレーションしてみます。

まず、一様乱数の平均値を 10000 個ほど生成します。

> vec <- mapply(function (n) return (mean(runif(10))), 1:10000)
> summary(vec)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2085 0.4391 0.5017 0.5010 0.5625 0.8141

最小値と最大値を考慮して、わかり易さ・やり易さのために [0:1] 区間を使うことにして、ヒストグラムを作ります。

> hist(vec, breaks=seq(0, 1, 0.01), probability=T, ylim=c(0, 5))

probability=TRUE を指定すると密度関数のように(おそらく、面積が 1 になるように)高さの単位を調節してくれます。
ちなみに ylim=c(0, 5) の 5 はヒストグラムを見て後から自分で調節しました。

0-1 一様分布の平均は 0.5、分散は 1/12 なので、n=10 のとき、\bar{X} は平均 0.5、分散 1/120 の正規分布に従います。

> par(new=T)
> plot(function(x) return (dnorm(x, mean=0.5, sd=sqrt(1/120))), xlim=c(0, 1), ylim=c(0, 5), xlab="", ylab="")

f:id:ti2236:20130510171903p:plain
それらしいグラフが表示されました。
ちなみに par はグラフィックスパラメータを取得/操作するための関数で、属性 new が TRUE のとき、新しい図が前の図に上書きされます。
new は図を描くたびに FALSE に戻るため、上書きするたびに設定してあげる必要があります。

分布の 100p パーセント点を求める

前回出てきた qDIST 系関数を使えば簡単です。

> qnorm(0.95) # 標準正規分布の95 パーセント点
[1] 1.644854
> qnorm(0.05) # 5 パーセント点
[1] -1.644854