【高级软件实习】蒙特卡洛模拟 | PRNG 伪随机数发生器 | LCG 线性同余算法 | 马特赛特旋转算法

写在前面：本篇博客将介绍经典的伪随机数生成算法，我们将重点讲解 LCG(线性同余发生器) 算法与马特赛特旋转算法，在此基础上顺带介绍 Python 的 random 模块。本篇博客还带有练习，无聊到喷水的练习，咳咳…… 学完前面的内容你就会了解到 Python 的 Random 模块的随机数生成的实现，是基于马特赛特旋转算法的，比如 random_uniform 函数。而本篇博客提供的练习会让你实现一个基于 LCG 算法的random_uniform，个人认为还是比较有意思的。练习题的环境为 Google Colaboratory（K80 GPU）Jupyter Notebook：https://colab.research.google.com/

Ⅰ. 蒙特卡洛模拟（Monte Carlo simulation）

0x00 引入：概念说明

引入：在 非确定性模型（deterministic model）的 概率模型（stochastic model）中，用分析法难以下手，这时我们就可以利用 “蒙特卡洛模拟” (或称蒙特卡洛仿真) 去解决问题。

概念：“蒙特卡洛模拟” 是一种反复使用随机采样技术进行模拟，并从总体概率分布计算出所需数值结果的一种计算算法。

命名的由来

蒙特卡洛模拟是在二战期间，在原子弹研制的项目中用来模拟裂变物质的中子随机扩散现象，由美国数学家冯·诺伊曼和乌拉姆研发的模拟算法。蒙特卡洛是在欧洲摩纳哥的一个城市，这个城市在当时是非常著名的一个赌城。因为赌博的本质是算概率，而蒙特卡洛模拟正是以概率为基础的一种方法，所以用赌城的名字为这种方法命名。

（摩纳哥公国 – 蒙特卡洛）

蒙特卡洛模拟可作用于数值积分、概率分布计算、基于概率的优化等领域。例如：

Simulation：为了模拟抛硬币，在

范围内取随机值：

将小于的视为
将大于的视为

Monte Carlo method：倒一盒硬币，数

和

的数字，由此获得

的概率。

Monte Carlo simulation：在

的范围内，反复提取随机值，由此获得

的概率。蒙特卡罗算法表示采样越多，越近似最优解。

百度百科：百科举了个容易理解的例子

举个例子，假如筐里有100个苹果，让我每次闭眼拿1个，挑出最大的。

于是我随机拿1个，再随机拿1个跟它比，留下大的，再随机拿1个……我每拿一次，留下的苹果都至少不比上次的小。拿的次数越多，挑出的苹果就越大，但我除非拿100次（即把筐中的苹果都拿了一遍），否则无法肯定是否挑出了最大的那一个。

这个挑苹果的算法，就属于蒙特卡罗算法。

告诉我们样本容量足够大，则最接近所要求解的概率。

0x01 蒙特卡洛仿真常规模式

统共分为如下四步：

STEP1：定义要进行采样的区域（Range）
STEP2：对定义的区域进行随机采样（Random sampling）
STEP3：对收集的样本进行确定型计算（Deterministic computation）
STEP4：统计结果并导出近似值（Approximate）

举个例子：两个简单的蒙特卡罗模拟的实例

① 掷两次骰子点数之和为
的概率：

STEP1：生成两个范围的随机数（掷两次骰子），并进行次随机运行模拟
STEP2：若生成的两个随机数之和为，则为，否则为，统计和的次数
STEP3：计算与总执行次数之比
STEP4：与数学计算的概率进行比较

② 求圆周率：

（撷取自维基百科）

利用以表现的圆去计算圆周率
该圆被包含于以表现的宽为 4 的正方形空间内。
范围限制在内
在此空间内提取个的随机数有序对（ordered pair）
在提取出的点中，统计圆内部点的个数。
利用圆内部的点与整体的点的个数比求出

Ⅱ. 伪随机数发生器（PRNG）

0x00 引入：随机数生成问题

为了给蒙特卡洛模拟生成随机数，我们需要置随机数 “种子”。

这里我们选择采用伪随机数发生器 (Pseudo Random Number Genarator) ，简称

。

伪随机数并非完全随机

伪随机数发生器用于在系统需要随机数的时候，通过一系列种子值计算出来的伪随机数。因为生成一个真正意义上的“随机数”对于计算机来说是不可能的，伪随机数也只是尽可能地接近其应具有的随机性，但是因为有“种子值”，所以伪随机数在一定程度上是可控可预测的。

本章我们将介绍以下几种经典的伪随机数生成算法：（我们将重点讲解前两种算法）

