蒙特卡洛方法原理

〖摘要〗贝叶斯统计需要在贝叶斯定理基础上，通过参数先验和数据似然对参数的后验概率分布作出推断。从推断精度上区分，贝叶斯推断方法大致包含精确推断和近似推断两大类，其中精确推断常见有变量消除法（Variable Elimination， VE）和信念传播法（Belief Propagation, BP）；而近似推断方法主要是蒙特卡洛法（Mente Carlo, MC）和变分近似推断法（Variational Inference,VI）。蒙特卡洛方法是一种以概率统计理论为指导的重要数值计算方法。它使用随机数来解决随机变量（或随机函数）期望值积分求解、仿真模拟等非常棘手的计算问题，特别适用于没有明确解析形式的复杂概率分布。蒙特卡洛方法中最为核心的部分是如何在给定一个复杂分布时，按照概率随机地、高效地获得样本，即采样方法问题。

〖原文〗蒙特卡洛方法数学基础、蒙特卡洛方法实践

1 引子

蒙特卡洛方法是一类通过随机采样来求解问题的算法的统称，通常要求解的问题是某随机事件的概率或某随机变量的期望。通过随机采样方法，以随机事件出现的频率估计其概率，并将其作为问题的解。

蒙特卡洛的基本做法：做大量重复试验，通过统计频率来近似概率（即概率密度越大的候选值采样频率越高），从而得到问题的解。

举一个简单例子来说明 “统计频率” 和 “概率” 之间的关系：

如 图 1 所示，一个矩形内有一个不规则图案，现在要求解不规则图形的面积。显然，矩形的面积可以简单计算为 $A=ab \times ac$ ，点位于不规则形状内的真实概率为 $p=A_{shape}/A$ 。现在重复往矩形范围内（均匀）随机地投射点，样本点有一定概率会落在不规则图形内。若假设做 $n$ 次试验，落到不规则图形内的次数为 $k$ ，频率为 $\frac{k}{n}$ 。 若样本数量足够大，则有：$p=\frac{A_{shape}}{A} \approx \frac{k}{n}$。根据大数定理：随着样本数增加，频率 $k/n$ 会收敛于概率 $p$ 。

如果不规则图形的面积很难求取，而恰好我们能够通过某种方式获得投射点的统计频率数据，则可以利用上述原理，通过矩形面积和统计频率来估计不规则图形的面积为 $A_{shape}=A\frac{k}{n}$ 。假设矩形面积为 $1$，投射了 $1000$ 次，有 $200$ 个点位于不规则形状内，则可以推算出不规则图形的面积为 $0.2$，并且投射的次数越大，计算出来的面积越精确。

该例子说明了蒙特卡洛方法的基本思路，虽然看起来很简单，但其实有着严格的数学基础作为依托。蒙特卡洛方法的特点是：求解的近似结果有误差，但重复试验次数越多，误差所引起的方差就会越低。

图 1：利用蒙特卡洛法估计不规则图形面积

2 采样的动机

2.1 现实案例

除了上面的例子外，生活中还有许多事情都很难准确估计，特别是当其涉及数量非常大的时候。例如：

• 计算一公斤罐子里面的糖豆数量。你需要一个接一个地数，但这需要很长时间。

• 计算一个国家成年人口的平均身高，需要测量每个人的身高，将这些数字相加，然后除以测量的总人数，这项任务可能需要很长时间。

取而代之的方法时：你可以从这些人口中抽取样本，并计算出其平均身高。这种采样方法可能无法给出非常准确的平均身高，但其结果通常是真实数字的很好近似。这里，我们实际上是在精度和速度之间做了平衡。不过有趣的是，近似值和精确值有时可能完全相同，尽管完全是偶然的。因此，我们可能想问一个问题，那就是两者之间到底有多大区别？事实上，随着样本大小的增加，近似值会收敛到精确值。换言之，随着样本量增加，近似结果和真实结果之间的误差越来越小。

上述案例就是一种典型的蒙特卡洛方法，体现了蒙特卡洛方法的基本理念：

