5 均匀分布随机数生成

设随机变量\(X\)有分布函数\(F(x)\)\(\{ X_i, i=1,2,\dots \}\)独立同分布\(F(x)\), 则\(\{ X_i, i=1,2,\dots \}\)的一次观测的数值 \(\{ x_i, i=1,2,\dots \}\)叫做分布\(F(x)\)的随机数序列, 简称随机数。 随机数是统计中一个重要的计算工具“随机模拟”的基本构成元素。

我们可以用物理方法得到一组真实的随机数, 比如,反复抛掷硬币、骰子、正二十面体骰子, 抽签、摇号,等等。 这些方法得到的随机数质量好, 但是数量不能满足随机模拟的需要。

另一种办法是预先生成大量的真实随机数存储起来, 进行随机模拟时读取存储的随机数。 这种方法的速度较低, 已经被取代了。

现在进行随机模拟的主流方法是使用计算机实时地产生随机数, 严格来说是伪随机数。 伪随机数是用计算机算法产生的序列 \(\{ x_i, i=1,2,\dots \}\), 其表现与真实的\(F(x)\)的独立同分布序列很难区分开来。 因为计算机算法的结果是固定的, 所以伪随机数不是真正的随机数; 但是,好的伪随机数序列与真实随机数序列表现相同, 很难区分,现在的计算机模拟都是使用伪随机数, 我们把伪随机数也叫做随机数。

需要某种分布的随机数时,一般先生成均匀分布随机数, 然后由均匀分布随机数再转换得到其它分布的随机数。 产生均匀分布随机数的算法叫做均匀分布随机数发生器

计算机中伪随机数序列是迭代生成的, 即\(x_n = g(x_{n-1}, x_{n-2}, \dots, x_{n-p}), n=1,2,\dots\), \(p\)是正整数,\(g\)是一个非线性函数。 如何找到这样的\(g\),使得产生的序列呈现出随机性? 首先要使结果有随机性。 比如,把一个大的整数平方后取中间的位数, 然后递推进行,叫做平方取中法。 这种方法均匀性差而且很快变成零, 所以已经不再使用。 还可以把序列的前两项相乘然后取中间的数字, 这种方法也有类似缺点。

现在常用的均匀分布随机数发生器有线性同余法、反馈位寄存器法 以及随机数发生器的组合。 均匀分布随机数发生器生成的是\(\{0, 1, \dots, M\}\)\(\{1, 2, \dots, M\}\) 上离散取值的离散均匀分布, 然后除以\(M\)\(M+1\)变成\([0,1]\)内的值当作 连续型均匀分布随机数,实际上是只取了有限个值。 因为取值个数有限, 根据算法\(x_n = g(x_{n-1}, x_{n-2}, \dots, x_{n-p})\)可知序列 一定在某个时间发生重复, 使得序列发生重复的间隔\(T\)叫做随机数发生器的周期。 好的随机数发生器可以保证\(M\)很大而且周期很长。

现代的统计计算编程语言如R、Python、Julia、MATLAB等都已经包含了性能很好的随机数发生器, 这里只是提供关于随机数发生器的一般知识, 一般的读者应该没有必要开发自己的随机数发生器, 这一部分用作教材时可以仅作介绍。

5.1 线性同余发生器(LCG)

5.1.1 同余

定义5.1 (同余) \(i,j\)为整数,\(M\)为正整数,若\(j-i\)\(M\)的倍数, 则称\(i\)\(j\)关于模\(M\)同余, 记为\(i \equiv j \pmod{M}\)。 否则称\(i\)\(j\)关于\(M\)不同余。

例5.1 \[\begin{aligned} 11 \equiv& 1 \pmod{10},& 1 \equiv& 11 \pmod{10},\\ -9 \equiv& 1 \pmod{10}. \end{aligned} \]

同余有如下性质:

  1. 对称性: \(i \equiv j \pmod{M} \Longleftrightarrow j \equiv i \pmod{M}\)
  2. 传递性:若\(i \equiv j \pmod{M}\), \(j \equiv k \pmod{M}\), 则\(i \equiv k \pmod{M}\)
  3. \(i_1 \equiv j_1 \pmod{M}\), \(i_2 \equiv j_2 \pmod{M}\), 则 \[\begin{aligned} i_1 \pm i_2 \equiv& j_1 \pm j_2 \pmod{M},& i_1 i_2 \equiv& j_1 j_2 \pmod{M}. \end{aligned}\]
  4. \(ik \equiv jk \pmod{M}\)(\(k\)为正整数), 则\(i \equiv j \pmod{\frac{M}{\mbox{gcd}(M,k)}}\), 其中\(\mbox{gcd}(M,k)\)表示\(M\)\(k\)的最大公约数。

