ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

乱数ストリームの作成と管理

RandStream クラスを使用すると、より柔軟に乱数ストリームを生成できます。たとえば、グローバル ストリームの状態に影響を与えずに、乱数値を生成することや、シミュレーションで使う乱数のソースを分けることや、MATLAB® 起動時の乱数発生器のアルゴリズムではなく、別のアルゴリズムを使用することが可能になります。RandStream コンストラクターを使用すると、独自のストリームを作成し、書き込み可能なプロパティを設定して乱数を生成できます。作成したストリームは、グローバル ストリームと同じ方法で制御することができます。また、グローバル ストリームを作成したストリームと置き換えることもできます。

ストリームを作成するには、以下のように RandStream コンストラクターを使用します。

myStream=RandStream('mlfg6331_64');
rand(myStream,1,5)

ans =

    0.6530    0.8147    0.7167    0.8615    0.0764

乱数ストリーム myStream はグローバル ストリームとは別に動作します。関数 rand、関数 randn、関数 randi は、グローバル ストリームから乱数を引き出し、myStream に適用した RandStream メソッド、関数 rand、関数 randn、関数 randi の結果に影響を与えません。

RandStream.setGlobalStream メソッドを使用すると、myStream をグローバル ストリームにできます。

RandStream.setGlobalStream(myStream)
RandStream.getGlobalStream

ans = 

mlfg6331_64 random stream (current global stream)
             Seed: 0
  NormalTransform: Ziggurat

RandStream.getGlobalStream==myStream

ans =

     1

サブストリーム

以前に実行したシミュレーション部分を返す場合を考えます。乱数ストリームは固定したチェックポイントにジャンプして制御することができ、Substream プロパティを使うと、複数のサブストリーム間で前後にジャンプすることができます。Substream プロパティを使用するには、サブストリームをサポートする発生器でストリームを作成します(発生器のアルゴリズムとプロパティのリストについては、乱数発生器の選択を参照してください)。

stream=RandStream('mlfg6331_64');
RandStream.setGlobalStream(stream)

Substream の初期値は 1 です。

stream.Substream

ans =

     1

サブストリームは逐次計算を行う場合に便利です。サブストリームは指定したチェックポイントに戻って、シミュレーションのすべてまたは一部分を作成しなおすことができます。たとえば以下のように、ループ内でサブストリームを使用できます。

for i=1:5
    stream.Substream=i;
    rand(1,i)
end

ans =
    0.6530

ans =
    0.3364    0.8265

ans =
    0.9539    0.6446    0.4913

ans =
    0.0244    0.5134    0.6305    0.6534

ans =
    0.3323    0.9296    0.5767    0.1233    0.6934

各サブストリームはループで反復した処理を再生成できます。たとえば 5 番目のサブストリームを返します。結果は上記 5 番目の出力と等価です。

stream.Substream=5;
rand(1,5)

ans =

    0.3323    0.9296    0.5767    0.1233    0.6934

乱数発生器の選択

MATLAB には発生器のアルゴリズムが 6 つ用意されています。次の表は、使用可能な発生器アルゴリズムの主要な性質と、発生器の選択に使うキーワードをまとめています。使用可能な発生器アルゴリズムのすべてを返すには、RandStream.list メソッドを使用します。

発生器アルゴリズム

キーワード発生器複数のストリームとサブストリームのサポート 高精度な場合の近似周期
mt19937arメルセンヌ・ツイスター (既定)なし2199371
dsfmt19937SIMD 指向高速メルセンヌ・ツイスター なし2199371
mcg16807乗算合同法発生器なし2312
mlfg6331_64乗法ラグ フィボナッチ発生器あり2124
mrg32k3a複数の再帰発生器の組み合わせあり2191
shr3cong線形合同法発生器を組み合わせたシフトレジスタ発生器なし264
swb2712修正桁下げ付き減算発生器なし21492

いくつかの発生器 (mcg16807shr3congswb2712) は MATLAB の前のバージョンと下位互換性があります。2 つの発生器 (mrg32k3amlfg6331_64) は明示的に並列乱数発生をサポートしています。残りの発生器 (mt19937ardsfmt19937) はシーケンシャルな用途のために主に設計されています。発生器がより高速に動作するか、より高精度な結果が得られるかは用途にも依存します。

