Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

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

RandStream クラスを使用すると、より柔軟に乱数ストリームを生成できます。このことはいくつかの理由で便利です。

  • グローバル ストリームの状態に影響を与えずに乱数値を生成できる。

  • シミュレーションにおいて乱数のソースを分けることができる。

  • MATLAB® ソフトウェアが起動時に使用するものとは異なる構成の発生器を使用できる。

RandStream オブジェクトを使用すると、独自のストリームを作成して、書き込み可能なプロパティを設定し、そのストリームを使用して乱数を生成できます。作成したストリームは、グローバル ストリームと同じ方法で制御することができます。また、グローバル ストリームを作成したストリームと置き換えることもできます。

ストリームを作成するには、関数 RandStream を使用します。

myStream = RandStream('mlfg6331_64');
rand(myStream,1,5)
ans =
    0.6986    0.7413    0.4239    0.6914    0.7255

乱数ストリーム myStream はグローバル ストリームとは別に動作します。myStream を最初の引数として関数 randrandnrandi、および randperm を呼び出した場合は、作成したストリームから値が引き出されます。myStream を指定せずに randrandnrandi、および randperm を呼び出した場合は、グローバル ストリームから値が引き出されます。

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

サブストリーム

サブストリームを使用して、ストリームから統計的に独立しているさまざまな結果を取得できます。乱数列での位置が正確には把握されないシードとは異なり、サブストリーム間の間隔は把握されているため、オーバーラップの可能性をなくすことができます。つまり、従来はシードを使用して行われていた多くのことを、サブストリームを使えば、より制御された方法で実行できます。また、サブストリームは並列ストリームよりも軽量なソリューションです。

サブストリームは、すばやく簡単に、異なるタイミングで同じコードから異なる結果を確実に取得する方法となります。RandStream オブジェクトの Substream プロパティを使用するには、サブストリームをサポートする発生器を使ってストリームを作成します (サブストリームをサポートする発生器のアルゴリズムとそのプロパティのリストは、乱数発生器の選択を参照してください)。たとえば、ループ内でいくつかの乱数を生成します。

myStream = RandStream('mlfg6331_64');
RandStream.setGlobalStream(myStream)
for i = 1:5
    myStream.Substream = i;
    z = rand(1,i)
end
z =
    0.6986

z =
    0.9230    0.2489

z =
    0.0261    0.2530    0.0737

z =
    0.3220    0.7405    0.1983    0.1052

z =
    0.2067    0.2417    0.9777    0.5970    0.4187

別のループ内で、反復 5 回の最初のセットとは独立した乱数を生成できます。

for i = 6:10
    myStream.Substream = i;
    z = rand(1,11-i)
end
z =
    0.2650    0.8229    0.2479    0.0247    0.4581

z =
    0.3963    0.7445    0.7734    0.9113

z =
    0.2758    0.3662    0.7979

z =
    0.6814    0.5150

z =
    0.5247

サブストリームは逐次計算を行う場合に便利です。サブストリームは指定したチェックポイントに戻って、シミュレーションのすべてまたは一部分を作成しなおすことができます。たとえば、ループ内の 6 番目のサブストリームに戻ることができます。結果には、上記の 6 番目の出力と同じ値が含まれています。

myStream.Substream = 6;
z = rand(1,5)
z =
    0.2650    0.8229    0.2479    0.0247    0.4581

乱数発生器の選択

MATLAB は、発生器アルゴリズムのさまざまなオプションを提供します。次の表は、使用可能な発生器アルゴリズムの主要な性質と、発生器の作成に使うキーワードをまとめています。使用可能な発生器アルゴリズムのすべてを返すには、RandStream.list メソッドを使用します。

キーワード発生器複数のストリームとサブストリームのサポート高精度な場合の近似周期
mt19937arメルセンヌ・ツイスター (MATLAB スタートアップ時に既定のストリームによって使用)なし219937-1
dsfmt19937SIMD 指向高速メルセンヌ・ツイスター なし219937-1
mcg16807乗算合同法発生器なし231-2
mlfg6331_64乗法ラグ フィボナッチ発生器あり2124 (長さが 272 の 251 個のストリーム)
mrg32k3a結合多重再帰発生器あり2191 (長さが 2127 の 263 個のストリーム)
philox4x32_10Philox 4x32 発生器、10 ラウンドあり2193 (長さが 2129 の 264 個のストリーム)
threefry4x64_20Threefry 4x64 発生器、20 ラウンドあり2514 (長さが 2258 の 2256 個のストリーム)
shr3cong線形合同法発生器を組み合わせたシフトレジスタ発生器なし264
swb2712修正桁下げ付き減算発生器なし21492

発生器 mcg16807shr3cong および swb2712 は、MATLAB の前のバージョンと下位互換性があります。mt19937ar および dsfmt19937 は主にシーケンシャルな用途のために設計されています。残りの発生器は明示的に並列乱数発生をサポートしています。

用途によっては、一部の発生器はより高速に動作するか、あるいはより高精度の値を返す場合があります。すべての疑似乱数発生器は決定性アルゴリズムに基づいており、また、すべての発生器はランダム性に関する十分に特定的な統計検定に合格しています。モンテ カルロ シミュレーションの結果をチェックする 1 つの方法は、2 つ以上の異なった発生器アルゴリズムを使ってシミュレーションを実行することです。MATLAB では発生器の選択ができるため、この方法が適用できます。異なった発生器を使った場合の結果の相違が、モンテ カルロ法のサンプリング誤差より大きくなることはまずありませんが、この種の検証によって特定の発生器のアルゴリズムに欠陥が見つかった例が文献にあります (例については、[13]を参照してください)。

