線形合同式法(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
%d bloggers like this: