PBRT学习笔记1 -- 蒙特卡洛积分

Monte Carlo Integration（蒙特卡洛积分）

Probability Review（概率论回顾）

一个随机变量（random variable）$X$是某个随机过程（random process）的取值，我们一般使用大写的希腊字母来表示一个随机变量。随机变量总是从某个域（domain）中抽取的，域可以是离散的或是连续的，我们可以应用函数$f$到随机变量$X$上产生新的随机变量$Y=f(X)$。

对于一个离散的随机变量$X$，某个事件发生的概率写作$p(X)$，这个概率被称为概率质量函数（probability mass function, PMF）。

两个独立事件的联合概率（joint probability）是他们概率的乘积，$p(X,Y)=p(X)p(Y)$。

对于相关随机变量（dependent random variables），一个的概率会影响另一个的概率，它们的联合概率建立在条件概率（conditional probability）之上，$p(X,Y)=p(X)p(Y|X)$。$p(Y|X)$是给定$X$的情况下$Y$的概率。

一个重要的随机变量是规范均匀随机变量（canonical uniform random variable），写作$\xi$，这个变量从域$[0,1)$中独立均匀取值。这个随机变量重要的原因有二。第一，这个随机变量的分布在软件中非常好实现，非常多库中都存在伪随机数生成器（pseudo-random number generator, PRNG）。第二，这个随机变量可用于映射到其他离散随机变量上，例如当${\sum }_{j=1}^{i-1}{p}_j\le \xi <{\sum }_{j=1}^i{p}_j$时，取值$X_i$。

累积分布函数（cumulative distribution function, CDF）$P(x)$是随机变量的取值小于或等于某个值$x$的概率： P ( x ) = Pr { X ≤ x } P(x)=\text{Pr}\{X\leq{}x\} P(x)=Pr{X≤x}

连续随机变量（Continuous random variables）是从连续的域中进行取值的随机变量，连续随机变量使用概率密度（Probability density function, PDF）描述随机变量取特定值的相对概率，这是对 PMF 的连续模拟。

连续随机变量的 PDF$p(x)$是它 CDF 的导数：$p(x)=\frac{\mathrm{d}P(x)}{\mathrm{d}x}$

对于均匀随机变量，$p(x)$是一个常数，例如$\xi$的 PDF 为：$p(x)=\{\begin{array}{ll}1& x\in [0,1)\\ 0& \text{otherwise}\end{array}$。PDF 的函数值始终非负，并且积分域上的总积分值为 1。

给定域上的一个区间$[a,b]$，对 PDF 进行积分可得出随机变量位于区间内的概率： Pr ⁡ { x ∈ [ a , b ] } = ∫ a b p ( x )   d x = P ( b ) − P ( a ) \Pr\{x\in[a,b]\}=\int_a^bp(x)\,\mathrm{d}x=P(b)-P(a) Pr{x∈[a,b]}=∫ab​p(x)dx=P(b)−P(a)

Expected Values（期望值）

使用一个随机变量$f(X)=Y$在一个域为$D$的函数$f$上，用分布$p(x)$进行采样，得到的样本均值（Sample Mean）由采样值$Y_1,Y_2,...,Y_n$的平均值给出：

$$\begin{array}{c}\bar{Y}=\frac{1}{n}\sum _{i=1}^n{Y}_i\end{array}$$

随机变量$Y$的期望，某个采样点$Y_i$的期望，样本均值$\bar{Y}$的期望，都由函数$f$和分布函数$p(x)$的积分给出：

$$\begin{array}{c}\mu =E[Y]=E[Y_i]=E[f(X_i)]=E[\bar{Y}]={\int }_{D}f(x)p(x)\mathrm{d}x\end{array}$$

期望有一些有用的性质：

$$\begin{array}{c}E[aX]=aE[X]\end{array}$$

$$\begin{array}{c}E[\sum _{i=1}^n{X}_i]=\sum _{i=1}^{n}E[X_i]\end{array}$$

Variance（方差）

方差（Variance），即随机变量或函数与其期望值的期望平方偏差，是表征随机变量收敛性的有效方法。

用随机变量$Y$采样在$p(x)$下的函数$f(x)$，它的方差为：

$$\begin{array}{c}V[Y]={\sigma }^2=\int f(x{)}^{2}p(x{)}^2\mathrm{d}x-{\mu }^2\end{array}$$

方差有以下性质：

$$\begin{array}{c}V[X]=E[(X-E[X]{)}^2]\end{array}$$

$$\begin{array}{c}V[X]=E[X^2]-E[X{]}^2\end{array}$$

如果一系列采样点$X_1,X_2,...,X_n$相互独立，则总和的方差就是各个随机变量方差的总和：

$$\begin{array}{c}V[\sum _{i=1}^n{X}_i]=\sum _{i=1}^{n}V[X_i]\end{array}$$

Monte Carlo: Basics

The Monte Carlo Estimator（蒙特卡洛估计器）

蒙特卡洛积分是一种数值积分方法（其他数值积分方法），使用蒙特卡洛估计器可以近似任意积分的值，假设我们要计算一维积分${\int }_b^{a}f(x)dx$，给定一组独立的均匀随机变量值$X_i\in [a,b]$，使用蒙特卡洛估计器$F_n$近似积分，有：

$$\begin{array}{c}{F}_n=\frac{b-a}{n}\sum _{i=1}^{n}f(X_i)\end{array}$$

这个近似结果的期望值$E[F_n]$等于其积分值，由于是均匀随机变量，因此对于任意$x$值，$p(x)=1/(b-a)$，根据(2), (3)可得：

$$\begin{array}{rl}E[F_n]& =E[\frac{b-a}{n}\sum _{i=1}^{n}f(X_i)]\\ & =\frac{b-a}{n}\sum _{i=1}^{n}E[f(X_i)]\\ & =\frac{b-a}{n}\sum _{i=1}^{n}E[X]\\ & =\frac{b-a}{n}\sum _{i=1}^n{\int }_a^{b}f(x)p(x)\mathrm{d}x\\ & ={\int }_a^{b}f(x)\mathrm{d}x\end{array}$$

我们可以非常直观的推广这个积分器到多维，考虑一个三维积分${\int }_{z_0}^{z_1}{\int }_{y_0}^{y_1}{\int }_{x_0}^{x_1}f(x,y,z)\mathrm{d}x\mathrm{d}y\mathrm{d}z$，使用均匀随机变量$X=(x_i,y_i,z_i)$采样，对于任意$X$值，$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)$。

