Visual Basicサンプル - 合同法による一様乱数の生成

乗算合同法による一様乱数

Ran0、Ran1、及びRan2と3種類のクラスがありますが、いずれも

NUMERICAL RECIPES in C　(日本語版)
William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
技術評論社
ISBN4-87408-560-1

に掲載されている関数を移植したもので、31bit整数の乱数を内部的に生成し、これを基に0～1のDouble型の一様乱数を生成するメソッドを備えています。(Double型で乱数が取り出されますが、内部的には高々31bitの精度しかありません。)

Ran0は、ParkとMillerが、乗算合同法でも線形合同法と同程度にランダムな乱数列を生成できることを示した例です。Schrageの技法により符号付き32bit整数型(VB.NetではInteger型、前のバージョンではLong型)でもオーバーフローしないように計算方法が工夫されています。正確な生成範囲は(0,1]で、周期は約2×10⁹です。
Ran1は、Ran0を改良し、予めある程度の乱数列を発生させて配列に格納しておき、引き続いて生成される乱数に基づいて配列の要素を選び出力し、要素そのものは置き換えていく方法で、格段にランダムさが向上しています。数千万個の乱数列でよければRan1が最も速度と統計的な性能のバランスに優れてるようです。生成範囲は(0,1)で、周期はRan0とほぼ同じです。
Ran2は、Ran0と同等の微妙に周期が異なる2つの乱数列を組み合わせ、各々の周期のほぼ2乗の周期を実現しています。合同法に基づく乱数生成法としては、最も優れているアルゴリズムのようで、実用上の要求をほぼ満足します。生成範囲は(0,1)で、周期は約2.3×10¹⁸です。

改定日	内　容
2007/06/08	コンパイルオプションによって、内部で保持されている31bit整数の乱数の上位23bitを、直接Single型変数の仮数部にセットできる処理も選べるように、サンプルを変更。ファイル名も変更。
2007/06/19	乱数をDouble型で返すように変更。ただし、精度は高々31bitのままである。初期値を省略して初期化を行った場合、コンパイルオプションによっては、内臓時計の値に基づいて初期値を算出するように変更。繰り返し同じ系列の乱数を生成できるように、初期値を返すプロパティーを追加。

応用にあたっての留意・注意事項

大方、出典のNumerical Recipesに書かれている内容の繰り返しになりますが、Visual Basicの特徴も配慮して、応用の際に留意あるいは注意して頂きたいことを掲げておきます。擬似乱数については、様々な誤解の流布や、個々の開発者の思い込みによって、誤用が多い分野ですので、参考にして頂ければ幸いです。

最初にも書きましたが、乱数は「予想不能」が良いとは限りません。統計的な振る舞いをシミュレーションする応用面で、かつ、様々なプラットフォームで共通に動かしたい場合は、どのプラットフォームでも、他の言語で記述しても、正確に同じ乱数列を再現できるアルゴリズムが求められます。その点で、ここに掲げる関数は、単純であり、内容を理解しやすく、移植性に優れたものです。
クラスの内部でInteger型変数で保持されている乱数の下位ビットを取り出して利用することはできません。下位ビットはランダムさが十分ではないためです。これは合同法の弱点であり最も注意すべき点です。一部のビットを取り出すならば、上位ビットから取り出さなければなりません。
複数の乱数を続けて発生させて「組」にして用いる場合には、Ran2かRan1を使ってください。Ran0は使えません。これは、Ran0は十分にランダムではなく、隣り合う乱数同士が強い相関を持つからです。例えば、xとyの座標値を乱数で与えて、平面状にランダムに点を分散させたい場合が該当します。
上述の2つの注意とも関連しますが、乱数を2つ発生させて、bit列同士を繋げて長いbit長の乱数として用いることは避けてください。Ran2やRan1であれば、2つの乱数の相関は十分に小さくてランダムと見なせるかもしれませんが、各々のbit列の下位ビットの方はあまりランダムではないかもしれません。どうしても、ということであれば、上位ビットからbit列を取り出して繋げるべきです。

インラインで実現できる32bit線形合同法(.Net 2.0のみ可能)

32bit符号無しの整数型(UInteger型あるいはUint32型)を用いて、かつ、整数のオーバーフローチェックを無効にした場合、以下のコードで、統計的な振る舞いは良くないですが非常に高速に乱数列を得られます。

Dim idum As UInteger ' = 何か初期値を代入する。0でもよい。
'何らかの処理のループ
    idum = 1664525UI * idum + 1013904223UI
    
    'idumの値を乱数として利用する。
    'idumの値はループ内の他のところで変更しないこと。
'ループの終わり

上述のidumの値を基にSingle型の浮動小数点数の乱数を得たい場合には、可読性には欠けますが、以下のような構造体のトリックを使って、やはりインラインで変換することができます。まずは、構造体SingleUnionの宣言を示します。

'構造体のメンバー変数としてSingle型とUInteger型(共に同じ
'4Byte長)をメモリ上の同じ位置を共有させることができる。
'CやC++で言うところの共用体(union)と同等である。

'この機能は.Net Frameworkの以下の名前空間で提供されている。
Imports System.Runtime.InteropServices


'属性を用いてメンバー変数のメモリ上の位置を細かく指定できる。
'詳しくは、.Net Frameworkのドキュメントを参照のこと。
<StructLayout(LayoutKind.Explicit)> _
Public Structure SingleUnion
    <FieldOffset(0)> Public S As Single
    <FieldOffset(0)> Public UI As UInteger
End Structure

この構造体を用いて、UInteger型の乱数値をSingle型に変換します。UInteger型でビット操作を行いSingle型の内部表現に変換しますが、この操作がそのままSingle型のメンバー変数の値に反映されます。なお、IEEEに準拠するSingle型にのみ適用可能ですのでご注意ください。

'変換操作を行うSingleUnionの変数を宣言する。
Dim r As SingleUnion


'上述のidumを算出する線形合同法の式
idum = 1664525UI * idum + 1013904223UI


'UInteger型のメンバー変数を9bit右シフトすることにより、idumの
'上位23bitがSingle型の仮数部に置かれる。さらに、指数部に1以上、
'2未満の数値を示す値をセットする。
r.UI = &H3F800000UI Or (idum >> 9)

'この時点で、r.Sは[1,2)の範囲の乱数になるので、1だけ減じて
'[0,1)に収まるようにする。
r.S -= 1.0!

'あとは、r.Sを利用する。

最初に掲げたサンプルのソースコードでコンパイルオプションを選択することにより、ここに示したものと同じ(ただし、31bitの整数を8bit右シフトする)方法を、最初に掲げたRan0、Ran1、及びRan2についても適用することができます。