2.1 Monte Carlo Basics
由于蒙特卡洛积分是基于随机性的算法，所以在本章的开始，我们先回顾一些概率学和统计学的相关知识。
2.1.1 Background and Probability Review
一个随机变量$X$是通过某种随机过程选择的值，我们通常用大写字母表示一个随机变量，对于某些特殊的随机变量我们会依照惯例使用希腊字母表示。随机变量总是从某个域中选取，该域可以是离散的，也可以连续的。对一个随机变量$X$应用函数$f$会生成一个新的随机变量$Y=f(X)$
例如，掷骰子的结果是一个从一组事件$X_i \in \lbrace 1,2,3,4,5,6 \rbrace $抽样的随机变量。每个事件的概率均等，且所有事件的概率总和为1。当随机变量对于其他所有潜在值具有的相同的概率时，我们称该随机变量是均分分布的。给出离散变量概率的函数$p(X)$被称为概率质量函数PMF，对于掷骰子这个例子，我们有$P(X)=\frac{1}{6}$。
如果一个随机变量的概率不会影响另一个随机变量的概率，那我们说这两个随机变量是相互独立的。此时，两个变量的联合概率joint probablity $p(X, Y)$由它们概率的乘积表示，即$P(X, Y)=p(X)p(Y)$。
对于相互依赖的随机变量，一个变量的概率会影响另一个变量的概率。考虑一个装有若干个黑球和若干个白球的袋子。如果我们从袋子里随机选出两个球，第二个球是白球的概率会受到第一个球颜色的影响，因为第一个球的选择改变了袋子里剩余某种类型球的数量。我们称第二个球的概率以第一个球的选择为条件。
在本书后续的内容中，随机变量的概率通常会以许多值为条件，比如，在选择采样光照的光源时，BVHLightSampler会考虑到着色点的位置与法线，因此光源的选择是以它们为条件的。
canonical uniform random variable， 这是一种特别重要的随机变量，记为$\zeta$，它的值域是$[0,1)$，并且每个值都有相同的概率。它具备两个重要的性质，首先对于计算机程序来说，$\zeta$很容易生成，此外，我们可以将$\zeta$映射到离散随机变量上，使用下面这个公式：
\[\sum_{j=1}^{i-1}p_j \le \zeta < \sum_{j=1}^{i}p_j\]其中，$p_j$是一系列满足$\sum_jp_j=1$的概率值。让我们来举一个例子进一步说明。在光照计算中，我们可以定义出每个光源的采样概率是：
\[p_i=\frac{\Phi _i}{\sum_j \Phi_j}\]其中，$\Phi _i$表示第$i$个光源的功率，$\sum_j \Phi_j$表示所有光源的总功率。这样一来，给定一个均匀随机变量$\zeta$，我们可以根据光源的功率来决定采样哪个光源，从而实现基于功率的光源采样
一个随机变量的累计分布函数CMF $P(X)$表示变量分布中该值小于或等于某个值的概率：
\[P(X) = Pr\lbrace X \le x \rbrace\]还是以掷骰子为例，$P(2)=\frac{1}{3}$表示点数小于等于2的概率为$\frac{1}{3}$
PDF描述了随机变量取某个特定值的相对概率，离散随机变量的PMF的一种连续性的形式。此外，PDF $p(x)$是随机变量CDF的导数。对于均匀的随机变量来说，PDF $p(x)$是一个常数。如果我们通过导数来理解，是因为均匀变量的CMF的斜率是固定的。我们通过公式将PDF与CMF的关系表示出来:
\[p(x) = \frac{dP(x)}{dx}\]PDF具有两个重要的性质：
- 非负性：$p(x) \ge 0$
- 归一化：$p(x)$在其整个定义域上的积分等于1，也就是总概率为1
要计算随机变量在给定区间$[a, b]$上的概率，需要在该区间上对PDF积分：
\[Pr\lbrace x \in [a, b] \rbrace = \int_a^b p(x)dx = P(b)-P(a)\]2.1.2 Expected Values
一个函数$f$的期望值$E_p[f(x)]$被定义为函数$f$在某个值分布$p(x)$下的平均值，数学表达式为： \(E_p[f(x)] = \int_Df(x)p(x)dx\)
其中：
- $D$是函数的定义域
- $f(x)$是需要计算期望值的函数
- $p(x)$是概率密度函数，描述了$x$在$D$上的分布情况
我们来举一个例子，要求余弦函数$cos(x)$在区间$[0, \pi]$上的期望值，其中$p(x)$是均匀分布的，也就是说$x$在$[0, \pi]$上的每个点的概率密度相同，即:
\[p(x) = \frac{1}{\pi}\]所以，余弦函数$cos(x)$在区间$[0, \pi]$上的期望值的计算如下：
\[E[cos(x)]=\int_0^{\pi}cos(x)\cdot \frac{1}{\pi}dx=\frac{1}{\pi}\int_0^{\pi}cos(x)\cdot dx = \frac{1}{\pi}(sin \pi-sin0) = 0\]此外，期望值有一些重要的形式，这些性质直接来源于期望值的定义：
线性性质： \(E[af(x)]=aE[f(x)]\)
加法性质： \(E\Big[\sum_{i=1}^nf(X_i)\Big]=\sum_{i=1}^nE[f(X_i)]\)
2.1.3 The Monte Carlo Estimator
现在我们可以来定义蒙特卡洛估值器了，它可以用来近似一个任意积分的值。假设我们想要计算一个一维积分$\int_a^bf(x)dx$，如果给定$X_i$是在区间$[a, b]$上均匀分布的独立随机变量，那么我们可以证明 估值器的期望值等于积分值，证明过程如下：
首先，蒙特卡洛估值器的定义如下：
\(F_n = \frac{b-a}{n}\sum_{i=1}^nf(X_i)\) 我们从估值器的期望值$E[F_n]$开始：
\[E[F_n] = E\Big[\frac{b-a}{n}\sum_{i=0}^nf(X_i) \Big]\]我们已经知道期望具有线性性质，所以可得：
\[E[F_n] = \frac{b-a}{n} \sum_{i=0}^nE[f(X_i) ]\]由于每个$X_i$是在区间$[a, b]$上均匀分布，因为它的概率密度函数$p(x)$为$1/(b-a)$，所以$f(X_i)$的期望值$E[f(X_i)]$为：
\[E[f(X_i)] = \int_a^b f(x)p(x)dx = \int_a^b f(x)\frac{1}{b-a}dx\]我们将期望值代入到估值器的期望值计算中：
\[\begin{aligned} E[F_n] & = \frac{b-a}{n} \sum_{i=0}^n \int_a^b f(x)\frac{1}{b-a}dx \\ & = \frac{1}{n} \sum_{i=0}^n \int_a^b f(x)dx \end{aligned}\]不难注意到，求和中的每一项都是相同的，所以我们进一步推导出：
\[\begin{aligned} E[F_n] & = \frac{1}{n} n \int_a^b f(x)dx \\ & = \int_a^b f(x)dx \end{aligned}\]在进一步深入之前，我们来总结一下：
- 蒙特卡洛估值器：通过在区间$[a. b]$上均匀采样并计算这些采样点的函数值平均值来近似积分
- 期望值的证明：利用期望值的线性性质与均匀分布的PDF，证明了估值器的期望值等于实际的积分值
接下来，我们就可以试着将估值器应用到多维或者复杂积分域上了，例如，下面是一个三维积分：
\[\int_{z_0}^{z_1} \int_{y_0}^{y_1} \int_{x_0}^{x_1} f(x, y, z) dxdydz\]如果我们依然采用刚才的思路，则我们需要从三维立方体$[x_0, x_1]\times[y_0, y_1]\times[z_0, z_1]$中均匀地采样$n$个采样点$X_i=(x_i, y_i, z_i)$。同样的，这些采样点的PDF也是常数：
\[p(X) = \frac{1}{(x_1-x_0)(y_1-y_0)(z_1-z_0)}\]对应地，估值器的定义如下：
\[F_n = \frac{(x_1-x_0)(y_1-y_0)(z_1-z_0)}{n}\sum_{i=1}^nf(X_i)\]但是，我们在这里做一个微小的调整，我们使用更通用的方式来选择采样点，也就是不再均匀采样，而是从特定的非均匀的概率密度函数$p(x)$中抽取样本，这样做的原因有如下几点：
- 提高采样效率：均匀随机采样虽然简单，但是在某些情况下效率不高。比如说，被积函数$f(x)$在某些区域变化剧烈而某些区域变换较小，那么在变化不大的区域内的均匀采样会浪费大量的采样点，从而增加近似的方差。
- 减少方差：从上一点中我们可以得出结论：从非均匀分布中抽样可以显著减少估计的方差。并且我们可以引入一个新的概念，重要性分布采样，它的基本思想是更多地在对结果贡献值大的区域上采样。从数学的角度上来说，如果$p(x)$与$f(x)$形状相似，则$\frac{f(x_i)}{p(x_i)}$会在所有$x$上接近一个常数，从而减少近似的方差。
- 更好地适应高维积分：在渲染领域中，有大量的高维积分问题，而在高维空间中，均匀采样的效率急剧下降，是因为高维空间的体积增长很快，均匀分布的采样点数量需要指数级地增长才能覆盖整个空间，而非均匀采样可以集中在对积分值贡献较大的区域，从而更高效地进行积分计算。
我们回到蒙特卡洛估值器上，如果随机变量$X_i$是从概率密度函数中选取的，则估值器的定义如下：
\[F_n=\frac{1}{n}\sum_{i=0}^n\frac{f(X_i)}{p(X_i)}\]需要注意的是，对于$p(x)$存在一个关键的限制条件：$p(x)$在$f(x) \ne 0$的所有$x$上必须为正，从而确保所有对积分有贡献的样本都能被采样到。
我们来继续证明估值器的期望值等于积分值：
\[\begin{aligned} E[F_n] & = E\Big[ \frac{1}{n}\sum_{i=1}^n \frac{f(X_i)}{p(X_i)}\Big] \\ & = \frac{1}{n}\sum_{i=1}^n \int_a^b\frac{f(x)}{p(x)}p(x)dx \\ & = \frac{1}{n}\sum_{i=1}^n \int_a^bf(x)dx \\ & = \int_a^bf(x)dx \end{aligned}\]然而，尽管理论上蒙特卡洛估值器的期望值等于积分值，但这是因为在无限多次采样的情况下，样本均值会收敛到期望值。但在实际应用中，由于样本数量有限，并且结果会依赖于随机采样，所有结果必然包含一定的随机误差。
2.1.4 Error in Monte Carlo Estimators
现在我们知道，蒙特卡洛估值器在理论上可以得到正确结果，在实际应用中必然存在误差后，我们就需要通过某种方式来评估一个蒙特卡洛估值器的收敛性。所以，我们在这里引入方差这个概念，它表示的是一个函数偏离其期望值的平方的期望值。对于一个估值器$F$，我们定义它的方差为：
\[V[F]=E[(F-E[F])^2]\]并由此得出：
\[V[aF]=a^2V[F]\]通过方差的定义和期望值的加法性质，我们还可以得到：
\[\begin{aligned} V[F] & = E[(F-E[F])^2] \\ & = E[F^2 - 2FE[F] + (E[F])^2] \\ & = E[F^2] - 2E[F]E[F] + E[(E[F])^2] \\ & = E[F^2] - 2(E[F])^2 + E[(E[F])^2 \\ \end{aligned}\]其中，$E[(E[F])^2]$是一个常数，任何常数的期望都等于常数自身，也就是等于$(E[F])^2$，代入到推导式中：
\[\begin{aligned} V[F] & = E[F^2] - 2(E[F])^2 + E[(E[F])^2 \\ & = E[F^2] - 2(E[F])^2 + (E[F])^2 \\ & = E[F^2] - (E[F])^2 \end{aligned}\]因此，我们可以得出结论：方差等于平方的期望减去期望的平方
如果估值器是独立随机变量的和（如蒙特卡洛估值器$F_n$），那么和的方差就等于是各个随机变量方差的和，即：
\[V\Big[ \sum_{i=1}^n X_i \Big] = \sum_{i=1}^nV[X_i]\]那么我们可以发现，方差与样本数量$n$呈反比关系，样本数量越大，方差越小。即方长随样本数量$n$的变化率为$O(\frac{1}{n})$，同时方差是误差的平方，因此在蒙特卡洛近似中，误差随样本数量$n$的变化率是$O(\frac{1}{\sqrt n})$
作为对比，标准的数值积分方法在一维中的收敛速度较快，但是随着被积函数的维度增加，性能则会呈指数级下降，而蒙特卡洛方法的收敛速度与维数无关，因此蒙特卡洛是高维积分问题中唯一实用的数值积分方法。
应用在渲染领域中，蒙特卡洛$O(\frac{1}{\sqrt n})$的特性会使得在初期增加样本数量的效果明显，而后期改进则会逐渐减慢。
所以，在对蒙特卡洛估值器的收敛速度有一定了解后，我们就可以比较两个估值器之间的优劣了。假设有两种估计器，第二种估计器的方差是第一种的一半，但计算一个估计值需要的时间是第一种的三倍。在这种情况下，第一种估计器更好，因为在第二种估计器计算一个估计值所需的时间内，第一种估计器可以增加三倍的样本数量，从而达到三倍的方差减少。
我们可以将估值器的效率抽象出一个概念，并定义为：
\[\epsilon[F]=\frac{1}{V[F]T[F]}\]其中，$V[F]$表示估值器的方法，而$T[F]$表示计算估值器所需要的时间。
并非所有估值器都可以计算出与积分相等的期望值，我们将这种估值器称为biased estimator。同时我们将估计值的期望值与实际积分值之间的差距定义为偏差，我们可以将偏差理解为估值器系统性地偏离真实值的程度，并且将偏差值定义为：
\[\beta = E[F] - \int(fx)dx\]通常，如果有偏差的估值器可以以更快的速度收敛，则我们认为该biased estimator仍然是可取的。
我们考虑下面这个例子：计算区间$[0, 1]$上均匀分布的随机数的平均值，使用蒙特卡洛近似，则有估值器定义如下：
\[\frac{1}{n}\sum_{i=1}^nX_i\]对应方差的阶数$O(\frac{1}{n})$。
或者，我们也可以使用有偏差的估值器：
\[\frac{1}{2}max(X_1, X_2,...,X_n)\]它的期望值为$0.5\frac{n}{n+1}$，与实际的期望值0.5不符，所以它确实是有偏差的估值器。它的方差的变化率为$O(n^{-2})$。尽管这个估值器是有偏差的，但是我们可以从它的方差阶数上得知，当样本数量$n$趋近于无限时，其误差趋近于0。我们将这种估值器称为一致估值器consistent estimator。
均方误差（Mean Squared Error，MSE）是衡量估计器质量的一个指标。它被定义为估计值和真实值之间差异的平方的期望：
\[MSE[F] = E\Big[ \Big( F-\int f(x)dx \Big)^2 \Big]\]对于无偏估计器，MSE等于方差；否则，MSE是方差与估计器偏差平方和的总和
在渲染中，我们通常无法直接计算一些估值器的方差与MSE，但是我们可以通过样本来近似这些值。样本方差可以通过一组独立的随机变量$X_i$来计算。假设样本均值（Samples mean）是这些随机变量的平均值：
\[\bar{X}=\frac{1}{n}\sum_{i=1}^nX_i\]则样本方差等于为：
\[s^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})^2\]在这里，我们使用$n-1$而非$n$是用到了贝塞尔矫正，以确保样本方差是方差的无偏估计。
样本方差本身是方差的一个估计，因此它自身也有方差。举个例子，如果一个随机变量有99.99%的概率取值为1，0.01%的概率取值为一百万，如果我们取十个随机样本并且它们的值都是1，那么样本方差会建议这个随机变量有零方差，尽管它的实际方差要大得多。