通过从总体中随机抽取样本，将大量人口的平均身高计算问题，转变成了少量采样人口的平均身高计算问题，其效果是用精度换取了时间效率。

人口身高这个案例与不规则图形面积求取案例有一个共同特点需要注意：案例中都隐含了一个特殊条件，那就是 由于机会均等，所以构成样本的元素都是以等概率（即均匀分布）随机选择的。但事实上，现实当中出现不平等（概率不同）条件的情况更为广泛（有些概率分布甚至存在多峰值、不规则、悬崖等非常复杂的情况，在多维随机变量场景中复杂度更高），并且是我们需要面对的主要情况。

2.2 数学解释

可以从数学上对上述案例做泛化解释。

首先，可以将身高视为一个随机数（当然，根据问题不同，它也可以是任何我们感兴趣的东西）。问题表述为：当对一个总体进行采样时，通过随机选取该总体中的个体并测量其身高，获得近似平均身高真实值的结果。注意：这里每个测量值都是一个随机数（不确定会抽到谁，因此身高数值不确定），而所有随机数之和是另一个随机数（抽取的元素不确定，因此平均身高数值也不确定）。

对于数学家来说，一个群体的身高会被称为随机变量，因为构成该群体的人的身高是随机变化的。我们通常用大写字母表示随机变量，例如：$X,Y$ 等。

在统计学中，组成总人口的个体也是随机的，通常用小写字母表示。例如，$x,y$ 等。 $x_2$ 表示人口中第二个人的身高（也是随机的）。所有这些 $x$ 构成的集合，被视为随机变量 $X$ 的一个可能结果。如果称 $X$ 为这个随机变量（人口高度），则可以用下面的伪数学公式来表达从样本中近似计算成人平均身高的概念：

$$\begin{align*} \text { Approximation }(\text { Average }(X))=\frac{1}{N} \sum_{n=1}^{N} x_{n}\\ \operatorname{Average}(X)=\frac{1}{M}\left(x_{1}+x_{2}+\ldots x_{M}\right) \end{align*}$$

可以理解为：随机变量 $X$ 均值（给定国家成年人口的身高）的近似值等于从该人口中随机选择的 $N$ 个成年人身高之和除以数字 $N$ （样本大小）。本质上，这是一种 蒙特卡罗近似，它通过从所有事物群体中随机选择 $N$ 个事物，并计算其某个属性的均值，用以近似总体上该属性的均值。在统计学中，随机变量 $X$ 的均值被称为期望值，记为 $\mathbb{E}(X)$。第二个公式为该国家成年人口平均高度的真实值计算公式，注意其中虽然每个 $x$ 都是随机的，但最终的平均值应是唯一的（为避免与样本数 $N$ 混淆，使用 $M$ 表示人口总数）。

注：

蒙特卡罗近似是一种常见的蒙特卡洛方法，它利用样本逼近随机变量的期望值，可以用以下公式在数学上进行定义：

$$\mathbb{E}(X) \approx \sum_{i=1}^N p_i x_i$$

其中： $p_i$ 为第 $i$ 个样本抽取值为 $x_i$ 的概率。当所有值的抽取概率相等时，上式可简化为：

$$\mathbb{E}(X) \approx \frac{1}{N} \sum_{n=1}^{N} x_{n}$$

数学符号 $\approx$ 意味着其右边的公式只给出了随机变量 $X$ 的期望 $\mathbb{E}(X)$ 真实值的一个“近似值”。而实际上，它只不过是若干随机数值（所有$x_n$ 的值）的平均值。

需要再次提醒：人口身高案例的情况比较特殊，由于人人平等（概率 $p$ 处处相等 ），因此期望值等于均值；而现实中存在大量机会不平等的情况（其中有些情况 $p$ 甚至复杂到没法用公式来描述），其期望值 $\mathbb{E}(X)$ 的计算会非常棘手。尤其是随着统计维度的增加（如：统计口径从原来的国民，扩展到不同性别、年龄、民族的国民），其期望值的计算更为复杂耗时。此时，采用蒙特卡洛近似可能是获得高精度期望值的唯一办法，因为随着抽取样本数量的增加，近似值会逐步收敛到真实值。

2.3 蒙特卡罗方法的主要应用

综上所述，蒙特卡罗方法是一种通过随机采样（或通过模拟随机变量）解决数学问题的数值方。所有蒙特卡洛方法都使用随机抽取的样本来计算给定问题的解。这些问题通常分为两大类：

• 模拟：当所求解的问题本身具有内在的随机性时，借助计算机的运算能力可以直接模拟这种随机的过程。例如在核物理研究中，分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，并将其作为反应堆设计的依据。

• 积分：数学中常用的计算，虽然计算积分的方法很多，但是随着分布复杂性和维度的增加，传统积分计算的代价会相当高。而蒙特卡洛积分可以通过随机采样方法实现积分计算，虽然没有达到积分真实解的最大收敛速度，但能够以 “更便宜的” 计算成本获得相当接近的结果。通常对于此类应用，所求解的问题可以被转化为某种随机分布的统计特征（或统计量），比如随机事件出现的概率、随机变量的期望值等。通过随机采样方法，用随机事件出现的频率（及样本出现的频率）来估计其真实概率，或者以样本的数字特征来估算随机变量的真实数字特征，并将其作为问题的解。

3 统计学说明

在统计学中，总体的某些特征可以通过该总体中所有元素来计算时，通常称其结果值为总体的参数，根据定义可知，参数为确定值。例如：总体均值是一个用于定义数量平均值的总体参数。另一方面，当使用样本来估计总体参数时，则称由样本得到的值是统计量。参数通常用希腊字母表示（例如 $\sigma$），而统计量通常使用罗马字母表示（例如：$\mathrm{S}$）。注意，由于可以从总体中抽取多个样本，因此统计量的值可能会因样本不同而变化。

计算均值或方差的方式取决于使用总体所有元素，还是只使用样本。下表有助于比较两者之间的差异：

总体（参数）	样本（统计量）
$\bar X \approx \mu$
参数是一个值，通常是未知的，用于表示特定的总体特征（例如：总体均值是用来表示一个数量平均值的参数）	统计量是从样本中计算出来的量。它用于提供有关总体中未知值的信息（例如：样本中数据的平均值可以提供有关总体均值的信息）
$N$ 指总体中所有元素的数量（例如：一个国家的总人口）	$n$ 是样本的数量（例如：我们抛掷硬币的次数）
均值： $\mu=\frac{\sum_{i=1}^{N} x_{i}}{N}$	$\bar{X}=\frac{\sum_{i=1}^{n} x_{i}}{n}$
方差： $\sigma^2 = \dfrac{\sum_{i=1}^N {(x_i - \mu)^2}}{N} = \dfrac{\sum_{i=1}^N x_i}{N} - \mu^2$	$S_n^2 = \dfrac{\sum_{i=1}^n (x_i - \bar X)^2 }{n} = \dfrac{\sum_{i=1}^n x_i}{n} - \bar X^2$

注意：统计量的方差用符号 $S_n^2$ 表示，以将其与总体方差 $\sigma^2$ 做区分。 $S_n^2$ 的公式很有帮助，因为从编程角度来看，可以在从总体中提取元素后使用该公式来计算方差，而不用等到掌握了所有数据。

下面举个具体的例子：

假设有一个 [0,20] 间的数字构成的总体，每个数字出现的次数在 [1,50] 之间，总体的大小为 $N=572$， 总体的平均值为 $\mu = 8.970280$， 总体的方差为 $\sigma ^2 =35.476387$ 。现在对该总体进行 $1000$ 次采样，每次采样的样本数量 $n$ 为 $2$ 到 $20$ 间的随机数，样本中每个元素均从总体的 $572$ 个元素中随机抽取。

由于样本中的元素是从总体中随机抽取的，因此每次采样我们都更新了用于计算样本均值和样本方差的变量。最后，我们得到了 $1000$ 个样本以及其对应的样本均值和方差（图 2）。样本数量小的采样用红色绘制，样本数量大的采样用蓝色绘制。图中的绿点表示总体均值和方差的真实值点。任何接近该点的样本都可以被视为提供了对总体均值和方差的良好估计，而样本离该点越远，估计越差。