证明: (1)、(2)和(3)的加减运算用同余定义很容易证明。 来证明(3)中乘法的结果。 由\(i_1 \equiv j_1 \pmod{M}\)\(i_2 \equiv j_2 \pmod{M}\) 的定义可知存在整数\(k_1\)\(k_2\)使得 \(j_1 - i_1 = k_1 M\), \(j_2 - i_2 = k_2 M\), 于是 \[ \begin{aligned} j_1 j_2 - i_1 i_2 =& j_1 j_2 - j_1 i_2 + j_1 i_2 - i_1 i_2 \\ =& j_1 ( j_2 - i_2) + i_2 (j_1 - i_1) = j_1 \cdot k_2 M + i_2 \cdot k_1 M \\ =& (j_1 k_2 + i_2 k_1) M \end{aligned} \] 即有\(i_1 i_2 \equiv j_1 j_2 \pmod{M}\)

再来证明(4)。 设\(ik - jk = nM\), \(n\)为整数, 则\(i - j = \frac{nM}{k}\)为整数。 设\(M = M_1 \cdot \mbox{gcd}(M, k)\), \(k = k_1 \cdot \mbox{gcd}(M, k)\), 则\(i - j = \frac{n M_1}{k_1}\)\(M_1\)\(k_1\)互素, 于是\(n\)\(k_1\)整除, 所以\(i \equiv j \pmod{M_1}\), 即\(i \equiv j \pmod{\frac{M}{\mbox{gcd}(M,k)}}\)

\(M\)是正整数,\(A\)为非负整数, \(A\)除以\(M\)的余数为 \(A \pmod{M} = A - [A/M]\times M\), 其中\([\cdot]\)表示取整。 显然\(0 \leq (A \pmod{M}) \leq M-1\)\(A\)\(A \pmod{M}\)关于模\(M\)同余。 在R中用x %% y表示\(x\)除以\(y\)的余数。

5.1.2 线性同余发生器

线性同余随机数发生器是利用求余运算的随机数发生器。 其递推公式为 \[ \begin{aligned} x_n = (a x_{n-1} + c) \pmod{M}, \ n=1,2,\dots \end{aligned} \] 这里等式右边的\((a x_{n-1} + c) \pmod{M}\)表示\(a x_{n-1} + c\)除以\(M\)的余数, 正整数\(M\)为除数,正整数\(a\)为乘数, 非负整数\(c\)为增量, 取某个非负整数\(x_0\)为初值可以向前递推, 递推只需要序列中前一项, 得到的序列\(\{ x_n \}\)为非负整数, \(0 \leq x_n < M\)。 最后, 令\(R_n = x_n / M\), 则\(R_n \in [0,1)\), 把\(\{ R_n \}\)作为均匀分布随机数序列。 这样的算法的基本思想是因为很大的整数前面的位数是重要的有效位数 而后面若干位有一定随机性。 如果取\(c=0\), 线性同余发生器称为乘同余发生器, 如果取\(c>0\), 线性同余发生器称为混合同余发生器

因为线性同余法的递推算法仅依赖于前一项, 序列元素取值只有\(M\)个可能取值, 所以产生的序列\(x_0, x_1, x_2, \ldots\)一定会重复。 若存在正整数\(n\)\(m\)使得 \(x_n = x_m\)(\(m<n\)), 则必有\(x_{n+k} = x_{m+k}, k=0,1,2,\dots\), 即\(x_{n}, x_{n+1}, x_{n+2}, \dots\) 重复了\(x_{m}, x_{m+1}, x_{m+2}, \dots\), 称这样的\(n-m\)的最小值\(T\)为此随机数发生器在初值\(x_0\)下的周期。 由序列取值的有限性可见\(T \leq M\)

例5.2 考虑线性同余发生器 \[ \begin{aligned} x_n = 7 x_{n-1} + 7 \pmod{10}. \end{aligned} \] 取初值\(x_0=7\), 数列为: \((7, 6, 9, 0, 7, 6, 9, 0, 7, \ldots)\), 周期为\(T=4 < M=10\)

例5.3 考虑线性同余发生器 \[ \begin{aligned} x_n = 5 x_{n-1} + 1 \pmod{10} \end{aligned} \] 取初值\(x_0=1\)。 数列为: \((1, 6, 1, 6, 1, \ldots)\), 显然周期\(T=2\)

例5.4 考虑线性同余发生器 \[ \begin{aligned} x_n = 5 x_{n-1} + 1 \pmod{8} \end{aligned} \] 取初值\(x_0=1\)。 数列为\((1, 6, 7, 4, 5, 2, 3, 0, 1, 6, 7, \ldots)\), \(T=8 = M\)达最大周期。

从例子发现\(T \leq M\)且有的线性同余发生器可以达到最大周期\(M\), 称为满周期。 满周期时,初值\(x_0\)取为\(0 \sim M-1\)间的任意数都是一样的, \(X_M=x_0\),序列从\(X_M\)开始重复。 如果发生器从某个初值不是满周期的, 那么它从任何初值出发都不是满周期的, 不同初值有可能得到不同序列。 比如随机数\(x_n = 7 x_{n-1} + 7 \pmod{10}\), 从不同初值出发的序列可能为 如\(7,6,9,0,7,\cdots\); \(1,4,5,2,1, \cdots\); \(3,8,3\)

不同的\(M, a, c\)选取方法得到不同的周期, 适当选取\(M, a, c\)才能使得产生的随机数序列 和真正的U\([0,1]\)随机数表现接近。

