[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
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 4 months ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 5 months ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 5 months ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 10 months ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 1 year ago
%d bloggers like this: