Skip to main content

4.2 蒙特卡洛积分

设想我们要计算一个一维函数的积分，例如∫abf(x)dx{\rm \int}_a^{b}f(x){\rm d}x，数值分析方法通常使用一些近似方法来计算积分，最简单的求积分的方法便是将被积函数沿作用域上划分成多个区域，然后计算这些区域面积的和。如果这些区域的划分在作用域上是均匀的，如图（1）所示，则对积分II的近似可以写成：

I=∫abf(x)dx≈∑i=1Nf(xi)(b−a)N \begin{aligned} I&={\rm \int}_{a}^{b}f(x){\rm d}x \\ &\approx \sum_{i=1}^{N}f(x_i) \cfrac{(b-a)}{N} \end{aligned}

（式1）

这里f(xi)f(x_i)是每个区间的中点位置的函数值，NN越大，则该近似越逼近积分II的真实结果。

图（1）：确定性的一维积分

当我们尝试将这种确定性的积分方法推广到d−d-维积分的时候，它要求将原始函数分为NdN^d个采样空间，显然这种方法不适合用于多维积分的计算。当然也存在一些其它的计算多维积分的方法, 但是其中在计算机图形学领域用的最多的还是蒙特卡洛方法，蒙特卡洛方法有很多种不同的用途，其中最重要的一种就是用于计算高维积分。

在推导具体的蒙特卡洛估计之前，我们首先需要了解一下蒙特卡洛方法和传统数值分析方法的区别。式（1）和式 [ref eq:mc-large-numbers]虽然看起来具有比较相似的思路：即它们都是计算NN个“采样点”的平均值，但是前者是一种传统的数值近似方法，它必须在每个维度上使用同样的近似结构（即对每个维度执行空间划分），因此它的采样的数量是和维度相关的；而大数定律的采样点来自于满足某个分布的一个随机数，这个随机数直接对一维或多维空间进行采样，这个采样的数量与函数的维度并无直接关系。

大数定律用于对期望表示的积分形式进行估计，即对积分∫−∞∞xf(x)dx{\rm \int}_{-\infty}^{\infty}xf(x){\rm d}x进行估计，但是我们通常需要求得一个任意函数的积分，假设函数g(x)g(x)的定义域为x∈Sx\in S(SS可以是一个多维空间，此时xx表示一个向量)，我们希望计算以下积分：

I=∫x∈Sg(x)dμ I={\rm \int}_{x\in S}g(x){\rm d}\mu

（式2）

根据前面的描述，给定任意一个实数函数ff以及一个随机变量x∼p(x)x\sim p(x)，我们可以使用以下的一个和的公式来近似f(x)f(x)的期望：

E[f(x)]=∫x∈Sf(x)p(x)dμ≈1N∑i=1Nf(xi) E[f(x)]={\rm \int}_{x\in S}f(x)p(x){\rm d}\mu\approx \cfrac{1}{N}\sum_{i=1}^{N}f(x_i)

（式3）

这个结果就是使用蒙特卡洛方法求积分的核心思想，因为期望表示的积分可以使用多个随机数采样值的和的形式来近似。对式（3）稍作变换， 使用g=fpg=fp 代替被积函数，则上式变为：

∫x∈Sg(x)dμ≈1N∑i=1Ng(xi)p(xi) {\rm \int}_{x\in S}g(x){\rm d}\mu\approx \cfrac{1}{N}\sum_{i=1}^{N} \cfrac{g(x_i)}{p(x_i)}

（式4）

\noindent 对于该公式，pp必须为正，我们定义II的估计（estimator）为INI_N，即：

IN=1N∑i=1Ng(xi)p(xi) I_N= \cfrac{1}{N}\sum_{i=1}^{N} \cfrac{g(x_i)}{p(x_i)}

（式5）

该估计的期望值为：

E[IN]=E[1N∑i=1Ng(xi)p(xi)]=1N∑i=1NE[g(xi)p(xi)]=1NN∫g(x)p(x)p(x)dx=∫g(x)dx=I \begin{aligned} E[I_N]&=E[ \cfrac{1}{N}\sum_{i=1}^{N} \cfrac{g(x_i)}{p(x_i)}]\\ &= \cfrac{1}{N}\sum_{i=1}^{N}E[ \cfrac{g(x_i)}{p(x_i)}]\\ &= \cfrac{1}{N}N{\rm \int} \cfrac{g(x)}{p(x)}p(x){\rm d}x\\ &={\rm \int} g(x){\rm d}x\\ &=I \end{aligned}

（式6）

由此可以看出该估计可以用于用作II的近似，该估计的方差为：

σ2=1N∫(g(x)p(x)−I)2p(x)dx \sigma^2= \cfrac{1}{N}{\rm \int}( \cfrac{g(x)}{p(x)}-I)^2p(x){\rm d}x

（式7）

由此可以看出，随着NN的增加，其方差随着NN线性降低，由于估计的误差正比于其标准差σ\sigma，因此估计的误差随着N\sqrt{N}线性降低。这就是一般蒙特卡洛方法的特点，实际上蒙特卡洛方法最大的问题就是估计逼近正确结果的速度非常慢，例如需要使用四倍多的采样点才可以仅仅减少一半的误差。

理论上，p(x)p(x)的选择可以是任意的，这也是蒙特卡洛方法的一大优点，因为通常很难生成与被积函数具有一致分布的随机数。但是从式（5）也可以看出，通过使g(xi)g(x_i)和p(xi)p(x_i)的比值尽可能地小也可以减少估计误差，在实践上这通常使p(x)p(x)的分布接近于g(x)g(x)，这称为本章后面会讨论的重要性采样的基础。

终上所述，蒙特卡洛方法是一种非常简单而强大，可以近似计算任意函数积分的数值方法。蒙特卡洛方法可以归纳为以下步骤：

  • 首先对一个满足某种概率分布的随机数进行采样。
  • 使用该采样值计算gip(xi) \cfrac{g_i}{p(x_i)}的值，这称为该样本的贡献值。
  • 最后对所有采样点计算的结果求平均值。

在这些步骤中，最困难的部分涉及怎样对一个任意分布函数进行采样，我们将在第（4.3）节讨论一些常用的采样方法；其次，蒙特卡洛方法的另一个重要内容是涉及怎样减少方差，第（4.5）节将会讨论这方面的内容。