どの発生器を選択すべきかは用途と一緒に考える必要があります。すべての疑似乱数発生器は決定性アルゴリズムに基づいているため、厳密な意味では真の乱数ではありません。よってモンテ カルロ シミュレーションの結果をチェックする 1 つの方法は、2 つ以上の別の発生器アルゴリズムを使ってシミュレーションを実行することです。MATLAB では発生器の選択ができるため、この方法が適用できます。他の発生器を使った場合の結果が、モンテ カルロ法のサンプリング誤差より大きい場合がまれにあります。これに関しては、特定の発生器アルゴリズムではこの種の検証に欠陥があることを説明した文献に例があります (例は、[13] を参照)。

発生器アルゴリズム

mcg16807

[14] で説明されている 32 ビットの乗算合同法発生器です。乗数は a=75、モジュロは m=2311 です。この発生器は 2312 の周期をもち、複数のストリームやサブストリームをサポートしていません。U(0,1) の各値は、32 ビット整数 1 つを使って発生器から作成されます。取り得る値はすべて (2311)1 の倍数で、厳密に区間 (0,1) 内にある値です。mcg16807 ストリームを既定で使用する randn アルゴリズムは、polar アルゴリズムです ([1]を参照)。メモ: この発生器は MATLAB Version 4 の初期に、関数 rand と関数 randn が使用していた発生器と同じであることに注意してください (rand('seed',s) または randn('seed',s) で実行されていました)。

mlfg6331_64

[10] で説明されている 64 ビットの乗法ラグ フィボナッチ発生器で、ラグは l=63k=31 です。この発生器は SPRNG パッケージに実装されている MLFG と同じです。近似周期は 2124 です。パラメーター化を通じて 261 個までの並列ストリームをサポートしており、それぞれ長さが 272 のサブストリーム 251 個をサポートしています。U(0,1) の各値は、64 ビット整数 1 つを使用して発生器から作成されます。取り得る値はすべて 253 の倍数で、厳密に区間 (0,1) 内にある値です。mlfg6331_64 ストリームを既定で使用する randn アルゴリズムは ziggurat アルゴリズム [7] ですが、mlfg6331_64 発生器も含みます。

mrg32k3a

[2]で説明されている 32 ビットの複数の再帰発生器の組み合わせです。この発生器は RngStreams パッケージに実装されている CMRG と同じです。周期は 2191 で、シーケンスの分割により、それぞれの長さが 2127 である 263 個までの並列ストリームをサポートしています。また、それぞれの長さが 276 のサブストリームを 251 個サポートしています。U(0,1) の各値は、32 ビット整数 2 つを使って発生器から作成されます。取り得る値はすべて 253 の倍数で、厳密に区間 (0,1) 内にある値です。mrg32k3a ストリームを既定で使用する randn アルゴリズムは ziggurat アルゴリズム [7] ですが、mrg32k3a 発生器も含みます。

mt19937ar

[11] で説明されているメルセンヌ・ツイスターで、周期は 2199371 です。U(0,1) の各値は 2 つの 32 ビット整数を使用して作成されます。取り得る値は、区間 (0,1) 内の 253 の倍数です。この発生器は複数のストリームやサブストリームをサポートしていません。mt19937ar ストリームを既定で使用する randn アルゴリズムは ziggurat アルゴリズム [7] ですが、mt19937ar 発生器も含みます。メモ: この発生器は MATLAB Version 7 の初期に、関数 rand が使用していた発生器と同じであることに注意してください (rand('twister',s) で実行されていました)。

dsfmt19937

[12] で説明されている倍精度 SIMD 指向高速メルセンヌ・ツイスターで、メルセンヌ・ツイスター アルゴリズムの高速実装です。周期は 2199371 で、取り得る値は区間 (0,1) 内の 252 の倍数です。この発生器は、ネイティブで [1,2) 内の倍精度値を生成し、これを変換して U(0,1) 値を作成します。この発生器は複数のストリームやサブストリームをサポートしていません。

shr3cong

