Monte Carlo 积分方法

重要性采样

Monte Carlo 积分方法是一种求解定积分的数值方法．假定需要计算的是下面的积分

$$I=\int_a^b f(x)dx.$$

先选取一个连续型随机变量 $X$，要求其概率密度函数 $p(x)$ 满足

$$\forall x\in[a,b],\quad p(x)>0.$$

然后令 $Y=g(X)=\frac{f(X)}{p(X)}$．注意到

$$\mathrm{E}\left[Y\right]=\mathrm{E}\left[\frac{f(X)}{p(X)}\right]=\int_{a}^{b} \frac{f\left(x\right)}{p\left(x\right)} p\left(x\right) d x=I,$$

我们可以将求解定积分 $I$ 的任务，直接转化为估计随机变量 $Y$ 的期望 $\mathrm{E}\left[Y\right]$.

估计 $\mathrm{E}\left[Y\right]$ 的一个常见做法，就是用 $Y$ 为总体进行简单随机抽样，然后以样本均值

$$\overline{Y}=\frac{1}{N}\sum_{i=1}^N Y_i$$

作为 $\mathrm{E}\left[Y\right]$ 的估计量．具体方案如下：

设 $X_i\sim p(x)(i=1,2,\cdots,n)$ 是独立同分布的 $n$ 个随机变量，则 $Y_i=\frac{f(X_i)}{p(X_i)}$ 也是独立同分布的 $n$ 个随机变量．令

$$I_N=\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}=\frac{1}{N}\sum_{i=1}^N Y_i,$$

以 $I_N$ 作为 $\mathrm{E}\left[Y\right]=I$ 的估计量．

统计学告诉我们：样本均值是总体期望的一致无偏估计量．这里我们可以再次证明一下这个结论．首先，$I_N$ 的期望为

$$\begin{aligned} \mathrm{E}\left[I_{N}\right]&=\mathrm{E}\left[\frac{1}{N}\sum_{i=1}^N Y_i\right]=\frac{1}{N}\mathrm{E}\left[\sum_{i=1}^N Y\right] = \mathrm{E}\left[Y\right] =I, \end{aligned}$$

这说明 $I_N$ 是无偏的．其次，当 $N\to\infty$ 时，$I_N$ 的方差

$$\begin{aligned} \mathrm{Var}\left[I_{N}\right]&=\mathrm{Var}\left[\frac{1}{N} \sum_{i=1}^{N} Y_i\right]\\ &=\frac{1}{N^2} \sum_{i=1}^{N} \mathrm{Var}\left[Y_i\right] \\ &=\frac{1}{N} \mathrm{Var}\left[Y\right] \to 0, \end{aligned}$$

这说明 $I_N$ 是一致的．证明完毕．

对于一个估计量，我们总希望减小它的方差．估计量 $I_{N}$ 的方差表达式为

$$\mathrm{Var}\left[I_{N}\right]=\frac{1}{N} \mathrm{Var}\left[\frac{f(X)}{p(X)}\right].$$

我们清楚的是：当 $p(x)=kf(x)$ 时，也即 $p(x)=\frac{1}{I}f(x)$ 时，方差会取得最小值 $0$．但是，求出 $I$ 是不可能的，因为这正是我们 Monte Carlo 积分方法的求解目标．我们所能做的，只是让 $p(x)$ 的形态尽量与 $f(x)$ 相接近．这种削减方差的手段叫作 重要性采样（importance sampling）．$f(x)$ 取值越大的区域，对积分的贡献会越多，也就越重要，需要更多的采样．

如果我们简单地让 $X$ 服从 $[a,b]$ 上的均匀分布．此时估计量为

$$I_N=\frac{b-a}{N}\sum_{i=1}^Nf(X_i).$$

直觉上看，这就是将积分区间等分成 $N$ 段，用 $N$ 个面积为 $\frac{b-a}{N}\cdot f(X_i)$ 的矩形来估计曲边梯形．如果我们希望在某个子区间采样数翻倍，那么这个子区间需要分成两半，将一个矩形换成宽度为 $\frac{b-a}{2N}$ 的两个矩形．推而广之就是采样概率越大的区域，单个样本的权重就越小．这也是为什么 $I_{N}$ 中要除以 $p(X_i)$ 的原因．

多重重要性采样

假定需要计算的是下面的积分

$$I=\int_a^b f(x)g(x)dx,$$

并且已知概率密度函数 $p_X(x)$ 和 $p_Y(x)$ 分别能较好地吻合 $f(x)$ 和 $g(x)$ 的形状，即它们分别能对

$$\int_a^b f(x)dx,\,\int_a^b g(x)dx$$

进行较好地重要性采样，那么该如何对 $I$ 进行重要性采样呢？

很容易想到的一种方法是用 $p_X(x)p_Y(x)$ 进行重要性采样，但是，有些时候为总体 $Z\sim p_X(x)p_Y(x)$ 生成样本 $Z_i$ 比较困难．这时，多重重要性采样（multiple importance sampling）不失为一种好的选择．

设 $X_i\sim p_X(x)$，$Y_i\sim p_Y(x)$，多重重要性采样使用了如下估计量

$$I_{n_{f},n_g}=\frac{1}{n_{f}} \sum_{i=1}^{n_{f}} \frac{f\left(X_{i}\right) g\left(X_{i}\right) w_{f}\left(X_{i}\right)}{p_{X}\left(X_{i}\right)}+\frac{1}{n_{g}} \sum_{j=1}^{n_{g}} \frac{f\left(Y_{j}\right) g\left(Y_{j}\right) w_{g}\left(Y_{j}\right)}{p_{Y}\left(Y_{j}\right)},$$

其中 $w_f(x),w_g(x)$ 是满足 $w_f(x)+w_g(x)=1$ 的权重函数．可以验证

$$\begin{aligned} \mathrm{E}\left[I_{n_{f},n_g}\right]&= \mathrm{E}\left[\frac{f\left(X\right) g\left(X\right) w_{f}\left(X\right)}{p_{X}\left(X\right)}\right]+\mathrm{E}\left[\frac{f\left(Y\right) g\left(Y\right) w_{g}\left(Y\right)}{p_{Y}\left(Y\right)}\right]\\ &=\int_a^b\frac{f\left(x\right) g\left(x\right) w_{f}\left(x\right)}{p_X(x)}p_X(x)dx+\int_a^b\frac{f\left(x\right) g\left(x\right) w_{g}\left(x\right)}{p_Y(x)}p_Y(x)dx\\ &=\int_a^bf\left(x\right) g\left(x\right) w_{f}\left(x\right)dx+\int_a^bf\left(x\right) g\left(x\right) w_{g}\left(x\right)dx\\ &=\int_a^bf\left(x\right) g\left(x\right) (w_{f}\left(x\right)+w_{g}\left(x\right))\:dx\\ &=\int_a^bf\left(x\right) g\left(x\right)dx, \end{aligned}$$

即 $I_{n_{f},n_{g}}$ 确实是 $I$ 的无偏估计量．实践中，常使用如下形式的启发式权重函数

$$w_{f}(x)=\frac{\left(n_{f} p_{X}(x)\right)^{\beta}}{\left(n_{f} p_{X}(x)\right)^{\beta}+\left(n_{g} p_{Y}(x)\right)^{\beta}},\quad w_{g}(x)=\frac{\left(n_{g} p_{Y}(x)\right)^{\beta}}{\left(n_{f} p_{X}(x)\right)^{\beta}+\left(n_{g} p_{Y}(x)\right)^{\beta}}.$$

当 $\beta=2$ 时，往往效果会较好．

光线传输方程

基本形式

之前我们已经见过了反射方程

$$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

如果我们将描述表面反射的双向反射分布函数 $f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)$ 替换为双向散射分布函数 $f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)$，并加入自发光辐照率 $L_e(\mathbf{p},\boldsymbol{\omega}_o)$，那么就得到了如下含自发光项的散射方程

$$L_o(\mathbf{p},\boldsymbol{\omega}_o)=L_e(\mathbf{p},\boldsymbol{\omega}_o)+\int_{S^2}f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\left|\cos\theta_i\right|d\Omega_i.$$

注意到在上式中，$L_e(\mathbf{p},\boldsymbol{\omega}_o)$ 和 $f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)$ 都是已知量．如果我们能将函数 $L_i$ 用函数 $L_o$ 表示出来，那么我们将得到一个只含 $L_o$ 的积分方程．为了实现这一点，我们需要引入 光线投射函数（ray-casting function）$t(\mathbf{p},\boldsymbol{\omega})$．假设在 $\mathbf{p}$ 点向 $\boldsymbol{\omega}$ 方向发出一条射线，射线第一次与场景中某个表面相交的位置就记作 $t(\mathbf{p},\boldsymbol{\omega})$．于是，我们得到

$$L_i(\mathbf{p},\boldsymbol{\omega}_i)=L_o(t(\mathbf{p},\boldsymbol{\omega}_i),-\boldsymbol{\omega}_i).$$

为了处理射线不与场景中任何表面相交的特殊情况，我们规定此时 $t(\mathbf{p},\boldsymbol{\omega})$ 取特殊值 $\Lambda$，并且对任意 $\boldsymbol{\omega}$，有 $L_o(\Lambda,\boldsymbol{\omega})=0$．

利用上面的关系式，并将 $L_o(\mathbf{p},\boldsymbol{\omega}_o)$ 简记为 $L(\mathbf{p},\boldsymbol{\omega}_o)$，由此导出的方程就称为 光线传输方程（light transport equation）

$$L(\mathbf{p},\boldsymbol{\omega}_o)=L_e(\mathbf{p},\boldsymbol{\omega}_o)+\int_{S^2}f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L(t(\mathbf{p},\boldsymbol{\omega}_i),-\boldsymbol{\omega}_i)\left|\cos\theta_i\right|d\Omega_i.$$

光线传输方程也叫做称渲染方程（rendering equation），它是关于函数 $L(\mathbf{p},\boldsymbol{\omega}_o)$ 的积分方程，描述了均衡状态下物体表面任意一点、沿任意方向的辐照率．实际上，我们可以也可以得到场景中任意一点、沿任意方向的辐照率，因为借助前面导出的关系式

$$L(\mathbf{p},\boldsymbol{\omega})=L(t(\mathbf{p},\boldsymbol{\omega}),-\boldsymbol{\omega}),$$

我们可以将场景中的辐照率转化为物体表面的辐照率．

表面形式

在上面的光线传输方程中，场景中物体表面的几何信息全部隐藏在了光线投射函数 $t(\mathbf{p},\boldsymbol{\omega})$ 里．为了让这些信息显式地表达出来，我们将导出表面形式的光线传输方程．

之前在对某点各个方向上的入射光做积分时，我们是在单位球上做积分．考虑到某点接收的入射光实际上来自于其他物体表面各点的出射光，我们也可以把单位球上的积分转化为其他物体表面上的积分．如果我们要用蒙特卡洛积分法估计该积分，那么我们的采样区域将从单位球转到其他物体的表面．

下面引入一些记号．如果从 $\mathbf{p}^\prime$ 点沿 $\boldsymbol{\omega}$ 方向出发的射线经过了 $\mathbf{p}$ 点，我们记

$$L(\mathbf{p}^\prime\to\mathbf{p}):=L(\mathbf{p}^\prime,\boldsymbol{\omega})=L\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert}\right).$$

如下图所示，

图1 - $f(\mathbf{p}^{\prime\prime}\to\mathbf{p}^\prime\to\mathbf{p})$ 示意图

