蒙特卡罗法
1. 模型原理
提出:
蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”的成员S.M.乌拉姆和J.冯·诺依曼首先提出。数学家冯·诺依曼用驰名世界的赌城一摩纳哥的Monte Carlo一来命名这种方法,为它蒙上了一层神秘色彩。在这之前,蒙特卡洛方法就已经存在。1777年,法国Buffon提出用投针试验的方法求圆周率,这被认为是蒙特卡罗方法的起源。
定义:
蒙特卡罗法又称统计模拟法,是一种随机模拟方法,以概率和统计理论方法为基础的一种计算方法,是使用随机数(或更常见伪随机数)来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系,用电子计算机实现统计模拟或抽样,以获得问题的近似解。为象征性地表明这一方法的概率统计特征,故借用赌城蒙特卡罗命名。
原理:
由大数定理可知,当样本容量足够大时,事件的发生频率即为其概率。
注意:蒙特卡洛不是一种算法,准确的来说只是一种思想,或者是一种方法,只要求解的问题与概率模型有关联,我们就可以采用这种方法,从数学建模的角度来看,蒙特卡洛是没有通用的代码的,每个问题对应的代码都是不同的。
2. 常见问题
2.1 圆周率问题
一个半径为1的圆,其外切正方形面积易知为4,若在正方形内随机撒大量的点,有些落在园内,有些落在圆外。
统计意义上:
圆内点数
总点数
=
圆面积
正方形面积
=
π
4
π
=
4
×
圆内点数
总点数
\frac{圆内点数}{总点数}=\frac{圆面积}{正方形面积}=\frac{\pi}{4}\\ \pi=4\times \frac{圆内点数}{总点数}
总点数圆内点数=正方形面积圆面积=4ππ=4×总点数圆内点数
求解的python代码:
import numpy as np
import matplotlib.pyplot as plt
# 参数初始化,投放10000个点,圆的半径为1,圆心为(1,1)
p = 10000 # 投放点的个数
r = 1 # 圆的半径
x0, y0 = 1, 1 # 圆心的坐标
n = 0 # 初始时还未投放点,有n个点落在圆内
# 设置绘图窗口
plt.figure()
plt.title("Monte Carlo Simulation of Pi")
plt.xlabel("x")
plt.ylabel("y")
# 保持绘制窗口,多次绘图
for i in range(p):
px = np.random.rand()*2 # 生成[0,2)之间的随机数作为点的横坐标
py = np.random.rand()*2 # 生成[0,2)之间的随机数作为点的纵坐标
# 判断点是否在圆内
if (px-x0)**2+(py-y0)**2 <= r**2:
plt.plot(px, py, '.', color='b') # 点在圆内,蓝色
n += 1
else:
plt.plot(px, py, '.', color='r') # 点在圆外,红色
plt.axis('equal') # 设置x,y轴的单位长度相等,保证圆不变形
plt.show()
# 计算pi的值
s = 4*n/p # s为圆的面积,p为正方形的面积,圆的面积为pi*r^2,正方形的面积为(2r)^2
pi_estimation = s
print("Estimated pi is", pi_estimation)
输出:
Estimated pi is 3.1416
![Figure_1](https://img-blog.csdnimg.cn/img_convert/9f6769c766f0ef43cd080790bd61a1b5.png)
2.2 三门问题
你参加一档电视节目,节目组提供了A、B、C三扇门,主持人告诉你,其中一扇门后边有辆汽车,其他两扇门后面是一头山羊,你可以选择一扇门打开获得门后的东西。假如你选择了B门,这时,主持人打开了C门,让你看到C门后是只山羊,然后问你要不要改选A门?(你想要汽车)
![]()
求解的python代码:
import numpy as np
n = 10000 # n代表蒙特卡罗模拟重复的次数
a = 0 # a代表不改变主意时能赢得汽车的次数
b = 0 # b代表改变主意时能赢得汽车的次数
# 重复n次蒙特卡罗模拟
for i in range(n):
door_car_in = np.random.randint(1, 4) # 生成1-3之间的随机数,代表汽车所在的门
door_selected = np.random.randint(1, 4) # 生成1-3之间的随机数,代表参赛者选择的门
# 分两种情况讨论
if door_car_in == door_selected: # 如果汽车在参赛者选择的门后
a += 1 # 不改变主意时能赢得汽车的次数加1
else: # 如果汽车不在参赛者选择的门后
b += 1 # 改变主意时能赢得汽车的次数加1
print("不改变主意时能赢得汽车的概率:", a/n)
print("改变主意时能赢得汽车的概率:", b/n)
输出:
不改变主意时能赢得汽车的概率: 0.3275
改变主意时能赢得汽车的概率: 0.6725