5.1.3 混合同余发生器(*)

线性同余发生器中\(c>0\)时称为混合同余发生器。 下列定理给出了混合同余发生器达到满周期的一个充分条件。

定理5.1 当下列三个条件都满足时,混合同余发生器可以达到满周期:

  1. \(c\)\(M\)互素;
  2. \(M\)的任一个素因子\(P\)\(a-1\)\(P\)整除;
  3. 如果\(4\)\(M\)的因子,则\(a-1\)\(4\)整除。

常取\(M = 2^L\)\(L\)为计算机中整数的尾数字长。 这时根据定理5.1的建议 可取\(a = 4 \alpha + 1\)\(c = 2 \beta + 1\)\(\alpha\)\(\beta\)为任意正整数, \(x_0\)为任意非负整数, 这样的混合同余发生器是满周期的, 周期为\(2^L\)

好的均匀分布随机数发生器应该周期足够长, 统计性质符合均匀分布, 序列之间独立性好。 把同余法生成的数列看成随机变量序列\(\{X_n\}\), 在满周期时,可认为\(X_n\)是从\(0 \sim M-1\)中随机等可能选取的,即 \[ \begin{aligned} P \{X_n = i \} = 1/M,\ i=0,1,\ldots,M-1 \end{aligned} \] 这时 \[ \begin{aligned} E X_n =& \sum_{i=0}^{M-1} i \cdot \frac{1}{M} = (M-1)/2, \\ \text{Var}(X_n) =& E X_n^2 - (E X_n)^2 = \sum_{i=0}^{M-1} i^2 \frac{1}{M} - \frac{(M-1)^2}{4} = \frac{1}{12}(M^2 - 1) \end{aligned} \] 于是当\(M\)很大时 \[ \begin{aligned} E R_n =& \frac{1}{2} - \frac{1}{2M} \approx \frac{1}{2}; \\ \text{Var}(R_n) =& \frac{1}{12} - \frac{1}{12M^2} \approx \frac{1}{12} \end{aligned} \] 可见\(M\)充分大时从一、二阶矩看生成的数列很接近于均匀分布。

随机数序列还需要有很好的随机性。 数列的排列不应该有规律, 序列中的两项不应该有相关性。

因为序列由确定性公式生成,所以不可能真正独立。 至少我们要求是序列自相关性弱。 对于满周期的混合同余发生器可以得序列中前后两项自相关系数的近似公式 \[ \begin{aligned} \rho(1) \approx \frac{1}{a} - \frac{6c}{aM}(1 - \frac{c}{M}) \end{aligned} \] 所以应该选\(a\)值大(但\(a<M\))。

例5.5 Kobayashi提出了如下的满周期\(2^{31}\)的混合同余发生器 \[\begin{align*} x_n = (314159269 x_{n-1} + 453806245) \pmod{2^{31}} \end{align*}\] 其周期较长,统计性质比较好。

Julia实现:

function rng_KS_fixed_seed(n=1, seed=0)
  y = Vector{Float64}(n)
  x::Int64 = seed
  for i in 1:n
    x = (314159269*x + 453806245) % 2147483648
    y[i] = x / 2147483648
  end
  y
end
## 测试:
y = rng_KS_fixed_seed(100)

上一版本需要用户输入种子,如果需要多次调用,就需要同时输出当前种子值。 利用“闭包”(closure)可以记住函数的当前状态,下一版本的实现可以自动记住当前种子:

function rng_KS(seed=0)
  local the_seed::Int64 = seed

  function rng(n=1)
    y = Vector{Float64}(n)
    x::Int64 = the_seed
    for i in 1:n
      x = (314159269*x + 453806245) % 2147483648
      y[i] = x / 2147483648
    end

    the_seed = x
    y
  end

  rng
end

测试:

myrng = rng_KS()
myrng(2)
## 2-element Array{Float64,1}:
##  0.21132  
##  0.0102247
myrng(2)
## 2-element Array{Float64,1}:
##  0.188331
##  0.665933

可以看出每次调用产生不同的序列。

5.1.4 乘同余法(*)

线性同余发生器中\(c=0\)时的生成方法称为乘同余法,或称积式发生器。 这时的递推公式为 \[ \begin{aligned} x_n = a x_{n-1} \pmod{M}, \ n=1,2,\dots \end{aligned} \] 问题是如何选\(M, a\)达到大周期且统计性质好。 显然各\(x_n\)都不能等于0; 如果某个\(x_k=0\), 则 \(x_{k+1} = x_{k+2} = \dots = 0\)。 乘同余法有可能进入0就不再周期变化,称为退化情况,这时周期为1。 所以乘同余法能够达到的最大周期是\(M-1\), 每个\(x_n\)都只在\(\{1,2,\dots, M-1 \}\)中取值。

定义5.2 (阶数) 设正整数\(a\)与正整数\(M\)互素, 称满足\(a^V \equiv 1 \pmod{M}\)的最小正整数\(V\)\(a\)对模\(M\) 的阶数(或次数),简称为\(a\)的阶数。

