Creating and Controlling a Random Number Stream

The RandStream class allows you to create a random number stream. This is useful for several reasons. For example, you might want to generate random values without affecting the state of the global stream. You might want separate sources of randomness in a simulation. Or you may need to use a different generator algorithm than the one MATLAB^® software uses at startup. With the RandStream constructor, you can create your own stream, set the writable properties, and use it to generate random numbers. You can control the stream you create the same way you control the global stream. You can even replace the global stream with the stream you create.

To create a stream, use the RandStream constructor.

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

ans =

    0.6530    0.8147    0.7167    0.8615    0.0764

The random stream myStream acts separately from the global stream. The functions rand, randn, and randi will continue to draw from the global stream, and will not affect the results of the RandStream methods rand, randn and randi applied to myStream.

You can make myStream the global stream using the RandStream.setGlobalStream method.

RandStream.setGlobalStream(myStream)
RandStream.getGlobalStream

ans = 

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

RandStream.getGlobalStream==myStream

ans =

     1

Substreams

You may want to return to a previous part of a simulation. A random stream can be controlled by having it jump to fixed checkpoints, called substreams. The Substream property allows you to jump back and forth among multiple substreams. To use the Substream property, create a stream using a generator that supports substreams. (See Choosing a Random Number Generator for a list of generator algorithms and their properties.)

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

The initial value of Substream is 1.

stream.Substream

ans =

     1

Substreams are useful in serial computation. Substreams can recreate all or part of a simulation by returning to a particular checkpoint in stream. For example, they can be used in loops.

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

Each of these substreams can reproduce its loop iteration. For example, you can return to the 5th substream. The result will return the same values as the 5th output above.

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

ans =

    0.3323    0.9296    0.5767    0.1233    0.6934

Choosing a Random Number Generator

MATLAB offers several generator algorithm options. The following table summarizes the key properties of the available generator algorithms and the keywords used to create them. To return a list of all the available generator algorithms, use the RandStream.list method.

Keyword	Generator	Multiple Stream and Substream Support	Approximate Period In Full Precision
mt19937ar	Mersenne twister (used by default stream at MATLAB startup)	No	2¹⁹⁹³⁷-1
dsfmt19937	SIMD-oriented fast Mersenne twister	No	2¹⁹⁹³⁷-1
mcg16807	Multiplicative congruential generator	No	2³¹-2
mlfg6331_64	Multiplicative lagged Fibonacci generator	Yes	2¹²⁴ (2⁵¹ streams of length 2⁷²)
mrg32k3a	Combined multiple recursive generator	Yes	2¹⁹¹ (2⁶³ streams of length 2¹²⁷)
philox4x32_10	Philox 4x32 generator with 10 rounds	Yes	2¹⁹³ (2⁶⁴ streams of length 2¹²⁹)
threefry4x64_20	Threefry 4x64 generator with 20 rounds	Yes	2⁵¹⁴ (2²⁵⁶ streams of length 2²⁵⁸)
shr3cong	Shift-register generator summed with linear congruential generator	No	2⁶⁴
swb2712	Modified subtract with borrow generator	No	2¹⁴⁹²

The generators mcg16807, shr3cong, and swb2712 provide for backwards compatibility with earlier versions of MATLAB. mt19937ar and dsfmt19937 are designed primarily for sequential applications. The remaining generators provide explicit support for parallel random number generation.

Depending on the application, some generators may be faster or return values with more precision. All pseudorandom number generators are based on deterministic algorithms, and all will fail a sufficiently specific statistical test for randomness. One way to check the results of a Monte Carlo simulation is to rerun the simulation with two or more different generator algorithms, and MATLAB software's choice of generators provide you with the means to do that. Although it is unlikely that your results will differ by more than Monte Carlo sampling error when using different generators, there are examples in the literature where this kind of validation has turned up flaws in a particular generator algorithm (see [13] for an example).

Generator Algorithms

