蒙特卡罗方法求解定积分

程锦华

【摘要】蒙特卡罗方法，又称随机模拟，是利用随机数进行统计试验以确定随机事件相应的概率与数学期望的方法。本文首先介绍均匀分布和强大数定律的有关内容，然后讨论用蒙特卡罗方法求定积分的理论基础，收敛速度以及在高维空间中的适用性。

【关键词】均匀分布  强大数定律  蒙特卡罗方法

【中图分类号】G633.6 【文献标识码】A 【文章编号】2095-3089（2018）42-0142-02

1.前言

蒙特卡罗（Monte Carlo）是摩纳哥国的世界著名赌城。第二次世界大战期间，美国原子弹“曼哈顿”计划的成员冯·诺依曼和乌拉姆对裂变物质中子的随机扩散进行模拟，并以蒙特卡罗来命名这种方法。传统的经验方法由于不能逼近真实的物理过程，所以很难得到令人满意的结果，而蒙特卡罗方法能够真实地模拟实际物理过程，故在解决实际问题时往往可以得到很圆满的结果。比如18世纪法国科学家蒲丰就提出了著名的用投针试验来估计圆周率的方法。本文也介绍了利用蒙特卡罗方法求解定积分的过程。

2.预备知识

2.1 均匀分布的随机变量

若随机变量的X的概率密度函数为：

P（x）=■，  a

则称X服从区间（a，b）上的均匀分布，记作X～U（a，b）。其分布函数为

F（x）=0，         x

例1：假设我们从（0，1）区间上取点，取得的点在数轴上的坐标记为随机变量X，则X服从均匀分布。那么X落在[■，■]内的概率是多少？

解：记随机变量X的概率密度函数为p（x），则

P（■≤x≤■）=■p（x）dx=■1dx=■

实际上（0，1）上的均匀分布和高中所学的几何概型息息相关，因为点落在[■，■]内的概率就是此区间的长度除以（0，1）区间的长度。

2.2 强大数定律

定义2.1 假设我们有概率空间（Ω，F，P），定义在其上的随机变量序列X1，X2，…，Xn…及X，满足P■Xn=X=1，则称序列{Xn}几乎处处收敛到X，记为Xn→Xa.s.（n→+∞）。

定理2.1（强大数定律）令X1，X2，…是两两独立同分布的随机变量，且期望存在，即EX1<∞。令EX1=μ，Sn=X1+X2+…+Xn。则当n→+∞时，■几乎处处收敛到μ。

下面，我们通过一个实验来解释强大数定律的直观含义。现在我们要抛一枚质地均匀的硬币n次，若出现正面向上则记为1，否则记为0。令Xi（i=1，2，…，n）表示第i次实验的结果。假设我们进行n次重复的伯努利试验，即可获得一组样本X1，X2，…，Xn。但在抽取样本前无法预知它们的数值，因此样本依然是随机的，具有如下的概率分布列：

表格1：随机变量X1的概率分布列

于是，我们可以证明X1的数学期望存在，且EX1=1/2，于是由强大数定律的结果可知■，即正面向上的次數所占的比例在n→+∞时的极限为1/2。

现在用C++语言进行模拟试验，得到抛硬币的结果如下，

表格2：抛硬币实验的结果

结果表明我们的抛硬币次数越多，正面向上的频率就会越接近1/2，和强大数定律的理论结果相符。

3.蒙特卡罗方法在定积分中的应用

有很多实际问题都需要计算定积分，比如在物理中讨论一些规则物体的质心、转动惯量的计算问题等。根据微积分基本定理，对于任意一个在区间[a，b]上可积的函数f（x），如果能够找到它的原函数F（x），则我们可以通过如下的牛顿-莱布尼茨公式求解定积分：

■f（x）dx=F（b）-F（a）

然而在实际使用这种求定积分方法的时候，往往会遇到很多困难，因为大量的函数，例如■，sinx2等找不到可以用初等函数表示的原函数;因此我们有必要研究求解定积分的其他方法。

3.1求解定积分

令f是[0，1]上的一个连续函数，我们想要计算■f（x）dx，传统的方法是将[0，1]区间划分为一些小区间，然后用梯形公式或者矩形公式进行近似求解。下面我们就介绍一种基于强大数定律的概率途径。令{X1，X2，…，Xn…}表示一列独立同分布服从[0，1]上均匀分布的随机变量，则由f的有界性可知EF（X1）<

+∞，于是应用强大数定律，可知■■■f（Xi）=Ef（X1）a.s.

而由数学期望的定义可知■f（x）dx=Ef（X1）

于是我们可以用来In（f）：=■■■■f（Xi）近似积分I（f）：=■f（x）dx

3.2收敛速度

首先，我们由期望的性质可得EIn（f）=E■■■■f（Xi）=■■■■Ef（Xi）=I（f）。其次，我们考虑这种估计式的均方误差，

E（In（f）-I（f））■=E■■f（Xi）-I（f）■

=■■E（f（Xi）-I（f））（f（Xj）-I（f））

=■E（f（X1）-I（f））■=■Var（f（X1））

于是我们可以得到In（f）-I（f）～■■

3.3蒙特卡罗方法的优点

（下转第145页）

（上接第142页）

当我们处理高维空间情形的时候，求解定积分的复杂度上升。如果我们采用划分小区间的方式，比如每个维度都分成n块，然后分块运用梯形公式或者矩形公式，那么在d维空间中就要有n■个小块，运算复杂度呈指数增长。而当我们采用蒙特卡罗方法求解定积分的时候，根据我们3.2节中的推导，我们发现收敛速度仅仅与模拟随机点的个数有关系，而与空间维度无关，于是在高维的情形下，蒙特卡罗方法几乎是唯一有用的途径。

3.4求定积分的变式

根据3.1中的结果，我们已经可以用蒙特卡罗方法求解在区间[0，1]上连续函数的积分，但是注意函数的连续性并不是必要的，我们只要求其可积性以及随机变量f（X1）的期望存在即可，其中X1服从[0，1]上的均匀分布。现在我们运用变量替换的方法推导其他形式积分的蒙特卡罗方法。

例2：假设g是定义在[a，b]上的连续函数，计算■g（x）dx。

解：作变量替换，令y=■，则dy=■■。我们有■g（x）dx=■g（a+（b-a）y）（b-a）dy=■h（y）dx，其中h（y）=g（a+（b-a）y）（b-a），那么新函数h是[0，1]上的连续函数，于是可以利用3.1中的结果进行近似。

例3：假设g是定义在[0，+∞）上的连续函数，计算■g（x）dx。

解：作变量替换y=■，则dy=-■=-y■dx■， 我们有■g（x）dx=■h（y）dx， 其中h（y）=■。

4.总结

本文首先介绍了均匀分布和强大数定律，并给出了相应的直观解释。然后在强大数定律的基础上推导出了蒙特卡罗方法，并研究了它的收敛速度和在高维空间中的适用性。

参考文献：

[1]茆诗松. 概率论与数理统计简明教程[M].高等教育出版社，2012.

[2]Brian W K， Dennis M R. C程序设计语言[M].机械工业出版社，2004.