创建和控制随机数流

通过 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 函数，它们将从您创建的流中获取。如果您调用 rand、randn、randi 和 randperm 而不使用 myStream，它们将从全局流中获取。

可以使用 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`	执行 10 轮的 Philox 4×32 生成器	是	2¹⁹³（2⁶⁴ 个流，长度为 2¹²⁹）
`threefry4x64_20`	执行 20 轮的 Threefry 4×64 生成器	是	2⁵¹⁴（2²⁵⁶ 个流，长度为 2²⁵⁸）
`shr3cong`	移位寄存器生成器与线性同余生成器求和	否	2⁶⁴
`swb2712`	修正的借位减法生成器	否	2¹⁴⁹²

生成器 mcg16807、shr3cong 和 swb2712 可与以前的 MATLAB 版本向后兼容。mt19937ar 和 dsfmt19937 主要是为序列应用程序设计的。其余的生成器明确支持并行随机数生成。

根据应用情况，某些生成器可能会更快或者返回精度更高的值。所有伪随机数生成器均基于确定性算法，并且所有生成器都通过足够具体的随机性统计学测试。检查蒙特卡罗模拟结果的一种方法是使用两种或更多种不同的生成器算法重新运行该仿真，MATLAB 中的生成器选项提供了执行该操作的方式。使用不同的生成器时，尽管结果不可能因多个蒙特卡罗采样误差而异，但资料中存在此类验证在特定的生成器算法中出现缺陷的示例。（有关示例，请参阅 [13]。）

生成器算法

mt19937ar: 梅森旋转算法（如 [11] 中所述）的周期为 $2^{19937} - 1$ ，并且每个 U(0,1) 值都是使用两个 32 位整数生成的。可能的值为 (0, 1) 区间内 $2^{- 53}$ 的倍数。此生成器不支持多个流或子流。默认情况下用于 mt19937ar 流的 randn 算法为 ziggurat 算法 [7]，但基于 mt19937ar 生成器。
注意
此生成器与 rand 函数从 MATLAB 版本 7 开始使用的生成器相同，通过 rand('twister',s) 激活。
dsfmt19937: 面向 SIMD 的双精度快速梅森旋转算法（如 [12] 中所述）是速度更快的梅森旋转算法实现。周期为 $2^{19937} - 1$ ，可能的值为 (0, 1) 区间内 $2^{- 52}$ 的倍数。生成器本身会生成 [1, 2) 内的双精度值，这些值经过变换生成 U(0,1) 值。此生成器不支持多个流或子流。
mcg16807: 32 位乘法同余生成器，如 [14] 中所述，其中乘数为 $a = 7^{5}$ ，模为 $m = 2^{31} - 1$ 。此生成器的周期为 $2^{31} - 2$ ，不支持多个流或子流。每个 U(0,1) 值都是使用单个 32 位整数通过该生成器创建的；可能的值为 ${(2^{31} - 1)}^{- 1}$ 的所有倍数（严格位于区间 (0, 1) 内）。对于 mcg16807 流，randn 使用的默认算法是极坐标算法（如 [1] 中所述）。
注意
此生成器与 rand 和 randn 函数从 MATLAB 版本 4 开始使用的生成器相同，通过 rand('seed',s) 或 randn('seed',s) 进行激活。
mlfg6331_64: 64 位的乘法滞后斐波那契生成器，如 [10] 中所述，其中滞后值为 $l = 63$ 、 $k = 31$ 。此生成器与 SPRNG 库中实现的 MLFG 类似。其周期约为 $2^{124}$ 。通过参数化，它最多支持 $2^{61}$ 个并行流以及 $2^{51}$ 个长度均为 $2^{72}$ 的子流。每个 U(0,1) 值都是使用一个 64 位整数通过该生成器创建的；可能的值为 $2^{- 64}$ 的所有倍数（严格位于区间 (0, 1) 内）。默认情况下用于 mlfg6331_64 流的 randn 算法为 ziggurat 算法 [7]，但基于 mlfg6331_64 生成器。
mrg32k3a: 32 位组合多递归生成器，如[2]中所述。此生成器类似于在 C 语言的 RngStreams 包中实现的 CMRG。其周期为 $2^{191}$ ，通过序列分割支持多达 $2^{63}$ 个并行流，每个流的长度为 $2^{127}$ 。它还支持 $2^{51}$ 个长度均为 $2^{76}$ 的子流。每个 U(0,1) 值都是使用两个 32 位整数通过该生成器创建的；可能的值为 $2^{- 53}$ 的倍数（严格位于区间 (0, 1) 内）。默认情况下用于 mrg32k3a 流的 randn 算法为 ziggurat 算法 [7]，但基于 mrg32k3a 生成器。
philox4x32_10: 执行 10 轮的 4×32 生成器，如 [15] 中所述。此生成器使用 Feistel 网络和整数乘法。该生成器专为高并行系统（如 GPU）的高性能而设计。它的周期为 2¹⁹³（2⁶⁴ 个流，长度为 2¹²⁹）。
threefry4x64_20: 执行 20 轮的 4×64 生成器，如 [15] 中所述。此生成器是 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 位整数通过该生成器创建的；可能的值为 $2^{- 32}$ 的所有倍数（严格位于区间 (0, 1) 内）。默认情况下用于 shr3cong 流的 randn 算法为早期形式的 ziggurat 算法 [9]，但基于 shr3cong 生成器。此生成器与 randn 函数从 MATLAB 版本 5 开始使用的生成器相同，通过 randn('state',s) 进行激活。
注意
[6] (1999) 中用到的 SHR3 生成器与 [7] (2000) 中不同。MATLAB 使用 [7] 中提供的最新版生成器。
swb2712: 修正的借位减法生成器，如[8]中所述。此生成器与加法滞后斐波那契生成器（其滞后值为 27 和 12）类似，但修改后具有约为 $2^{1492}$ 的更长周期。此生成器本身会使用双精度类型来创建 U(0,1) 值，并且开区间 (0, 1) 中的所有值都是可能的。默认情况下用于 swb2712 流的 randn 算法为 ziggurat 算法 [7]，但基于 swb2712 生成器。
注意
此生成器与 rand 函数从 MATLAB 版本 5 开始使用的生成器相同，通过 rand('state',s) 进行激活。

转换算法

Inversion: 通过对均匀的随机变量应用标准正态逆累积分布函数来计算正态随机变量。每个正态值恰好使用一个均匀值。
Polar: 极坐标抑制算法，如[1]中所述。每个正态值平均约使用 1.27 个均匀值。
Ziggurat: ziggurat 算法，如[7]中所述。每个正态值平均约使用 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

您可以分别使用 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