将其推广到非均匀采样的随机变量$X_i$，取自 PDF$p(x)$，有：

$$\begin{array}{c}{F}_n=\frac{1}{n}\sum _{i=1}^n\frac{f(X_i)}{p(X_i)}\end{array}$$

唯一的限制是对于所有$|f(x)|>0$的 x 值，PDF 都是非零的。类似于均匀情况，这个近似结果的期望值$E[F_n]$也等于其积分值。

$$\begin{array}{rl}E[F_n]& =E[\frac{1}{n}\sum _{i=1}^n\frac{f(X_i)}{p(X_i)}]\\ & =\frac{1}{n}\sum _{i=1}^n{\int }_a^b\frac{f(x)}{p(x)}p(x)\mathrm{d}x\\ & =\frac{1}{n}\sum _{i=1}^n{\int }_a^{b}f(x)\mathrm{d}x\\ & ={\int }_a^{b}f(x)\mathrm{d}x\end{array}$$

蒙特卡洛的优点在于样本数量可以任意选取，不管被积函数的维度如何，这是蒙特卡洛相对于传统确定性求积技术的另一个重要优势，传统确定性求积技术通常需要大量在维度上呈指数的样本。

Error in Monte Carlo Estimators（蒙特卡洛估计器误差）

证明蒙特卡罗估计器收敛于正确答案并不足以证明它的使用，它的收敛速度也很重要，收敛速度可以由方差进行评判。

假设估计器$F$需要估计积分${\int }_b^{a}f(x)\mathrm{d}x$，考虑使用概率密度为$p(x)$的随机变量$X$来采样函数$f(x)$，计算$V[F]$，有：

$$\begin{array}{rl}V[F]={\hat{\sigma }}^2& =V[\frac{1}{n}\sum _{i=1}^n\frac{f(X_i)}{p(X_i)}]\\ & =\frac{1}{n^2}V[\sum _{i=1}^n\frac{f(X_i)}{p(X_i)}]\\ & =\frac{1}{n^2}\sum _{i=1}^{n}V[Y_i]\\ & =\frac{1}{n}V[Y]\end{array}$$

这里我们设使用单样本估计积分的方差为$V[Y]$，可以推断出方差随着样本数量线性递减。