线性同余发生器（Linear congruential generator）
马特赛特旋转演算法（Mersenne Twister）
冯·诺伊曼平方取中法（John von Neumann mid-square method）
乘法取中法（mid-product method）
常数乘子法（constant multiplier method）
斐波那契法（Fibonacci Method）
时滞斐波那契法（Lagged Fibonacci method）
移位法（Shift Method）

0x01 线性同余发生器（LCG）

线性同余发生器（Linear Congruential Generator），简称

。

它能产生具有不连续计算的伪随机序列的分段线性方程，生成器由循环关系定义如下：

【高级软件实习】蒙特卡洛模拟 | PRNG 伪随机数发生器 | LCG 线性同余算法 | 马特赛特旋转算法 | Python Random 模块由下列参数唯一决定：

，，，【高级软件实习】蒙特卡洛模拟 | PRNG 伪随机数发生器 | LCG 线性同余算法 | 马特赛特旋转算法 | Python Random 模块

为保持

的最大周期

的充分必要条件如下：

与互质（coprime）
的所有质因数都能整除
若是的倍数则也是的倍数

其中参数 a, C, m 比较敏感，它们直接影响了伪随机数产生的质量，所以参数的取值非常重要！

算法优点：

能以较小的内存去生成随机数，因此适用于嵌入式开发环境。

有可观的计算速度，且在不需要考虑相关性的情况下也能使用，实现起来非常简单。

算法缺点：

如果知道了的参数和最后生成的随机数，就可以预测其后生成的所有因数，因此从密码学角度来看，我们不能认为它是安全的随机数生成器。

0x02 马特赛特旋转算法（MersenneTwister）

“马特赛特旋转算法是当今生成随机数质量最好的算法”

算法实现的伪随机数很不错，但是周期不够长，很容易被 “不好的人” 推算出随机数种子。两个小日子过得不错的学者 —— 松本和西村，在 1997 年又研究提出了新的伪随机数算法：

马特赛特旋转算法（MersenneTwister），简称

。是基于有限二进制字段上的矩阵线性递归。马特赛特旋转算法可迅速产生高质量的伪随机数，修正了古老随机数产生算法的很多缺陷。该算法中用到了梅森素数，所以又称为梅森旋转算算法。（但翻译大部分选用了 “马特赛特” 作为翻译，所以这里我们还是称之为马特赛特旋转算法）

补充：梅森素数

梅森素数的意思是形如的素数。利用梅森素数的性质可以设计出周期长度为梅森素数长度的随机数周期。比如目前Python、C++11等语言当中用的随机数计算包都是用的这种算法。目前常用的版本周期是，这是一个巨大的天文数字。

即使你从来没有听过这个算法，但你一定用过 Python 中的 random 模块：

Python 中的 random 模块就是采用了马特赛特旋转算法实现随机数的，不仅 Python 用，PHP、Perl 等流行编程语言内置的

也都是采用该算法实现，可见马特赛特旋转算法真是个备受大佬青睐的算法，它的地位可见一斑！

算法的三个阶段：

第一阶段：获得基础的马特赛特旋转链
第二阶段：对于旋转链进行旋转算法
第三阶段：对于旋转算法所得的结果进行处理

算法优点：

马特赛特旋转能更快的生成高质量的伪随机数，随机数的反复周期取决于马特赛特的因数。

在 MT19937 上跑，周期最长可达到，具有 623 维均匀分布的性质，真的很大！（使用 624 个整型，623*32=19936！）它的序列关联比较小，能通过多种随机性测试，甚至包括 Hdamard Test 量子计算！

只需位运算就能实现的算法，速度异常迅速。

算法缺点：

该算法因其必须分配能容纳个数字的空间，对容量有限的嵌入式环境来说是一大缺陷。

从密码学的角度来看，马特赛特旋转算法仍然是不安全的，当你知道随机数的特性时，仅利用有限的随机数（624个）就可以知道当前生成器所处状态并预测出即将出现的随机数。

0x03 冯·诺伊曼平方取中法（Mid-square Method）

平方取中法（Mid-square Method）是最早出现的伪随机数发生器，最早由冯诺依曼提出。

平方取中法的意思如同其名：

平方：是将一个 2s 位十进制随机数开平方后得到的一个 4s 位数（不足 4s 时高位补 0）。
取中：去头截尾，取中间 2s 位数作为一个新的随机数。

并将此数规范化（化为小于1的2s位数值），即第一个

上的随机数，以此类推，重复以上过程就能得到一系列伪随机数。其递推公式如下：

产生伪随机数数列：

平方取中法的优点：

易于实现，内存占用少，计算简单。

平方取中法的缺点：

存对小数目偏倚现象，均匀性差，对初始数据依赖大，很难说明取什么样的种子值可以保证有足够长的周期。

很容易出现重复元素的段循环，而且容易退化为一常数，甚至退化为零，如果某个元素退化为0，那么后面所有元素都将会是0。

平方取中法的变形和推广：