蒙特卡洛方法的名字来源于摩纳哥的一个城市蒙特卡洛,该城市以赌博业闻名,而蒙特卡洛方法正是以概率为基础的方法。
蒙特卡洛方法(Monte Carlo method),也称统计模拟方法,是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明,而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。
当所求解问题是某种随机事件出现的概率,或者是某个随机变量的期望值时,通过某种“实验”的方法,以这种事件出现的频率估计这一随机事件的概率,或者得到这个随机变量的某些数字特征,并将其作为问题的解。有一个例子可以使你比较直观地了解蒙特卡洛方法:
假设我们要计算一个不规则图形的面积,那么图形的不规则程度和分析性计算(比如,积分)的复杂程度是成正比的。蒙特卡洛方法是怎么计算的呢?假想你有一袋豆子,把豆子均匀地朝这个图形上撒,然后数这个图形之中有多少颗豆子,这个豆子的数目就是图形的面积。当你的豆子越小,撒的越多的时候,结果就越精确。在这里我们要假定豆子都在一个平面上,相互之间没有重叠。
根据上面的思想,我们来求算圆周率。假设有半径为 R 的圆,其面积为 πR^2,则其外接正方形边长为 2R,面积为 4R^2 。假设向这个由正方形和圆组成的图形中随机扔豆子,扔的次数为 N ,落入圆中的次数为 n ,则根据面积比例,有:
πR^2 / 4R ^ 2 ≈ n / N
注意这里是约等于。根据概率论,当 N 越大时,则上面左右两边的差值会越小,这样,圆周率π,可以使用下式来逼近:
π ≈ 4 * n / N
我们使用 Python 语言将上面算法实现:
from __future__ import division
import random
import time
num = 1
for j in range(2, 100000):
startT = time.clock()
# 落入圆内计数
counter = 0
# 往正方形中扔了 10 的 j 次方次
for i in range(10 ** j):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
# 落入了圆内
if x**2 + y**2 < 1:
counter = counter + 1
endT = time.clock()
print ( num)
print ( 'pi:{0}'.format(4 * (counter / 10 ** j)))
print ( 'Time:{0}'.format(endT - startT))
print ('')
num += 1
尽管可以使用上面代码来计算π,但是实际的计算速度是很慢的(当然还可以使用扩展库或 ctypes 进行优化)。而且,计算机产生的随机数只能精确到某位数,并不能产生任何实数(例如无理数等等),即使取10的9次方个随机点时,其结果也仅在前4位与圆周率吻合。