由于方差是平方误差，因此蒙特卡罗估计中的误差（标准差）仅以$O(n^{-1/2})$的速率下降。虽然标准正交技术（standard quadrature techniques or Gaussian quadrature）在一维上收敛速度更快，但随着被积函数维数的增加，它们的性能会呈指数级下降，而蒙特卡罗的收敛速度与维数无关，使蒙特卡罗成为唯一实用的高维积分数值积分算法。

$O(n^{-1/2})$的性质使得在渲染场景，对每个 pixel 进行采样时，样本数较少时，蒙特卡洛增加样本数对画面的改善是非常迅速的，但随着样本数的增加，误差消除的速度会逐渐增加。

方差的线性递减性还可以很好的用于指导比较不同的蒙特卡洛积分器，考虑两个估计器，其中第二个估计器在同等样本数下的方差是第一个的一半，但需要三倍的时间来计算估计;两者中哪一个更好?在这种情况下，首选第一种方法，因为在第二个估计器的消耗的时间内，第一个可以采集三倍的样本，在这种情况下，它会实现三倍的方差减少。由此可以给估计器$F$封装一个新的属性 效率（efficiency）$ϵ$。有：

$$\begin{array}{c}ϵ[F]=\frac{1}{V[F]T[F]}\end{array}$$

其中，$V[F]$表示方差，$T[F]$表示计算到此方差值需要的时间。

此外，并不是所有的估计器的期望值都等于积分的值，这种估计器被称为是有偏（biased）的，差值表示为：

$$\begin{array}{c}\beta =E[F]-\int f(x)\mathrm{d}x\end{array}$$

如果有偏估计器能够比无偏估计器更快地接近正确的结果，那么有偏估计器仍然是可取的。

例如 Kalos and Whitlock 提出的一种计算$[0,1]$上均匀分布$X_i\sim p$的平均值的估计器。第一种估计器为$\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}\ne 0.5$，是有偏的，但方差下降速度为$O(n^{-2})$。

由此可见第二种估计器更加优秀。第二种估计器还有一个有用的性质是当采样点数量趋于无限时，误差会降为 0，这种估计器被称为是 consistent 的。

与方差密切相关的是均方误差（mean squared error, MSE），它被定义为估计器与真实值的平方差的期望：

$$\begin{array}{c}\text{MSE}[F]=E[(F-\int f(x)\mathrm{d}x{)}^2]\end{array}$$

对于无偏估计器，MSE 等于方差；反之是方差与估计器的偏差的平方之和。

使用 MSE 的概念，我们就可以通过样本均值（sample mean）来计算样本方差（sample variance），写作：

$$\begin{array}{c}\frac{1}{n-1}\sum _{i=1}^n(X_i-\bar{X}{)}^2\end{array}$$

除数为$n-1$而不是$n$的原因可以查阅 Bessel’s correction。

样本方差是基于给定的样本而不是真实概率分布的，因此它只是方差的估计值，本身就具有误差。例如对于一个 99.99% 取 1，0.01% 取 1e10 的随机变量，样本方差可能由于取得的样本值全是 1 而等于 0，但方差并不为 0。

如果一个积分的准确估计$\tilde{F}\approx \frac{1}{n}{\sum }_{i=1}^n(f(X_i)-\tilde{F}{)}^2$可以被计算出来（例如使用极大的采样数），那么 MSE 可以被估计为：

$$\begin{array}{c}\text{MSE}[F]\approx \frac{1}{n}\sum _{i=1}^n(f(X_i)-\tilde{F}{)}^2\end{array}$$

Sampling Strategies for Improving Efficiency

Stratified Sampling（分层采样）

分层采样的思想是以更好的方式确保采样点在域上的均匀性。

典型的分层采样的域是 s 维的单位 cube $[0,1{]}^s$，其他域可以通过合适的变换$T:[0,1{]}^s\to \Lambda$实现，通过选择不同的映射函数$T$，变换后的采样可以是不同的概率密度函数（若概率密度函数$p(x)$太过复杂，将难以找到与之对应的映射函数$T$），下文的分层采样所描述的情况使用 s 维的单位 cube 作为域，且概率均匀。

分层采样将积分域$\Lambda$细分成$n$个非重叠的区域${\Lambda }_1,{\Lambda }_2,...,{\Lambda }_n$，每个区域被称为一个 stratum，并且所有 stratum 的并集为原始域，有${\bigcup }_{i=1}^n{\Lambda }_i=\Lambda$。

假设总估计器$F^{\prime }$需要从${\Lambda }_i$中采样$n_i$个点，stratum${\Lambda }_i$的概率密度函数为$p_i$（不同于$p$，$p_i$定义在${\Lambda }_i$上，且在${\Lambda }_i$上的积分为 1），在每个 stratum 上使用蒙特卡洛估计可得：$F_i=\frac{1}{n_i}{\sum }_{i=1}^{n_i}\frac{f(X_{i,j})}{p_i(X_{i,j})}$，再通过$F^{\prime }={\sum }_{i=0}^n{v}_i{F}_i$计算得到总体估计，其中$v_i$是 stratum$i$的体积分数（fractional volume），取值为$(0,1]$。

Stratum $i$上被积函数的真值为：${\mu }_i=E([f(X_{i,j})])=\frac{1}{v_i}{\int }_{{\Lambda }_i}f(x)\mathrm{d}x$，方差为：${\sigma }_i^2=\frac{1}{v_i}{\int }_{{\Lambda }_i}(f(x)-{\mu }_i{)}^2\mathrm{d}x$，由上文线性递减证明可知，当在${\Lambda }_i$上采样$n_i$点时，$V[F_i]$将会趋近${\sigma }_i^2/n_i$。因此总估计器的方差为：

$$\begin{array}{rl}V[F^{\prime }]& =E[\sum v_i{F}_i]\\ & =\sum V[v_i{E}_i]\\ & =\sum v_i^{2}V[E_i]\\ & =\sum \frac{v_i^2{\sigma }_i^2}{n_i}\end{array}$$

我们合理假设${\Lambda }_i$上的采样数$n_i$与体积分数$v_i$成正比，有$n_i=v_{i}n$，则总估计器的方差为：

$$\begin{array}{c}V[F^{\prime }]=\frac{1}{n}\sum v_i{\sigma }_i^2\end{array}$$

将$V[F^{\prime }]$和不进行分层采样的估计器的方差$V[F]$进行比较，我们可以注意到，非分层采样点相当于根据$v_i$随机选取了一个 stratum${\Lambda }_I$， 再在${\Lambda }_I$上采样点$X$，在这种情况下，$X$是有条件地在$I$上选择的，因此可以使用条件概率表示方差：

$$\begin{array}{c}V[F]=\frac{1}{n}[\sum v_i{\sigma }_i^2+\sum v_i({\mu }_i-Q{)}^2]\end{array}$$

$Q$是整个$\Lambda$域上函数$f$的平均值（具体推导需要条件方差的性质，详见 Veach 1997）

从这个公式中我们可以看出，相比$V[F^{\prime }]$，$V[F]$多出了一项$\frac{1}{n}\sum v_i({\mu }_i-Q{)}^2$，同时，无论如何取 stratum，$V[F]$的值始终是不变的，因此增加$\frac{1}{n}\sum v_i({\mu }_i-Q{)}^2$就相当于在减小$V[F^{\prime }]$，所以分层采样的任务就在于增加不同域${\Lambda }_i$下${\mu }_i$与$Q$的差值。因此就算我们不知道函数$f$，也可以通过创建更加紧凑（compact）的 strate 来减小方差。

分层抽样的主要缺点是它和高斯型积分一样，在提高域的维度的时候，需要进行的分层数会指数上升，但幸运的是，独立地对某些维度进行分层是可行的，完成独立分层后，随机地将来自不同维度的样本关联起来即可得到最终结果，同时，在选择分层维度时，应该对那些对被积函数值的影响相关度最高的维度进行分层。

Importance Sampling（重要性采样）

在蒙特卡洛方法中，重要性采样是一种强大的用于减小方差的技术，它的基本思想在于，输入随机变量的某些值对估计参数的影响比其他值更大。如果通过更频繁的采样来强调这些“重要”的值，那么估计器的方差就可以减小，因此，重要性采样的基本方法是选择一个更偏向于重要值的分布。

假设估计器$F$需要估计积分$\int f(x)\mathrm{d}x$，采样使用的 PDF 为$p(x)$。若$p(x)\propto f(x)$，即$p(x)=cf(x)$（隐含了$c=\frac{1}{\int f(x)\mathrm{d}x}$），那么估计器的方差即为：

$$\begin{array}{c}V[F]=V[\frac{1}{n}\sum _{i=1}^n\frac{f(X_i)}{p(X_i)}]=V[\frac{1}{n}\cdot \sum _{i=1}^n\frac{1}{c}]=0\end{array}$$