我们来证明阶数存在。 考虑乘同余发生器\(x_n = a x_{n-1}\), \(n=1,2,\dots\), 取\(x_0=1\)。 由同余的传递性可知 \(x_n \equiv a^n x_0 \pmod{M}\)\(x_n \equiv a^n \pmod{M}\), 而\(0 \leq x_n \leq M-1\)所以\(x_n = a^n \pmod{M}\)。 因为\(a\)\(M\)互素所以\(x_n \neq 0\), 序列\(\{ x_n \}\)只在\(1,\dots,M-1\)中取值, 必存在\(0 \leq i < j < i + M\)使得 \(x_j = x_i\), 其中\(x_i \equiv a^i \pmod{M}\)\(x_j \equiv a^j \pmod{M}\), 于是\(a^i \equiv a^j \pmod{M}\), 由同余的性质(4)有 \(1 \equiv a^{j-i} \pmod{\frac{M}{\mbox{gcd}(a^i,M)}}\), 其中\(\mbox{gcd}(a^i,M)=1\)所以 \(a^{j-i} \equiv 1 \pmod{M}\)。 取\(V = j-i\)\(V\)为正整数且\(a^V \equiv 1 \pmod{M}\)

Lemma 5.1 \(a\)\(M\)互素,初值\(x_0\)\(M\)互素, 则乘同余法的周期为\(a\)对模\(M\)的阶数\(V\)

证明: 由同余的传递性可知\(x_V \equiv a^V x_0 \pmod{M}\), 而\(a^V \equiv 1 \pmod{M}\)所以 \[ \begin{aligned} a^V x_0 \equiv x_0 \pmod{M} \end{aligned} \] 于是 \[ \begin{aligned} x_V \equiv x_0 \pmod{M}. \end{aligned} \] 这时必有\(x_V = x_0\), 所以乘同余法的周期\(T \leq V\)

如果\(T < V\), 则存在\(0 \leq i < j < i+V\)使得\(x_j = x_i\), 因\(x_j \equiv a^j x_0 \pmod{M}\), \(x_i \equiv a^i x_0 \pmod{M}\), 由同余的传递性可知\(a^j x_0 \equiv a^i x_0 \pmod{M}\), 由同余的性质(4)有 \[ \begin{aligned} a^j \equiv a^i \pmod{\frac{M}{\mbox{gcd}(M,x_0)}} \end{aligned} \]\(\mbox{gcd}(M,x_0)=1\)所以 \(a^j \equiv a^i \pmod{M}\), 由同余的性质(4)可知 \[ \begin{aligned} a^{j-i} \equiv 1 \pmod{\frac{M}{\mbox{gcd}(a^i, M)}} \end{aligned} \] 其中\(\mbox{gcd}(a^i, M)=1\)所以 \(a^{j-i} \equiv 1 \pmod{M}\)\(0<j-i<V\), 与\(V\)是满足\(a^V \equiv 1 \pmod{M}\)的最小正整数矛盾。证毕。

※※※※※

对乘同余法, 当\(M\)\(a\)互素且\(M\)\(x_0\)互素时,\(x_n = a^n x_0 \pmod{M}\), 易见\(x_n\)\(M\)互素,序列不会退化为0.

例5.6 考虑如下乘同余发生器 \[\begin{align*} x_n =& (8\alpha + 5) x_{n-1} \pmod{2^L}\\ x_0 =& 4b+1 \end{align*}\] (其中\(\alpha\), \(b\)为非负整数)。 再考虑如下的混合同余发生器 \[\begin{align*} x_n^* =& (8\alpha+5) x_{n-1}^* + (2\alpha+1) \pmod{2^{L-2}} \\ x_0^* =& b \end{align*}\]\(x_n = 4 x_n^* + 1, n=0,1,2,\dots\)

用归纳法。当\(n=0\)\[ x_0^*=b, \ x_0 = 4b+1 \] 所以\(x_0 = 4 x_0^* + 1\)成立。 设\(x_n = 4 x_n^* + 1\)成立,由 \[ x_n = (8\alpha+5) x_n \pmod{2^L}, \] 其中 \[\begin{aligned} (8\alpha+5) x_n =& (8\alpha+5) (4 x_n^* + 1) \\ =& 4[ (8\alpha+5) x_n^* + (2\alpha+1)] + 1 \end{aligned} \] 于是 \[ \begin{aligned} (8\alpha+5) x_n \pmod{M} =& 4 \{ [ (8\alpha+5) x_n^* + (2\alpha+1)] \pmod{2^{L-2}} \} + 1 \\ =& 4 x_{n+1}^* + 1 \end{aligned} \] 得证。

这个例子中,由定理5.1可知 \(\{ x_n^* \}\)是满周期的, \(\{ x_n \}\)\(\{x_n^* \}\)一一对应, 所以\(\{ x_n \}\)周期为\(2^{L-2}\)。 这里的乘同余发生器与一个混合同余发生器等价。

※※※※※

对乘同余法得到的随机序列\(\{ X_n \}\), 当\(M\)充分大时可以类似得到 \[ \begin{aligned} E R_n \approx& \frac{1}{2}, & \text{Var}(R_n) \approx& \frac{1}{12}. \end{aligned} \] 为了使得前后两项的相关系数较小, 应取\(a\)较大且使得\(a\)的二进制表示排列无规律。

