均匀分布随机数的生成算法简介

关键词：均匀分布随机数，线性同余法，反馈位寄存器法，组合随机数发生器

一、引言

\qquad设随机变量 X X 的分布函数为F(X)F(X)， {Xi,i=1,2,⋯} \{X_i,i=1,2,\cdots\} 独立同分布 F(X) F(X) ，则 {Xi,i=1,2,⋯} \{X_i,i=1,2,\cdots\} 的一次观察值 {x1,x2,x3,⋯} \{x_1,x_2,x_3,\cdots\} 称为分布 F(X) F(X) 随机数序列，简称随机数。随机数是随机模拟的基本构成元素，其质量的优劣将直接影响模拟研究的成败。人们本能的可以通过物理实验产生一些常见分布的随机数，如可以通过反复的投硬币来产生二项分布的随机数，可以通过反复投骰子来产生多项分布的随机数，通过对排队的观察记录来产生泊松分布的随机数，该方法产生的随机数质量好，但是数量有限，而且成本很高。随后人们尝试预先生成大量的真实随机数存储起来，进行随机模拟时读取存储的随机数，但是这种方法的速度较低，如今已经被取代。现在随机模拟采用的随机数通常是由计算机按确定的递推公式实时地产生的伪随机数，其在一定程度上体现随机性。好的伪随机数序列与真实随机数序列在统计检验上表现几乎相同，很难区分，因此伪随机数通常也被称为随机数。计算机最容易产生的随机数是均匀分布随机数，产生这些随机数的发生器主要有线性同余发生器、反馈位移寄存器发生器以及组合发生器。

二、线性同余发生器

\qquad线性同余发生器是由Lehmer在1951年提出的，其是利用数论中的同余运算来产生随机数，其递推公式为：

xn=(axn−1+c)(modM);a>0,c≥0,M>0n=1,2,⋯

x_n = (ax_{n-1}+c)\pmod M; \quad a \gt 0,\quad c \ge 0, \quad M \gt 0 \quad n = 1,2,\cdots 其中 c>0 c\gt0 时，称为混合同余发生器， c=0 c=0 时，称为乘同余发生器。

1. 混合同余发生器

\qquad 为了对混合同余发生器有更加直观的认识，我们列举一些线性同余发生器的例子，如:

xn=(7xn−1+7)(mod10)

x_n = (7x_{n-1}+7)\pmod {10} 初始值取 x0=7 x_0=7 ，生成的序列为 {7,6,9,0,7,6,9,0,7,6,⋯} \{7,6,9,0,7,6,9,0,7,6,\cdots\} ，周期 T=4<M=10 T=4 \lt M=10；

xn=(5xn−1+1)(mod10)

x_n = (5x_{n-1}+1)\pmod {10} 初始值取 x0=1 x_0=1 ，生成的序列为 {1,6,1,6,1,6,⋯} \{1,6,1,6,1,6,\cdots\} ，周期 T=2<M=10 T=2 \lt M=10；

xn=(5xn−1+1)(mod8)

x_n = (5x_{n-1}+1)\pmod {8} 初始值取 x0=1 x_0=1 ，生成的序列为 {1,6,7,4,5,2,3,0,1,6,7,⋯} \{1,6,7,4,5,2,3,0,1,6, 7,\cdots\} ，周期 T=8=M=8 T=8 = M=8。
\qquad 从上述的例子中,我们发现，不同的 a,c,M a,\quad c,\quad M取值，序列周期也不相同。若 T<M {T \lt M} 时，则序列的取值种类不能遍历 0∼M−1 0 \sim M-1，也即在 0∼M−1 0 \sim M-1 中取值不是均匀等可能的，因此需要选取使得 T=M T=M （也即满周期）的取值。那么什么样的取值会使产生的这样的序列呢？满足以下三个条件的参数取值，能够使得序列达到满周期: (1) (1) c c 与 MM 互素， (2) (2) 对于M的任意素因子 P P，a−1a-1能被 P P整除，(3)(3) 如果 4 4 是 MM 的因子，则 a−1 a-1 能被 4 4 整除。具体的证明可参考该文献。由于在实际的应用过程的中，使用得随机数的数量很大，因此周期（MM）也需要尽可能的大，同时为了利用计算机整数溢出yuanl原理简化计算， M M 的取值通常为 2L2^L（ L L 为计算机整数的尾字长数）。由条件（1）（1）可知，在 M=2L M = 2^L 时， c c 只能为通常的取值为奇数 c=2β+1c = 2\beta + 1，由条件 (2) (2) 和 (3) (3) 可知， 4α+1 4\alpha + 1 是 a a 合理的取值。同时，由序列中前后两项自相关系数的近似公式：

