
3.1 蒙特卡洛法求定积分
蒙特卡洛(Monte Carlo)法是一类随机算法的统称。它是20世纪40年代中期由于科学技术的发展,尤其是电子计算机的发明,而被提出并发扬光大的一种以概率统计理论为基础的数值计算方法。它的核心思想就是使用随机数(或更准确地说是伪随机数)来解决一些复杂的计算问题。现今,蒙特卡洛法已经在诸多领域展现出了超强的能力。本节,我们将通过蒙特卡洛法最为常见的一种应用——求解定积分,来演示这类算法的核心思想。
3.1.1 无意识统计学家法则
作为一个预备知识,先来介绍一下无意识统计学家法则(Law of the Unconscious Statistician,LOTUS)。在概率论与统计学中,如果知道随机变量X的概率分布,但是并不显式地知道函数g(X)的分布,那么LOTUS就是一个可以用来计算关于随机变量X的函数g(X)之期望的定理。该法则的具体形式依赖于随机变量X之概率分布的描述形式。
如果随机变量X的分布是离散的,而且我们知道它的PMF是fX,但不知道fg(X),那么g(X)的期望是

其中和式是在取遍X的所有可能之值x后求得。
如果随机变量X的分布是连续的,而且我们知道它的PDF是fX,但不知道fg(X),那么g(X)的期望是

简而言之,已知随机变量X的概率分布,但不知道g(X)的分布,此时用LOTUS公式能计算出函数g(X)的数学期望。其实就是在计算期望时,用已知的X的PDF(或PMF)代替未知的g(X)的PDF(或PMF)。
3.1.2 投点法

图3-1 投点法求定积分
投点法是讲解蒙特卡洛法基本思想的一个最基础也最直观的实例。这个方法也常常被用来求圆周率π。现在我们用它来求函数的定积分。如图3-1所示,有一个函数f(x),若要求它从a到b的定积分,其实就是求曲线下方的面积。
可以用一个比较容易算得面积的矩型罩在函数的积分区间上(假设其面积为Area)。然后随机地向这个矩形框里面投点,其中落在函数f(x)下方的点为菱形,其他点为三角形。然后统计菱形点的数量占所有点(菱形+三角形)数量的比例为r,那么就可以据此估算出函数f(x)从a到b的定积分为Area×r。
注意由蒙特卡洛法得出的值并不是一个精确值,而是一个近似值。而且当投点的数量越来越大时,这个近似值也越接近真实值。
3.1.3 期望法
下面来重点介绍利用蒙特卡洛法求定积分的第二种方法——期望法,有时也称为平均值法。
任取一组相互独立、同分布的随机变量{Xi},Xi在[a,b]上服从分布律fX,也就是说fX是随机变量X的PDF(或PMF)。令,则g∗(Xi)也是一组独立同分布的随机变量,而且因为g∗(x)是关于x的函数,所以根据LOTUS可得

由强大数定理

若选

则依概率1收敛到I。平均值法就用
作为I的近似值。
假设要计算的积分有如下形式

其中,被积函数g(x)在区间[a,b]上可积。任意选择一个有简便办法可以进行抽样的概率密度函数fX(x),使其满足下列条件:
(1)当g(x)≠0时,fX(x)≠0,a≤x≤b;
(2)。
如果记

那么原积分式可以写成

因而求积分的步骤是:
(1)产生服从分布律fX的随机变量Xi,i=1,2,…,N;
(2)计算均值

并用它作为I的近似值,即。
如果a,b为有限值,那么fX可取作为均匀分布

此时原来的积分式变为

因而求积分的步骤是:
(1)产生[a,b]上的均匀分布随机变量Xi,i=1,2,…,N;
(2)计算均值

并用它作为I的近似值,即。
最后来看一下平均值法的直观解释。注意积分的几何意义就是[a,b]区间曲线下方的面积,如图3-2所示。
当在[a,b]随机取一点x时,它对应的函数值就是f(x),然后便可以用f(x)·(b-a)来粗略估计曲线下方的面积(也就是积分),如图3-3所示,当然这种估计(或近似)是非常粗略的。

图3-2 积分的几何意义

图3-3 对积分值进行粗略估计
于是我们想到在[a,b]随机取一系列点xi时(xi满足均匀分布),然后把估算出来的面积取平均来作为积分估计的一个更好的近似值,如图3-4所示。可以想象,如果这样的采样点越来越多,那么对于这个积分的估计也就越来越接近。

图3-4 对积分值进行估计
按照上面这个思路,得到积分公式为

其中,就是均匀分布的PMF。这跟之前推导出来的蒙特卡洛积分公式是一致的。