5.1.5 素数模乘同余法(*)

若取\(M\)为小于\(2^L\)的最大素数, 选取适当的\(a\)可以达到\(T = M-1\)周期, 这样的发生器叫素数模乘同余发生器

定义5.3 (素元) 设正整数\(a\)和素数\(M\)互素, 若\(a\)对模\(M\)的阶数\(V\)满足\(V = M-1\), 则称\(a\)\(M\)的素元(或原根)。

乘同余发生器中取\(a\)\(M\)的素元可以达到最大周期\(M-1\)

例5.7 \(M=7\)是素数,取\(a=3\)\(M\)互素。 考虑素数模乘同余发生器 \[\begin{aligned} x_n = a^n \pmod{M}, \ n=1,2,\dots, \ x_0=1, \end{aligned} \] 序列为\(1, 3, 2, 6, 4, 5, 1, 3, \ldots\)。 周期达到\(M-1=6\), 所以\(a=3\)\(M=7\)的素元。

可以验证\(a=5\)也是\(M=7\)的素元, 序列为\(1, 5, 4, 6, 2, 3, 1, 5, \ldots\)\(1, 2, 4, 6\)不是\(M=7\)的素元。 此例说明素元可以不唯一。

在选取素数模乘同余发生器的参数\(M\)\(a\)的时候, 可以取\(M\)为小于\(2^L\)的最大素数(\(L\)为计算机整数表示的尾数位数), 取\(a\)\(M\)的素元,这样保证周期\(T = M-1\)\(a\)应尽可能大,\(a\)的二进制表示尽可能无规律。 这样在一个周期内可以得到\(1,2,\ldots,M-1\)的一个排列。 初值\(x_0\)可以取\(1,2,\ldots,M-1\)中任一个。

例5.8 Lewis-Goodman-Miller(1969)的素数模乘同余发生器为 \[ \begin{aligned} x_n = a x_{n-1} \pmod{2^{31} - 1}, \ n=1,2,\dots, \ \text{$x_0$ is arbitrary positive integer} \end{aligned} \] \(a\)取如下四个值之一: \((7^5=16807, 397204094, 764261123, 630360016)\)

其中\(a=16807\)的版本有许多应用, 但是其高维表现很差, 现在应该避免使用。

Julia实现,实现一个用户指定种子的版本:

function rng_LGM_fixed_seed(n=1, seed=11111, form=4)
  mults = [16807, 397204094, 764261123, 630360016]
  a = mults[form]
  y = Vector{Float64}(n)
  x::Int64 = seed
  for i in 1:n
    x = (a*x) % 2147483647
    y[i] = x / 2147483647
  end
  y
end
## 测试:
rng_LGM_fixed_seed(5, 112233)
## 5-element Array{Float64,1}:
##  0.230227
##  0.556961
##  0.770818
##  0.692613
## 0.289898

注意,在设计线性同余法程序时, 如果使用程序语言中的整数型或无符号整数型来存储 \(x_n\),则在计算\(a x_{n-1}\)的乘法时可能发生溢出。 对\(M = 2^L\)情形, 如果程序语言支持溢出而不作为错误, 则溢出恰好完成了求除以\(2^L\)的余数的运算。 对于\(M < 2^L\)的情形也可以巧妙设计以利用溢出。 对于较小的\(a\)\(M\)(如\(a\)\(M\)都小于\(2^{L/2}\)时)不会溢出。 如果使用双精度实数型来保存\(x_n\)则不需要考虑溢出问题。

例5.9 设计素数模乘同余发生器算法, 要求不产生溢出。

\(a\)为正整数,\(M\)为素数, 令\(b = [M/a]\), \(c = M \pmod{a}\), 则\(M = a b + c\), 设\(a<b\)。 乘同余递推中的除法满足 \[\begin{align*} \frac{a x_{n-1}}{M} =& \frac{a x_{n-1}}{a b + c} \leq \frac{x_{n-1}}{b} \end{align*}\]\[\begin{align*} &\frac{x_{n-1}}{b} - \frac{a x_{n-1}}{a b + c} < \frac{a x_{n-1}}{a b} - \frac{a x_{n-1}}{a b + a} = \frac{x_{n-1}}{b(b+1)} \\ <& \frac{M}{b(b+1)} = \frac{ab + c}{b(b+1)} < \frac{ab + a}{b(b+1)} = \frac{a}{b} < 1 \end{align*}\]\[\begin{align*} \frac{a x_{n-1}}{M} \leq \frac{x_{n-1}}{b} < \frac{a x_{n-1}}{M} + 1 \end{align*}\]\(k_0 = \left[ \frac{a x_{n-1}}{M} \right]\), \(k_1 = \left[ \frac{x_{n-1}}{b} \right]\), 则\(k_0 \leq k_1 \leq k_0+1\), 计算\(k_0\)不用计算乘法所以不会溢出, \(k_1=k_0\)\(k_0+1\)。 于是 \[\begin{align*} x_n =& a x_{n-1} - k_0 M = a x_{n-1} - k_1 M + (k_1 - k_0) M \\ =& a x_{n-1} - k_1 (ab+c) + (k_1 - k_0) M \\ =& a(x_{n-1} - k_1 b) - k_1 c + (k_1 - k_0) M \\ =& a( x_{n-1} \pmod{b}) - k_1 c + (k_1 - k_0) M \end{align*}\] 其中\(k_1 - k_0\)等于0或1, 因为\(x_n > 0\), 所以如果\(a( x_{n-1} \pmod{b}) - k_1 c < 0\)\(k_1 - k_0 = 1\)