如果 $t(\mathbf{p}^\prime,\boldsymbol{\omega}_o)=\mathbf{p}$，$t(\mathbf{p}^\prime,\boldsymbol{\omega}_i)=\mathbf{p}^{\prime\prime}$，我们记

$$f(\mathbf{p}^{\prime\prime}\to\mathbf{p}^\prime\to\mathbf{p}):=f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)=f\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert},\frac{\mathbf{p}^{\prime\prime}-\mathbf{p}^\prime}{\Vert\mathbf{p}^{\prime\prime}-\mathbf{p}^\prime\Vert}\right).$$

设 $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ 是可视性函数，当 $\mathbf{p}$ 与 $\mathbf{p}^{\prime}$ 两点之间没有遮挡、相互可见时，即

$$t\left(\mathbf{p},\frac{\mathbf{p}^\prime-\mathbf{p}}{\Vert\mathbf{p}^\prime-\mathbf{p}\Vert}\right)=\mathbf{p}^\prime,\quad t\left(\mathbf{p}^\prime,\frac{\mathbf{p}-\mathbf{p}^\prime}{\Vert\mathbf{p}-\mathbf{p}^\prime\Vert}\right)=\mathbf{p}$$

时，$V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ 取值为 $1$；否则 $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ 取值为 $0$．

设几何耦合项（geometric coupling term）为可视性函数与积分区域转换映射的 Jacobi 行列式的乘积，即

$$G\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)=V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right) \frac{|\cos \theta|\left|\cos \theta^{\prime}\right|}{\left\|\mathbf{p}-\mathbf{p}^{\prime}\right\|^{2}},$$

其中 $\theta$ 为 $\mathbf{p}$ 点处法向量与 $\mathbf{p}^\prime-\mathbf{p}$ 的夹角， $\theta^\prime$ 为 $\mathbf{p}^\prime$ 点处法向量与 $\mathbf{p}-\mathbf{p}^\prime$ 的夹角．

当 $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)=1$ 时，光线传输方程可以重写成场景中所有物体表面 $A$ 上的积分

$$L\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right)=L_{\mathrm{e}}\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right)+\int_{A} f\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime} \rightarrow \mathbf{p}\right) L\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime}\right) G\left(\mathbf{p}^{\prime \prime} \leftrightarrow \mathbf{p}^{\prime}\right)\: d A\left(\mathbf{p}^{\prime \prime}\right).$$

这就是表面形式的光线传输方程．

路径积分形式

表面形式的光线传输方程给出了 $L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)$ 关于 $i$ 的递推关系

$$\begin{aligned} &L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)\\ =\;&L_{\mathrm{e}}\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)+\int_{A} f\left(\mathbf{p}_{i+2} \rightarrow \mathbf{p}_{i+1}\rightarrow \mathbf{p}_i\right) L\left(\mathbf{p}_{i+2}\rightarrow \mathbf{p}_{i+1}\right) G\left(\mathbf{p}_{i+2}\leftrightarrow \mathbf{p}_{i+1}\right)\: d A\left(\mathbf{p}_{i+2}\right). \end{aligned}$$

从 $i=0$ 出发，先将等式右边的 $L\left(\mathbf{p}_{2}\rightarrow\mathbf{p}_{1}\right)$ 用 $L\left(\mathbf{p}_{3}\rightarrow \mathbf{p}_{2}\right)$ 表示，再将先将等式右边的 $L\left(\mathbf{p}_{3}\rightarrow \mathbf{p}_{2}\right)$ 用 $L\left(\mathbf{p}_{4}\rightarrow \mathbf{p}_{3}\right)$ 表示，以此类推，最终将得到如下等式

