python编程通过蒙特卡洛法计算定积分详解

想当初，考研的时候要是知道有这么个好东西，计算定积分。。。开玩笑，那时候计算定积分根本没有这么简单的。但这确实给我打开了一种思路，用编程语言去解决更多更复杂的数学问题。下面进入正题。

如上图所示，计算区间[a b]上f(x)的积分即求曲线与X轴围成红色区域的面积。下面使用蒙特卡洛法计算区间[2 3]上的定积分：∫(x2+4*x*sin(x))dx

# -*- coding: utf-8 -*-import numpy as npimport matplotlib.pyplot as pltdef f(x):  return x**2 + 4*x*np.sin(x) def intf(x):   return x**3/3.0+4.0*np.sin(x) - 4.0*x*np.cos(x)a = 2;  b = 3; # use N draws N= 10000X = np.random.uniform(low=a, high=b, size=N) # N values uniformly drawn from a to b Y =f(X)  # CALCULATE THE f(x) # 蒙特卡洛法计算定积分：面积=宽度*平均高度Imc= (b-a) * np.sum(Y)/ N;exactval=intf(b)-intf(a)print "Monte Carlo estimation=",Imc, "Exact number=", intf(b)-intf(a)# --How does the accuracy depends on the number of points(samples)? Lets try the same 1-D integral # The Monte Carlo methods yield approximate answers whose accuracy depends on the number of draws.Imc=np.zeros(1000)Na = np.linspace(0,1000,1000)exactval= intf(b)-intf(a)for N in np.arange(0,1000):  X = np.random.uniform(low=a, high=b, size=N) # N values uniformly drawn from a to b   Y =f(X)  # CALCULATE THE f(x)   Imc[N]= (b-a) * np.sum(Y)/ N;   plt.plot(Na[10:],np.sqrt((Imc[10:]-exactval)**2), alpha=0.7)plt.plot(Na[10:], 1/np.sqrt(Na[10:]), 'r')plt.xlabel("N")plt.ylabel("sqrt((Imc-ExactValue)$^2$)")plt.show()

>>>

Monte Carlo estimation= 11.8181144118 Exact number= 11.8113589251

从上图可以看出，随着采样点数的增加，计算误差逐渐减小。想要提高模拟结果的精确度有两个途径：其一是增加试验次数N；其二是降低方差σ2. 增加试验次数势必使解题所用计算机的总时间增加，要想以此来达到提高精度之目的显然是不合适的。下面来介绍重要抽样法来减小方差，提高积分计算的精度。

重要性抽样法的特点在于，它不是从给定的过程的概率分布抽样，而是从修改的概率分布抽样，使对模拟结果有重要作用的事件更多出现，从而提高抽样效率，减少花费在对模拟结果无关紧要的事件上的计算时间。比如在区间[a b]上求g(x)的积分，若采用均匀抽样，在函数值g(x)比较小的区间内产生的抽样点跟函数值较大处区间内产生的抽样点的数目接近，显然抽样效率不高，可以将抽样概率密度函数改为f(x)，使f(x)与g(x)的形状相近，就可以保证对积分计算贡献较大的抽样值出现的机会大于贡献小的抽样值，即可以将积分运算改写为：