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

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 を最初の引数として関数 rand、randn、randi、および randperm を呼び出した場合は、作成したストリームから値が引き出されます。myStream を指定せずに rand、randn、randi、および 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`	メルセンヌ・ツイスター	なし	2¹⁹⁹³⁷-1
`dsfmt19937`	SIMD 指向高速メルセンヌ・ツイスター	なし	2¹⁹⁹³⁷-1
`mcg16807`	乗算合同法発生器	なし	2³¹-2
`mlfg6331_64`	乗法ラグフィボナッチ発生器	あり	2¹²⁴ (長さが 2⁷² の 2⁵¹ 個のストリーム)
`mrg32k3a`	結合多重再帰発生器	あり	2¹⁹¹ (長さが 2¹²⁷ の 2⁶³ 個のストリーム)
`philox4x32_10`	Philox 4x32 発生器、10 ラウンド	あり	2¹⁹³ (長さが 2¹²⁹ の 2⁶⁴ 個のストリーム)
`threefry4x64_20`	Threefry 4x64 発生器、20 ラウンド	あり	2⁵¹⁴ (長さが 2²⁵⁸ の 2²⁵⁶ 個のストリーム)
`shr3cong`	線形合同法発生器を組み合わせたシフトレジスタ発生器	なし	2⁶⁴
`swb2712`	修正桁下げ付き減算発生器	なし	2¹⁴⁹²

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

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

発生器アルゴリズム

mt19937ar: [11]で説明されているメルセンヌ・ツイスターで、周期は $2^{19937} - 1$ です。U(0,1) の各値は 2 つの 32 ビット整数を使用して作成されます。取り得る値は、区間 (0, 1) 内にある $2^{- 53}$ の倍数です。この発生器は複数のストリームやサブストリームをサポートしていません。既定で mt19937ar ストリームに使用される randn アルゴリズムは ziggurat アルゴリズム[7]ですが、mt19937ar 発生器も含みます。
メモ
この発生器は MATLAB Version 7 以降に関数 rand で使用されていたものと同一です。rand('twister',s) により起動していました。
dsfmt19937: [12]で説明されているこの倍精度 SIMD 指向高速メルセンヌ・ツイスターは、メルセンヌ・ツイスターアルゴリズムの高速実装です。周期は $2^{19937} - 1$ で、取り得る値は区間 (0, 1) 内にある $2^{- 52}$ の倍数です。この発生器は、ネイティブで [1, 2) 内の倍精度値を生成し、これを変換して U(0,1) 値を作成します。この発生器は複数のストリームやサブストリームをサポートしていません。
mcg16807: [14]で説明されている 32 ビットの乗算合同法発生器で、乗数は $a = 7^{5}$ 、モジュロは $m = 2^{31} - 1$ です。この発生器は $2^{31} - 2$ の周期をもち、複数のストリームやサブストリームをサポートしていません。U(0,1) の各値は、32 ビット整数 1 つを使って発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある ${(2^{31} - 1)}^{- 1}$ の倍数です。mcg16807 ストリームの場合、randn で使用される既定のアルゴリズムは polar アルゴリズムです ([1]で説明)。
メモ
この発生器は MATLAB Version 4 以降に関数 rand と関数 randn で使用されていたものと同一です。rand('seed',s) または randn('seed',s) により起動していました。
mlfg6331_64: [10]で説明されている 64 ビットの乗法ラグフィボナッチ発生器で、ラグは $l = 63$ 、 $k = 31$ です。この発生器は SPRNG ライブラリに実装されている MLFG と類似しています。周期はおよそ $2^{124}$ です。パラメーター化を通じて $2^{61}$ 個までの並列ストリームをサポートしており、それぞれ長さが $2^{72}$ のサブストリーム $2^{51}$ 個をサポートしています。U(0,1) の各値は、64 ビット整数 1 つを使用して発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある $2^{- 64}$ の倍数です。既定で mlfg6331_64 ストリームに使用される randn アルゴリズムは ziggurat アルゴリズム[7]ですが、mlfg6331_64 発生器も含みます。
mrg32k3a: [2]で説明されている 32 ビットの結合多重再帰発生器です。この発生器は C の RngStreams パッケージに実装されている CMRG と同様です。周期は $2^{191}$ で、シーケンスの分割により、それぞれ長さが $2^{127}$ である $2^{63}$ 個までの並列ストリームをサポートしています。また、それぞれ長さが $2^{76}$ のサブストリームを $2^{51}$ 個サポートしています。U(0,1) の各値は、32 ビット整数 2 つを使って発生器から作成されます。取り得る値は、厳密に区間 (0, 1) 内にある $2^{- 53}$ の倍数です。既定で mrg32k3a ストリームに使用される randn アルゴリズムは ziggurat アルゴリズム[7]ですが、mrg32k3a 発生器も含みます。
philox4x32_10: [15]で説明されている 4x32 発生器、10 ラウンド。この発生器は Feistel ネットワークおよび整数乗算を使用しています。この発生器は特に、GPU のような並列性の高いシステムにおいて高度なパフォーマンスを発揮するよう設計されています。これは 2¹⁹³ (長さが 2¹²⁹ の 2⁶⁴ 個のストリーム) の周期をもちます。
threefry4x64_20: [15]で説明されている 4x64 発生器、20 ラウンド。この発生器は、Skein ハッシュ機能の Threefish ブロック暗号の非暗号化適応です。これは 2⁵¹⁴ (長さが 2²⁵⁸ の 2²⁵⁶ 個のストリーム) の周期をもちます。
shr3cong: 線形合同法発生器を組み合わせた Marsaglia の SHR3 シフトレジスタ発生器で、乗数は $a = 69069$ 、加数は $b = 1234567$ 、モジュラスは $2^{- 32}$ です。SHR3 は $u = u (I + L^{13}) (I + R^{17}) (I + L^{5})$ で定義された 3 シフトレジスタ発生器で、 $I$ は恒等作用素、 $L$ は左シフト演算子、R は右シフト演算子です。組み合わせた発生器 (SHR3 については[7]で説明) は周期がおよそ $2^{64}$ になります。この発生器は複数のストリームやサブストリームをサポートしていません。U(0,1) の各値は、32 ビット整数 1 つを使用して発生器から作成されます。取り得る値はすべて、厳密に区間 (0, 1) 内にある $2^{- 32}$ の倍数です。既定で shr3cong ストリームに使用される randn アルゴリズムは ziggurat アルゴリズム[9]の初期の形式ですが、shr3cong 発生器も含みます。この発生器は MATLAB Version 5 以降に関数 randn で使用されていたものと同一です。randn('state',s) により起動していました。
メモ
[6](1999) で使用されている SHR3 発生器は、[7](2000) で使用されている SHR3 発生器とは異なります。MATLAB では、[7]に示されている最新バージョンの発生器を使用します。
swb2712: [8] で説明されている修正桁下がり付き減算発生器です。この発生器はラグが 27 と 12 である加算遅れフィボナッチ発生器と似ていますが、はるかに長いおよそ $2^{1492}$ の周期をもつように修正されています。発生器は 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

rand、randi、randn、および 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.

参考

rng | RandStream