$$\begin{aligned} L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)&=L_{\mathrm{e}} \left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) \\ &\quad+\int_{A} L_{\mathrm{e}}\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1}\right) f\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) G\left(\mathbf{p}_{2} \leftrightarrow \mathbf{p}_{1}\right)\: d A\left(\mathbf{p}_{2}\right) \\ &\quad+\int_{A} \int_{A} L_{\mathrm{e}}\left(\mathbf{p}_{3} \rightarrow \mathbf{p}_{2}\right) f\left(\mathbf{p}_{3} \rightarrow \mathbf{p}_{2} \rightarrow \mathbf{p}_{1}\right) G\left(\mathbf{p}_{3} \leftrightarrow \mathbf{p}_{2}\right)\\ &\quad\quad\times f\left(\mathbf{p}_{2} \rightarrow \mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right) G\left(\mathbf{p}_{2} \leftrightarrow \mathbf{p}_{1}\right) \:d A\left(\mathbf{p}_{3}\right) d A\left(\mathbf{p}_{2}\right)\\ &\quad+\cdots \end{aligned}$$

记

$$T\left(\overline{\mathbf{p}}_{n}\right):=\prod_{i=1}^{n-1} f\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_{i} \rightarrow \mathbf{p}_{i-1}\right) G\left(\mathbf{p}_{i+1} \leftrightarrow \mathbf{p}_{i}\right)$$

为路径 $\overline{\mathbf{p}}_{n}=(\mathbf{p}_n,\mathbf{p}_{n-1},\cdots,\mathbf{p}_{0})$ 的吞吐量（throughput），又令

我们可以将 $L\left(\mathbf{p}_{1}\to\mathbf{p}_{0}\right)$ 写成级数形式

$$L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right).$$

这被称作路径积分形式的光线传输方程．注意到等式右端不再包含 $L$，因此，我们实际上以级数形式显式地写出了光线传输方程的解．$P\left(\overline{\mathbf{p}}_{1}\right)$ 表示从光源直接射向照相机的光，$P\left(\overline{\mathbf{p}}_{2}\right)$ 表示经过一次散射的直接光，$P\left(\overline{\mathbf{p}}_{3}\right),P\left(\overline{\mathbf{p}}_{4}\right),\cdots$ 表示经过多次散射的间接光．

路径追踪算法

路径追踪算法（path tracing algorithm）的目标是估计

$$L\left(\mathbf{p}_{1} \rightarrow \mathbf{p}_{0}\right)=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right).$$

的值．需要解决的问题有两个：第一个是如何估计无穷求和，第二个是如何估计积分项 $P\left(\overline{\mathbf{p}}_{n}\right)$．

俄罗斯轮盘赌

无穷求和的困难可以由所谓的俄罗斯轮盘赌（Russian roulette）解决．

假设我们有办法得到任意一个 $P\left(\overline{\mathbf{p}}_{n}\right)$ 的无偏估计量 $\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)$，那么我们可以立刻写出 $\sum\limits_{n=1}^{\infty} {P} \left(\overline{\mathbf{p}}_{n}\right)$ 的一个无偏估计量

$$\sum_{n=1}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right).$$

但问题是在现实中，我们只能获得有限项的 $\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)$．俄罗斯轮盘赌的想法，就是以一定概率 $q$ 丢弃第 $N$ 项之后的所有项之和 $\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)$．而在没有被丢弃的情况下，放大 $\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)$ 作为补偿．

对于前 $2$ 项 $P\left(\overline{\mathbf{p}}_{1}\right),P\left(\overline{\mathbf{p}}_{2}\right)$，我们并不采用俄罗斯轮盘赌，而是将它们直接计算出来．从第 $3$ 项开始，我们执行俄罗斯轮盘赌．设 $\eta_1$ 服从 0-1 分布，

$$\mathrm{P}(\eta_1=k) = q^{k}(1-q)^{1-k},\quad k=0,1,$$

当 $\eta_1$ 取 $1$ 时，我们丢弃之后的项，只取前两项之和

$$\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right),$$

而当 $\eta_1$ 取 $0$ 时，我们将后面的项乘上一个缩放因子 $\frac{1}{1-q}$，即得

$$\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1}{1-q}\sum_{n=3}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right).$$

