Python蒙特卡洛模拟 | LCG 算法 | 马特赛特旋转算法 | Random 模块

Posted 2022-12-08 柠檬叶子C

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

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

---

Ⅰ. 蒙特卡洛模拟（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），简称 。

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

 由下列参数唯一决定：

 ，，  ，   

📜 为保持   的最大周期  的充分必要条件如下：

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

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

✅ 算法优点：

•  能以较小的内存去生成随机数，因此适用于嵌入式开发环境。
• 有可观的计算速度，且在不需要考虑相关性的情况下也能使用 ，实现起来非常简单。

❌ 算法缺点：

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

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

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

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

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

补充：梅森素数

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

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

import 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。

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

• 乘法取中法（mid-product method）： 

若要产生具有10进制 2s 位的伪随机数序列，任意选取两个初始随机数 ，递推公式：

产生伪随机数数列：

• 常数乘子法（constant multiplier method）：

递推公式：

产生伪随机数数列：

0x04 斐波那契法（Fibonacci Method）

遗憾的是，该方法并不是一个好的伪随机数发生器，该方法基于 Fibonacci 数列，其递推公式为：

     

该方法有   两个初始种子，最大的优点就是计算周期快，且达到满周期。此发生器没有乘法运算，产生速度很快。缺点是序列中的数可能会重复出现，独立性较差，且有让人难以容忍的不居中现象。

📜 斐波那契法的变形和推广（了解）：

• 时滞斐波那契法（Lagging Fibonacci Method）：

时滞斐波那契生成器，简称 。时滞斐波那契生成器的理论相当复杂，理论也不够充分到能指导如何选择  与  。生成器的初始化也非常敏感，其递推公式为：

其中，新项由两个老项计算生成。 通常是 2 的幂（）， 经常是  或 。其中 ★ 运算符表示一般的二元运算符，这可以是加法、减法、乘法或者位运算异或。

0x05 移位法（Shift Method）

计算机善于进行移位等逻辑运算，移位法就是利用计算机的这个特点去实现的。移位法是关于平方取中法的另一种改进形式，它是移位运算于指令相加运算相结合的迭代过程。方法如下：

取一个初始值 ，使它左右移位，分别为  ，然后指令相加得到 ，再对  进行上述过程，如此重复下去，就能产生随机数序列了。递推公式如下（适用于32位机器）：

产生伪随机数数列：

移位法运算速度快，但是对初始值的依赖性很大，一般初始值不能取的太小，初始值选的不好会使伪随机数序列长度较短。同时独立性交叉，随后随机数序列周期  与计算机的字长有关。

Ⅲ. Python 的随机数生成函数（Random 模块）

0x00 引入：Random 模块介绍

我们刚才讲解马特塞特旋转算法时提到了，Python 中的 Random 模块就是采用该算法实现的。

其生成  位精度的 float，周期为  。

提供了可提取分布的函数：Uniform，Normal，lognormal，negative exponential，gamma, beta 

 📚 使用前需引入 random 模块：

import random

0x01 random - 生成 0.0 ~ 1.0 间的随机数

random.random()   # 在 0.1 与 1 之间的实数中生成随机数

📚 作用：随机生成一个  范围的实数。

💬 代码演示：生成实数  之间的随机数：

print(random.random())

🚩 运行结果如下：

0x02  uniform - 生成指定范围的随机数

random.uniform(a, b)   # 在 a 和 b 之间生成随机数

📚 作用：随机生成一个指定范围的实数。

💬 代码演示：生成一个 100~300 之间的数

print(random.uniform(10, 30))

🚩 运行结果如下：

0x03  randint - 生成指定范围的整数随机数

random.randint(a, b)      # 在 a 和 b 之间生成整数随机数

📚 作用：随机生成一个指定范围的整数。

💬 代码演示：生成一个100~300之间的整数

print(random.randint(100, 300))

🚩 运行结果如下：

0x04  choice - 在样本集中随机选择

random.choice(sample_set)    # 在样本集 sample_set 中随机选择一个样本

📚 作用：在样本集中随机选择一个样本

💬 代码演示：从样本集  中随机选择一个样本

L = [1, -3, 5, -2, 6, 8, -3, 0]  
print(random.choice(L))

🚩 运行结果如下：

0x05 sample - 在样本集中随机选择多次

random.sample(sample_set, n)    # 在样本集 sample_set 中随机选择 n 个样本

📚 作用：执行 n 次 "在样本集中随机选择一个样本" 。

💬 代码演示：从样本集  中随机选择多次

L2 = [1, 'CSDN', 5, 'bilibili', 6, -3, 0]
print(random.sample(L2, 3))   # 在样本集L2中选择3次
print(random.sample(L2, 6))   # 在样本集L2中选择6次
print(random.sample(L2, 3))   # 在样本集L2中选择3次

🚩 运行结果如下：

Ⅳ. 练习

练习1：绘制对数直方图（利用LCG实现）📈

• 在 Colaboratory 中使用 Python 实现。
• 利用线性同余发生器（LCG）实现。
• 生成  范围内的随机数。
• 使用实现的  将生成随机数的过程重复多次，并在如下次数时生成随机数的分布：
  • 分别在第  次、第  次、第 次和第  次生成。
• 使用 Python 的 Matplotlib 模块，以生成直方图的形式展现分布，如下所示：

import random
import matplotlib.pyplot