発生器アルゴリズム

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) 値を作成します。この発生器は複数のストリームやサブストリームをサポートしていません。

mcg16807

[14]で説明されている 32 ビットの乗算合同法発生器で、乗数は a=75、モジュロは m=2311 です。この発生器は 2312 の周期をもち、複数のストリームやサブストリームをサポートしていません。U(0,1) の各値は、32 ビット整数 1 つを使って発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある (2311)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 つを使用して発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある 264 の倍数です。既定で mlfg6331_64 ストリームに使用される randn アルゴリズムは ziggurat アルゴリズム[7]ですが、mlfg6331_64 発生器も含みます。

mrg32k3a

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

philox4x32_10

[15]で説明されている 4x32 発生器、10 ラウンド。この発生器は Feistel ネットワークおよび整数乗算を使用しています。この発生器は特に、GPU のような並列性の高いシステムにおいて高度なパフォーマンスを発揮するよう設計されています。これは 2193 (長さが 2129 の 264 個のストリーム) の周期をもちます。

threefry4x64_20

[15]で説明されている 4x64 発生器、20 ラウンド。この発生器は、Skein ハッシュ機能の Threefish ブロック暗号の非暗号化適応です。これは 2514 (長さが 2258 の 2256 個のストリーム) の周期をもちます。

shr3cong

線形合同法発生器を組み合わせた Marsaglia の 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) の各値は、32 ビット整数 1 つを使用して発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある 232 の倍数です。既定で 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 の一様な値が、平均して通常の値ごとに使われます。

ストリームの設定

乱数ストリーム s は、その動作を制御するプロパティをもっています。プロパティへのアクセスやその変更を行うには、構文 p = s.Property および s.Property = p を使用します。

たとえば、randn を使用する際に、正規分布の疑似乱数値を生成するための変換アルゴリズムを設定できます。既定の Ziggurat 変換アルゴリズムを使用して、正規分布の疑似乱数値を生成します。

s1 = RandStream('mt19937ar');
s1.NormalTransform
ans = 'Ziggurat'
r1 = randn(s1,1,10);

Polar 変換アルゴリズムを使用するようストリームを設定して、正規分布の疑似乱数値を生成します。

s1.NormalTransform = 'Polar'
s1 =
mt19937ar random stream
             Seed: 0
  NormalTransform: Polar
r2 = randn(s1,1,10);

rand を使用して一様分布の乱数を生成する場合は、対称な疑似乱数値 (一様な値の場合、1 から通常値を減算したもの) を生成するようストリームを設定することもできます。

ストリーム s から一様分布の 6 つの乱数を生成します。

s2 = RandStream('mt19937ar');
r1 = rand(s2,1,6)
r1 =
    0.8147    0.9058    0.1270    0.9134    0.6324    0.0975

ストリームの初期状態を復元します。Antithetic プロパティを true に設定し、もう一度 6 つの乱数を作成します。これら 6 つの乱数が、前回生成した乱数を 1 から減算したものと等しいことを確認します。

reset(s2)
s2.Antithetic = true;
r2 = rand(s2,1,6)
r2 =
    0.1853    0.0942    0.8730    0.0866    0.3676    0.9025
isequal(r1,1 - r2)
ans =
   1

ストリームのプロパティを 1 つずつ設定する代わりに、A = get(s)set(s,A) をそれぞれ使用することで、ストリーム s のすべてのプロパティを保存および復元できます。たとえば、ストリーム s2 のすべてのプロパティを、ストリーム s1 と同じになるように設定します。

A = get(s1)
A =
               Type: 'mt19937ar'
         NumStreams: 1
        StreamIndex: 1
          Substream: 1
               Seed: 0
              State: [625x1 uint32]
    NormalTransform: 'Polar'
         Antithetic: 0
      FullPrecision: 1
set(s2,A)
               Type: 'mt19937ar'
         NumStreams: 1
        StreamIndex: 1
          Substream: 1
               Seed: 0
              State: [625x1 uint32]
    NormalTransform: 'Polar'
         Antithetic: 0
      FullPrecision: 1

関数 get と関数 set によってストリームの設定全体を保存および復元できるため、後で結果を厳密に再現することができます。

出力を再現するための、乱数発生器の状態の復元

State プロパティには、乱数発生器の内部状態が設定されます。乱数を生成した後で結果を再現するために、シミュレーションの特定のポイントにおけるグローバル ストリームの状態を保存できます。

RandStream.getGlobalStream を使用してグローバル ストリームのハンドル、すなわち、rand が乱数の生成元とする現在のグローバル ストリームを返します。グローバル ストリームの状態を保存します。

globalStream = RandStream.getGlobalStream;
myState = globalStream.State;

myState を使うと、globalStream の状態を保存し、実行結果を再現することができます。

A = rand(1,100);
globalStream.State = myState;
B = rand(1,100);
isequal(A,B)
ans = logical
   1

randrandirandn、および randperm はグローバル ストリームにアクセスします。これらすべての関数は同じストリームにアクセスするため、呼び出した関数の生成値は前に呼び出した関数の影響を受けます。

globalStream.State = myState;
A = rand(1,100);
globalStream.State = myState;
C = randi(100);
B = rand(1,100);
isequal(A,B)
ans = logical
   0

関数 reset を使用して、ストリームを初期設定にリセットすることもできます。

reset(globalStream)
A = rand(1,100);
reset(globalStream)
B = rand(1,100);
isequal(A,B)
ans = logical
   1

参照

[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 https://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 https://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 https://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.

[15] Salmon, J. K., M. A. Moraes, R. O. Dror, and D. E. Shaw. "Parallel Random Numbers: As Easy As 1, 2, 3." In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC11). New York, NY: ACM, 2011.

参考

|

関連するトピック