图 2： 样本统计量分布（红色代表该次采样为 2 个样本，蓝色代表该次采样有 20 个样本）。

观察 图 2 可以得出什么结果？

首先，可以沿横坐标 ( $x$ ) 和纵坐标 ( $y$ ) 理解图形。任何接近总体均值线（垂直青色线）的样本都可以被视为提供了对总体均值的良好估计。靠近总体方差线（水平青色线）的任何样本都可以很好地估计总体方差。然而，真正“好”的样本是接近两条线的样本。我们在 图 2 中用黄色形状说明了这一点，显示了绿点周围的一组样本，这些样本可以被认为是“好”样本，因为其统计量与相应参数相当接近。黄色形状内的大多数样本为蓝色，但并非都为蓝色（存在小部分红色或微红色样本）。这表明，样本量越大越有可能得出收敛于总体参数的估计值。但请注意，此处并没有说样本量越大，估计就一定越好，而只是好的可能性更大了。这捕捉到了统计学和近似概念之间的差异：统计学中样本给出的值总是有可能（即使很小）与总体参数完全相同。无论多么不可能，无论样本量有多大（或小），这种机会都存在。但与规模大得多的样本相比，规模非常小的样本更有可能提供偏离的估计。事实上，根据大数定律，当 $n$ 接近无穷大时，样本平均值在概率上收敛并且“几乎”肯定会达到总体的期望值 $μ$。

图 3：随着 n 增加，样本平均值在概率上收敛并且“几乎”肯定会达到期望值 μ（红色 = 16 个样本，蓝色 = 32 个样本）。

我们不能肯定它会收敛，但可以说随着样本量增加，它收敛的概率会越来越高。图 3 显示了随着样本量增加，样本的“质量”会发生什么变化。样本仍然分布在总体均值和方差周围，但样本聚簇范围缩小了，这表明该图中估计的质量确实比 图 1 中的质量要好。换句话说，随着 $n$ 增加，样本的聚簇范围缩小（表明估计质量的整体改进）。也可以说随着 $n$ 的增加，统计量越来越聚集在某个值上。该值是总体均值，即样本均值随着 $n$ 接近无穷大而接近总体均值。

$$\begin{array}{l} \lim_{n \rightarrow \infty} \text{ Pr }(|\bar X - \mu| < \epsilon) = 1 \\ \bar X \rightarrow \mu \end{array}$$

图 4 为总体的分布，该总体不按照任何特定概率分布（尤其不是正态分布）。由于该分布是离散且有限的，因此该总体具有明确定义的均值和方差。现在要做的是从这个总体中抽取一个大小为 $n$ 的样本，计算它的均值和方差，并重复实验 $1000$ 次。样本均值四舍五入到最接近的整数值，以便转换为 $0$ 到 $20$ 之间的任何整数值。在该过程结束时，我们将统计均值为 $0、1、2、...、20$ 的样本计数，图 5 显示了计数结果。值得注意的是，该分布遵循正态分布。此处需非常清楚两者的区别：图 4 是一个卡片数值的总体分布，而 图 5 是一个统计量的分布。另外，图 5 并不是完美的正态分布，其与完美的正态分布（红色曲线）之间存在一些差异。图 4 和 图 5 说明：即使总体分布是任意的，多次采样产生的统计量的分布却收敛于正态分布。

图 4 总体的分布。纵轴表示总体中值为 0 到 20 之间任意值的元素数量。由于随机性，该分布是任意的。

图 5：样本（统计量）的分布。纵轴代表样本均值等于 0 到 20 之间任意值的样本数量。该分布不是任意的，而是服从正态分布。

