蒙特卡洛方法求π值的可视化
什么是蒙特卡洛
蒙特卡络不是一个人名,而是一个地名,因摩纳哥著名的赌场而得名,而该方法的提出者是大名鼎鼎的数学家冯·诺伊曼(现代计算机之父)。
蒙特卡洛(Monte Carlo)方法,又称为随机抽样或统计试验方法,是以概率和统计理论方法为基础的一种计算方法,本质是使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。它将所求解的问题同一定的概率模型相联系,以获得问题的近似解。
这里要说明的一点是蒙特卡络方法是一种基于概率方法的统称,包含蒙特卡洛算法、模拟、过程、搜索树等。而且使用可能得到的数值并非是最终最精确的那个,会有一定的误差,而误差的大小与模拟的样本大小直接相关,模拟样本越大,误差越小,它是一种概率数值的逼近。对于简单问题来说,蒙特卡洛是个“笨”办法。但对许多问题来说,它往往是个有效,有时甚至是唯一可行的方法。
用蒙特卡洛方法求圆周率π值
考虑到π和圆形是紧密相关的,所以我首先构造一个边长为1的正方形和半径为1的圆形。为什么选正方形呢,因为正方形的边长相等,又和圆形半径类似,而且正方形在某个角度可以看做是平面坐标系,方便后面计算。
然后使两个图案相交,重叠部分是一个$\frac{1}{4}$个单位圆部分(将这个重叠部分简称为P(A))。
接着,在边长为1的正方形区域内随机落点,统计落在重叠部分P(A)区域内点的所占比例。这个比例其实很好计算,因为随着随机落点数量的增加,每个随机点的落点概率会越来越接近四分之一圆面积与正方形面积之比(概率比值即为$\frac{π}{4}$),而这个概率就是随机落点到重叠部分P(A)的比例。
π值求解的推导公式如下:
$P(A)=\frac{四分之一圆面积}{正方形面积}=\frac{\frac{1}{4}πr^2}{r^2}=\frac{π}{4}$
$P(A) = \frac{π}{4}$
$π = P(A)\times4$
考虑到用文字难以形象的表达清楚计算π的过程,下面用我的手绘稿来讲解下原理。
注:由于在我的手稿中以平面坐标系来看,将$\frac{1}{4}$圆形放在第二象限,所以x横坐标取值就为负数。(当然无论实际设定在哪个象限,计算出来的结果都是一样的。)
为了更好的呈现这个落点的过程,我做成了动态可视化,就像下面这样。
可以从上面的动态模拟随机落点看到,使用计算机不停的模拟随机落点,由于圆的半径是1,所以小于等于1单位的随机点都标记为红色,其余则是绿色。最终红色部分很明显现形成了之前提到的$\frac{1}{4}$个单位圆部分。
代码部分如下:
import numpy as np
import pandas as pd
import math
import random
import matplotlib.pyplot as plt
random.seed(123)
# 总随机数量
total = 50000
# 区域内的数量
in_count = 0
x, y = [], []
# 红点x,y坐标
x_red, y_red = [], []
# 绿点x,y坐标
x_green, y_green = [], []
for i in range(total):
# random()方法默认为0,1边界。另外,由于演示的可视化为第二象限,所以x坐标取值为负数。
x.append(-random.random())
y.append(random.random())
# 也可以直接计算(x**2 + y**2)**0.5,因为开根号就等于0.5次方
distance = math.sqrt((x[i]**2 + y[i]**2))
if distance <= 1:
in_count += 1
x_red.append(x[i])
y_red.append(y[i])
else:
x_green.append(x[i])
y_green.append(y[i])
fig, ax = plt.subplots(figsize=(6, 6))
plt.scatter(x_red, y_red, c='red', edgecolor='black')
plt.scatter(x_green, y_green, c='green', edgecolor='black')
print(f'π: {4*in_count/total}')
上面的代码会输出最终的可视化图形,并且输出π值为: 3.13232,可见这个数值并非很精准,因为只设定了最大随机落点50000次,如果随着数量的增大就越会逼近最优解。
最后,要说的是蒙特卡洛方法并没有什么高深的理论支撑,如果一定要说有理论也就只有概率论或统计学中的大数定律了。蒙特卡洛的基本原理简单描述是先大量模拟,然后计算一个事件发生的次数,再通过这个发生次数除以总模拟次数,得到想要的结果。蒙特卡洛方法当然也可以运用在很多领域,如金融,工程,物理,生物医学等等。