避免溢出的素数模发生器算法
设置\(M, a, b, c, x_0\)的值
for(\(n\) in \(1:N\)){
\(\qquad\) \(k_1 \leftarrow \mbox{floor}(x_{n-1} / b)\)
\(\qquad\) \(x_n \leftarrow a (x_{n-1} - k_1 b) - k_1 c\)
\(\qquad\) if \((x_n < 0)\ x_n \leftarrow x_n + M\)
}

对例5.8的素数模发生器, 取\(a=16807\)可以按上述算法计算。

上述的算法使用了伪代码来描述, 伪代码的控制结构仿照R语言的控制结构。

一些常见的素数模乘同余发生器参数:

  • \(m=2^{31}-1\)\(a=16807\)(Lewis, Goodman, Miller), \(39373\)(L’Ecuyer), \(742938285\), \(950706376\), \(1226874159\)(Fishman and Moore);
  • \(m=2147483399\), \(a=40692\)(L’Ecuyer);
  • \(m=2147483563\), \(a=40014\)(L’Ecuyer)。

5.2 FSR发生器(*)

线性同余法的周期不可能超过\(2^L\)(\(L\)为整数型尾数的位数),
而且作为多维随机数相关性大,分布不均匀。 基于Tausworthe(1965)文章的FSR方法是一种全新的做法, 对这些方面有改善。

FSR(反馈位移寄存器法)按照某种递推法则生成一系列二进制数 \(\alpha_1, \alpha_2, \ldots, \alpha_k, \ldots\), 其中\(\alpha_k\)由前面的若干个\(\{\alpha_i\}\)线性组合并求除以2的余数产生: \[\begin{align} \alpha_k = (c_p \alpha_{k-p} + c_{p-1} \alpha_{k-p+1} + \dots + c_1 \alpha_{k-1}) \pmod{2}, \ k=1,2,\dots \tag{5.1} \end{align}\] 线性组合系数\(\{c_i\}\)只取0, 1, 这样的递推可以利用程序语言中的整数二进制运算快速实现。 给定初值\((\alpha_{-p+1}, \alpha_{-p+2}, \dots, \alpha_{0})\)向前递推, 得到\(\{\alpha_k, k=1,2,\dots\}\) 序列后依次截取长度为\(L\)的二进制位组合成整数\(x_n\)\(R_n=x_n / 2^L\)。 不同的组合系数和初值选择可以得到不同的随机数发生器, 巧妙设计可以得到很长的周期, 作为多维均匀随机数性质较好。

FSR算法中系数\((c_1, c_2, \dots, c_p)\)如果仅有两个为1, 比如\(c_p=c_{p-q}=1\)(\(1<q<p\)), 则算法变成 \[ \begin{aligned} \alpha_k =& (\alpha_{k-p} + \alpha_{k-p+q}) \pmod{2} \\ =& \begin{cases} 0 & \text{若} \alpha_{k-p} = \alpha_{k-p+q} \\ 1 & \text{若} \alpha_{k-p} \neq \alpha_{k-p+q} \end{cases} \end{aligned} \]\(\oplus\)表示二进制异或运算, 则 \[ \begin{aligned} \alpha_k = \alpha_{k-p} \oplus \alpha_{k-p+q}, \ k=1,2,\dots. \end{aligned} \] 比如,取\(p=98, q=27\)

设计FSR计算机程序时, 直接对包含\(M\)个二进制位的非负整数\(\{x_n\}\)的数列用异或递推更方便。 递推公式为 \[ \begin{aligned} x_n =& x_{n-p} \oplus x_{n-p+q}, \quad (1<q<p), \ n=1,2,\dots \\ R_n =& x_n / 2^M \end{aligned} \] 这需要由\(p\)\(M\)位二进制非负整数作为种子(初值)。 这种算法只需要异或运算,不受计算机字长限制, 适当选取\(p,q\)后周期可以达到\(2^p-1\), 作为多维随机数的性质可以很好, 需要预先研究得到种子表而不能随便取初值。

5.3 组合发生器法(*)

随机数设计中比较困难的是独立性和多维的分布。 可以考虑把两个或若干个发生器组合利用, 可以比单个发生器由更长的周期和更好的随机性。

例5.10 Wichmann和Hill(1982)(见 Wichmann and Hill (1982))设计了如下的线性组合发生器。

利用三个16位运算的素数模乘同余发生器: \[\begin{eqnarray*} U_n &=& 171 U_{n-1} \pmod{30269} \\ V_n &=& 172 V_{n-1} \pmod{30307} \\ W_n &=& 170 W_{n-1} \pmod{30323} \end{eqnarray*}\] 作线性组合并求余: \[\begin{align*} R_n = (U_n / 30269 + V_n / 30307 + W_n / 30323) \pmod{1} \end{align*}\] 这个组合发生器的周期约有\(7\times 10^{12}\)长, 超过\(2^{32} \approx 4\times 10^9\)

Julia实现(需要用户输入种子):

function rng_WH_fixed_seed(n, seed1=11, seed2=13, seed3=17)
  y = Vector{Float64}(n)
  u::Int64 = seed1
  v::Int64 = seed2
  w::Int64 = seed3
  for i in 1:n
    u = (171*u) % 30269
    v = (172*v) % 30307
    w = (170*w) % 30323
    y[i] = (u/30269 + v/30307 + w/30323) % 1.0
  end
  y
end
## 测试:
rng_WH_fixed_seed(5)
## 5-element Array{Float64,1}:
##  0.231228
##  0.518513
##  0.153344
##  0.502289
##  0.875749

下面是重复调用时自动记忆上一次结束时种子的随机数发生器制造机:

function make_rng_WH(seed=[1011,2013,3017])
  local the_seed = seed

  function rng(n)
    y = Vector{Float64}(n)
    u::Int64 = the_seed[1]
    v::Int64 = the_seed[2]
    w::Int64 = the_seed[3]
    for i in 1:n
      u = (171*u) % 30269
      v = (172*v) % 30307
      w = (170*w) % 30323
      y[i] = (u/30269 + v/30307 + w/30323) % 1.0
    end

    the_seed = [u, v, w]
    y
  end

  rng
end

例5.11 MacLaren和Marsaglia(1965)设计了组合同余法, 组合两个同余发生器,一个用来“搅乱”次序。

设有两个同余发生器A和B。 用A产生\(m\)个随机数(如\(m=128\)), 存放在数组\(T=(t_1, t_2, \ldots, t_m)\). 需要产生\(x_n\)时, 从B中生成一个随机下标\(j \in \{1,2,\ldots,m\}\), 取\(x_n = t_j\), 但从A再生成一个新随机数\(y\)代替\(T\)中的\(t_j\), 如此重复。

这样组合可以增强随机性,加大周期(可超过\(2^L\))。 也可以只使用一个发生器,用\(x_{n-1}\)来选择下标。

在R软件中, 用runif(n)产生\(n\)个U(0,1)均匀分布的随机数。 R软件提供了若干种随机数发生器, 可以用RNGkind函数切换。 在使用随机数进行随机模拟研究时, 如果希望模拟研究的结果可重复, 就需要在模拟开始时设置固定的随机数种子。 虽然不同的随机数发生器种子的形式有所不同, 在R中, 总是可以用函数set.seed(seed)来设置种子, 其中seed是一个序号,实际的种子由这个序号决定。

5.4 随机数的检验(*)

文献中已经有许多随机数生成算法, 统计软件中也已经包含了许多公认的可靠的随机数发生器。 但是,我们要解决一个自己的模拟问题时, 还是要反复确认所用的随机数在自己的问题中是好的, 没有明显缺陷的。 比如,一个经典的称为RANDU的随机数发生器用在一维积分时效果很好, 但连续三个一组作为三维均匀分布则很不均匀。

为了验证随机数的效果, 最好是找一个和自己的问题很接近但是有已知答案的问题, 用随机模拟计算并考察其精度。

