線形合同式法(Linear Congruence Method)による疑似乱数生成

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

R バージョン

congruence_method <- function(a, c, m, r) {
  return (function() {
    r <<- (a * r + c) %% m
    return (r / m)
  })
}

meth <- congruence_method(1229, 351750, 1664501, 1.0)

for (i in 1:10) {
  print (meth())
}

★実行結果★

$ Rscript  rand.R
[1] 0.2120630
[1] 0.8366958
[1] 0.5104659
[1] 0.573915
[1] 0.5527921
[1] 0.5928131
[1] 0.7786261
[1] 0.1428590
[1] 0.7850737
[1] 0.06693237
”<<-” の使い方は “An Introduction to R” の ”10.7 Scope” にテクニカルな説明がある。

Python バージョン

def congruence_method(a, c, m, r):
    while True:
        r = (a * r +c) % m
        yield r / m

gen = congruence_method(1229, 351750, 1664501, 1.)
for i in range(10):
    print gen.next()

★実行結果★

$ python rand.py
0.212062954603
0.836695802526
0.51046589939
0.57391494508
0.552792098052
0.592813101344
0.778626146815
0.142859031025
0.785073724798
0.0669323719241

ヒストグラムで一様性を確認

1万回実行した結果を確認するには、R から次のコマンドを実行する。

hist(replicate(10000, meth()))

出力結果

Advertisements
Tagged with: , , , ,
Posted in algorithm, python, 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 6 months ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 6 months ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 6 months ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 1 year 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: