乱数の生成について

Pocket

こんにちは。データ解析部の西岡(敏)と申します。 広告配信の効果分析とアルゴリズム実装を担当しております。

今回は、システムで利用する乱数の生成について書きたいと思います。

乱数について

過去、コンピュータの無い時代には、サイコロを振ってその数字を記録したり、コインを投げてその裏表を調べたりすることで乱数を得ることができました。 ただし、こういった乱数は大量のデータを必要とする実験では不便であるため、現在はそのような実験を行う際にはコンピュータを利用して擬似的に乱数を生成しています。

こういった形でコンピュータを用いて生成する乱数のことを、前述のサイコロの乱数などとは区別し、算術乱数や擬似乱数と呼びます。

プログラムで使われる乱数

前述のように、疑似乱数はコンピュータから生成を行うため、乱数が本来持っているべき予測不可能性を持ちません。 ですが、出来る限りランダムであることが担保され、非常に長い周期を持てば、シミュレーションなどに利用するのに差し支えはないということになります。

疑似乱数は主にコンピュータ上で生成・使用され、そのために生成される乱数の質や周期の長さと同等に、コンピュータ上で生成される速度、プログラミングの簡素さ、使用する記憶領域の大きさが乱数生成器の評価のポイントとなっています。

以下では、システムで広く使われている代表的なアルゴリズムである、 線形合同法と、メルセンヌ・ツイスタの説明を行っていきます。

線形合同法

概要

1948年頃にレーマーが提案した方法で、今日でも広く使われている手法になります。 BASICやPascal, UNIXなどでこちらの方法が乱数生成に使われています。

基本的な算出方法は、 適当なX(0), 乗数a, 加数b, 法(modulus)Mを指定し、

X(n) = (a * X(n-1) + b) % M として、計算を行うものになります。

このとき、 bがMと互いに素 a-1がMを割り切るすべての素数の倍数 Mが4の倍数 であれば

a-1も4の倍数

という条件が、X(n)の周期が最大値Mとなる必要十分条件になります。

pythonによるプログラミング

実際にpythonで線形合同法を用いて乱数の生成を行ってみます。

線形合同法の問題点

上記のように、簡単に生成ができる一方で、この手法には 周期の限界が2^32にあること、下位ビットには周期的な偏りがあること(上記の例でも、下一桁が奇数→偶数の繰り返しになっている) などの問題点があります。 そのため、下位ビット部分を切り捨てて実装を行っている場合もあるそうです。

メルセンヌ・ツイスタ

概要

続いては、今日様々なシステムで利用されているメルセンヌ・ツイスタに関して説明を行います。 メルセンヌ・ツイスタはM系列を松本眞氏と西村拓士氏が発展させたものであり、疑似乱数生成器として現在最も高い評価を得ています。

メルセンヌ・ツイスタの説明

メルセンヌ・ツイスタの説明に当たって、まずM系列の説明を行います。

M系列

M系列とは、Xn = Xn-p + Xn-qで表現できる1ビットの数列になります。

例えば、 X0 = 1, X1 = 0, X2 = 0, p = 3, q = 1 としたとき、

X3 = X0 + X2 = 1 + 0 = 1
X4 = X1 + X3 = 0 + 1 = 1
X5 = X2 + X4 = 0 + 1 = 1
X6 = X3 + X5 = 1 + 1 = 0

という計算によって算出できます。

TGFSR

M系列を発展させたものに、 Twisted General Feedback Shift Register(TGFSR)があり、こちらは

Xn+p = Xn+q + Xn *A (p>q>0)

を用いて表現することができます。 ここで、Aはω✕ωの行列、Xはω次元の行ベクトルになります。 また、このときωは32, 64など、生成する乱数のビット数となります。

メルセンヌ・ツイスタ

さらにTGFSRを

Xn+p = Xn+q + Xn+1 *B + Xn *C (p>q>0)

として改良したものがメルセンヌ・ツイスタになります。

特に、こちらのこちらの最大周期が2^19937 – 1のときにMT19937と呼ばれ、 通常のメルセンヌ・ツイスタはこのMT19937を指します。

具体的な計算方法としては、 init_genrand関数とgenrand_int関数を準備し、 init_genrand関数を実行した後に、生成する乱数分genrand_int関数を実行するという流れになります。 それぞれの関数は以下の形式となります。

  • init_genrand 関数
    • 0 番地に種を代入する。
    • i 番地の値を得るために、i-1 番地の値を使用する。シフト、EXOR、乗算、加算の順に演算を行う。
    • 623 番地まで計算を繰り返し、値を代入する。
  • genrand_int 関数 a. (624n+1)回目の場合(n=0,1,2,…)
    • 0 番地から 227 番地まで AND、OR、EXOR 演算を行う。
    • 228 番地から 623 番地まで、別の数値を使って AND、OR、EXOR 演算を行う。
    • 再び別の数値を使って、AND、OR、EXOR 演算を行い、y に代入する。 b. 毎回共通の部分
    • y にシフト、AND、EXOR 演算を行い、生成された乱数を返す。

pythonによるプログラミング

実際にpythonで線形合同法を用いて乱数の生成を行ってみます。

以上のような形で、乱数の生成を行うことができました。 このメルセンヌ・ツイスタに関しては、線形合同法であったような、 周期が短かったり、値によっては下一桁が偶数, 奇数と交互に出現するといった問題点が発生せず、 非常に優れた擬似乱数生成器として、様々なプログラミング言語に幅広く利用されています。

まとめ

このように乱数と言っても複数の種類があり、それぞれの乱数には特徴があります。 現在は、高い精度の検証をする場合などは、線形合同法で作成されたよりもメルセンヌ・ツイスタの方が良いと言われています。 もし何かしらの実験を行っているときなどに、こういった乱数に着目すると、意外な解決策が見いだせるかもしれません。

参考

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/TEACH/ichimura-sho-koen.pdf http://www.soi.wide.ad.jp/class/20010000/slides/03/sfc2002.pdf http://random.ism.ac.jp/info01/random_number_generation/node6.html http://www.utp.or.jp/book/b300859.html http://www001.upp.so-net.ne.jp/isaku/rand.html https://narusejun.com/archives/5/

Windows Mixed RealityのadstirSDK研究開発について
HadoopからSparkへの移行