对均匀分布随机数发生器产生的序列\(\{ R_n, n=1,2,\dots,N \}\) 可以进行各种各样的检验以确认其均匀性与独立性。 下面列举一些检验的想法:

  • 对随机数列\(\{ R_n, n=1,2,\dots,N \}\)计算 \[\begin{aligned} \bar R = \frac{1}{N} \sum_{n=1}^N R_n, \ \overline{R^2} = \frac{1}{N} \sum_{n=1}^N R_n^2, \ S^2 = \frac{1}{N}\sum_{n=1}^N (R_n - \frac{1}{2})^2 \end{aligned}\]\(N \to \infty\)\(\bar R\)\(\overline{R^2}\)\(S^2\)均渐近服从正态分布, 可以用Z检验法检验这三个统计量与理论期望值的偏离程度。

  • \([0,1]\)等分成\(k\)段,用拟合优度卡方检验法检验\(\{ R_n, n=1,2,\dots,N \}\) 落在每一段的取值概率是否近似为\(1/k\)

  • 用Kolmogorov-Smirnov检验法进行拟合优度检验, 看\(\{ R_n, n=1,2,\dots,N \}\)是否与U[0,1]分布相符。

  • \(\{ R_n, n=1,2,\dots,N \}\)\(d\)个组合在一起成为\(\mathbb R^d\)向量, 把超立方体\([0,1]^d\)每一维均分为\(k\)份,得到\(k^d\)个子集, 用卡方检验法检验组合得到的\(\mathbb R^d\)向量落在每个子集的概率是否近似为\(k^{-d}\)

  • \(\{ R_n, n=1,2,\dots, N \}\)看作时间序列样本,计算其样本自相关函数列 \(\{ \hat\rho_j, j=1,2,\dots \}\), 在\(\{ R_n, n=1,2,\dots \}\)独立同分布情况下 \(\{ \sqrt{N}\hat\rho_j, j=1,2,\dots,N \}\)应该渐近服从独立的标准正态分布, 可以据此进行白噪声检验。

  • \(\{ R_n \}\)离散化为\(y_n = \text{floor}(k R_n), n=1,2,\dots, N\), 令\(\xi_n = y_n\), \(\eta_n = y_{n+b}\)(\(b\)为正整数), \(n=1,2,\dots, N-b\), 用列联表检验法检验\(\xi_n\)\(\eta_n\)的独立性。

  • 游程检验。把序列\(\{ R_n, n=1,2,\dots,N \}\)按顺序切分为不同段, 每一段中的数值都是递增的,这样的一段叫做一个上升游程, 如\((0.855)\), \((0.108, 0.226)\), \((0.032,0.123)\), \((0.055, 0.545, 0.642, 0.870)\), \(\ldots\)。 在相邻的两个上升游程中前一个游程的最后一个值大于后一个游程的第一个值。 在\(\{ R_n, n=1,2,\dots, \}\)独立同分布条件下游程长度的理论分布可以得知, 然后可以比较实际游程长度与独立同分布条件下期望长度。

  • 扑克检验。把\(\{ R_n \}\)离散化为\(y_n = \text{floor}(8 R_n), n=1,2,\dots, N\), 然后每连续的8个作为一组, 计数每组内不同数字的个数(\(1\sim8\)个)。 在\(\{ R_n \}\)独立同均匀分布条件下每组内不同数字个数的理论分布概率可以计算出来, 然后用卡方检验法检验实际观测与理论概率是否相符。

  • 配套检验。\(\{ R_n \}\)离散化为\(y_n = \text{floor}(k R_n)\)(\(k\)为正整数), \(n=1,2,\dots, N\), 然后顺序抽取\(\{ y_n, n=1, 2,\dots,N \}\)直到\(\{ y_n \}\)的可能取值 \(\{0,1,\dots,k-1\}\)都出现过为止, 记录需要抽取的\(\{ y_n \}\)的个数\(L\), 反复抽取并记录配齐数字需要抽取的值的个数\(l_j, j=1,2,\dots\)。 在\(\{ R_n \}\)独立同U(0,1)分布条件下这样的\(L\)分布可以得到, 可以计算\(\{ l_j \}\)的平均值并用渐近正态分布检验观测均值与理论均值的差异大小, 或直接用卡方检验法比较\(\{ l_j \}\)的样本频数与理论期望值。

  • 正负连检验。令\(y_n = R_n - \frac12\), 把连续的\(\{ y_n \}\)的正值分为一段, 把连续的\(\{ y_n \}\)的负值分为一段, 每段叫做一个“连”, 连长\(L\)的分布概率为\(P(L=k) = 2^{-k}, k=1,2,\dots\), 可以用卡方检验法检验\(L\)的分布; 总连数\(T\)满足\(ET = \frac{n+1}{2}\), \(\text{Var}(T) = \frac{n-1}{4}\), 可以用Z检验法检验\(T\)的值与理论期望的差距。

  • 升降连检验。计算\(y_n=R_n - R_{n-1}\), \(n=2,3,\dots,N\), 把连续的正值的\(\{ y_n \}\)分为一段叫做一个上升连, 把连续的负值的\(\{ y_n \}\)分为一段叫做一个下降连, 可以用卡方检验法比较连的长度与\(\{ R_n \}\)独立同分布假设下的理论分布, 或用Z检验法比较总连数与理论期望值的差距。

习题

习题1

对线性同余发生器 \[ \begin{aligned} x_n = (a x_{n-1} + c) \pmod{M}, \ n=1,2,\dots \end{aligned} \] 证明若从初值\(x_0\)出发周期等于\(M\), 则以\(0 \sim M-1\)中任何一个整数作为初值时周期都等于\(M\); 若从某初值\(x_0\)出发周期小于\(M\), 则以\(0 \sim M-1\)中任何一个整数作为初值时周期都小于\(M\)

习题2

编写R程序,对某个均匀分布随机数发生器产生的序列做均匀性的卡方检验。

习题3

设随机变量\(X\)密度函数\(p(x)\)和分布函数\(F(x)\)已知, 编写R程序,对\(X\)的某个随机数发生器产生的序列做拟合优度卡方检验。

习题4

设随机变量\(X\)的分布密度为 \[\begin{aligned} p(x) = \frac12 + x, \ 0 \leq x \leq 1, \end{aligned}\] 分别给出用逆变换法、舍选法、复合法生成\(X\)随机数的算法, 并比较三种方法的效率。

参考文献

Wichmann, B. A., and I. D. Hill. 1982. “Algorithm as 183: An Efficient and Portable Pseudo-Random Number Generator.” Applied Statistics 31: 188–190. Remarks: 34, 198 and 35, 89.