引入——圆周率的计算
如何计算?由
知可以由圆面积求得
那么如何求圆面积?有一种近似的方法:往单位圆的外切正方形里投掷点,落在圆内的点数与总点数的比值近似等于圆面积与正方形面积的比值,也即,由此可以近似计算圆周率。
其实这里就用到了蒙塔卡罗的方法
蒙塔卡罗法
又称统计模拟法,是一种随机模拟方法,以概率和统计理论方法为基础的一种计算方法,是使用随机数来解决很多计算问题的方法。将所求解问题同一定的概率模型相联系,用计算机模拟,获得近似解。
原理:大数定理。
蒙特卡洛并不是一种算法或者模型,而是一种思想,因此没有统一的模板。
蒙特卡洛法运用实例的Python代码
圆周率计算
from numpy import random
def cal_pai_mc(n=1000000):
r = 1.0
a, b = (0.0, 0.0)
x_neg, x_pos = a - r, a + r
y_neg, y_pos = b - r, b + r
m = 0
for i in range(0, n + 1):
x = random.uniform(x_neg, x_pos)
y = random.uniform(y_neg, y_pos)
if x ** 2 + y ** 2 <= 1.0:
m += 1
return (m / float(n)) * 4
print(cal_pai_mc())
计算定积分
若要求函数f(x)从a到b的定积分,我们可以用一个比较容易算得面积的矩型包围在函数的积分区间上(假设其面积为Area),定积分值其实就是求曲线下方的面积。随机地向这个矩形框里面投点,统计落在函数f(x)下方的点数量占所有点数量的比例为P,那么就可以据此估算出函数f(x)从a到b的定积分为Area×P。此处我们将a和b设为0和1,函数
from numpy import random
def cal_integral_mc(n=1000000):
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
m = 0
for i in range(0, n + 1):
x = random.uniform(x_min, x_max)
y = random.uniform(y_min, y_max)
# x*x > y 表示该点位于曲线的下面。
if x * x > y:
m += 1
# 所求的积分值即为曲线下方的面积与正方形面积的比
return m / float(n)
print(cal_integral_mc(1000000))
计算函数极值
此处在区间[-2,2]上随机生成一个数,求出其对应的y,找出其中最大值认为是函数在[-2,2]上的极大值。
from numpy import random
import numpy as np
def cal_extremum_mc(n=1000000):
y_max = 0.0
x_min, x_max = -2.0, 2.0
y = lambda x: 200 * np.sin(x) * np.exp(-0.05 * x) # 匿名函数
for i in range(0, n + 1):
x0 = random.uniform(x_min, x_max)
if y(x0) > y_max:
y_max = y(x0)
x_max = x0
return y_max, x_max
print(cal_extremum_mc(1000000))