線形合同法発生器を組み合わせた Marsaglia's SHR3 シフトレジスタ発生器で、乗数は a=69069、加数は b=1234567、モジュラスは 232 です。SHR3 は u=u(I+L13)(I+R17)(I+L5) で定義された 3 シフトレジスタ発生器で、I は恒等作用素、L は左シフト演算子、R は右シフト演算子です。組み合わせた発生器 (SHR3 については[7] で説明) は近似周期が 264 になります。この発生器は複数のストリームやサブストリームをサポートしていません。U(0,1) の各値は、1 つの 32 ビット整数を使って発生器から作成されます。取り得る値はすべて 232 の倍数で、厳密に区間 (0,1) 内にある値です。shr3cong ストリームを既定で使用する randn アルゴリズムは ziggurat アルゴリズム [9] の初期の形式ですが、shr3cong 発生器も含みます。この発生器は MATLAB Version 5の初期に、関数 randn が使用していた発生器と同じであることに注意してください (randn('state',s) で実行されていました)。

    メモ:   [6](1999) で使用されている SHR3 発生器は、[7](2000) で使用されている SHR3 発生器とは異なります。MATLAB では、[7]に示されている最新バージョンの発生器を使用します。

swb2712

[8] で説明されている修正桁下がり付き減算発生器です。この発生器はラグが 27 と 12 である加算遅れフィボナッチ発生器と同様ですが、より長い近似周期 21492 をもつように修正されています。発生器は U(0,1) 値の作成を倍精度で行い、開区間 (0,1) の値が有効になります。swb2712 ストリームを既定で使用する randn アルゴリズムは ziggurat アルゴリズム [7] ですが、swb2712 発生器も含みます。メモ: この発生器は MATLAB Version 5の初期に、関数 rand が使用していた発生器と同じであることに注意してください (rand('state',s) で実行されていました)。

変換アルゴリズム

Inversion

標準正規分布の累積分布逆関数を一様なランダム変量に適用して、通常のランダム変量を計算します。厳密に一様な値が、通常の値ごとに使われます。

Polar

[1] で説明されている polar rejection アルゴリズムです。近似的に 1.27 の一様な値が、平均して通常の値ごとに使われます。

Ziggurat

[7] で説明されている ziggurat アルゴリズムです。近似的に 2.02 の一様な値が、平均して通常の値ごとに使われます。

参照

[1] Devroye, L. Non-Uniform Random Variate Generation, Springer-Verlag, 1986.

[2] L’Ecuyer, P. “Good Parameter Sets for Combined Multiple Recursive Random Number Generators”, Operations Research, 47(1): 159–164. 1999.

[3] L'Ecuyer, P. and S. Côté. “Implementing A Random Number Package with Splitting Facilities”, ACM Transactions on Mathematical Software, 17: 98–111. 1991.

[4] L'Ecuyer, P. and R. Simard. “TestU01: A C Library for Empirical Testing of Random Number Generators,” ACM Transactions on Mathematical Software, 33(4): Article 22. 2007.

[5] L'Ecuyer, P., R. Simard, E. J. Chen, and W. D. Kelton. “An Objected-Oriented Random-Number Package with Many Long Streams and Substreams.” Operations Research, 50(6):1073–1075. 2002.

[6] Marsaglia, G. “Random numbers for C: The END?” Usenet posting to sci.stat.math. 1999. Available online at http://groups.google.com/group/sci.crypt/browse_thread/
thread/ca8682a4658a124d/
.

[7] Marsaglia G., and W. W. Tsang. “The ziggurat method for generating random variables.” Journal of Statistical Software, 5:1–7. 2000. Available online at http://www.jstatsoft.org/v05/i08.

[8] Marsaglia, G., and A. Zaman. “A new class of random number generators.” Annals of Applied Probability 1(3):462–480. 1991.

[9] Marsaglia, G., and W. W. Tsang. “A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions.” SIAM J.Sci.Stat.Comput. 5(2):349–359. 1984.

[10] Mascagni, M., and A. Srinivasan. “Parameterizing Parallel Multiplicative Lagged-Fibonacci Generators.” Parallel Computing, 30: 899–916. 2004.

[11] Matsumoto, M., and T. Nishimura.“Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator.” ACM Transactions on Modeling and Computer Simulation, 8(1):3–30. 1998.

[12] Matsumoto, M., and M. Saito.“A PRNG Specialized in Double Precision Floating Point Numbers Using an Affine Transition.” Monte Carlo and Quasi-Monte Carlo Methods 2008, 10.1007/978-3-642-04107-5_38. 2009.

[13] Moler, C.B. Numerical Computing with MATLAB. SIAM, 2004. Available online at http://www.mathworks.com/moler

[14] Park, S.K., and K.W. Miller. “Random Number Generators: Good Ones Are Hard to Find.” Communications of the ACM, 31(10):1192–1201. 1998.

この情報は役に立ちましたか?