结合这两种情况，我们的统计量应该写成

$$\widehat{Q}_3=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\sum_{n=3}^{\infty}\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)$$

注意到

$$\mathrm{E}\left[\frac{1_{\eta_1=0}}{1-q}\right]=\frac{\mathrm{P}(\eta_1=0)}{1-q}=\frac{1-q}{1-q}=1,$$

我们有

$$\mathrm{E}\left[ \sum_{n=1}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right]=\mathrm{E}\left[\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\sum_{n=3}^{\infty}\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right],$$

即 $\widehat{Q}_3$ 也是 $\sum\limits_{n=1}^{\infty}{P}\left(\overline{\mathbf{p}}_{n}\right)$ 的无偏估计量．

设 $\eta_2$ 与 $\eta_1$ 独立同分布，我们还可以用同样的方法继续构造出

$$\widehat{Q}_4=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_2=0}}{1-q}\sum_{n=4}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right).$$

显然，$\mathrm{E}[\widehat{Q}_4]=\mathrm{E}[\widehat{Q}_3]$，故 $\widehat{Q}_4$ 也是一个无偏估计量．

一般地，若 $\eta_1,\eta_2,\cdots$ 独立且都服从成功率为 $q$ 的 0-1 分布，则

$$\begin{aligned} \widehat{Q}_{N+2}&=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_2=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{4}\right)+\cdots+\frac{1_{\eta_N=0}}{1-q}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right)\right)\\ &=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\frac{1_{\eta_1=0}}{1-q}\widehat{P}\left(\overline{\mathbf{p}}_{3}\right)+\frac{1_{\eta_1,\eta_2=0}}{(1-q)^2}\widehat{P}\left(\overline{\mathbf{p}}_{4}\right)+\cdots+\frac{1_{\eta_1,\eta_2,\cdots,\eta_N=0}}{(1-q)^N}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) \end{aligned}$$

是 $\sum\limits_{n=1}^{\infty}{P}\left(\overline{\mathbf{p}}_{n}\right)$ 的无偏估计量．令 $\widehat{Q}=\lim\limits_{N\to\infty}\widehat{Q}_{N}$，则 $\widehat{Q}$ 也是无偏估计量．

设

$$\tau=\inf\{\:i\mid\eta_i=1\:\},$$

即第 $\tau+2$ 项及其后的所有项都被丢弃，我们可以将 $\widehat{Q}$ 写作

$$\widehat{Q}=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\tau+1}\frac{1}{(1-q)^{n-2}} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)=\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)1_{\tau>n-2}.$$

不难求出，$\tau$ 服从几何分布，其概率质量函数为

$$\mathrm{P}(\tau=k)=(1-q)^{k-1} q,\quad k=1,2, \cdots.$$

于是，我们可以用计算验证 $\widehat{Q}$ 的无偏性

$$\begin{aligned} \mathrm{E}\left[\widehat{Q}\right]&=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\mathrm{P}(\tau>n-2)\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\sum_{k=n-1}^{\infty}(1-q)^{k-1}q\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)(1-q)^{n-2}\\ &=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right). \end{aligned}$$

注意到 $\mathrm{P}(\tau<\infty)=1$，即 $\widehat{Q}$ 为有限项的概率为 $1$，我们完全可以在现实中求出 $\widehat{Q}$ 的估计值．

伪代码

路径追踪算法的伪代码如下：

shade(p, wo):
    # Contribution from the light source.
    uniformly sample the light at x’ (pdf_light = 1 / A);
    shoot a ray from p to x’;
    L_dir = 0.0;
    if the ray is not blocked in the middle:
        L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light;

    # Contribution from other reflectors.
    L_indir = 0.0;
    test Russian Roulette with probability P_RR;
    if passing Russian Roulette test:
        uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi);
        trace a ray r(p, wi);
        if ray r hit a non-emitting object at q:
            L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR;
    
    return L_dir + L_indir;