根据上述对比分析，样本的均值和方差也可以被视为随机变量（因为每个样本都可能具有不同值），因此像任何其他随机变量一样，可以研究样本均值和方差的概率分布（如图 5 所示）。换句话说，我们并不是研究事物的某个属性值（如：特定国家人口的身高）如何分布，而是通过从该事物总体中抽取的样本来估计该属性值的某些统计量（如：平均身高），并看看这些统计量的值在样本中如何分布。在统计学中，样本（或统计量）分布称为采样分布。同总体分布类似，采样分布也可以用概率分布（或概率密度）来表示。它显示了给定总体和样本数量情况下，某些统计量的样本分布情况。

某个统计量的采样分布是将该统计量视为随机变量后的概率分布，该统计量的每次取值均来自于一个大小为 $n$ 的随机样本。换句话说，均值的采样分布是样本均值的分布、方差的采样分布是样本方差的分布。

请记住，统计量是一些随机变量的函数，因此它本身就是一个随机变量。与所有其他随机变量一样，统计量也具有分布，而该分布就是采样分布。根据概率分布的定义，采样分布反映了“统计量可能的值以及取用这些值的可能性”。

采样分布也是一种概率分布，也具有期望、方差等参数，并沿用期望和方差的计算公式。

$$\mu_{\bar X} = E(\bar X) = \sum_{i=1}^n p_i \bar X_i,$$

其中 $p_i$ 是与样本均值 $\bar X_i$ 相关的概率，这些概率根据采样分布而分布。当以相同概率（即均匀分布）抽取样本时，概率为 $\frac{1}{n}$，其中 $n$ 代表样本数量（注意不是单个样本中样本点的数量）。

注意：不能将 $μ_{\bar X}$ 视为与均值完全相同的事物（因为随着样本量接近无穷大，该值没有定义），但它可以被解释或视为均值。这种均值统计量采样分布的预期值，可被视为样本均值（有限数量）的加权平均，其权重是每个可能值的概率。

可以证明，该采样分布的期望值 $μ_{\bar X}$ 等于总体的均值 $μ$ （证明略）。在统计学中，当统计量的期望值（如： $μ_{\bar X}$ ）等于总体参数（如：$\mu$ ）时，统计量本身被称为该参数的无偏估计量。

如果现在增加样本数量 $n$（注：此处并非指单一样本中的样本点的数量，而是指总样本数），那么样本集更有可能对总体参数做出准确估计，并非所有样本都接近参数值，但通常样本数量越大，相应统计量采样分布的均值就越接近总体参数值。新的采样分布与 图 5 一样具有正态分布的形状，但 图 6 比 图 5 更为紧凑。随着样本数量 $n$ 的增加，越来越接近于使用总体所有元素来估计参数值。因此从逻辑上讲，随着样本数量的增加，样本均值这个统计量平均下来更有可能接近总体均值。这个实验显示了类似于在 图 2 和 图 3 中看到的东西，即样本量（此处指单一样本中的样本点数量）越大，越有可能收敛到真实平均值。可以通过观察 $μ_{\bar X}$ （及其标准偏差 $σ_{\bar X}$ ）的采样分布的预期值随样本数量增加而变化的规律轻松证实这一点。图 6 中红色曲线为 $1000$ 个样本，而绿色曲线对应于 $10000$ 个样本。

图 6： 在增加样本数量同时，采样分布变得更接近完美的正态分布 $N(μ,0)$。

有关样本、样本集合、统计量、统计量的均值等概念，很容易造成理解的混乱，图 7 给出了一个形象的解释。

图 7： 统计量、统计量的均值及其关系

4 蒙特卡洛积分

如 第 2 节 所述，蒙特卡洛方法主要应用于解决模拟和积分两种问题，本节重点介绍其中积分问题的解决思路和原理（注：模拟部分不在本系列文章范围内）。

4.1 最简单的定积分

基本的蒙特卡罗估计原理是：想象我们要对一个从 $a$ 到 $b$ 的一维函数 $f(x)$ 进行积分，例如：

$$F = \int_a^b f(x)\;dx$$

函数 $f(x)$ 的积分可解释为计算函数曲线下方的面积。该想法如图 7 所示。现在假设我们选择一个随机值，比如范围 $[a,b]$ 中的 $x$，在 $x$ 处评估函数 $f(x)$ 并将结果乘以 $(b-a)$。图 7 中红色区域显示了该结果：它是一个高度为 $f(x)$、宽度为 $(b-a)$ 的矩形。某种程度上，可以将该矩形区域面积近似作为相应曲线的面积。再如 图 7 下图所示，如果评估 $x1$ 处的函数值，我们会大大低估这个区域的面积。如果我们在 $x2$ 处评估函数，我们就会高估面积。但随着不断在 $a$ 和 $b$ 之间的随机点做评估，并将所有矩形面积相加并求平均值，会发现数字越来越接近积分的实际结果。事实上，随着计算中使用样本数量的增加，平均面积实际上会收敛到真实的积分“面积”。如 图 7 下图所示，在四个不同位置评估了 $x$ 处的函数值，并分别乘以相应的 $(b-a)$ 得到各自面积，然后求面积和与平均值，最终的面积平均值可被认为是实际积分的近似值。

图 8：可以在 x 处评估曲线值 f(x)，结果乘以 (b - a)。这定义了一个矩形，可以将其视为积分的粗略的近似。

可以用以下公式形式化该想法：

$$\langle F^N\rangle = (b-a) \dfrac{1}{N } \sum_{i=0}^{N-1} f(X_i).$$

其中 $N$ 是此近似中使用的样本数量。在数学符号中，$⟨S⟩$ 表示 $S$ 中所有元素的平均值，$⟨F^N⟩$ 指 使用 $N$个样本的得到的 $F$ 近似值。它等价于样本均值符号 $\bar X_n$ 。该公式被称为基本蒙特卡罗估计量。区间 $[a,b]$ 中的随机点 $X_i$ 可以通过 $[0,1]$ 区间的均匀随机生成器结果与 $(b-a)$ 相乘来获得，即：$X_i = a + \xi (b - a)$ ，其中 $ξ$ 均匀分布在 0 和 1 之间。 结果 $X_i$ 的概率密度函数 $PDF$ 是 ${1}/{(b - a)}$。（注：对于离散型变量，由于随机数是等概率产生的，不同取值的概率只需将 1 除以可能的结果总数即可。但本例函数为连续的，所以采用概率密度函数表示，将 1 除以区间 $[a,b]$ ）。

重要的是要注意：

$$Pr(\lim_{ N \to \infty} \langle F^N \rangle = F ) = 1.$$

大数定律告诉我们，随着 $N$ 接近无穷大，蒙特卡罗近似会大概率收敛到正确答案（即概率为 1）。

还要注意 $⟨F^N⟩$ 是一个随机变量，因为它是由随机数的总和组成。现在可以证明 $⟨F^N⟩$ 的期望值等于 $F$：

$$\begin{array}{l} E[\langle F^N \rangle] & = & E \left[ (b-a) \dfrac{1}{N } \sum_{i=0}^{N-1} f(x_i)\right],\\ & = & (b-a)\dfrac{1}{N } E[f(x)],\\ & = &(b-a)\dfrac{1}{N} \sum_{i=0}^{N-1} \int_a^b f(x)pdf(x)\:dx\\ & = & \dfrac{1}{N} \sum_{i=0}^{N-1} \int_a^b f(x)\:dx,\\ &=& \int_a^b f(x)\:dx,\\ &=&F\\ \end{array}$$

请记住，$pdf$ 等于 $1/(b-a)$ ，因此它抵消了积分符号右侧的 $(b-a)$ 项（第 3 行）。

4.2 泛化到任意概率密度函数

前述用于蒙特卡洛估计器的公式是非常基本的，因为它只有在随机变量 $X$ 的概率密度函数相等时才有效（即当 $X$ 是随机变量时，其在 [a,b] 区间内取任意值的概率相等，但现实中可能存在很多取值概率不等的情况） 。事实上可以将蒙特卡洛积分扩展到使用任意概率密度函数的随机变量。更通用的公式是：

$$\langle F^N \rangle = \dfrac{1}{N} \sum_{i=0}^{N-1} \dfrac{f(X_i)}{pdf(X_i)}.$$

这是更广义形式的蒙特卡洛估计量。 需要明确的是，分母中的 $pdf$ 就是指随机变量 $X$ 的 $PDF$ 。

与上一节的基本蒙特卡洛估计器一样，为了确保此公式有效，需要检查此估计器是否具有正确的期望值：

$$\begin{align*} E[\langle F^N \rangle ] & = E \left[ \dfrac{1}{N } \sum_{i=0}^{N-1} \dfrac{f(X_i)}{pdf X_i)} \right],\\ & = \dfrac{1}{N} \sum_{i=0}^{N-1} E\left[ \dfrac{f(X_i)}{pdf(X_i) }\right],\\ & = \dfrac{1}{N} \sum_{i=0}^{N-1} \int_\Omega \dfrac{f(x)}{pdf(x)} pdf(x)\;dx, \\ & = \dfrac{1}{N} \sum_{i=0}^{N-1} \int_\omega f(x) \; dx, \\ & = F \end{align*}$$

这里的技巧是第三行对第二行 $E \left[ f(X_i) \right]$ 替换：

$$E[f(x)] = \int_{\Omega} f(x) pdf(x) \;dx.$$

这是迄今为止所有内容中最重要的结果，也是后续蒙特卡洛所有算法的支柱。

让我们直观地了解上式中为什么需要将 $f(x)$ 除以 $pdf(x)$ 。$pdf$ 给出了随机变量 $X$ 获得某个值 $x_i$ 的概率。当从任意概率密度函数中抽取样本时，样本不是均匀分布的：在概率密度高的地方应当生成更多样本，在概率密度低的地方应当生成更少样本（见下图）。但在蒙特卡罗积分中，样本需要均匀分布。如果在函数的某个区域生成高密度样本（因为该区域的概率密度很高），则蒙特卡罗积分结果将明显有偏差。把 $f(x)$ 除以 $pdf(x)$ 将抵消这种影响，而实际上当概率密度很高时，将 $f(x)$ 除以 $pdf(x)$ 会降低这些样本在总和中的“权重”。因此，本质上是通过通过降低样本的贡献来补偿较高的样本密度。另一方面，当 $pdf$ 较低时，生成的样本较少，但如果将 $f(x)$ 除以一个较低的值（例如 $1$ 除以 $0.1$），则这些样本的贡献就会放大，通过赋予它们更多权重来补偿较小的样本密度。这就是 $f(x)$ 除以 $pdf(x)$ 的作用。它对样本的贡献进行加权以补偿从任意概率密度生成的样本分布不均匀这一事实：它缩小在较高密度区域生成的样本，并按比例放大在较低密度区域中生成的样本。

关于期望值的定义

当知道随机变量 $X$ 的分布后（注：大写字母 $X$ 表示向量，即多维随机变量），如何计算与其相关的另外一个随机变量 $f(X)$ 的期望值？

随机变量 $f(X)$ 的期望值可以定义为：

$$E[f(X)] = \int f(X) P_X(X) \;dX$$

其中 $P_X$ 是随机变量 $X$ 的概率分布。

在期望值概念基础上，可以延伸出很多有用的概念，如 方差 和标准差 的定义。例如，方差可以采用两种等价的方式定义：

$$\begin{array}{l} Var(X)& =& E[(X-E[X])^2],\\ &=&E[X^2] - E[X]^2. \end{array}$$

4.3 蒙特卡罗积分的性质

蒙特卡罗估计只不过是样本均值，只是我们将总体替换为一个实值的任意函数。因此，蒙特卡罗估计和样本均值具有相同的属性：

• 随着样本量接近无穷大，蒙特卡罗估计收敛到函数 $f(x)$ 的期望值。这是一个非常重要的属性。它表示，$N$ 越高，越有可能收敛到正确答案 $(F)$。

• 蒙特卡洛估计量估计量是无偏和一致的。

