Creating and Controlling a Random Number Stream

The RandStream class allows you to create a random number stream. This is useful for several reasons:

You can generate random values without affecting the state of the global stream.
You can separate sources of randomness in a simulation.
You can use a generator that is configured differently than the one MATLAB^® software uses at startup.

With a RandStream object, you can create your own stream, set the writable properties, and use the stream 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 function.

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

ans =
    0.6986    0.7413    0.4239    0.6914    0.7255

The random stream myStream acts separately from the global stream. If you call the rand, randn, randi, and randperm functions with myStream as the first argument, they draw from the stream you created. If you call rand, randn, randi, and randperm without myStream, they draw from the global stream.

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 can use substreams to get different results that are statistically independent from a stream. Unlike seeds, where the locations along the sequence of random numbers are not exactly known, the spacing between substreams is known, so any chance of overlap can be eliminated. In short, substreams are a more-controlled way to do many of the same things that seeds have traditionally been used for. Substreams are also a more lightweight solution than parallel streams.

Substreams provide a quick and easy way to ensure that you get different results from the same code at different times. To use the Substream property of a RandStream object, create a stream using a generator that supports substreams. For a list of generator algorithms that support substreams and their properties, see the table in the next section. For example, generate several random numbers in a loop.

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

In another loop, you can generate random values that are independent from the first set of 5 iterations.

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

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, you can return to the 6th substream in the loop. The result contains the same values as the 6th output above.

myStream.Substream = 6;
z = rand(1,5)

z =
    0.2650    0.8229    0.2479    0.0247    0.4581

Choosing a Random Number Generator

MATLAB offers several generator algorithm options. The 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	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 backward 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 might be faster or return values with more precision. All pseudorandom number generators are based on deterministic algorithms, and all generators pass 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 the choice of generators in MATLAB provides you with the means to do that. Although it is unlikely that your results will differ by more than the 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). For mcg16807 streams, the default algorithm used by randn 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 library. 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 in C. It has a period of $2^{191}$ and supports up to $2^{63}$ parallel streams through 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. The generator 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 it 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

When creating a random number stream with RandStream, you can specify the NormalTransform name-value argument using one of the following algorithms: 'Ziggurat', 'Polar', or 'Inversion'. This algorithm transforms a uniform generator into a Gaussian one, enabling the generation of normally distributed numbers using randn. See [16] for details on the 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.

Configuring a Stream

A random number stream s has properties that control its behavior. To access or change a property, use the syntax p = s.Property and s.Property = p.

For example, you can configure the transformation algorithm to generate normally distributed pseudorandom values when using randn. Generate normally distributed pseudorandom values using the default Ziggurat transformation algorithm.

s1 = RandStream('mt19937ar');
s1.NormalTransform

ans = 'Ziggurat'

r1 = randn(s1,1,10);

Configure the stream to use the Polar transformation algorithm to generate normally distributed pseudorandom values.

s1.NormalTransform = 'Polar'

s1 =
mt19937ar random stream
             Seed: 0
  NormalTransform: Polar

r2 = randn(s1,1,10);

When generating random numbers with uniform distribution using rand, you can also configure the stream to generate antithetic pseudorandom values, that is, the usual values subtracted from 1 for uniform values.

Create 6 random numbers with uniform distribution from the stream s.

s2 = RandStream('mt19937ar');
r1 = rand(s2,1,6)

r1 =
    0.8147    0.9058    0.1270    0.9134    0.6324    0.0975

Restore the initial state of the stream. Create another 6 random numbers with the Antithetic property set to true. Check that these 6 random numbers are equal to the previously generated random numbers subtracted from 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

Instead of setting the properties of a stream one-by-one, you can save and restore all properties of a stream s by using A = get(s) and set(s,A), respectively. For example, configure all properties of the stream s2 to be the same as the stream 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

The get and set functions enable you to save and restore a stream's entire configuration so that results are exactly reproducible later on.

Restore State of Random Number Generator to Reproduce Output

The State property is the internal state of the random number generator. You can save the state of the global stream at a certain point in your simulation when generating random numbers to reproduce the results later on.

Use RandStream.getGlobalStream to return a handle to the global stream, that is, the current global stream that rand generates random numbers from. Save the state of the global stream.

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

Using myState, you can restore the state of globalStream and reproduce previous results.

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

ans = logical
   1

rand, randi, randn, and randperm access the global stream. Since all of these functions access the same underlying stream, a call to one affects the values produced by the others at subsequent calls.

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

ans = logical
   0

You can also reset a stream to its initial settings using the reset function.

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

ans = logical
   1

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.

[16] Thomas, D. B., P. H. W. Leong, and J. D. Villasenor. “Gaussian Random Number Generators.”ACM Computing Surveys, 39(4): 11-es. 2007.