我们再回过头来看一下  的循环关系定义：

我们刚才提到过，LCG 算法的 a, C, m 参数比较敏感，它们直接影响了伪随机数产生的质量，所以这里的取值是很重要的！这里介绍一种业内常用的参数取值：

a = 25214903917
C = 11
m = 2**48

这些值都是经过专家精密计算得到的，并非是信手拈来，随随便便取的！

参考代码：（仅供参考）

import matplotlib.pyplot as plt
from tqdm import tqdm # optional
import random


# my LCG
def LCG(x, a, c, m):
  while True:
      x = (a * x + c) % m
      yield x


num_iterations = 100    # 次数
random_integers = []

def random_uniform_sample(n, space, seed=0):
  # Random seed
  a = 214013
  C = 2531011
  m = 2**32
  
  get = LCG(seed, a, C, m)
  lower, upper = space[0], space[1]

  # create the random value
  for i in range(n):
    value = (upper - lower) * (next(get) / (m - 1)) + lower
    random_integers.append(value)   

  return random_integers     # return the result


X = random_uniform_sample(num_iterations, [0, 99])
# print(X)

fig = plt.figure()
plt.hist(X)
plt.title(f"Check Uniform Distribution of num_iterations iterations")
plt.xlabel("Numbers")
plt.ylabel("N of each number")
plt.xlim([0, 99])
plt.show()

🚩 运行结果：100次

🚩 运行结果：1000次

 

🚩 运行结果：10000次 

🚩 运行结果：100000次

 

练习2：检验 LCG 和 马特赛特旋转 中谁生成的分布更接近 🔍

利用  LCG（刚才练习1中自己实现的） 和 MT算法（Python random 模块)  查看分布， 查看哪种随机数生成器更符合均匀分布。

• 利用 Python random 模块的 randint（1,6）得到随机的骰子点数
• 记录第一个掷骰子的点数和第二个掷骰子的点数
• 点数之和为 8 时，视为 “命中” 。
• 使用蒙特卡洛算法，求两次掷骰子时点数之和为 8 的概率，确认误差率 20%
• 按下列形式输出结果（颜色可加可不加，代码达到效果即可）：

参考代码：

# 模拟掷骰子
def roll_dice():  
  roll = random.randint(1, 6)   # 按题目要求，使用randint在1~6范围取值
  return roll                   # 将结果返回

# 记录区
dice_tries1 = []
dice_tries2 = []
num_iterations = 100    # 次数
hits = 0   # 命中数

# 投掷次数
for i in range(num_iterations):
    dice_tries1.append(roll_dice())  # 投掷100次 存入数组
    dice_tries2.append(roll_dice())  # 投掷100次 存入数组

print("="*100)  

print(colored("* 红色表示两个骰子点数之和为8*", 'red'))   # 打印标题

change = 0
for j in range(num_iterations):
  if (change == 5):
      print()
      change = 0
  if (dice_tries1[j] + dice_tries2[j] == 8):
    print(colored("try %2s : %s %s" % (j, dice_tries1[j], dice_tries2[j]), "blue"), end = " ")
    change += 1
    hits += 1
  else:
    print("try %2s : %s %s" % (j, dice_tries1[j], dice_tries2[j]), end = " ")
    change += 1


print(colored("\\n实际值 : 0.138889", "red"))
print(colored(f"计算值 : round(hits / num_iterations,6)", "red"))
print(colored(f"误差率 : abs(hits / num_iterations - 5/36) / (5/36) * 100 %", "red"))
print("="*100)

练习3：骰子游戏胜利预测 🎲

玩家在进行骰子游戏，初始拥有费用为1000美元。
掷两个骰子，如果骰子的数字相同，得到 4 美元，骰子的数字如果不同，将失去 1美元
使用蒙特卡罗算法，根据游戏执行次数计算玩家的平均速率，确认平均余额。

 参考代码：

import matplotlib.pyplot as plt
import random

def roll_dice():
  A = random.randint(1, 6)
  B = random.randint(1, 6)
  return A == B

# Inputs
num_simulations = 10000
max_num_rolls = 1000
bet = 1

win_probability = []
end_balance = []


plt.figure(figsize=(12,6))

for i in range(num_simulations):  # 投一万次骰子
  arr = [1000]
  rolls = [0]
  wins = 0

  while (rolls[-1] < max_num_rolls):
    if (roll_dice() == True):
      arr.append(arr[-1] + 4 * bet)
      wins += 1
    else:
      arr.append(arr[-1] - bet)

    rolls.append(rolls[-1] + 1)

  win_probability.append(wins / rolls[-1])
  end_balance.append(arr[-1])
  plt.plot(rolls, arr)



plt.title("Monte Carlo Dice Game [" + str(num_simulations) + " players simulations]")
plt.xlabel("Roll Number")
plt.ylabel("Balance [$]")
plt.xlim([0, max_num_rolls])

…… 待更新

📌 [ 笔者 ]   王亦优
📃 [ 更新 ]   2022.11.11
❌ [ 勘误 ]   /* 暂无 */
📜 [ 声明 ]   由于作者水平有限，本文有错误和不准确之处在所难免，
              本人也很想知道这些错误，恳望读者批评指正！

📜 参考资料  Microsoft. MSDN(Microsoft Developer Network)[EB/OL]. []. . 百度百科[EB/OL]. []. https://baike.baidu.com/.