Loading... # 前言 前一篇有讲过蒙特卡罗方法需要大量的模拟抽样获得结果，在计算机领域我们通常使用确定性算法生成的随机数来达到我们的目的。那么我们有必要知道这样的随机数是否能够满足我们的模拟需要 # 随机数发生器 ## 确定性算法 如果有在代码里使用过随机数的同学应该知道，正常情况下我们直接输出一个随机数的值，实验3次，3次得到的随机数都是不相同的；然而如果我们在输出随机数前设置一个随机数种子，那么我们就会得到相同的的随机数结果，如下图Python程序所示。 ```python import random for i in range(3): random.seed(996) print('第{}次获取随机数的结果为{:.6f}'.format(i+1, random.random())) # 代码结果 # 第1次获取随机数的结果为0.915454 # 第2次获取随机数的结果为0.915454 # 第3次获取随机数的结果为0.915454 ``` 通过设置随机数种子能够获取到相同的随机数结果，是因为计算机生成随机数的底层算法，在给定特定输入的情况下，总是会产生相同的输出[<sup>[1]</sup>](#ref1)。当我们确定了随机数种子后，其实也就确定了这个算法的初始状态，根据这个初始状态，我们能够确定此后一系列的状态序列的值。 > 可以想象，我们通过初始状态$x_{0}$，通过算法可以唯一确定一个输出$x_{1}$，而$x_{1}$又是下一个算法过程的初始状态，因此可以确定出$x_{2}$。以此类推，也就是给定初始$x_{0}$，可以确定出一组状态序列。 ## 随机数发生器 一般来说，我们把生成随机数的确定性算法称之为随机数发生器，按照前面的介绍，其实我们也不难定义出随机数发生器的基本数学结构，即 $$ \begin{aligned} x_{i+1} &= f(x_{i}) \end{aligned} $$ 其中$x_{i}$是我们的第$i$个随机数，$f$是我们的随机数发生器，**特别的，我们规定$x_{0}$为我们的初值，也称之为种子(seed)**。 为什么说随机数发生器产生的是伪随机数而不是真正的随机数呢？这是因为随机数发生器生成的随机数是存在周期的，在给定了随机数生成器后，我们可以根据算法推演出其随机数周期。 假设随机数的周期为$M$，第$M+1$随机数的值和第1个随机数的值必定是相等的，即$x_{M+1}=x_{1}$。那么，我们如何保证我们的伪随机数具备随机数那样的统计特征呢？一般来说需要做到以下两点： - 通过随机数检验，比如均匀性检验、独立性检验等。这些检验主要是保证生成的序列是在统计上满足我们的需求的。 - 有较长的循环周期。循环周期长，我们生成的序列就越不容易被检验出规律，其随机数性质也就越好 在统计推断中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布$U(0,1)$的随机数。通常情况下，我们会用另一种算法$g$，通过$u_{t+1}=g(x_{t+1})$，将产生的随机数映射到$[0, 1]$区间，综合起来即： $$ \begin{aligned} x_{i+1} &= f(x_{i}) \\ u_{t+1} &= g(x_{t+1}) \end{aligned} $$ ## 常见的随机数生成器：线性同余生成器 常见的一种随机数生成器为线性同余生成器，这种生成器的基本形式为 $$ \begin{aligned} x_{i+1} &= ax_{i}+c \quad mod \quad m \end{aligned} $$ 其中$a, c, m$为三个常数，且满足$a, c, m \in Z$ 假如我们设$a=3, c=0, m=5$，指定种子为996，生成10个随机数，得到的结果为[3, 4, 2, 1, 3, 4, 2, 1, 3, 4]。从结果可以看到这个生成器的周期为4，这样的随机数序列是没法直接应用的。 如果我们尝试调整上面的三个参数，比如$a=10000, c=0, m=9145$，可以得到[1095, 3435, 1380, 195, 2115, 6760, 160, 8770, 8595, 5290]。仅看这10个数的结果，是不是感觉这10个数似乎看起来有那么一点“随机”了呢？这就是为什么我们需要找到好的随机数算法，需要达到我们前面提到的两个条件的原因。 > 可以证明，在给定参数$m$的情况下，该生成器的周期与$a$和$x_{0}$有关 线性同余生成器是最简单的一种伪随机数算法，常见的算法还有 - 冯诺依曼平方取中法（midsquare method） - 滞后斐波那契发生器（lagged Fibonacci generator） - ... 在Python中，我们使用的random模块的随机数发生器为Mersenne Twister算法，对随机数发生器感兴趣的同学可以自行百度/Google。 ## 简要介绍：将伪随机数映射到$[0, 1]$区间（标准化） 由于使用时大多数情况下都希望得到基于均匀分布$U(0,1)$的随机数，因此需要一些方式将上述的随机数序列进行区间映射，这类操作通常叫做标准化。其中最常用的方法就是最大最小值法，即利用一个序列的最大、最小值将随机数进行缩放，用数学语言描述如下： $$ \begin{aligned} u_{t+1} &= \frac{x_{t+1}-x_{min}}{x_{max}-x_{min}} \end{aligned} $$ 接着前面的例子，假设我们能够得到的随机数序列是[3, 4, 2, 1, 3, 4, 2, 1, 3, 4]，那么映射后得到序列[0.67, 1.0, 0.33, 0.0, 0.67, 1.0, 0.33, 0.0, 0.67, 1.0] 当然还有很多方法，具体可以自行百度。这些方法不是唯一的，只要合理即可，比如前例还可以使用更简单的方式进行标准化： $$ \begin{aligned} u_{t+1} &= \frac{x_{t+1}}{m} \end{aligned} $$ ## 逆变换方法（Inverse Transform Method） 由于实际研究中我们遇到的通常并不是均匀分布，而是各种各样的概率分布，因此对我们来说，获得$U(0,1)$的随机数并不是重点，更重要的是需要通过这些均匀分布的随机数，生成我们需要的分布的随机数。为了达到我们的目的，需要介绍如下定理 > 反函数定义补充：设随机变量$X$的累积分布函数为$F_{X}(x)$，那么定义$F_{X}^{-1}(y)=inf\{x:F_{X}(x) \geq y \}$为$F_{X}(x)$的反函数。 **定理**： 如果随机变量$U$满足$U \sim (0, 1)$，那么随机变量$X=F_{X}^{-1}(U)$的累积分布函数为$F_{X}(x)$ 要想证明这一点，只需要证明随机变量$X$的分布函数$Pr(X \leq x)$将$X$用$F_{X}^{-1}(U)$替换后仍然能得到原来的分布即可。 证明： $$ \begin{aligned} Pr(X \leq x) = Pr(F_{X}^{-1}(U) \leq x) = Pr(U \leq F_{X}(x))=F_{U}(F(X))=F_{X}(x) \end{aligned} $$ 这个定律的意义在于，只要随机变量$X$的累积分布函数可逆，那么**我们也就可以用服从均匀分布去构建随机变量$X$的分布。** > 累积分布函数可逆的条件为该函数严格单调递增 整体步骤如下： - 对于随机变量$X$，利用均匀分布的随机变量和反函数，构建出$F_{X}^{-1}(U)$ - 生成$U(0,1)$的随机数，并代入$F_{X}^{-1}(U)$进行计算 假设$F_{X}(x) = 1-e^{-\frac{x}{\lambda}}(x > 0)$，利用反函数计算，可以得出 $$ \begin{aligned} X = F_{X}^{-1}(U) = -\lambda \log (1-U) \quad U \in (0, 1) \end{aligned} $$ 此时我们生成均匀分布随机数，通过上述方法，自然就可以得到服从随机变量$X$的分布的随机数了。 > 这个方法非常重要，后面我们还是会再重复提到的 # 附录 ## 线性同余生成器代码（Python） ```python # 线性同余生成器代码 def linear_congruence_generator(iter_num: int, a: int, c: int, m: int, x0: int=996): cnt = 0 while True: if cnt < iter_num: if cnt == 0: x = (a * x0 + c ) % m else: x = (a * x + c ) % m yield x cnt += 1 else: raise StopIteration # 最大最小值变换 def min_max_transformer(random_num_list: list): x_max = max(random_num_list) x_min = min(random_num_list) uniform_list = [round((x - x_min) / (x_max - x_min), 2) for x in random_num_list] return uniform_list x0 = 996 a = 3 c = 0 m = 5 iter_num = 10 # 生成10个随机数 f = linear_congruence_generator(iter_num, a, c, m, x0) random_num_list = [next(f) for i in range(iter_num)] uniform_random_num_list = min_max_transformer(random_num_list) print(uniform_random_num_list) ``` # 参考文献 <div id="ref1"></div> [1] Deterministic Algorithm. (n.d.). https://en.wikipedia.org/wiki/Deterministic_algorithm 最后修改：2022 年 06 月 04 日 © 允许规范转载 打赏 赞赏作者 赞 0 如果觉得我的文章对你有用，请随意赞赏