(正規)分布の確率密度関数を乱数から推定する

すごい基本的なお話な気がするけど、数回は書き直しているのでメモっておく。ここでは

  • 正規分布に従う乱数を生成して、それがある幅(dx)の中に入っているか否かをカウントして確率にする
  • 正規分布確率密度関数に幅(dx)をかけて確率を計算する

の2つの方法を比較して、それらが一致することを確認する。コードは以下。

#乱数のシードの指定
set.seed(100)
#確率密度関数が一定とみなせるような適当な幅
dx <- 0.1
#標準偏差3の乱数生成
x <- rnorm(10^6, sd=3)
#5 ± dxに入っている|ないを計算
hit <- sapply(x, function(z)abs(z-5) < dx)
#確率を計算
sum(hit)/length(hit)
#幅×確率密度関数
(2*dx)*dnorm(5, sd=3)

最後の部分だけ実行結果を載せると

> sum(hit)/length(hit)
[1] 0.006558
> (2*dx)*dnorm(5, sd=3)
[1] 0.006631809

となって、乱数から生成された確率が、正規分布確率密度関数で計算したものとほぼ同じ値になっていることがわかる。