创建和控制随机数流
通过 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 | 梅森旋转 | 否 | 219937-1 |
dsfmt19937 | 面向 SIMD 的快速梅森旋转算法 | 否 | 219937-1 |
mcg16807 | 乘法同余生成器 | 否 | 231-2 |
mlfg6331_64 | 乘法滞后斐波那契生成器 | 是 | 2124(251 个流,长度为 272) |
mrg32k3a | 组合多递归生成器 | 是 | 2191(263 个流,长度为 2127) |
philox4x32_10 | 执行 10 轮的 Philox 4×32 生成器 | 是 | 2193(264 个流,长度为 2129) |
threefry4x64_20 | 执行 20 轮的 Threefry 4×64 生成器 | 是 | 2514(2256 个流,长度为 2258) |
shr3cong | 移位寄存器生成器与线性同余生成器求和 | 否 | 264 |
swb2712 | 修正的借位减法生成器 | 否 | 21492 |
生成器 mcg16807
、shr3cong
和 swb2712
可与以前的 MATLAB 版本向后兼容。mt19937ar
和 dsfmt19937
主要是为序列应用程序设计的。其余的生成器明确支持并行随机数生成。
根据应用情况,某些生成器可能会更快或者返回精度更高的值。所有伪随机数生成器均基于确定性算法,并且所有生成器都通过足够具体的随机性统计学测试。检查蒙特卡罗模拟结果的一种方法是使用两种或更多种不同的生成器算法重新运行该仿真,MATLAB 中的生成器选项提供了执行该操作的方式。使用不同的生成器时,尽管结果不可能因多个蒙特卡罗采样误差而异,但资料中存在此类验证在特定的生成器算法中出现缺陷的示例。(有关示例,请参阅 [13]。)
生成器算法
mt19937ar
梅森旋转算法(如 [11] 中所述)的周期为 ,并且每个 U(0,1) 值都是使用两个 32 位整数生成的。可能的值为 (0, 1) 区间内 的倍数。此生成器不支持多个流或子流。默认情况下用于
mt19937ar
流的randn
算法为 ziggurat 算法 [7],但基于mt19937ar
生成器。注意
此生成器与
rand
函数从 MATLAB 版本 7 开始使用的生成器相同,通过rand('twister',s)
激活。dsfmt19937
面向 SIMD 的双精度快速梅森旋转算法(如 [12] 中所述)是速度更快的梅森旋转算法实现。周期为 ,可能的值为 (0, 1) 区间内 的倍数。生成器本身会生成 [1, 2) 内的双精度值,这些值经过变换生成 U(0,1) 值。此生成器不支持多个流或子流。
-
mcg16807
32 位乘法同余生成器,如 [14] 中所述,其中乘数为 ,模为 。此生成器的周期为 ,不支持多个流或子流。每个 U(0,1) 值都是使用单个 32 位整数通过该生成器创建的;可能的值为 的所有倍数(严格位于区间 (0, 1) 内)。对于
mcg16807
流,randn
使用的默认算法是极坐标算法(如 [1] 中所述)。注意
此生成器与
rand
和randn
函数从 MATLAB 版本 4 开始使用的生成器相同,通过rand('seed',s)
或randn('seed',s)
进行激活。mlfg6331_64
64 位的乘法滞后斐波那契生成器,如 [10] 中所述,其中滞后值为 、。此生成器与 SPRNG 库中实现的 MLFG 类似。其周期约为 。通过参数化,它最多支持 个并行流以及 个长度均为 的子流。每个 U(0,1) 值都是使用一个 64 位整数通过该生成器创建的;可能的值为 的所有倍数(严格位于区间 (0, 1) 内)。默认情况下用于
mlfg6331_64
流的randn
算法为 ziggurat 算法 [7],但基于mlfg6331_64
生成器。mrg32k3a
32 位组合多递归生成器,如[2]中所述。此生成器类似于在 C 语言的 RngStreams 包中实现的 CMRG。其周期为 ,通过序列分割支持多达 个并行流,每个流的长度为 。它还支持 个长度均为 的子流。每个 U(0,1) 值都是使用两个 32 位整数通过该生成器创建的;可能的值为 的倍数(严格位于区间 (0, 1) 内)。默认情况下用于
mrg32k3a
流的randn
算法为 ziggurat 算法 [7],但基于mrg32k3a
生成器。philox4x32_10
执行 10 轮的 4×32 生成器,如 [15] 中所述。此生成器使用 Feistel 网络和整数乘法。该生成器专为高并行系统(如 GPU)的高性能而设计。它的周期为 2193(264 个流,长度为 2129)。
threefry4x64_20
执行 20 轮的 4×64 生成器,如 [15] 中所述。此生成器是 Skein 散列函数中的 Threefish 块加密的非加密改编。它的周期为 2514(2256 个流,长度为 2258)。
shr3cong
Marsaglia 的 SHR3 移位寄存器生成器与线性同余生成器求和,其中乘数为 ,加数为 ,模为 。SHR3 是一个定义为 的 3 移位寄存器生成器,其中 为单位运算符, 为向左移位运算符,R 为向右移位运算符。此组合生成器(SHR3 部分在 [7] 中有述)的周期约为 。此生成器不支持多个流或子流。每个 U(0,1) 值都是使用一个 32 位整数通过该生成器创建的;可能的值为 的所有倍数(严格位于区间 (0, 1) 内)。默认情况下用于
shr3cong
流的randn
算法为早期形式的 ziggurat 算法 [9],但基于shr3cong
生成器。此生成器与randn
函数从 MATLAB 版本 5 开始使用的生成器相同,通过randn('state',s)
进行激活。swb2712
修正的借位减法生成器,如[8]中所述。此生成器与加法滞后斐波那契生成器(其滞后值为 27 和 12)类似,但修改后具有约为 的更长周期。此生成器本身会使用双精度类型来创建 U(0,1) 值,并且开区间 (0, 1) 中的所有值都是可能的。默认情况下用于
swb2712
流的randn
算法为 ziggurat 算法 [7],但基于swb2712
生成器。注意
此生成器与
rand
函数从 MATLAB 版本 5 开始使用的生成器相同,通过rand('state',s)
进行激活。
转换算法
配置流
随机数流 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.