蒙特卡洛估算定积分的Python实现

蒙特卡洛方法计算定积分

主要介绍随即投点法

我们要估算y=x^2在到1上的定积分

首先我们知道，该定积分的计算出来是1/3，也就是0.3333333……

下面我们利用蒙特卡洛的随机投点方法估算上述定积分mtkl以及两者的误差errorate:

# -*- coding: cp936 -*-

import numpy as np

"""

利用蒙特卡洛随机投点方法估算y=x^2的定积分

"""

def f(x):

return x**2

a=sum(np.where(y

mtkl=a/10000.0

s=1/3.0

errorate=abs(s-mtkl)

print "蒙特卡洛的估计值为：",mtkl

print "误差值为:",errorate

运行上述代码，得到结果如下：

蒙特卡洛的估计值为： 0.3316

事实上，随着投点的个数的增加，误差值将会不断降低，下面我们来验证一下：

首先封装代码为可以调用的函数 ，然后再画一下投点的散点图，修改后如下：

# -*- coding: cp936 -*-

import numpy as np

import matplotlib.pyplot as plt

"""

利用蒙特卡洛随机投点方法估算y=x^2的定积分

"""

deff(x):

return x**2

defnum(N):

a=sum(np.where(y

mtkl=a*1.0/N

s=1/3.0

errorate=abs(s-mtkl)

print"蒙特卡洛的估计值为：",mtkl

print"误差值为:",errorate

plt.scatter(x,y)

plt.show()

下面开始调用：

>>> num(10)

蒙特卡洛的估计值为： 0.3

>>> num(100)

蒙特卡洛的估计值为： 0.38

>>> num(1000)

蒙特卡洛的估计值为： 0.339

>>> num(100000000)

这次没图了，N太大了，因为算这个电脑卡爆了，可以脑补一下，密密麻麻的那种，肉眼看上去几乎覆盖了整个画布。

再次实验多个值之后会发现，误差值随着N的增加而减小。

寝室要关灯了，今天先到这儿吧。

参考：https://www.cnblogs.com/hhh5460/p/6713287.html