第3章〓随机数生成器
3.1本章引言
随机数(RNs)是任何蒙特卡罗模拟的重要组成部分。任何蒙特卡罗模拟的质量都取决于所使用随机数的质量(或随机性)。如果随机数均匀分布，则实现了高度的随机性。因此，需要设计出一种方法，生成随机的数字序列，在重复之前有很长的周期，并且不需要耗费大量资源。
计算机上随机数生成器的早期实现可以追溯到冯·诺依曼在与曼哈顿计划(1941—1945年)相关的蒙特卡罗模拟中使用了随机数生成器。从那时起，许多生成随机数的方法［15,40,5859,6264,6970］被开发出来。LEcuyer指出，一般的计算机用户通常认为计算机生成均匀随机数的问题已经得到解决。尽管在开发“好”生成器方面取得了重大进展，但仍有“坏”或“不合适”的生成器，并产生了不良结果。因此，任何蒙特卡罗用户都应该清楚这些问题及其局限性。
完整起见，本章介绍了实验和算法技术，并讨论了随机性检验的选择，通过实例证明了随机数生成器的随机性和周期高度依赖几个参数的“正确”选择。这些例子表明，参数的微小变化会对估计的随机数产生重大影响。
本章将回顾用于确定随机数的实验和算法，检查随机数生成器的行为，回顾几个随机性检验，并详细说明参数对随机数序列长度和生成随机数随机性的影响。
3.2随机数生成方法
产生随机数的常见方法有两种： ①实验； ②算法。
1. 实验(查表，在线)
实验(物理过程)用于生成一系列随机数，这些随机数以表的形式保存在计算机内存中。例如： ①抛硬币或掷骰子; ②从瓮中抽球，即抽彩票; ③旋转第2章中介绍的标记圆盘; ④在飞镖游戏中从中心测量飞镖的位置。Frigerio和Clark［36］提出了一种方法，记录放射性物质在给定时间间隔内的衰变次数，例如，计数20 ms，如果数字是奇数，则记录0位(bit)，否则记录1位，然后形成31位的数字，这个过程在1 h内产生了不到6000个数字。显然，在蒙特卡罗模拟中使用这种方法的速度是相当慢的。(阿贡国家实验室(ANL)程序中心有一盘包含2.5×106个随机数的磁带。)另一种技术是监控计算机部件(如计算机内存)以生成随机数。注意，这是在执行实际模拟时完成的。
2. 算法(或确定性方法)
算法用于生成随机数。这种方法由于其确定性的本质，通常被称为伪随机数生成器(PRNG)，所生成的数字被称为伪随机数(PRNs)。
每种方法都有其优缺点，这直接影响到对方法的使用。选择“正确”方法有如下6个重要因素。
（1） 随机性： 随机数应该是真正随机的； 也就是说，它们应该均匀分布。在实验方法中，如果过程遵循均匀分布，则实现随机性。在算法方法中，生成的序列必须满足几个统计检验。
（2） 再现性: 为了测试模拟算法或进行敏感性/扰动研究，有必要能够再现随机数序列。
（3） 随机数序列的长度: 为了对实际工程问题进行蒙特卡罗模拟，需要数百万个随机数; 因此，生成器必须能够产生大量的随机数。
（4） 计算机内存: 生成器不应占用过多的计算机内存。(请注意，这个问题可能对于下一代计算机而言并不重要。)
（5） 生成时间: 生成随机数序列所需的时间(工程师/分析师)不是很重要，可以是几天或者几个月。
（6） 计算机时间: 生成数字所需的计算机时间应该明显短于实际模拟的时间。
表3.1基于上述因素对实验随机数生成器和算法随机数生成器进行了比较。
表3.1实验随机数生成器与算法随机数生成器的比较
要素
实验
查表法在线法算法
随机性
好好
待检验
再现性
可以不可以
可以
周期*
有限无限
有限
计算机内存
大小
小
生成时间**
长无
无
计算机时间
短长
短
* 随机数序列的最大长度。
** 生成随机数所需时间。
在实践中，算法方法是首选，主要是因为它的序列是可重复的，并且只需要最小的努力(计算机资源，工程师时间)。在使用任何算法生成器时，都必须进行一系列随机性检验，并确定/测量生成器重复其序列的周期，即随机性的丧失。
下面几节将讨论不同的算法随机(伪随机)数生成器和相关的随机性检验方法。
3.3伪随机数生成器
一个好的伪随机数生成器必须产生一个在(0,1)上均匀分布的随机数序列。此外，它应该有一个很长的周期，并且必须通过一系列的随机性检验。本节将介绍常用的生成器，包括: ①同余生成器; ②多重递归生成器。
3.3.1同余生成器
同余生成器是德里克·亨利·莱默（D.H.Lehmer）提出的一种整数生成器。它采用:
xk+1=(axk+b)，mod M,b<M(3.1)
其中，x0(种子)、a、b、M是给定的整数； M是计算机表示的最大整数，例如，在32位的二进制机器上，最大的无符号整数是232-1，最大的有符号整数是231-1。模函数决定了(α=axk+b)除以M的余数，例如，35对16取模等于3。注意，同余一词源于αM与xk+1M具有相同余数，即xk+1在模 M下与α同余。值得注意的是，式（3.1）被称为线性同余生成器，如果b=0，则称为乘法同余生成器。
这里讨论的方法将给出一个在［0,M-1］上的随机整数x，为了将生成的随机整数转换为［0,1)上的随机数，使用关系η=xM-1。
用一个简单的例子来讨论同余生成器的细节和相关问题是有指导意义的。考虑由式（3.2）给出的线性同余生成器［20］：
xk+1=(5xk+1),mod 16(3.2)
假设种子为1，则随机数序列可以表示为一个随机数循环(顺时针)，如图3.1所示。
图3.1由式（3.2）给出的生成器随机循环示意
注意： 等高线给出了随机数值，径向索引指的是序列中随机数的顺序
图3.1中的循环展示了一个生成器的预期行为。
（1） 这个序列的周期是16，即等效于模量M=24，最大值为24-1。
（2） 选择任何其他数(小于M=16)作为种子将导致上述序列的周期性移位。
检验乘子(a)对上述生成器周期的影响是有指导意义的。表3.2给出了选择不同乘子（3~15）时生成器的周期。
表3.2乘子(a)对伪随机数生成器周期的影响(式（3.2）)
参数
值
乘子（a）
3
4
5
6
7
8
9
10
11
12
13
14
15
周期
8
2
16
4
4
2
16
4
8
2
16
4
2
表3.2表明，如果乘子(a)等于5、9或13，则得到完整周期为16，否则得到部分周期为2、4或8。更具体地说，图3.2给出了乘子9(实线)和乘子13(虚线)的随机数序列。
图3.2乘子9(实线)和乘子13(虚线)的随机数序列
图3.3给出了导致得到部分周期的乘子8（实线）和乘子12（虚线）的随机数序列。
图3.3乘子8(实线)和乘子12(虚线)的随机数序列
正如预期的那样，在具有完整周期的情况下，随机数形成具有不同值的形状，而所有部分周期的情况都呈现重复的形状。接下来，检验常数(b)的影响是有指导意义的，对于乘子为5的情况，本节分析了常数(b)对生成器周期的影响。表3.3比较了不同b值（2~8）时的生成器的周期。
表3.3表明，奇数b能得到一个完整的周期，而偶数b只能得到部分周期。
表3.3常数(b)对PRNG周期的影响(见式（3.2）)
参数
值
常数（b）
2
3
4
5
6
7
8
周期
8
16
2
16
8
16
4
考虑到上述观察结果，对于任意n，生成器xk+1=(5xk+1)mod 2n应该实现一个完整的周期。图3.4和图3.5分别给出了n=5和n=6时的随机数序列。
图3.4M=25时的式（3.2）PRNG的随机数序列
图3.5M=26时的式（3.2）PRNG的随机数序列
由图3.4和图3.5可知，正如期望的那样，M=32和M=64分别对应生成32个和64个随机数，即全周期。
为了实现线性同余生成器的全周期，可以使用Hull和Dobell［53］定理的推论，该定理设置了不同参数的性质，包括乘子、常数和模。表3.4给出了这些属性。
表3.4实现全周期的线性同余生成器参数的性质
参数
值
注释
乘子(a)
4N+1
N>0
常数(b)
奇数
模(M)
2k
k>1
现在，设置常数参数(b)为0来考虑一个乘法同余生成器，即
xk+1=(axk)mod 16(3.3)
假设种子为1，表3.5给出了不同乘子（2~15）的不同随机数序列的周期。
表3.5乘子(a)对PRNG周期的影响(见式（3.3）)
参数
值
乘子（a）
2
3
4
5
6
7
8
9
10
11
12
13
14
15
周期
—
4
—
4
—
2
—
2
—
4
—
4
—
2
图3.6显示了乘法生成器(见式（3.3）)对乘子3和乘子11生成的随机数序列。这演示了生成器的部分周期，因为每种情况都会导致多边形通过重复种子而终止。
图3.6乘子3(实线)和乘子11(虚线)的随机数序列，式（3.3）PRNG
同理可以得到乘子分别为5和13的结果。这些结果表明，若N≥0，则当a=8N+3或a=8N+5时，乘法同余生成器的周期为2k-2，其中k对应于模的指数，如16=24。现在，考虑一个具有奇数乘数(16339)和偶数模2k的生成器，使用表3.6中给出的算法确定在任意幂N下的周期。
表3.6线性同余随机数生成器的一个算法
算法
注释
readt,*,multi,const,imod读取种子、乘子、常量和模数
ixx(1)= x0初始化随机数数组
do i=2,imod
iab=multi * ixx(i-1)+const
ixx(i)=mod(iab,imod)
if (ixx(i).lt.0)then
ixx(i)=ixr(i)+imod
end if
循环生成随机数，
计算axi-1+b
用“mod”函数得到新的随机数
检查随机数是否为正
if (ixx(i).eq.xo)then
iperiod=i-1
go to 10
else
iperiod = i
end if
end do
10 continue
确定周期
注意，mod函数可以用式（3.4）代替：
ixx(i)=iab-INT(iab/mod)imod(3.4)
注意，该算法的实现依赖特定系统的整数运算知识。(有关这方面的讨论，见附录A。)
如果考虑2k的模(mod)的不同幂k，可以证明，每个序列的周期是模的25%。
为了实现乘法同余生成器的大周期(大于模量的25%)，已经证明［20］，如果考虑以下几点，可以实现M-1的周期:
（1） M是质数。
（2） 乘子a是M的基本因子。
Park和Miller［79］证明了模数为231-1，乘子为16807时，周期为M-1。关于最佳随机数生成器更进一步的信息，读者可以参考以下文献： Bratley、Fox和Schrage［15］、Park和Miller［79］、Zeeb和Burns［121］、LEcuyer［63］、Gentle［39］及Marsaglia［69］。
总之，具有合理周期的同余生成器对于大多数物理模拟可以产生令人满意的结果，因为物理系统通过将相同的随机数应用于不同的物理过程来引入随机性。然而，对于解决数学问题，这种说法并不总是正确的。
3.3.2多重递归生成器
一组伪随机数生成器称为多重递归生成器［39］，可以表示为
xk+1=(a0xk+a1xk-1+…+ajxk-j+b),mod M(3.5)
其通过选择j+1个随机数(可能来自更简单的生成器)来初始化生成器。生成器的长度和随机性(或统计属性)取决于aj、b和M的值。这类生成器的一个特例是斐波那契生成器。
斐波那契生成器是一个浮点数生成器。它的特点是通过前两个数字的组合(差、和或积)计算出一个新数字。如式（3.6）所示［55］:
xk=xk-17-xk-5(3.6)
是一个17和5滞后的斐波那契生成器，所生成的序列取决于初始xk的选择，如17。由于上面的斐波那契公式xk+1依赖前面的17个值，所以它的周期(p)可以很大，即p=(217-1)2n，其中n是xi的小数部分的位数。例如,对于32位的浮点运算，n=24,因此，p约为241或1012。由于预期周期大，因此对于一些大型问题，斐波那契生成器是一个很好的选择。超级计算机上更大的问题使用了滞后97和33的斐波那契生成器。
要启动斐波那契生成器，必须生成初始随机数，如17个数字。一种方法是以二进制形式表示每个初始随机数:
r=r12+r222+…+rm2m,对于其中m≤n(尾数)(3.7)
其中，ri是二进制数，即0或1。需要一个更简单的同余生成器来设置每个位ri为0或1。例如，我们可以根据整数同余生成器的输出是大于0还是小于0，将ri设置为0或1。因此，斐波那契公式的质量确实取决于初始数字的质量，或者更简单的整数生成器的质量。表3.7给出了推荐的斐波那契生成器列表。
表3.7推荐的斐波那契生成器列表
生成器
期 望 周 期
xk=xk-17-xk-5
(217-1)×224=2.2×1012
xk=xk-31-xk-13
(231-1)×224=3.6×1016
xk=xk-97-xk-33
(297-1)×224=2.7×1036
请注意，当使用32位整数时，同余生成器可以拥有的最大周期是232，或4.3×109，这比表3.7中的任何一个建议都要小得多至少是其1500。
3.4随机性检验
如果伪随机数生成器通过了一系列随机性检验，那么它就是一个可以被接受的生成器。许多随机性检验［59,68］通过考察各种参数来检验随机数序列的独立性和一致性。本节的其余部分将讨论一组检验，这些检验根据随机数的位数和整数来检查随机数［77］。
3.4.1χ2检验
χ2检验测量样本与假定概率分布(假设)之间的偏差。χ2检验的公式由式（3.8）给出：
χ2=∑ni=1(Ni-Npi)2Npi(3.8)
其中，{p1,p2，…，pn}是与N个事件相关的假设概率集合，这些事件属于N个类别，观测频率为N1,N2，…，Nn。请注意，与假设分布相比，该检验一次检查整个抽样分布，并且在这个意义上比检查样本均值、样本方差等更普遍。对于较大的N值，随机变量χ2近似服从自由度为n-1的χ2分布密度函数。
3.4.1.1χ2分布
如果随机变量ω=χ2具有χ2分布［19］，且遵循如下概率密度函数：
fm(w)dw=wm2-12m2Γm2e-w2dw(3.9)
其中，m为正整数，表示自由度； Γ为伽马函数，ω>0。一般来说，我们感兴趣的是找到ω=χ2小于给定值χ20的概率，但可用的χ2分布表通常给出χ2≥χ20的概率。因此，有必要使用概率的补集，即
P(χ2≤χ20)=1-P(χ2≥χ20)(3.10)
其中，
P（χ2≥χ20）=∫∞χ20dwfm（w）（3.11）
值得注意的是，随着m的增大，χ2分布趋于正态分布，分布的均值和方差分别等于m和2m。
3.4.1.2χ2检验的流程
根据式（3.8）得到卡方(χ2)值，然后将这些值与χ2分布表中给出的确定值进行比较，如表3.8所示。
表3.8一个χ2分布表的例子(P(χ2≥χ20))
自由度
0.990
0.950
0.050
0.010
0.001
1
0.000
0.004
3.840
6.640
10.830
2
0.020
0.103
5.990
9.210
13.820
3
0.115
0.352
7.820
11.350
16.270
4
0.297
0.711
9.490
13.280
18.470
5
0.554
1.145
11.070
15.090
20.520
6
0.872
1.635
12.590
16.810
22.460
7
1.239
2.167
14.070
18.480
24.320