ρ(1)≈1a−6caM(1−cM)

\rho(1) \approx \frac{1}{a}-\frac{6c}{aM}(1-\frac{c}{M})可知， a a 应该尽可能大，但是应小于 MM。

2. 乘同余发生器

\qquad c=0 c=0 的线性同余发生器为乘同余发生器，其递推公式转变成：

xn=axn−1(modM);a>0,M>0n=1,2,⋯

x_n = ax_{n-1}\pmod M; \quad a \gt 0,\quad M \gt 0 \quad n = 1,2,\cdots从公式可知，由于 xi≠0 x_i \neq 0 ，所以其产生的序列不能达到满周期（ M M），其可能的最大周期为 M−1M-1。那么 a a 和 MM 的取什么样的值能够使得序列的周期足够的接近甚至达到 M−1 M-1 呢？由于当 (M,a)≠1 (M,a)\neq 1 或 (M,x0)≠1 (M,x_0)\neq 1 时，发生器产生的序列会退化到 xi=0 x_i = 0 的状态，也即周期为 1 1，研究的意义不大，在以下的讨论中不予考虑。通过同余的传递性，很容易就能证明序列的周期为满足aV≡1(modM)a^V\equiv 1 \pmod M的最小整数 V V，也即 aa 对 M M 的阶数。因此取合适的值使得 V=M−1V = M-1 ，序列就能达到可能的最大周期。
\qquad Hutchinson提出了一种取值方式： M M 为小于2L2^L 的最大素数，取 a a 为 MM 的素元。这种取值方式形成的发生器也称为素数模乘同余发生器，该类型的发生器是目前使用最广泛的随机数发生器。素数模乘同余发生器面临两个问题：（1）由于 M M 不是 2L2^L 的形式，不能利用计算机的溢出原理来减少除法运算，（2）如何求得 M M 的素元 aa。1969年Payne，Rabang和Bogyo提出的“模拟式除法”基本解决了问题（1），问题（2）涉及很多数论的知识，其仍有很多值得探究的价值。除了素数模乘同余发生器，前人的研究已经给出了其他两种种类型的取值，但是周期都与 M−1 M-1 相差较大：(1) 当 M M 为2L（L≥4）2^L（L\ge 4）， x0 x_0 为奇数， a≡3或5(mod8) a \equiv 3 或 5 \pmod 8 时，序列的周期为 2L−2 2^L-2 ；(2)当 M M 为10s（s≥5）10^s（s \ge 5）， x0 x_0 不是2或5的倍数时， a a 取合适的值能够使得序列的周期达到 5∗10s−2 5 * 10^s-2 。

三、反馈位移寄存器发生器

\qquad 线性同余发生器虽然应用普遍，优点众多，但是其有两大缺陷：（1）随机数的序列周期很难大于 2L 2^L，（2）生成多维均匀分布随机数时，相关性大。因此，1965年Tausworthe提出了一种通过对寄存器位移直接在存储单元产生随机数的方法，使用该方法的发生器被称为反馈位移寄存器发生器，其递推公式为：

αk=(cpαk−p+cp−1αk−p+1+⋯+c1αk−1)(mod2)

\alpha _k=(c_p \alpha _{k-p} + c_{p-1} \alpha _{k-p+1} + \cdots + c_1 \alpha _{k-1})\pmod 2其中 p p 是给定的正整数，ci∈{0,1}c_i \in \{0,1\} 。反馈位移寄存器发生器产生的随机数序列在最佳的情况下，周期能够达到 2p−1 2^{p-1} ，不受计算机字长的限制，且自相关系数近似为零；当 mL≤p mL\le p，且 mL mL 与 2p−1 2^p-1 互素时，可构成维的均匀随机数序列。

四、组合发生器

\qquad 组合发生器是采用二个或多个随机数发生器以某种方式进行组合而形成的发生器，其产生的随机数有更好的随机性，更长的周期，统计性质更优。最著名的组合发生器是由MacLaren和Marsaglia设计的组合线性同余发生器，其原理是通过第二个线性同余发生器来扰乱第一个线性同余发生器产生的随机数，从而得到最终的结果。Gebhardt也证明了这种组合发生器的随机性增强，周期增大的性质。

五、总结与不足

\qquad 本文简要的梳理了产生随机数的常用的算法，对其原理给出了简单的说明，没有给出过于详细的证明，有兴趣的可以留言私信我。下一篇文章我们将介绍非均匀分布随机数是如何生存。