方差降为了 0，但这是不可能的，毕竟我们想要拟合的目标即是$\int f(x)\mathrm{d}x=\frac{1}{c}$，但这给了我们一个思路，即我们需要尽可能找到与$f(x)$成比例的概率密度函数$p(x)$。

Multiple Importance Sampling（多重重要性采样）

我们经常会遇到两个或多个函数积的积分，例如$\int f_a(x)f_b(x)\mathrm{d}x$。在这种情况下，我们通常可以为单个函数单独推导出各自的重要性抽样策略，尽管不是与其乘积相似的策略。这种情况在涉及 light transport 的积分中尤其常见，例如在计算 light transport equation 时 BSDF，incident radiance 和余弦因子的乘积。

假设我们现在有采样分布$p_a$和$p_b$以及对应的采样$f_a$和$f_b$的方法，为使用蒙特卡洛积分，我们也许可以直接应用$p_a$作为 PDF，采样随机变量$X$，这时估计器的采样值为$\frac{f(X)}{p_a(X)}=\frac{f_a(X)f_b(X)}{p_a(X)}=c{f}_b(X)$，由于我们假设$p_a(X)\propto f_a(X)$，且$p_a$的积分为 1，所以$c$是$f_a(X)$的积分值，由这个方程可以看出，估计器的值将与$f_b$直接相关，估计器的方差会与$f_b$的方差成正比，可能导致方程非常大，同样的，使用$p_b$作为 PDF 也会导致相同的问题。

不幸的是，从每个分布中取一些样本并对两个估计量取平均值的解决方案并没有好到哪里去。因为方差是可加性的，一旦方差进入一个估计器，我们就不能通过将它添加到另一个低方差估计器中来消除它。

多重重要性采样(Multiple Importance Sampling, MIS)通过一种易于实现的方差减少技术解决了这个问题。基本思想是，当估计一个积分时，我们应该从多个抽样分布中抽取样本，希望至少有一个样本与被积函数的形状匹配得很好，即使我们不知道这是哪一个。多重重要性采样甚至鼓励专门的采样例程只考虑不寻常的特殊情况，因为它们减少了这些情况发生时的方差，通常成本相对较小。

我们现在有采样分布$p_a$和$p_b$，从这两个分布中各采一个样本，$X\sim p_a$，$Y\sim p_b$，MIS 蒙特卡洛积分器写作：

$$\begin{array}{c}{w}_a(X)\frac{f(X)}{p_a(X)}+w_b(Y)\frac{f(Y)}{p_b(Y)}\end{array}$$

$w_a,w_b$是确保估计器的期望值等于$f(x)$积分的权值函数。

更一般地，给定$n$个采样点，其中每个分布$p_i$采样$n_i$次，$X_{i,j}$是第 i 个分布下第 j 个采样点，MIS 积分器写作：

$$\begin{array}{c}{F}_n=\sum _{i=1}^n\frac{1}{n_i}\sum _{j=1}^{n_i}{w}_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})}\end{array}$$

使积分器无偏的方法是当$f(x)\ne 0$时，都有${\sum }_{i=1}^n{w}_i(x)=1$，并且$p_i(x)=0$时$w_i(x)=0$。

如果将各个权值函数设置为$w_i(x)=1/n$，MIS 就退化成了取平均值的解决方案，这并不能降低方差。我们更希望权值函数$w_i$在对应的$p_i$与 $f$匹配较好时取相对大的值，反之，取相对小的值，这样才能有效的降低方差。

在实践中，一种较好的方法被称为平衡启发式法（balance heuristic），它考虑了采样点$X_{i,j}$在所有采样分布下的概率，而单是在$p_i$下的概率。

$$\begin{array}{c}{w}_i(x)=\frac{n_i{p}_i(x)}{{\sum }_j{n}_j{p}_j(x)}\end{array}$$

各采样一个样本的(21)积分器可写作：

$$\begin{array}{c}\frac{f(X)}{p_a(X)+p_b(X)}+\frac{f(Y)}{p_b(Y)+p_b(Y)}\end{array}$$

这样一来，在$p_a$上采样的点$X$就算在$p_a$上的 PDF 较小，只要$p_b(X)$比较大，也能有效的降低这个点的权重，这种方法有效的降低了最终估计器的方差。

在实践中，指数启发式法（power heuristic）一般会更大程度的降低方差，使用指数$\beta$，有指数启发式赋权法：

$$\begin{array}{c}{w}_i(x)=\frac{(n_i{p}_i(x){)}^{\beta }}{{\sum }_j(n_j{p}_j(x){)}^{\beta }}\end{array}$$

即使不从所有分布中抽样，也可以应用 MIS。这种方法被称为单样本模型（single sample model）。对于被积函数$f(x)$，如果选取采样策略$p_i$的概率是$q_i$，从$p_i$采样得到样本$X$，单样本估计器写作：

$$\begin{array}{c}\frac{w_i(X)}{q_i}\frac{f(X)}{p_i(X)}\end{array}$$

这个估计器也是积分的无偏估计，对于单样本模型，可以证明平衡启发式算法是最优的。

多重重要性抽样的一个缺点是，如果其中一种采样策略$p_i$与被积函数$f$非常匹配，MIS 可能会略微增加方差。对于渲染应用来说，MIS 几乎总是值得的，因为它在可能存在高方差的情况下提供了方差减少。

MIS Compensation

在使用 MIS 时，若一块区域的被积函数$f$很小，但对应所有采样分布函数$p_i$都不为 0，则可能导致这块区域被过采样，同时$f$较大的区域会欠采样。MIS Compensation 可以在一定程度上解决这个问题，减小方差。

MIS Compensation 是基于锐化一个或多个(但不是全部)概率分布的思想——例如，通过调整它们，使它们在以前具有低概率的区域具有零概率。例如，一个新的抽样分布可以定义为：

$$\begin{array}{c}{p}^{\prime }(x)=\frac{\max (0,p(x)-\delta )}{\int \max (0,p(x)-\delta )\mathrm{d}x}\end{array}$$

这种技术特别容易应用于表化抽样分布（tabularized sampling distributions）的情况，它在采样环境贴图光源时效果很好。

Russian Roulette

俄罗斯轮盘赌方法跳过对最终结果贡献较小的样本，从而提高蒙特卡洛效率。在渲染中，我们经常使用形如$\frac{f(X)v(X)}{p(X)}$的估计器，其中，被积函数由一些容易评估的函数$f(X)$（eg. 表面如何散射光）和较为昂贵的函数$v(X)$（eg. 光线是否可见的二元函数）组成，在这种情况下，大部分计算费用都与$v$相关。

若$f(X)=0$，我们可以直接跳过计算$v(X)$，但若当$f(X)$较小时我们也跳过计算$v(X)$，这么做只会增加 bias，Russian Roulette 即用于解决这个问题，基本思路在于，跳过一些$v(X)$的计算，并增加不跳过时$v(X)$的权值。

为应用 Russian Roulette，我们先要根据采样点$X$计算出一个终止概率（termination probability）$q$，这个值可以通过很多方法获得，例如在上面的例子中计算一个与$f(X)$成反比的$q$。基于$q$，新估计器需要一个常数$c$（$c$一般等于 0），然后在 $1-q$概率下，估计器进行计算，并且附有权值$1/(1-q)$，反之跳过计算，直接使用$c$作为结果，新估计器为：

$$\begin{array}{c}{F}^{\prime }=\{\begin{array}{ll}\frac{F-qc}{1-q}& \xi >q\\ c& \text{otherwise}\end{array}\end{array}$$

新估计器的期望值为$E[F^{\prime }]=(1-q)(\frac{E[F]-qc}{1-q})+qc=E[F]$，是无偏估计法。

Russian Roulette 会提高估计器的方差，但更关键的是，它跳过一些低贡献的采样点的计算，从而提高了运行效率。

Splitting

Russian roulette 通过减少采样数来提高效率，而 Splitting 通过在多维积分上增加采样数来提高效率。考虑积分${\int }_A{\int }_{B}f(x,y)\mathrm{d}x\mathrm{d}y$，使用标准重要性采样估计器，我们需要在独立分布$X_i\sim p_x$，$Y_i\sim p_y$上独立采样$n$个点并计算$\frac{1}{n}{\sum }_{i=1}^n\frac{f(X_i,Y_i)}{p_x(X_i)p_y(Y_i)}$。而 Splitting 采取的策略是在一个$A$上的采样点上采样多个$B$上的值，对于每一个$X_i$采样点，都采样$m$个$Y_{i,j}$，从而将估计器改写成：

$$\begin{array}{c}\frac{1}{n}\sum _{i=1}^n\frac{1}{m}\sum _{j=1}^m\frac{f(X_i,Y_{i,j})}{p_x(X_i)p_y(Y_{i,j})}\end{array}$$

这样一来，计算$nm$个采样点就不需要分别采样$X,Y$各$nm$次，$X$只需要采样$n$次即可。例如在进行光追计算时，从 pixel 出发的光线打到一个半球表面后可以分裂成多条，相当于后续多条光线都共用了第一条光线的信息，从而提高了效率。

Sampling Using the Inversion Method

为了评估蒙特卡洛估计器，我们需要从给定的概率密度函数中随机抽取样本，有很多方法可以做到这一点，在 rendering 一种重要的方法被称为 Inversion Method。Inversion Method 使用反转 CDF 的方法，将$[0,1)$上的均匀样本映射到给定的一维概率分布上，当使用分布良好的样本时，Inversion Method 会非常高效。BSDFs，光源，相机，散射介质等的分布都可以使用 Inversion Method 实现。

Discrete Case

利用规范均匀随机变量，我们可以给出一种从一组离散概率中采样的算法。假设我们有一个有四种可能结果的 process，每种结果产生的概率为$p_1,p_2,p_3,p_4$，有$\sum p_i=1$，那么对应的 PMF 如下图所示：

离散概率的 CDF 的形式为$P_i={\sum }_{j=1}^i{p}_j$，以图像表示如下图所示：

我们可以在$P_{i-1}<\xi <P_i$时，选取$i-1$作为选择的结果，这种方式可以被解释成在反转 CDF，因此得名 Inversion Method。使用图形化解释，可以将$\xi$看做是在$[0,1]$的纵轴上进行随机采样，由于$\xi$击中任何特定柱块的概率正好等于该柱块的高度，因此结果是正确的。

Continuous Case

连续情况下的 Inversion Method 与离散情况下很相像，只是把 PMF 函数换成了 PDF，CDF 变成对 PDF 积分。

假设我们从 PDF$p(x)$中采样$X_i$，需要通过以下几步：

1. 积分 PDF 找到与之对应的 CDF，$P(x)={\int }_{-\infty }^{x}p(x^{\prime })\mathrm{d}{x}^{\prime }$
2. 得到规范均匀随机变量$\xi$
3. 生成符合$\xi =P(X)$的采样点，即$X=P^{-1}(\xi )$

Sampling a Linear Function

设被积函数$f(x)=(1-x)a+xb$定义在$[0,1]$上，且$a,b\ge 0$，PDF$p(x)\propto f(x)$，有$p(x)=\frac{2f(x)}{a+b}$，

积分 PDF 得到 CDF 函数为二次函数$P(x)=\frac{x(a(2-x)+bx)}{a+b}$，通过$\xi =P(X)$，反转得到$X=\frac{a-\sqrt{(1-\xi )a^2+\xi b^2}}{a-b}=\frac{\xi (a+b)}{a+\sqrt{(1-\xi )a^2+\xi b^2}}$

更多例子

Transforming between Distributions

使用 Inversion Method 可以通过$\xi$生成 1 维的分布函数，但一个更加一般的问题是如何从一个任意的分布函数转换成其他分布，理解这种转换，将使我们能够推导出多维采样算法。

假设我们有一个从$p(x)$采样的随机变量$X$，CDF 为$P(x)$。现在给定一个方程$f(x)$，有$f(x)=y$，考虑一个新的随机变量$Y=f(X)$。在这种情况下，$f(x)$必须是一个一对一的转换（如果有多个$x$值映射到同个$y$值，则不可能明确地描述这个$y$值的概率密度），同时，还需要保证$f(x)$的导数严格大于 0 或小于 0（严格递增或递减），这意味着对于任意给定的$x$值，都有 Pr { Y ≤ f ( x ) } = Pr { X ≤ x } \text{Pr}\{Y\leq{}f(x)\}=\text{Pr}\{X\leq{}x\} Pr{Y≤f(x)}=Pr{X≤x}，由 CDF 的定义可知：$P_f(y)=P_f(f(x))=P(x)$。

CDF 之间的关系直接导致它们 PDF 之间的关系，假设$f$是严格递增的，则有微分：

$$\begin{array}{c}{p}_f(y)\frac{\mathrm{d}f}{\mathrm{d}x}=p(x)\end{array}$$

或者写作：

$$\begin{array}{c}{p}_f(y)=(\frac{\mathrm{d}f}{\mathrm{d}x}{)}^{-1}p(x)\end{array}$$

更一般的，由于$f$的导数严格大于 0 或小于 0，因此有：

$$\begin{array}{c}{p}_f(y)=|\frac{\mathrm{d}f}{\mathrm{d}x}{|}^{-1}p(x)\end{array}$$

例如随机变量$X$的概率分布为$p(x)=2x,x\in [0,1]$，$f(x)=\sin x$，为计算随机变量$Y=f(X)$，我们先计算出$\mathrm{d}f/\mathrm{d}x=\cos x$，然后计算出$p_f(y)=\frac{p(x)}{|\cos x|}=\frac{2x}{\cos x}=\frac{2\arcsin y}{\sqrt{1-y^2}}$。

这个过程看起来是反向的，因为在实际应用中通常我们已知 PDF，进而从而计算转换。

假设我们用服从$p(x)$分布的随机变量$X$和服从$p_f(y)$分布的随机变量$Y$，如何计算转换方程，我们需要做的是利用 CDF 的相等，$P_f(y)=P(x)$计算出：

$$\begin{array}{c}f(x)=P_f^{-1}(P(x))\end{array}$$

这个方法是 inversion method 的推广。

Transformation in Multiple Dimensions

将一维的随机变量推广到多维，假设我们有一个 d 维的随机变量$X$服从概率分布$p(x)$，有随机变量$Y$服从$Y=T(X)$，$T$为双射函数，这种情况下，$X$和$Y$的概率密度函数的关系为：

$$\begin{array}{c}{p}_T(y)=p_T(T(x))=\frac{p(x)}{|J_T(x)|}\end{array}$$

其中$|J_T(x)|$是$T$的雅可比行列式的绝对值，$T$的雅可比矩阵为

$$J_T(x)=\left[\begin{array}{ccc}\frac{\partial T_1}{\partial x_1}& \cdots & \frac{\partial T_1}{\partial x_d}\\ ⋮& \ddots & ⋮\\ \frac{\partial T_d}{\partial x_1}& \cdots & \frac{\partial T_d}{\partial x_d}\end{array}\right]$$

例如我们想关联极坐标的采样$X=(r,\theta )$和笛卡尔坐标的采样$Y=(x,y)$，已知$x=r\cos \theta$，$y=r\sin \theta$和$X$的概率密度函数$p(r,\theta )$，计算$p(x,y)$。

我们先计算这个变换的雅可比矩阵$J_T$，有

$$J_T(x)=\left[\begin{array}{cc}\frac{\partial x}{\partial r}& \frac{\partial x}{\partial \theta }\\ \frac{\partial y}{\partial r}& \frac{\partial y}{\partial \theta }\end{array}\right]=\left[\begin{array}{cc}\cos \theta & -r\sin \theta \\ \sin \theta & r\cos \theta \end{array}\right]$$

$|J_T|=r$，因此$p(x,y)=p(r,\theta )/r$。

若想从笛卡尔坐标系上的采样转换到极坐标下的采样，有：

$$\begin{array}{c}p(r,\theta )=rp(x,y)\end{array}$$

在三维情况下通过方向坐标求球坐标，有$|J_T|=r^2\sin \theta$，对应的概率密度函数为：

$$\begin{array}{c}p(r,\theta ,\varphi )=r^2\sin \theta p(x,y,z)\end{array}$$

Sampling with Multidimensional Transformation

假设我们有一个二维密度函数$p(x,y)$并期望采样$(X,Y)$，如果密度是独立的，你们我们可以将其表示成一维密度的乘积，$p(x,y)=p_x(x)p_y(y)$。

更一般的，当考虑$(X,Y)$相关时，我们可以先将一个维度积分掉，从而计算出边际密度函数（marginal density function）$p(x)$，有：

$$\begin{array}{c}p(x)=\int p(x,y)\mathrm{d}y\end{array}$$

$p(x)$可以被看做是某个特定$x$在所有可能的$y$值上的平均密度，由此我们可以采样$X\sim p(x)$，然后使用条件概率密度函数采样$Y\sim p(y|x)$，有：

$$\begin{array}{c}p(y|x)=\frac{p(x,y)}{\int p(x,y)\mathrm{d}y}\end{array}$$

采样更高维的分别也可以使用同样的策略。

Reference

https://www.pbr-book.org/4ed/Monte_Carlo_Integration

https://www.math.arizona.edu/~tgk/mc/book_chap6.pdf

PS

本文主要是对PBRT内容的翻译总结，欢迎交流~