[R]疑似乱数生成(Poisson分布)

キーポイント確率統計 (理工系数学のキーポイント)」という本の「ポイント8確率過程とシミュレーション」 式(8.37)(P。139) で説明されている 疑似Poisson乱数の生成を R で実装。実際に手を動かしたほうが理解も深まるしね、、、

一様乱数法

poisson_random_number <- function(lambda) {
  base <- exp(lambda)
  n <- 0
  repeat {
    base <- base * runif(1)
    if (base <= 1) {
      return (n)
    }
    n <- n + 1
  }
}

wikipedia の”Poisson distribution:Generating Poisson-distributed random variables”によると、このアルゴリズムは Knuth が開発したらしい。

ヒストグラムで分布を確認
ビルトインのPoisson分布と手製Poisson分布の違いを10000個生成してヒストグラムで比較

まとめて実行
乱数の生成からグラフの作成まで。

$ cat pseudo-poisson-distribution.R
#!/usr/bin/Rscript

print ('built-in pisson distribution')
rpois(10, 2)

poisson_random_number <- function(lambda) {
  base <- exp(lambda)
  n <- 0
  repeat {
    base <- base * runif(1)
    if (base <= 1) {
      return (n)
    }
    n <- n + 1
  }
}

print ('knuth poisson distribution generation')
for (i in 1:10) {
  print (poisson_random_number(2))
}

# make graph:
png(file="poisson-distribution.png")
par(mfrow=c(2, 1))
hist(rpois(10000, 2),                            main='built-in poisson distribution', br=100)
hist(replicate(10000, poisson_random_number(2)), main='Knuth poisson distribution',    br=100)
dev.off()
# make graph:END

★実行結果★

$ ./pseudo-poisson-distribution.R
[1] "built-in pisson distribution"
 [1] 3 3 2 2 2 2 2 3 6 1
[1] "knuth poisson distribution generation"
[1] 2
[1] 1
[1] 0
[1] 3
[1] 4
[1] 1
[1] 1
[1] 3
[1] 3
[1] 3
null device
          1
Advertisements
Tagged with: , , , ,
Posted in algorithm, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
%d bloggers like this: