设随机变量\(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 (同余) 设\(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)、(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\)的余数。

线性同余随机数发生器是利用求余运算的随机数发生器。 其递推公式为 \[ \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\)(\(m0\)时称为混合同余发生器。 下列定理给出了混合同余发生器达到满周期的一个充分条件。

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

常取\(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