蒙特卡洛估算定积分的Python实现--平均值法

蒙特卡洛估算定积分的第二种方法

本文首发于个人微信公众号“我将在南极找寻你”

平均值法
前几天，我们利用蒙特卡洛的随机投点法实现了 y=x^2在 0到 1上的定积分的估算（传送门），今天，我们介绍另一种蒙特卡洛估算定积分的方法：平均值法

原理：辛钦大数定律

这里的大数定律，指的是伯努利大数定律，即

绝对值里面的前者为平均观测值，后者为概率（由二点分布的的公式得到），这就是说，当n很大时，后者（概率）可以由前者（平均观测值）近似代替。而后者又可以写成期望的形式，连续型随机变量求期望就是求对应的定积分，我们由此来估算定积分

具体的推导证明不再叙述（因为我也没看太懂），原理有点复杂，但是我们用代码来实现起来还是挺容易的。

直接上代码，还是以y=x^2为例

# -*- coding: cp936 -*-

import time

import numpy as np

def f(x):

         return x**2

def n(N):

         start_time=time.time()

         a=1/3.0

         x=np.random.uniform(0,1,N)#随机生成N个0-1之间的一维点

         c=f(x)#把x的值代入f（x）计算

         s=0

         for i in c:

                  s=s+i

         j=s*1.00/N

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

         print "与真实值之间的误差为：",abs(a-j)

         end_time=time.time()

         clap=end_time-start_time

         result=round(clap,3)

         print "程序运行耗时： ",result,"秒"

下面进行效果测试：

>>> n(1)

蒙特卡洛估计值为：  0.878686797506

与真实值之间的误差为： 0.545353464173

程序运行耗时：  0.079 秒

>>> n(10)

蒙特卡洛估计值为：  0.206862514025

与真实值之间的误差为： 0.126470819309

程序运行耗时：  0.118 秒

>>> n(100)

蒙特卡洛估计值为：  0.344701364962

与真实值之间的误差为： 0.0113680316288

程序运行耗时：  0.073 秒

>>> n(1000)

蒙特卡洛估计值为：  0.329302957955

与真实值之间的误差为： 0.00403037537882

程序运行耗时：  0.095 秒

>>> n(10000)

蒙特卡洛估计值为：  0.334215109

与真实值之间的误差为： 0.000881775667086

程序运行耗时：  0.082 秒

>>> n(100000)

蒙特卡洛估计值为：  0.333468173775

与真实值之间的误差为： 0.000134840441518

程序运行耗时：  0.141 秒

>>> n(1000000)

蒙特卡洛估计值为：  0.333231069323

与真实值之间的误差为： 0.000102264010195

程序运行耗时：  0.78 秒

>>> n(10000000)

蒙特卡洛估计值为：  0.333309418901

与真实值之间的误差为： 2.39144322688e-05

程序运行耗时：  8.005 秒

>>> n(100000000)

蒙特卡洛估计值为：  0.333349619401

与真实值之间的误差为： 1.62860680264e-05

程序运行耗时：  80.876 秒

当N取100000000时，误差已经非常小了，只是计算时间变成了80多秒，有点耗时，但，毕竟，欲速则不达，对吧。

从开学到现在感觉还没有睡醒，头脑昏昏沉沉的，错误难免，还望指出，睡觉，晚安。