mt19937ar: The Mersenne Twister, as described in [11], has period $2^{19937} - 1$ and each U(0,1) value is created using two 32-bit integers. The possible values are multiples of $2^{- 53}$ in the interval (0,1). This generator does not support multiple streams or substreams. The randn algorithm used by default for mt19937ar streams is the ziggurat algorithm [7], but with the mt19937ar generator underneath. Note: This generator is identical to the one used by the rand function beginning in MATLAB Version 7, activated by rand('twister',s).
dsfmt19937: The Double precision SIMD-oriented Fast Mersenne Twister, as described in [12], is a faster implementation of the Mersenne Twister algorithm. The period is $2^{19937} - 1$ and the possible values are multiples of $2^{- 52}$ in the interval (0,1). The generator produces double precision values in [1,2) natively, which are transformed to create U(0,1) values. This generator does not support multiple streams or substreams.
mcg16807: A 32-bit multiplicative congruential generator, as described in [14], with multiplier $a = 7^{5}$ , modulo $m = 2^{31} - 1$ . This generator has a period of $2^{31} - 2$ and does not support multiple streams or substreams. Each U(0,1) value is created using a single 32-bit integer from the generator; the possible values are all multiples of ${(2^{31} - 1)}^{- 1}$ strictly within the interval (0,1). The randn algorithm used by default for mcg16807 streams is the polar algorithm (described in [1]). Note: This generator is identical to the one used beginning in MATLAB Version 4 by both the rand and randn functions, activated using rand('seed',s) or randn('seed',s).
mlfg6331_64: A 64-bit multiplicative lagged Fibonacci generator, as described in [10], with lags $l = 63$ , $k = 31$ . This generator is similar to the MLFG implemented in the SPRNG package. It has a period of approximately $2^{124}$ . It supports up to $2^{61}$ parallel streams, via parameterization, and $2^{51}$ substreams each of length $2^{72}$ . Each U(0,1) value is created using one 64-bit integer from the generator; the possible values are all multiples of $2^{- 64}$ strictly within the interval (0,1). The randn algorithm used by default for mlfg6331_64 streams is the ziggurat algorithm [7], but with the mlfg6331_64 generator underneath.
mrg32k3a: A 32-bit combined multiple recursive generator, as described in [2]. This generator is similar to the CMRG implemented in the RngStreams package. It has a period of $2^{191}$ and supports up to $2^{63}$ parallel streams via sequence splitting, each of length $2^{127}$ . It also supports $2^{51}$ substreams, each of length $2^{76}$ . Each U(0,1) value is created using two 32-bit integers from the generator; the possible values are multiples of $2^{- 53}$ strictly within the interval (0,1). The randn algorithm used by default for mrg32k3a streams is the ziggurat algorithm [7], but with the mrg32k3a generator underneath.
philox4x32_10: A 4x32 generator with 10 rounds as described in [15]. This generator uses a Feistel network and integer multiplication, and is specifically designed for high performance in highly parallel systems such as GPUs. It has a period of 2¹⁹³ (2⁶⁴ streams of length 2¹²⁹).
threefry4x64_20: A 4x64 generator with 20 rounds as described in [15]. This generator is a non-cryptographic adaptation of the Threefish block cipher from the Skein Hash Function. It has a period of 2⁵¹⁴ (2²⁵⁶ streams of length 2²⁵⁸).
shr3cong: Marsaglia's SHR3 shift-register generator summed with a linear congruential generator with multiplier $a = 69069$ , addend $b = 1234567$ , and modulus $2^{- 32}$ . SHR3 is a 3-shift-register generator defined as $u = u (I + L^{13}) (I + R^{17}) (I + L^{5})$ , where $I$ is the identity operator, $L$ is the left shift operator, and R is the right shift operator. The combined generator (the SHR3 part is described in [7]) has a period of approximately $2^{64}$ . This generator does not support multiple streams or substreams. Each U(0,1) value is created using one 32-bit integer from the generator; the possible values are all multiples of $2^{- 32}$ strictly within the interval (0,1). The randn algorithm used by default for shr3cong streams is the earlier form of the ziggurat algorithm [9], but with the shr3cong generator underneath. This generator is identical to the one used by the randn function beginning in MATLAB Version 5, activated using randn('state',s).
Note
The SHR3 generator used in [6] (1999) differs from the one used in [7] (2000). MATLAB uses the most recent version of the generator, presented in [7].
swb2712: A modified Subtract-with-Borrow generator, as described in [8]. This generator is similar to an additive lagged Fibonacci generator with lags 27 and 12, but is modified to have a much longer period of approximately $2^{1492}$ . The generator works natively in double precision to create U(0,1) values, and all values in the open interval (0,1) are possible. The randn algorithm used by default for swb2712 streams is the ziggurat algorithm [7], but with the swb2712 generator underneath. Note: This generator is identical to the one used by the rand function beginning in MATLAB Version 5, activated using rand('state',s).

Transformation Algorithms

Inversion: Computes a normal random variate by applying the standard normal inverse cumulative distribution function to a uniform random variate. Exactly one uniform value is consumed per normal value.
Polar: The polar rejection algorithm, as described in [1]. Approximately 1.27 uniform values are consumed per normal value, on average.
Ziggurat: The ziggurat algorithm, as described in [7]. Approximately 2.02 uniform values are consumed per normal value, on average.

References

[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.