数值模拟
数值模拟
用计算机模拟自然现象，达到研究现实问题的目的
Monte Carlo simulation
蒙特卡罗模拟
- 数值模拟的一种
- 利用随机数进行数值模拟，以获得问题的近似解
蒙特卡罗算法的思想，最早出现于1777年，法国数学家布丰使用投针法计算圆周率，算法实际诞生于1946年
19世纪50年代，摩纳哥两个小镇独立，皇室陷入破产边缘，卡洛琳王妃提出建造一个赌场创收，以城市命名：蒙特卡罗赌场创立
赌博（掷骰子、轮盘赌）是现代统计学和概率论的起源
轮盘赌的传奇故事
蒙特卡罗模拟方法实际诞生于曼哈顿计划，1946年
塔尼斯拉夫·乌拉姆打牌时候想计算52张牌中抽出同花顺的概率，
他发现通过计算无法完成，想到如果随机模拟抽牌几百次后查看同花顺出现的次数，就可以近似得出概率的方法
他把想法告诉了冯诺依曼，后者在ENIAC上编程实现了算法
为了保密，乌拉姆和和冯诺依曼的同事，尼古拉斯·梅特罗波利斯提议将这种算法命名为 蒙特卡罗算法
（一方面因为随机概率的关系。另一方面因为乌拉姆的叔叔经常在那里输钱）
使用这种概率型算法可以解决（近似解决）大量确定性算法无法解决的问题
蒙特卡罗模拟被大量应用于曼哈顿计划中的各种计算和模拟，
随后被大规模应用于各种领域
蒙特卡罗算法分支
- 蒙特卡罗积分（用于计算高维积分）
- 蒙特卡罗定位
- 蒙特卡罗搜索树（AlphaGo）
- 元启发算法（带有随机性的算法）等
蒙特卡罗模拟的三个主要步骤
- 构造或描述概率过程
- 实现从已知概率分布抽样
- 建立各种估计量
蒙特卡诺模拟的优缺点：
- 优点：简单、快速。省去了反复的数学推导演算过程，一般人也可以做实验
- 缺点：计算量大（越多越精确）；如果随机数实际上不随机（有难以察觉的相关性），会导致模拟失败
案例：圆周率pi的计算
圆周率公式：pi = C/d = 周长/直径
祖冲之使用“割圆法”计算圆周率
使用蒙特卡罗模拟求圆周率
正方形内部有一个相切的圆，半径为r，
正方形面积/圆面积 = 2r*2r/pi*r*r = 4r^2/pi*r^2 = 4/pi
在正方形内部随机生成一万个点，计算它们到中心点的距离，以判断是否落在圆内部
image.png
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle # 画圆
生成数据
n = 10000 # 生成多少随机数
r = 1 # 圆半径
a, b = 0, 0 # 原点坐标轴
xmin = a - r
xmax = a + r
ymin = b - r
ymax = b + r
xmin, xmax, ymin, ymax
# 带起始结束值的均匀分布随机数
x = np.random.uniform(xmin, xmax, n)
y = np.random.uniform(ymin, ymax, n)
x
array([ 0.82045309, 0.62834121, 0.56162509, ..., 0.31759717,
-0.72173351, 0.42144282])
可视化
fig = plt.figure(figsize=(8, 8)) # 父图
axis = fig.add_subplot(1,1,1) # 子图
axis.scatter(x, y, s=1) # 随机散点
axis.scatter(0, 0, s=100, color='r') # 原点
axis.plot([0, 0],[0, 1], color='black') # 半径
# x/y轴范围
axis.set_xlim([-1, 1])
axis.set_ylim([-1, 1])
# 网格
axis.grid(linewidth=0.2, alpha=0.5)
# 画圆
circle = Circle(xy=(a, b), radius=r, alpha=0.2, color='r')
axis.add_patch(circle)
<matplotlib.patches.Circle at 0x9295e48>
output_6_1.png
计算哪些点在圆内
# x、y坐标到中心坐标的距离，小于等于r，就在圆内
# 勾股定理
d = ((x - a) ** 2 + (y - b) ** 2) ** 0.5
d
res = np.sum(np.where(d <= r, 1, 0)) # 圆内点为1，圆外点为0，求和得出圆内点的个数
res
7891
求圆周率pi
# 正方形面积/圆面积 = 4r^2 / pi * r^2 = 4/pi
# 面积=点个数，所以：4/pi = 正方形面积/圆面积
# pi = 4 * 圆内点个数 / 正方形内点个数
pi = 4 * res / n
pi
3.1564
接近圆周率，增加随机数个数更接近（预防死机）