• 函数的收敛速度与其方差 $σ^2$ 成正比。蒙特卡洛估计量的方差为 $σ^2/n$ ， $\sigma[\langle F^N \rangle] \propto { 1 \over \sqrt{N} }$ （符号 $∝$ 表示“正比于”）。实践中，这意味着需要四倍样本才能将估计误差减少一半。这种收敛速度与积分维数无关，这与用于求解积分的其他确定性方法（例如黎曼和）不同。

• 中心极限定理精确了估计的渐近分布性质。换句话说，$\bar X_n$ 服从正态分布。

• 随机变量 $X_i$ 是随机的，这意味着样本均值 $\bar X_n / ⟨F^N⟩$ 是随机的，每次运行估计时，都会得到不同的 $\bar X_n / ⟨F^N⟩$ 值。

4.4 蒙特卡罗积分相对于确定性方法的优势

蒙特卡洛估计器背后的想法很简单，但它只是随着 1940 年代后期计算机技术的出现而兴起。多次评估函数并对结果求平均值是一项任务，计算机可以比人类速度快无数倍。能够有效运行这些模拟，有助于解决众多科学领域中的大量重要和复杂的问题。它也可以应用于其他领域，例如金融、计算机图形学。

但您可能想知道为什么我们会对这种技术感兴趣。事实上，计算积分的方法很多确定性方法（如：黎曼和），有些也很简单。但随着积分维数增加，这些确定性方法的使用成本越来越高。事实上，它们受到维数灾难的影响，随着积分维数增加，收敛速度呈指数级恶化（它们需要 $N^d$ 个样本在 $d$ 维上做积分）。

$$\int \cdots \int_D f(x_1, x_2, \cdots, x_n) dx_1 \cdots dx_n$$

蒙特卡罗积分原理可以很容易地扩展到更高维度，且其收敛速度与维度数无关。当有时必须对具有许多自变量的函数做积分（或多重积分）时，蒙特卡洛积分更适用。当然，对于固定数量样本，近似质量会随维数的增加而降低，但仍然可以保证以固定成本（样本数 $N$）获得结果。

注意：蒙特卡罗估计的关键是能够生成并使用随机数序列，我们使用这些随机数序列来评估函数 $f(x)$ 在所需区间 $[a,b]$ 内 x 的“随机”值。因此，如果您希望使用蒙特卡罗方法，掌握生成随机数的艺术是非常重要的。

4.5 改进蒙特卡罗积分：方差缩减

蒙特卡洛估计器适用于简单情况，但有时我们会希望解决更复杂的实际问题。如前所述，蒙特卡洛积分的问题在于其收敛速度，尽管它是很小的常数（ $O \sqrt N$）。出于此原因，很多研究都致力于开发减少误差（或方差）以加快收敛的技术，专业称呼为方差缩减（Variance Reduction）。例如：重要性采样就是其中一个例子。

4.6 小结

蒙特卡罗积分方法的优点：

• 简单，且很好地适应多维积分

• 公正、一致

• 并行性质，并行计算机的每个处理器都可以分配进行随机试验的任务。

蒙特卡罗积分方法的缺点：

• 收敛速度慢
• 很难评估函数 $f(x)$ 的方差，很难知道近似值的误差是多少

[1] Matt Pharr, Wenzel Jakob, and Greg Humphreys. Physically based rendering: From theory to implementation. Morgan Kaufmann, 2016.

[2] 概率论与数理统计，浙江大学第 4 版。

[3] Scratchapixel 2.0

[4] 解析 Monte-Carlo 算法（基本原理，理论基础，应用实践）

[5] Monte Carlo 数学原理

[6] 无偏估计量

[7] Philip Dutre, Philippe Bekaert, and Kavita Bala. Advanced global illumination. CRC Press, 2018.

[8] Bruce Walter, et al. “Microfacet models for refraction through rough surfaces.” Proceedings of the 18th Eurographics conference on Rendering Techniques. Eurographics Association, 2007.

[9] 低差异序列（一）- 常见序列的定义及性质

[10] 低差异序列（二）- 高效实现以及应用

[11] wikipedia, Low-discrepancy sequence

[12] wikipedia, Van der Corput sequence