はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

トップ > 乱数 > 乱数発生方法と線形合同法のメモ

乱数発生方法と線形合同法のメモ

乱数 math

乱数発生方法と線形合同法のメモ

メルセンヌ・ツイスタの前知識として，乱数発生方法と線形合同法のメモ

乱数発生方法

乱数発生法：「一様かつ独立に次々とでたらめな数を発生させる方法」

例：

  • サイコロ：{1, 2, 3, 4, 5, 6}の出目を発生
  • コイン投げ：表と裏で{0, 1}の要素を発生
  • ルーレット：{0, 1, . . . , 36}の要素を発生
  • くじ引き

これらの乱数の並びを乱数列

物理乱数

乱数を、どうやって生成するか？

“決定的な動作しかしないコンピュータでは、乱数は生成できない”

コンピュータを使わない物理乱数発生器:

  • 物理雑音から拾ってくる: もっとも素朴な方法
  • 例：サイコロ、熱雑音、株価の変動
  • 問題点: コスト、スピード、再現不能性

参考： 講演 あなたの使っている乱数、大丈夫？-危ない標準乱数と、メルセンヌ・ツイスター開発秘話,p.6,2014

物理乱数発生器の問題点：再現不能性

再現不能性：同じ乱数列を再現するのに、それらをすべて記録しておかないとならないこと．

再現性が必要となる場合： - 追試 - 最適化

例：核シミュレーションでは乱数は何兆個も消費される

⇒ それらをすべて記録しておくのは効率が悪い。

参考： 講演 あなたの使っている乱数、大丈夫？-危ない標準乱数と、メルセンヌ・ツイスター開発秘話,p.7,2014

疑似乱数

計算で乱数に見えるもの（疑似乱数）をつくろう.

擬似乱数発生法：漸化式を用いて、乱数のように見える数列を生成する方法．

擬似乱数のメリット:

  • 漸化式と初期シードを記録しておけば、誰でも同じ数列を再現できる
  • 高速で低コスト

問題点：

  • 「乱数と呼んでいいのか」

    • 擬似乱数の創始者von Neumann「漸化式で乱数を作るのはある種の罪」“Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin”

  • 周期性が存在する:周期Nをもっている生成器は，循環することになる． 周期6でも100でも$O(10^{6001})$でも常に成立する．

疑似乱数の出力器を疑似乱数生成器（Pseudo random number generator, PRNG）

=> PRNGってどんなもの？

PRNG

  1. 乱数のseedを与える
  2. PRNG: seed値に対して処理
  3. 疑似乱数を生成（出力）

PRNGの種類

PRNGの主なもの

線形合同法(LCG)

線形合同法(LCG, Lehmer 1960)

  • ある整数$X_{1}$を初期シードとして選ぶ。
  • 次の例のような漸化式により，$X_2, X_3, \cdots$ を次々に生成
X_{n+1}=\left(A \times X_{n}+ C \right)\bmod M , (n=1,2,\cdots)

$A,C,M$は任意の定数 - 乗数 $A$ 及び法 $M$ は正の整数であり， - 加数 $C$ は非負の整数

メリット： - 周期や分布を求めるアルゴリズムがある

デメリット： - M系列法などと比較すると周期が短い - 周期を大きくするにはM を大きくするしかなく、掛け算や割り算が必要で生成が遅くなる - 線形合同法を「 mod 偶数」で使った場合,必ず奇数と偶数が交互に出る。

参考リンク：

実装例：

線形合同法は、70年代から80年代にかけて，ANSI-Cなどの標準擬似乱数rand()であった．いまでも教科書にのっていて、広く使われている。

  • この数列の周期は、初期シードの選び方によらず$232$
  • 現代のパソコンは数分で$2^{32}$個の乱数を使ってしまう
  • 生成される数列はかなり乱数っぽく見えるが，数千万個の出力を使うと、非乱数性が現れてくる

例 A=3, C=5, M=13

X_{n+1}=\left(3X_{n}+ 5 \right)\bmod 13 , (n=1,2,\cdots)

3回で周期がくる

A = 3
C = 5
M = 13
seed = 8

x = seed
for i in range(0, 10):
    x = (A * x + C )%M
    print(i, x)

0 3
1 1
2 8
3 3
4 1
5 8
6 3
7 1
8 8
9 3

例：A=3, C=1, M=8

X_{n+1}=\left(3X_{n}+ 1 \right)\bmod 8 , (n=1,2,\cdots)

周期:4回

Mを偶数にしているので，奇数偶数の繰り返しになる．

A = 3
C = 1
M = 8
seed = 6

x = seed
for i in range(0, 15):
    x = (A * x + C )%M
    print(i, x)

0 3
1 2
2 7
3 6
4 3
5 2
6 7
7 6
8 3
9 2
10 7
11 6
12 3
13 2
14 7

例：32bit LCG, 64bit LCG

32bitLCG

周期：$2^{32}$ 10桁

X_{n+1} = (X_{n} * 0x41C64E6D + 0x6073) \bmod 0xffffffff 
A = 0x41C64E6D #32bit
C = 0x6073
M = 0xFFFFFFFF
seed = 0x9FF1E41D

x = seed
for i in range(0, 15):
    x = (A * x + C )%M
    print(i, "{:x}".format(x))

0 c205a7cd
1 5cf10427
2 36ee37fe
3 d30ea4e3
4 b445f4f3
5 67af1e71
6 d7e08c7f
7 761c2201
8 8875cb73
9 af26a7fb
10 2ec6de87
11 f1e8fc92
12 a63bf30d
13 7f58d0fd
14 49e4a4f1

64bit LCG

周期 64bit LCG: $2^{64}$ 20桁

X_{n+1} = (X_{n} * 0x5D588B656C078965 + 0x269EC3) \bmod 0xffffffffffffffff 
A = 0x5D588B656C078965 #64bit
C = 0x269EC3
M = 0xFFFFFFFFFFFFFFFF
seed = 0x9FF1E41D

x = seed
for i in range(0, 15):
    x = (A * x + C )%M
    print(i, "{:x}".format(x))

0 a7c1802240c65550
1 906abd90cb41218f
2 9cffec733a4f9ba5
3 e02466f6c5c3bbe6
4 99ae6e868cbe7e84
5 6348eedc16af88aa
6 b849c06757738304
7 210bb2d539308fe
8 2f39ac81c5f6d193
9 970ed5826793ee1f
10 e7c4ad8749c63f96
11 f118a921a16c1398
12 cbb3580b72efef56
13 220f50a14755a6a
14 cfa6129ab294cc45

パラメータ(A,C,M)の決め方

JIS9031「乱数生成及びランダム化の手順 Procedure for random number generation and randomization」 5.3.2で示されている．

$ A，M $及び$ C $の値を勝手に定めたのでは良い疑似乱数列が得られないことが多いので，次のような基準に従って定めることが望ましい。

  • 法$ M $は，得られる数列の周期の上界となるので，できるだけ大きくするのがよい。1 語が 32ビットからなる計算機を使用する場合には，$M＝2^{32}$又は$M＝2^{31}−1$ とすることが望ましい。

  • 乗$A$ としては，スペクトル検定で良好な成績を修めたものを使うことが望ましい。

  • 加数$ C $ の選び方については，特に厳しい基準はない。ただし，これを 0 とするか，又は正の整数とするかによって，生成される数列の周期が異なる場合がある。

    • $ M $ が2の整数べき（冪）乗の場合，$C＝0$ とすると，周期は m/2 以下である。$ C $ を奇数，$A$ を 4で除して 1 余る数，周期は $ M $ となる。