2.1 蒙特卡罗：基础知识（Monte Carlo: Basics）
由于蒙特卡罗积分基于随机性（randomization），我们将以对概率（probability）和统计学（statistics）概念的简要回顾开始本章，这些概念为该方法提供了基础。这样做将使我们能够介绍基本的蒙特卡罗算法以及评估其误差的数学工具。
2.1.1 背景与概率回顾（Background and Probability Review）
我们将首先定义一些术语并回顾概率的基本概念。我们假设读者已经熟悉基本的概率概念；想要更全面了解这一主题的读者应参考谢尔顿·罗斯的 《概率模型导论》(Introduction to Probability Models)（2002）。
随机变量（random variable） \( X \) 是由某种随机过程选择的值。我们通常使用大写字母来表示随机变量，但对于一些表示特殊随机变量的希腊字母则另作例外。随机变量总是从某个定义域中抽取，这个定义域可以是离散的（例如，一个固定的有限可能集合），也可以是连续的（例如，实数集 \( \mathbf{R} \) ）。对随机变量 \( X \) 应用一个函数 \( f \) 会产生一个新的随机变量 \( Y = f(X) \) 。
例如，掷骰子的结果是从事件集合 \( X_i \in \{1,2,3,4,5,6\} \) 中抽取的离散随机变量。每个事件的概率为 \( p_i = \frac{1}{6} \) ，而概率的总和 \( \sum p_i \) 必然为一。像这样的随机变量，对于其所有潜在值具有相同的概率，称为 均匀的（uniform） 。给出离散随机变量概率的函数 \( p(X) \) 称为 概率质量函数（probability mass function）（PMF），因此在这种情况下我们可以等效地写作 \( p(X) = \frac{1}{6} \) 。
两个随机变量是 独立的（independent） ，如果一个的概率不影响另一个的概率。在这种情况下，两个随机变量的 联合概率（joint probability） \( p(X, Y) \) 由它们的概率的乘积给出：
\[ p(X,Y) = p(X)p(Y) \]
例如，表示骰子六个面的随机样本的两个随机变量是独立的。
对于 相依（dependent） 随机变量，一个的概率会影响另一个的概率。想象一个装有若干个黑球和若干个白球的袋子。如果我们从袋子中随机选择两个球，第二个球是白色的概率会受到第一个球颜色的影响，因为第一个球的选择改变了袋子中剩下的某种类型球的数量。我们将说第二个球的概率是以第一个球的选择 为条件的（conditioned） 。在这种情况下，选择两个球 \( X \) 和 \( Y \) 的联合概率为
\[ p(X,Y) = p(X)p(Y|X) \]
其中 \( p(Y|X) \) 是在给定 \( X \) 的值时 \( Y \) 的 条件概率（conditional probability）。
在后面的内容中，随机变量的概率通常会依赖于多个值；例如，在选择光源以进行照明采样时，第 12.6.3 节中的 BVHLightSampler 考虑了接收点的三维位置及其表面法线，因此光源的选择依赖于这些因素。然而，在随机变量依赖于多个变量且列举这些变量会使符号变得模糊的情况下，我们通常会省略这些变量。
一个特别重要的随机变量是 标准均匀随机变量（canonical uniform random variable） ，我们将其写作 \( \xi \) 。该变量在其定义域 \( [0,1) \) 中独立地以均匀概率取所有值。这个特定变量重要的原因有两个。首先，在软件中生成具有这种分布的变量很简单——大多数运行时库都有一个伪随机数生成器（pseudo-random number generator），正是用来实现这一点的。†（虽然蒙特卡罗理论是基于使用真正的随机数，但在实际应用中，一个编写良好的伪随机数生成器（PRNG）已经足够。pbrt 使用了一个质量非常高的伪随机数生成器，它返回的伪随机数序列在效果上与真正的随机数一样“随机”。真正的随机数可以通过测量诸如原子衰变或大气噪声等随机现象获得，对于那些认为伪随机数生成器不可接受的人来说，可以从诸如 www.random.org 之类的来源获取。）其次，我们可以将标准均匀随机变量 \( \xi \) 映射到一个离散随机变量，选择 \( X_i \), 若
\[ \sum_{j=1}^{i-1}{p_i} \leq \xi < \sum_{j=1}^{i}{p_j} \]
对于照明应用，我们可能希望根据场景中每个光源相对于所有光源的总功率 \( \Phi_i \) 来定义从每个光源采样照明的概率：
\[ p_i = \frac{\Phi_i}{\sum_j{\Phi_j}} \]
注意，这些 \( p_i \) 值的总和也为 1。给定这样的每个光源的概率，\( \xi \) 可以用来选择一个光源以采样照明。
随机变量的 累积分布函数（cumulative distribution function）（CDF）\( P(x) \) 是指从该变量的分布中取出的值小于或等于某个值 \( x \) 的概率（"随机变量 \( X \) 小于或等于 \( x \) 的概率"）：
\[ P(x) = \text{Pr}\{ X \leq x \} \]
对于骰子的例子，\( P(2) = \frac{1}{3} \) ，因为六种可能性中有两个小于或等于 2。
连续随机变量（Continuous random variables） 在连续域的范围内取值（例如，实数、单位球面上的方向或场景中形状的表面）。除了 \( \xi \) ，另一个连续随机变量的例子是取值于 0 到 2 之间的实数的随机变量，其中其取特定值 \( x \) 的概率与值 \( 2 - x \) 成正比：该随机变量取值接近 0 的可能性是取值接近 1 的两倍，依此类推。
概率密度函数（probability density function） （PDF）形式化了这一概念：它描述了随机变量取特定值的相对概率，是概率质量函数（PMF）的连续模拟（continuous analog）。PDF \( p(x) \) 是随机变量累积分布函数（CDF）的导数，
\[ p(x) = \frac{\text{d}P(x)}{\text{d}x} \]
对于均匀随机变量， \( p(x) \) 是一个常数；这是均匀性的直接结果。对于 \( \xi \) 我们有
\[ p(x) = \begin {cases} 1 &\text{x} \in [0, 1) \\ 0 &\text{otherwise} \end {cases} \]
PDF 必然是非负的，并且在其定义域上总是积分为 1。请注意，它们在某一点 \( x \) 的值 不 一定小于 1。
在定义域内给定一个区间 \( [a, b] \) ，对 PDF 进行积分可以得到随机变量位于该区间内的概率：
\[ \text{Pr}\{ x \in [a, b] \} = \int_{a}^{b}{p(x) \text{d}x} = P(b) - P(a) \]
这直接源于微积分的第一基本定理和 PDF 的定义。
2.1.2 期望值（Expected Values）
函数 \( f \) 的 期望值（expected value） \( E_p[f(x)] \) 被定义为在其定义域 \( D \) 上某些值 \( p(x) \) 的分布下函数的平均值。它被定义为
\[ E_p[f(x)] = \int_{D}{f(x)p(x)\text{d}x} \]
作为一个例子，考虑在 \( 0 \) 和 \( \pi \) 之间寻找余弦函数的期望值，其中 \( p \) 是均匀分布的。因为 PDF \( p(x) \) 必须在该区间内积分为 1， \( p(x) = 1 / \pi \) ，所以†（当计算均匀分布的期望值时，我们将从 \( E_p \) 中删除下标 \( p \) ）
\[ E[\cos{x}] = \int_{0}^{\pi} \frac{\cos{x}}{\pi} \text{d}x = \frac{1}{\pi}(\sin{\pi} - \sin{0}) = 0 \]
这正是预期的结果。（考虑 \( \cos{x} \) 在 \( [0,\pi] \) 区间的图形以了解原因。）
期望值具有一些有用的属性，这些属性源于其定义：
\[ E[af(x)] = aE[f(x)] \]
\[ E\left[\sum_{i=1}^{n}{f(X_i)}\right] = \sum_{i=1}^{n}{E[f(X_i)]} \]
我们将在后面章节的推导中反复使用这些性质。
2.1.3 蒙特卡罗估计（The Monte Carlo Estimator）
我们现在可以定义蒙特卡罗估计，它用于近似任意积分的值。假设我们想要评估一个一维积分 \( \int_{a}^{b}{f(x)\text{d}x} \) 。给定一组独立均匀随机变量 \( X_i \in [a,b] \) ，蒙特卡罗估计表示估计量(estimator)的期望值为
\[ F_n = \frac{b-a}{n} \sum_{i=1}^{n}{f(X_i)} \]
\( E[F_n] \) 等于积分。这个事实可以通过几个步骤来证明。首先，注意到与随机变量 \( X_i \) 对应的 PDF \( p(x) \) 必须等于 \( 1/(b-a) \) ，因为 \( p \) 不仅必须是一个常数，还必须在域 \( [a,b] \) 上积分为 1。然后，利用方程 (2.4) 和 (2.5) 的性质进行代数运算，可以得出：
\[ \begin{align} E[F_n] &= E\left[\frac{b-a}{n}\sum_{i=1}^{n}{f(X_i)}\right] \\ &= \frac{b-a}{n}\sum_{i=1}^{n}{E[f(X_i)]} \\ &= \frac{b-a}{n}\sum_{i=1}^{n}{\int_{a}^{b}{f(x)p(x)\text{d}x}} \\ &= \frac{1}{n}\sum_{i=1}^{n}{\int_{a}^{b}{f(x)\text{d}x}} \\ &= \int_{a}^{b}{f(x)\text{d}x} \end{align} \]
将此估计量推广到多个维度或复杂积分域是简单的：从均匀多维的 PDF 中取 \( n \) 个独立样本 \( X_i \) ，并以相同的方式应用估计量。例如，考虑三维积分
\[ \int_{z_0}^{z_1} \int_{y_0}^{y_1} \int_{x_0}^{x_1} f(x,y,z) \text{d}x \text{d}y \text{d}z \]
如果样本 \( X_i = (x_i,y_i,z_i) \) 是从立方体 \( [x_0,x_1] \times [y_0,y_1] \times [z_0,z_1] \) 中均匀选择的，则 PDF \( p(X) \) 是常数值
\[ \frac{1}{(x_1-x_0)} \frac{1}{(y_1-y_0)} \frac{1}{(z_1-z_0)} \]
估计量是
\[ \frac{(x_1-x_0)(y_1-y_0)(z_1-z_0)}{n} \sum_{i=1}^{n} f(X_i) \]
对均匀随机变量的限制可以通过小的推广来放宽。这是一个重要的步骤，因为仔细选择抽取样本的 PDF 可以引出减少蒙特卡罗方法中的误差的关键技术，这将在第 2.2.2 节中介绍。如果随机变量 \( X_i \) 是从 PDF \( p(x) \) 中抽取的，则估计量
\[ F_n = \frac{1}{n} \sum_{i=1}^{n} \frac{f(X_i)}{p(X_i)} \]
可以用来估计积分。对 \( p(x) \) 的唯一限制是对于 \( |f(x)| > 0 \) 所有的 \( x \) 必须是非零的。
同样不难看出，这个估计量的期望值是所需的函数 \( f \) 的积分：
\[ \begin{align} E[F_n] &= E \left[ \frac{1}{n} \sum_{i=1}^{n} \frac{f(X_i)}{p(X_i)} \right] \\ &= \frac{1}{n} \sum_{i=1}^{n} \int_{a}^{b} \frac{f(x)}{p(x)} p(x) \text{d}x \\ &= \frac{1}{n} \sum_{i=1}^{n} \int_{a}^{b} f(x) \text{d}x \\ &= \int_{a}^{b} f(x) \text{d}x \end{align} \]
我们现在可以理解在 RandomWalkIntegrator 的实现中系数 \( 1/(4\pi) \) 的含义：在面积为 \( 4\pi \) 的单位球面上均匀采样方向。由于 PDF 在采样域上是归一化的，因此它必须具有常数值 \( 1/(4\pi) \) 。当应用方程（2.7）的估计量时，该值出现在除数中。
通过蒙特卡罗方法，可以任意选择样本数量 \( n \) ，而不受被积函数维度的限制。这是蒙特卡罗方法相对于传统确定性求积技术的另一个重要优势，后者通常需要的样本数量在维度上是指数级的。
2.1.4 蒙特卡罗估计量的误差（Error in Monte Carlo Estimators）
显示蒙特卡罗估计量收敛到正确答案并不足以证明其使用的合理性；其收敛速度也很重要。 方差（Variance） ，即一个函数与其期望值之间的期望平方偏差，是表征蒙特卡罗估计量收敛性的一种有用方法。估计量 \( F \) 的方差定义为
\[ V[F] = E[(F - E[F])^2] \]
由此可得
\[ V[aF] = a^2 V[F] \]
该性质和方程（2.5）提供了方差的另一种表达式：
\[ V[F] = E[F^2] - E[F]^2 \]
因此，方差是平方的期望值减去期望值的平方。
如果估计量是独立随机变量的和（如蒙特卡罗估计量 \( Fn \) ），那么和的方差是各个随机变量方差的总和：
\[ V \left[ \sum_{i=1}^{n} X_i \right] = \sum_{i=1}^{n} V[X_i] \]
从方程（2.10）可以很容易地证明，方差随着样本数量 \( n \) 的增加而线性减少 。由于方差是平方误差，因此蒙特卡罗估计中的误差仅以样本数量的 \( O(n^{-1/2}) \) 的速率下降。尽管标准求积技术在一维中以更快的速率收敛，但随着被积函数维度的增加，它们的性能会呈指数级恶化，而蒙特卡罗的收敛速率与维度无关，使得蒙特卡罗成为高维积分的唯一实用数值积分算法。
蒙特卡罗误差减少的 \( O(n^{-1/2}) \) 特性在观察在所有像素中逐步采集额外样本的场景的渐进渲染时显而易见。因为这相对较少的额外工作，在前几个样本中，当样本数量翻倍时，图像迅速改善。随后，一旦采集了数十或数百个样本，每增加一个样本的翻倍所需的时间就会大大增加，图像中剩余的误差消失也需要很长时间。
随着样本数量的增加，方差线性减少使得比较不同的蒙特卡罗估计量变得容易。考虑两个估计量，其中第二个的方差是第一个的一半，但计算一个估计所需的时间是第一个的三倍；那么这两者中哪个更好？在这种情况下，第一个是更可取的：它可以在第二个所消耗的时间内获取三倍的样本，这样就能实现 \( 3\times \) 的方差减少。这个概念可以用估计量 \( F \) 的 效率（efficiency） 来概括，效率定义为
\[ \epsilon[F] = \frac{1}{V[F]T[F]} \]
其中 \( V[F] \) 是其方差，\( T[F] \) 是计算其值的运行时间。
并非所有积分的估计量都有期望值等于积分的值。这些估计量被称为 有偏的（biased），其中差值
\[ \beta = E[F] - \int f(x)\text{d}x \]
是偏差的量。如果有偏估计量能够比无偏估计量更快地接近正确结果，它们仍然可能是可取的。Kalos 和 Whitlock（1986 年，第 36-37 页）给出了以下例子：考虑计算均匀分布 \( X_i \sim p \) 在区间 0 到 1 上的均值估计的问题。可以使用以下估计量：
\[ \frac{1}{n} \sum_{i=1}^{n} X_i \]
或者可以使用有偏估计量
\[ \frac{1}{2} \max(X_1,X_2,...,X_n) \]
第一个估计量是无偏的，但其方差的阶数为 \( O(n^{-1}) \) 。第二个估计量的期望值是
\[ 0.5 \frac{n}{n+1} \neq 0.5 \]
因此它是有偏的，尽管它的方差是要好得多的 \( O(n^{-2}) \) 。这个估计量具有一个有用的特性，即当样本数量 \( n \) 趋向于无穷大时，它的误差趋向于 0；这样的估计量是 一致的（consistent） 。†（作为一个技术性说明，一个具有无限方差的估计量可能是无偏的，但不一致。然而，这样的估计量通常不会出现在渲染中。） 在 pbrt 中使用的大多数蒙特卡罗估计量是无偏的，唯一的例外是 SPPMIntegrator ，它实现了光子映射算法。
均方误差（mean squared error）（MSE）与方差密切相关，其定义为估计量与真实值之间平方差的期望
\[ MSE[F] = E \left[ (F - \int f(x) \text{d}x)^2 \right] \]
对于无偏估计量，均方误差 MSE 等于方差；否则，它是方差与估计量的平方偏差之和。
可以以封闭形式计算一些简单估计量的方差和均方误差，但对于大多数在渲染中感兴趣的估计量，这并不可行。然而，能够量化这些值仍然是有用的。为此，可以使用一组独立随机变量 \( X_i \) 计算 样本方差（sample variance） 。方程 (2.8) 指出了一种计算一组 \( n \) 个随机变量 \( X_i \) 的样本方差的方法。如果 样本均值（sample mean） 被计算为它们的平均值 \( \overline{X} = (1/n) \sum X_i \) ，那么样本方差为
\[ \frac{1}{n-1} \sum_{i=1}^{n}(X_i - \overline{X})^2 \]
对 \( n-1 \) 而不是 \( n \) 的除法是 贝塞尔校正（Bessel’s correction） ，确保样本方差是方差的无偏估计。（另见 B.2.11 节，其中介绍了一种数值稳定的样本方差计算方法。）
样本方差本身就是方差的一个估计，因此它本身也有方差。例如，考虑一个随机变量，其值在 99.99% 的时间里为 1，而在 0.01% 的时间里为 一百万 。如果我们对它进行了十次随机抽样，所有样本的值都是 1，那么样本方差会表明这个随机变量的方差为零，尽管它的实际方差要高得多。
如果可以计算出积分 \( \tilde{F} \approx \int f(x)\text{d}x \) 的准确估计（例如，使用大量样本），那么均方误差可以通过以下方式估计
\[ MSE[F] \approx \frac{1}{n} \sum_{i=1}^{n}(f(X_i) - \tilde{F})^2 \]
在 pbrt 的发行版中提供的 imgtool 工具程序可以通过其 diff 选项计算图像相对于参考图像的均方误差（MSE）。