# C13. 蒙特卡洛积分 Monte Carlo Integeration
[TOC]
许多积分方程没有解析解，需要采用数值方法。尽管标准的数值积分方法如梯形积分和高斯求积能高效求解低维平滑积分，他们的收敛速度在高维和不连续积分时会变得糟糕。
蒙特卡罗数值积分法是解决这一问题的一种方法。他们利用随机性来计算积分，其收敛速度与被积函数的维数无关。
合理地使用随机性已经彻底改变了算法设计领域。随机算法可大致地分为两类
- Las Vegas：使用随机性，但最终结果相同。
- Monte Carlo：依赖于随机数，给出不同的结果，但平均上结果正确。
蒙特卡洛的一个非常有用的性质是，只要有能力计算被积函数 $f(x)$ 在域内任意点的值，就能估计积分 $\int f(x)\mathrm{d}x$。
渲染中遇到的许多积分很难或者不可能直接计算。而这些问题可以用蒙特卡洛积分来解决。
蒙特卡洛的主要缺点：如果用 n 个样本来估计积分，算法收敛到正确结果的速率是 $O(n^{-1/2})$。换句话说，如果要减少一半的误差，需要 4 倍的样本。
## 13.1 背景与概率回顾 Background and Probability Review
随机变量 $X$ 是一个从某个随机过程选择的值。一个函数作用在随机变量 $X$ 上得到新的随机变量 $Y=f(X)$。
采样离散随机变量 $X$，其满足 $p_i=P(X=X_i)$，我们可以用一个连续均匀分布的随机变量 $\xi\in[0,1)$，然后将其映射到离散随机变量上，选择 $X_i$ 如果
$$
\sum _ { j = 1 } ^ { i - 1 } p _ { j } < \xi \leq \sum _ { j = 1 } ^ { i } p _ { j }
$$
对于光照任务，概率基于功率，为
$$
p _ { i } = \frac { \Phi _ { i } } { \sum _ { j } \Phi _ { j } }
$$
随机变量的累积分布函数 cumulative distribution function (CDF) $P(x)$ 定义如下
$$
P ( x ) = \operatorname { Pr } \{ X \leq x \}
$$
### 13.1.1 连续随机变量 Continuous Random Variables
一个特别重要的随机变量是标准均匀随机变量 canonical uniform random variable，记为 $\xi$，在 $[0,1)$ 上等概率取值。
重要性体现在两方面
- 容易从软件层面实现这个随机变量
- 可以用该随机变量和适当的变换来生成任意分布的样本
概率分布函数 probability density function PDF 描述随机变量取特定值的相对概率，定义为
$$
p ( x ) = \frac { \mathrm { d } P ( x ) } { \mathrm { d } x }
$$
对于 $\xi$，我们有
$$
p ( x ) = \left\{ \begin{array} { l l } { 1 } & { x \in [ 0,1 ) } \\ { 0 } & { \text { otherwise } } \end{array} \right.
$$
PDFs 要求非负和定义域上积分值为 1。
区间概率
$$
P ( x \in [ a , b ] ) = \int _ { a } ^ { b } p ( x ) \mathrm { d } x
$$
### 13.1.2 期望和方差 Expected Values and Variance
函数 $f$ 的期望值 $E_p[f (x)]$ 定义为函数在其定义域上的某个值的分布 $p(x)$ 上的平均值。域 D 上的期望值定义为
$$
E _ { p } [ f ( x ) ] = \int _ { D } f ( x ) p ( x ) \mathrm { d } x
$$
函数的方差是函数与期望值的平方差的期望，即
$$
V [ f ( x ) ] = E \left[ ( f ( x ) - E [ f ( x ) ] ) ^ { 2 } \right]
$$
期望满足线性性
$$
\begin{aligned} E [ a f ( x ) ] & = a E [ f ( x ) ] \\ E \left[ \sum _ { i } f \left( X _ { i } \right) \right] & = \sum _ { i } E \left[ f \left( X _ { i } \right) \right] \end{aligned}
$$
方差满足
$$
\begin{aligned}
V [ a f ( x ) ] &= a ^ { 2 } V [ f ( x ) ]\\
\sum _ { i } V \left[ f \left( X _ { i } \right) \right] &= V \left[ \sum _ { i } f \left( X _ { i } \right) \right]\\
V [ f ( x ) ] &= E \left[ ( f ( x ) ) ^ { 2 } \right] - E [ f ( x ) ] ^ { 2 }\\
\end{aligned}
$$
## 13.2 蒙特卡洛估计 the Monte Carlo Estimator
假设我们要计算 1D 积分 $\int _ { a } ^ { b } f ( x ) \mathrm { d } x$。给定一组均匀随机变量 $X_i\in[a,b]$，蒙特卡洛估计表明估计式
$$
F _ { N } = \frac { b - a } { N } \sum _ { i = 1 } ^ { N } f \left( X _ { i } \right)
$$
的期望值 $E[F_N]$ 等于积分值。
> **推导** 
>
> 已知随机变量 $X_i$ 的 PDF $p(x)=1/(b-a)$，可推得
> $$
> \begin{aligned} E \left[ F _ { N } \right] & = E \left[ \frac { b - a } { N } \sum _ { i = 1 } ^ { N } f \left( X _ { i } \right) \right] \\ & = \frac { b - a } { N } \sum _ { i = 1 } ^ { N } E \left[ f \left( X _ { i } \right) \right] \\ & = \frac { b - a } { N } \sum _ { i = 1 } ^ { N } \int _ { a } ^ { b } f ( 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{aligned}
> $$
对均匀随机变量的限制可以用一个小的泛化来放宽。这是非常重要的一步，因为在蒙特卡罗中，仔细选择抽取样本的PDF是减少方差的一项重要技术。
如果随机变量 $X_i$ 服从某任意 PDF $p(x)$，则估计式
$$
F _ { N } = \frac { 1 } { N } \sum _ { i = 1 } ^ { N } \frac { f \left( X _ { i } \right) } { p \left( X _ { i } \right) }
$$
的期望 $E \left[ F _ { N } \right]$ 是积分值，要求当 $f(x)\ne 0$ 时 $p(x)>0$。
> **推导** 
> $$
> \begin{aligned} E \left[ F _ { N } \right] & = E \left[ \frac { 1 } { N } \sum _ { i = 1 } ^ { N } \frac { f \left( X _ { i } \right) } { p \left( X _ { i } \right) } \right] \\ & = \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{aligned}
> $$
## 13.3 采样随机变量 Samping Random Variables
### 13.3.1 反演方法 the Inversion Method
> 正确的叫法应该是[逆变换采样 inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling)
反演方法使用一个或多个均匀随机变量，并将它们映射成服从期望分布的随机变量。
采样服从 PDF $p(x)$ 的 $X_i$ 的步骤
1. 计算 CDF $P(x)=\int_{-\infty}^x p(x')dx'$ 
2. 计算逆函数 $p^{-1}(x)$ 
3. 获取均匀分布的随机数 $\xi$ 
4. 计算 $X_i=P^{-1}(\xi)$ 
> **证明** 
>
> 我们有均匀分布的随机变量 $U\sim \text{Unif}[0,1]​$，还有一个严格单调递增的概率累计函数 CDF $F_X​$。
>
> 我们想要用 $U​$ 来生成 $X​$，期望用某严格单调变换 $T:[0,1]\mapsto \mathbb{R}​$ 来实现，满足 $T(U)\stackrel{d}{=}X​$。
>
> 我们想要求出 T，如下推导
> $$
> \begin{aligned}
> F_X(x) &= Pr(X\le x)\\
> &= Pr(T(U)\le x)\\
> &= Pr(U\le T^{-1}(x))\\
> &= T^{-1}(x)\\
> \end{aligned}
> $$
> 因此有
> $$
> T(u)=F_X^{-1}(u)
> $$
> 搞定
#### 示例：分段常数 1D 函数 Example: Piecewise-Constant 1D Functions
![1557142569174](assets/1557142569174.png)
每段长度 $\Delta=1/N$。积分值为
$$
c = \int _ { 0 } ^ { 1 } f ( x ) \mathrm { d } x = \sum _ { i = 0 } ^ { N - 1 } \Delta v _ { i } = \sum _ { i = 0 } ^ { N - 1 } \frac { v _ { i } } { N }
$$
这样 PDF $p(x)=f(x)/c$。因此 CDF $P(x)$ 为
$$
\begin{aligned}
P \left( x _ { 0 } \right) & = 0 \\
P \left( x _ { 1 } \right) & = \int _ { x _ { 0 } } ^ { x _ { 1 } } p ( x ) \mathrm { d } x = \frac { v _ { 0 } } { N c } = P \left( x _ { 0 } \right) + \frac { v _ { 0 } } { N c } \\
P \left( x _ { 2 } \right) & = \int _ { x _ { 0 } } ^ { x _ { 2 } } p ( x ) \mathrm { d } x = \int _ { x _ { 0 } } ^ { x _ { 1 } } p ( x ) \mathrm { d } x + \int _ { x _ { 1 } } ^ { x _ { 2 } } p ( x ) \mathrm { d } x = P \left( x _ { 1 } \right) + \frac { v _ { 1 } } { N c } \\
P \left( x _ { i } \right) & = P \left( x _ { i - 1 } \right) + \frac { v _ { i - 1 } } { N c }
\end{aligned}
$$
在两点 $x_i$ 和 $x_{i+1}$ 间，CDF 以斜率 $v_i/c$ 线性增加。
![1557142777793](assets/1557142777793.png)
由此我们可以利用反演方法去采样。
### 13.3.2 拒绝方法 the Rejection Method
对于有些函数 $f(x)$，难以积分它们来得到 PDFs，或者难以求逆 CDFs。拒绝方法只要求我们知道如何采样服从 PDF $p(x)$ 的随机变量，以及要求存在 $c$，使得 $f(x)<c\ p(x)$。算法如下
![1557133611082](assets/1557133611082.png)
拒绝方法并没有用在 `pbrt` 的蒙特卡洛算法里，一般更喜欢采样直接采样的方法。但这也是一个值得注意的重要方法。
## *13.4 大都市采样 Metropolis Sampling
## 13.5 分布变换 Transforming between Distributions
给定一个服从 PDF $p_x(x)$ 的随机变量 $X_i$，我们想知道随机变量 $Y_i=y(X_i)$ 的分布。函数 $y(x)$ 必须是单射的。即 y 的导数必须严格大于 0（或严格小于 0），意味着
$$
Pr\{Y\le y(x)\}=Pr\{X\le x\}
$$
即
$$
P_y(y)=P_y(y(x))=P_x(x)
$$
由 PDF 的定义可得微分关系为
$$
p_y(y)=\left|\frac{\mathrm{d}y}{\mathrm{d}x}\right|^{-1}p_x(x)
$$
已知 $X$ 服从 PDF $p_x(x)$，$Y$ 服从 PDF $p_y(y)$，要求累计分布函数相等，即 $P_y(y) = P_x(x)$，则
$$
y(x)=P_y^{-1}(P_x(x))
$$
如果 $X$ 是 $[0,1)$ 上均匀分布的随机变量，则 $P_x(x)=x$，这样我们就得到了 13.3.1 节给出的步骤。
### 13.5.1 多维变换 Transformation in Multiple Dimensions
假设我们有一个 n 维随机变量 $X$，其概率密度函数为 $p_x(x)$，令 $Y=T(X)$，且 $T$ 是一个双射，这样有
$$
p _ { y } ( y ) = p _ { y } ( T ( x ) ) = \frac { p _ { x } ( x ) } { \left| J _ { T } ( x ) \right| }
$$
其中 $T ( x ) = \left( T _ { 1 } ( x ) , \ldots , T _ { n } ( x ) \right)$，$|J_T|$ 是 $T$ 的雅克比矩阵的行列式的绝对值，雅克比矩阵如下
$$
\left( \begin{array} { c c c } { \partial T _ { 1 } / \partial x _ { 1 } } & { \cdots } & { \partial T _ { 1 } / \partial x _ { n } } \\ { \vdots } & { \ddots } & { \vdots } \\ { \partial T _ { n } / \partial x _ { 1 } } & { \cdots } & { \partial T _ { n } / \partial x _ { n } } \end{array} \right)
$$
### 13.5.2 极坐标 Polar Coordinates
极坐标变换为
$$
\begin{array} { l } { x = r \cos \theta } \\ { y = r \sin \theta } \end{array}
$$
假设我们从 $p(r,\theta)$ 中采样，想知道 $p(x,y)$，雅克比矩阵为
$$
J _ { T } = \left( \begin{array} { c c } { \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} { c c } { \cos \theta } & { - r \sin \theta } \\ { \sin \theta } & { r \cos \theta } \end{array} \right)
$$
其行列式为 $r \left( \cos ^ { 2 } \theta + \sin ^ { 2 } \theta \right) = r$，因此 $p ( x , y ) = p ( r , \theta ) / r$。当然，这个与我们期望的相反了——通常我们从笛卡尔坐标采样然后变换它到极坐标中。这样我们有
$$
p ( r , \theta ) = r\ p ( x , y )
$$
### 13.5.3 球坐标 Spherical Coordinates
球坐标变换为
$$
\begin{aligned} x & = r \sin \theta \cos \phi \\ y & = r \sin \theta \sin \phi \\ z & = r \cos \theta \end{aligned}
$$
可求得 $\left| J _ { T } \right| = r ^ { 2 } \sin \theta$，相应的概率密度函数为
$$
p ( r , \theta , \phi ) = r ^ { 2 } \sin \theta p ( x , y , z )
$$
## 13.6 多维变换的二维采样 2D Sampling with Multidimensional Transformations
假设我们要采样服从 $p(x,y)$ 的 $(X,Y)$。有时多为目的是可分的 $p ( x , y ) = p _ { x } ( x ) p _ { y } ( y )$，这样可以分别采样。
给定一个 2D 密度函数，边际密度函数 $p(x)$ 为
$$
p ( x ) = \int p ( x , y ) \mathrm { d } y
$$
条件密度函数 $p ( y | x )$ 为
$$
p ( y | x ) = \frac { p ( x , y ) } { p ( x ) }
$$
从联合分布中进行 2D 采样的基本思想是首先计算边际密度，以分离一个特定的变量，然后使用标准 1D 采样技术从该密度中采样。一旦得到该采样，就可以计算给定值的条件密度函数并从该分布中采样，同样使用标准的 1D 采样技术。
### 13.6.1 半球均匀采样 Uniformly Sampling a Hemisphere
我们要在半球上关于立体角均匀采样，则有 $p(\omega)=c$。
根据 PDF 的性质，有
$$
\begin{aligned}
\int _ { \mathcal { H } ^ { 2 } } p ( \omega ) \mathrm { d } \omega &= 1\\
c \int _ { \mathcal { H } ^ { 2 } } \mathrm { d } \omega &= 1\\
c &= \frac { 1 } { 2 \pi }\\
\end{aligned}
$$
则 $p ( \omega ) = 1 / ( 2 \pi )$ 或 $p ( \theta , \phi ) = \sin \theta / ( 2 \pi )$。
注意，这个密度函数是可分的，但下边用边际和条件密度来说明多维采样技术。
首先采样 $\theta$，我们需要边际密度函数 $p(\theta)$ 
$$
p ( \theta ) = \int _ { 0 } ^ { 2 \pi } p ( \theta , \phi ) \mathrm { d } \phi = \int _ { 0 } ^ { 2 \pi } \frac { \sin \theta } { 2 \pi } \mathrm { d } \phi = \sin \theta
$$
然后计算 $\phi$ 的条件密度函数
$$
p ( \phi | \theta ) = \frac { p ( \theta , \phi ) } { p ( \theta ) } = \frac { 1 } { 2 \pi }
$$
现在我们使用 1D 反演技术来各自采样
$$
\begin{aligned} P ( \theta ) & = \int _ { 0 } ^ { \theta } \sin \theta ^ { \prime } \mathrm { d } \theta ^ { \prime } = 1 - \cos \theta \\ P ( \phi | \theta ) & = \int _ { 0 } ^ { \phi } \frac { 1 } { 2 \pi } \mathrm { d } \phi ^ { \prime } = \frac { \phi } { 2 \pi } \end{aligned}
$$
求逆，并用 $\xi$ 去替代 $1-\xi$ 
$$
\begin{aligned} \theta & = \cos ^ { - 1 } \xi _ { 1 } \\ \phi & = 2 \pi \xi _ { 2 } \end{aligned}
$$
再将他们转换为笛卡尔坐标
$$
\begin{aligned} x & = \sin \theta \cos \phi = \cos \left( 2 \pi \xi _ { 2 } \right) \sqrt { 1 - \xi _ { 1 } ^ { 2 } } \\ y & = \sin \theta \sin \phi = \sin \left( 2 \pi \xi _ { 2 } \right) \sqrt { 1 - \xi _ { 1 } ^ { 2 } } \\ z & = \cos \theta = \xi _ { 1 } \end{aligned}
$$
另外我们需要计算 PDF，要注意是哪个 PDF（如关于立体角 $\omega$，关于 $(\theta,\phi)$）。对于所有的方向采样，我们需要的都是关于立体角的 PDF。所以结果为 $1/(2\pi)$。
整个球面上采样也类似
$$
\begin{aligned} x & = \cos \left( 2 \pi \xi _ { 2 } \right) \sqrt { 1 - z ^ { 2 } } = \cos \left( 2 \pi \xi _ { 2 } \right) 2 \sqrt { \xi _ { 1 } \left( 1 - \xi _ { 1 } \right) } \\ y & = \sin \left( 2 \pi \xi _ { 2 } \right) \sqrt { 1 - z ^ { 2 } } = \sin \left( 2 \pi \xi _ { 2 } \right) 2 \sqrt { \xi _ { 1 } \left( 1 - \xi _ { 1 } \right) } \\ z & = 1 - 2 \xi _ { 1 } \end{aligned}
$$
PDF 为 $1/4\pi$。
### 13.6.2 单位圆盘上采样 Samping a Unit Disk
直觉上的 $r = \xi _ { 1 } , \theta = 2 \pi \xi _ { 2 }$ 是不均匀采样。
![1557138995302](assets/1557138995302.png)
关于面积均匀采样，因此 $p(x,y)=c$，通过标准化可知，$p(x,y)=1/\pi$。将其转换为极坐标，有 $p ( r , \theta ) = r / \pi$。计算下边际和条件概率分布
$$
\begin{aligned} p ( r ) & = \int _ { 0 } ^ { 2 \pi } p ( r , \theta ) \mathrm { d } \theta = 2 r \\ p ( \theta | r ) & = \frac { p ( r , \theta ) } { p ( r ) } = \frac { 1 } { 2 \pi } \end{aligned}
$$
计算 $P(r)$，$P^{-1}(r)$，$P(\theta)$ 和 $P^{-1}(\theta)$，得到
$$
\begin{aligned} r & = \sqrt { \xi _ { 1 } } \\ \theta & = 2 \pi \xi _ { 2 } \end{aligned}
$$
![1557139229425](assets/1557139229425.png)
尽管这个映射解决了均匀的问题，但它扭曲了区域（单位正方形被拉长和/或压扁到圆盘上）。
![1557139773236](assets/1557139773236.png)
一个更好的方法是“同心”映射，将单位正方形变成单位圆。它将 $[-1,1]^2$ 上的点映射到单位圆盘上通过均匀的映射同心正方形到同心圆。
![1557140043318](assets/1557140043318.png)
对于上图阴影部分，有
$$
\begin{aligned} r & = x \\ \theta & = \frac { y } { x } \frac { \pi } { 4 } \end{aligned}
$$
![1557140764959](assets/1557140764959.png)
不同部分对应的公式不同。
> 实现代码
>
> ```c++
>     // Map uniform random numbers to $[-1,1]^2$
>     Point2f uOffset = 2.f * u - Vector2f(1, 1);
> 
>     // Handle degeneracy at the origin
>     if (uOffset.x == 0 && uOffset.y == 0) return Point2f(0, 0);
> 
>     // Apply concentric mapping to point
>     Float theta, r;
>     if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
>         r = uOffset.x;
>         theta = PiOver4 * (uOffset.y / uOffset.x);
>     } else {
>         r = uOffset.y;
>         theta = PiOver2 - PiOver4 * (uOffset.x / uOffset.y);
>     }
>     return r * Point2f(std::cos(theta), std::sin(theta));
> ```
### 13.6.3 余弦加权半球采样 Cosine-Weighted Hemisphere Sampling
要求 $p ( \omega ) \propto \cos \theta$，标准化
$$
\begin{aligned}
\int _ { \mathcal { H } ^ { 2 } } p ( \omega ) \mathrm { d } \omega & = 1 \\
\int _ { 0 } ^ { 2 \pi } \int _ { 0 } ^ { \frac { \pi } { 2 } } c \cos \theta \sin \theta \mathrm { d } \theta \mathrm { d } \phi & = 1 \\
c \  2 \pi \int _ { 0 } ^ { \pi / 2 } \cos \theta \sin \theta \mathrm { d } \theta & = 1 \\
c & = \frac { 1 } { \pi }
\end{aligned}
$$
则
$$
p ( \theta , \phi ) = \frac { 1 } { \pi } \cos \theta \sin \theta
$$
我们可以计算边际和条件密度，可以使用 Malley 方法来生成 cosine-weighted 点。该方法的思想是，我们均匀的在单位圆盘上采样，然后将其投影到半球面上即可，得到的方向分布就是余弦分布的。
![1557141707836](assets/1557141707836.png)
> 证明
>
> 记 $(r,\phi)$ 为圆盘上的极坐标，有 $p(r,\phi)=r/\pi$。投影后 $\sin\theta=r$。为了完成变换 $( r , \phi ) = ( \sin \theta , \phi ) \rightarrow ( \theta , \phi )$，我们需要雅克比行列式
> $$
> \left| J _ { T } \right| = \left| \begin{array} { c c } { \cos \theta } & { 0 } \\ { 0 } & { 1 } \end{array} \right| = \cos \theta
> $$
> 因此
> $$
> p ( \theta , \phi ) = \left| J _ { T } \right| p ( r , \phi ) = \cos \theta \frac { r } { \pi } = ( \cos \theta \sin \theta ) / \pi
> $$
这个方法不要求具体的圆盘采样方法，我们可以使用之前的同心映射或者简单的 $( r , \theta ) = \left( \sqrt { \xi _ { 1 } } , 2 \pi \xi _ { 2 } \right)$。
相应的 PDF 为 $\cos\theta/\pi$。
### 13.6.4 圆锥采样 Sampling a Cone
### 13.6.5 三角形采样 Sampling a Triangle
### *13.6.6 相机采样 Sampling Cameras
### 13.6.7 分段常数 2D 分布 Piecewise-Constant 2D Distributions
考虑一个 2D 函数 $f(u,v)$，它由$n_u\times n_v$ 个值的 $f[u_i,v_j]$ 定义，其中 $u _ { i } \in \left[ 0,1 , \ldots , n _ { u } - 1 \right] , v _ { j } \in \left[ 0,1 , \ldots , n _ { v } - 1 \right]$，且 $f[u_i,v_j]$ 给出了 $f(x)$ 在区间 $\left[ i / n _ { u } , ( i + 1 ) / n _ { u } \right) \times \left[ j / n _ { v } , ( j + 1 ) / n _ { v } \right)$ 上的常数值。
用 $( \tilde { u } , \tilde { v } )$ 记作 $(u_i,v_j)$ 相应的离散值，有 $\tilde { u } = \left\lfloor n _ { u } u \right\rfloor$ 和 $\tilde { v } = \left\lfloor n _ { v } v \right\rfloor$，则有
$$
f ( u , v ) = f [ \tilde { u } , \tilde { v } ]
$$
积分值为
$$
I _ { f } = \iint f ( u , v ) \mathrm { d } u \mathrm { d } v = \frac { 1 } { n _ { u } n _ { v } } \sum _ { i = 0 } ^ { n _ { u } - 1 } \sum _ { j = 0 } ^ { n _ { v } - 1 } f \left[ u _ { i } , v _ { j } \right]
$$
得到 PDF
$$
p ( u , v ) = \frac { f ( u , v ) } { \iint f ( u , v ) \mathrm { d } u \mathrm { d } v } = \frac { f [ \tilde { u } , \tilde { v } ] } { 1 / \left( n _ { u } n _ { v } \right) \sum _ { i } \sum _ { j } f \left[ u _ { i } , v _ { j } \right] }
$$
边际密度分布 $p(v)$ 为
$$
p ( v ) = \int p ( u , v ) \mathrm { d } u = \frac { \left( 1 / n _ { u } \right) \sum _ { i } f \left[ u _ { i } , \tilde { v } \right] } { I _ { f } }
$$
这是一个分段常数 1D 函数，13.3.1 节有说明采样方法。
条件密度分布 $p(u|v)$ 为
$$
p ( u | v ) = \frac { p ( u , v ) } { p ( v ) } = \frac { f [ \tilde { u } , \tilde { v } ] / I _ { f } } { p [ \tilde { v } ] }
$$
给定一个特定的 $\tilde{v}$，$p[\tilde{u}|\tilde{v}]$ 是一个分段常数 1D 函数，总共有 $n_v$ 个这样的 1D 条件密度。
采样的过程分 2 步，先采样 v，然后选出特定的条件分布，再采样 u。两个 pdf 相乘得到联合的 pdf。
时间复杂度为 $O(\log n_u \log n_v)$，空间复杂度为 $O(n_u n_v)$。
> 可以用 Alias Method 算法将时间复杂度减为 $O(1)$，空间复杂度还是 $O(n_u n_v)$。
## 13.7 俄罗斯轮盘赌和分裂 Rissian Roulette and Spitting
估计式 F 的效率 efficiency 定义为
$$
\epsilon [ F ] = \frac { 1 } { V [ F ] T [ F ] }
$$
其中 $V[F]$ 是方差，$T[F]$ 是运行时间。
俄罗斯轮盘赌解决了样本评估昂贵但对最终结果贡献较小的问题。分裂是一种可以在被积函数的重要维度上放置更多样本的技术。
考虑一个关于俄罗斯轮盘赌的例子，计算直接光照积分
$$
L _ { \mathrm { o } } \left( \mathrm { p } , \omega _ { \mathrm { o } } \right) = \int _ { \delta ^ { 2 } } f _ { \mathrm { r } } \left( \mathrm { p } , \omega _ { \mathrm { o } } , \omega _ { \mathrm { i } } \right) L _ { \mathrm { d } } \left( \mathrm { p } , \omega _ { \mathrm { i } } \right) \left| \cos \theta _ { i } \right| \mathrm { d } \omega _ { \mathrm { i } }
$$
只采样两次
$$
\frac { 1 } { 2 } \sum _ { i = 1 } ^ { 2 } \frac { f _ { \mathrm { r } } \left( \mathrm { p } , \omega _ { \mathrm { o } } , \omega _ { i } \right) L _ { \mathrm { d } } \left( \mathrm { p } , \omega _ { i } \right) \left| \cos \theta _ { i } \right| } { p \left( \omega _ { i } \right) }
$$
开销主要花在计算可见性上。
对于所有 $f _ { \mathrm { r } } \left( \mathrm { p } , \omega _ { \mathrm { o } } , \omega _ { i } \right) = 0$ 的方向 $\omega_i$，我们应该跳过可见性计算。俄罗斯轮盘赌让我们当积分值很小时跳过可见性计算，如 $f _ { \mathrm { r } } \left( \mathrm { p } , \omega _ { \mathrm { o } } , \omega _ { i } \right)$ 很小，或者 $|\cos\theta_i|$ 很小。
为了用俄罗斯轮盘赌，我们选择停止概率 q，这个值可以用很多方式来抽取。以概率 q 终止计算，并用常数 c 替代（经常 $c=0$）。以概率 $1-q$ 继续计算，但是加权 $1/(1-q)$。
如此得到一个新的估计式
$$
F ^ { \prime } = \left\{ \begin{array} { l l } { \frac { F - q c } { 1 - q } } & { \xi > q } \\ { c } & { \text { otherwise } } \end{array} \right.
$$
期望为
$$
E \left[ F ^ { \prime } \right] = ( 1 - q ) \left( \frac { E [ F ] - q c } { 1 - q } \right) + q c = E [ F ]
$$
俄罗斯轮盘赌从不减少方差，只有 $c=F$ 的时候方差不变。
### 13.7.1 分裂 Spiltting
忽略像素过滤，计算只考虑直接光的积分
$$
\int _ { A } \int _ { 8 ^ { 2 } } L _ { \mathrm { d } } ( x , y , \omega ) \mathrm { d } x \mathrm { d } y \mathrm { d } \omega
$$
计算时，需要一个射线找到可见点，然后在来一条射线判断该点与光源的可见性。
如果有 100 个图像样本，则需要 200 条射线。可能样本数过多了。可以考虑分裂，采样 N 个，一个图像样本对应 M 个光源样本。这样估计式为
$$
\frac { 1 } { N } \frac { 1 } { M } \sum _ { i = 1 } ^ { N } \sum _ { j = 1 } ^ { M } \frac { L \left( x _ { i } , y _ { i } , \omega _ { i , j } \right) } { p \left( x _ { i } , y _ { i } \right) p \left( \omega _ { i , j } \right) }
$$
取 N = 5，M = 20。则总共只要 105 条光线。
## 13.8 仔细的样本布置 Careful Sample Placement
分层抽样的工作原理是将积分域细分为 n 个互不重叠的区域 $\Lambda _ { 1 } , \Lambda _ { 2 } , \dots , \Lambda _ { n }$。每个区域称为一个层 stratum。要求完全覆盖原本的区域
$$
\bigcup _ { i = 1 } ^ { n } \Lambda _ { i } = \Lambda
$$
在单层 $\Lambda_i$ 中，蒙特卡洛估计为
$$
F _ { i } = \frac { 1 } { n _ { i } } \sum _ { j = 1 } ^ { n _ { i } } \frac { f \left( X _ { i , j } \right) } { p _ { i } \left( X _ { i , j } \right) }
$$
总的估计为
$$
F=\sum _ { i = 1 } ^ { n } v _ { i } F _ { i }
$$
其中 $v_i$ 是层 i 的体积占比（$v_i\in(0,1]$）。
层 i 的积分真值为
$$
\mu _ { i } = E \left[ f \left( X _ { i , j } \right) \right] = \frac { 1 } { v _ { i } } \int _ { \Lambda _ { i } } f ( x ) \mathrm { d } x
$$
方差为
$$
\sigma _ { i } ^ { 2 } = \frac { 1 } { v _ { i } } \int _ { \Lambda _ { i } } \left( f ( x ) - \mu _ { i } \right) ^ { 2 } \mathrm { d } x
$$
层 i 的方差为 $\sigma^2_i/n_i$。则总的估计式的方差为
$$
\begin{aligned} V [ F ] & = V \left[ \sum v _ { i } F _ { i } \right] \\ & = \sum V \left[ v _ { i } F _ { i } \right] \\ & = \sum v _ { i } ^ { 2 } V \left[ F _ { i } \right] \\ & = \sum \frac { v _ { i } ^ { 2 } \sigma _ { i } ^ { 2 } } { n _ { i } } \end{aligned}
$$
如果样本数 $n_i$ 正比于体积 $v_i$，我们有 $n_i=v_i N$，这样方差变为
$$
V \left[ F _ { N } \right] = \frac { 1 } { N } \sum v _ { i } \sigma _ { i } ^ { 2 }
$$
未分层的方差为
$$
\begin{aligned} V [ F ] & = E _ { x } V _ { i } F + V _ { x } E _ { i } F \\ & = \frac { 1 } { N } \left[ \sum v _ { i } \sigma _ { i } ^ { 2 } + \sum v _ { i } \left( \mu _ { i } - Q \right) \right] \end{aligned}
$$
> See Veach (1997) for a derivation of this result. 
右边部分非负，说明分层永远不会增加方差。只有当每层的均值都为总均值时为 0。为了让分层抽样工作的更好，我们会去最大化右边的和，所以我们最好让层的均值尽可能的不同。所以在不知道 $f$ 的情况下使用致密的层是可取的。如果层很大，他们包含更多的变化，使得层均值更接近与 Q。
> 示例
>
> ![1557156779823](assets/1557156779823.png)
>
> 左图随机采样，右图分层采样
>
> 左图比右图方差更大，噪声更多
### 13.8.2 准蒙特卡洛 Quasi Monte Carlo
准蒙特卡洛的关键就是用低差异采样替代伪随机数。
### 13.8.3 Warpping Samples and Distortion
当使用分层采样或低差异采样时，`pbrt` 会变换样本。这个过程隐式的要求是保持样本的分层性质。
## 13.9 偏差 Bias
引入偏差 bias 也是降低方差的一种方法：有时计算一个由偏差的估计式方差更小。
估计式是无偏的 unbiased，当它的期望值等于真值，否则，偏差为
$$
\beta = E [ F ] - \int f ( x ) \mathrm { d } x
$$
> 求在 $[0,1)$ 上均匀分布的 $X_i$ 的均值
>
> 一种估计式为
> $$
> \frac { 1 } { N } \sum _ { i = 1 } ^ { N } X _ { i }
> $$
> 它是无偏的，方差$O(N^{-1})$。
>
> 另一种估计式为
> $$
> \frac { 1 } { 2 } \max \left( X _ { 1 } , X _ { 2 } , \ldots , X _ { N } \right)
> $$
> 它是有偏的
> $$
> 0.5 \frac { N } { N + 1 } \neq 0.5
> $$
> 但方差 $O(N^{-2})$。
## 13.10 重要性采样 Importance Samping
重要性采样是很强力的方差减少技术，它利用蒙特卡洛估计
$$
F _ { N } = \frac { 1 } { N } \sum _ { i = 1 } ^ { N } \frac { f \left( X _ { i } \right) } { p \left( X _ { i } \right) }
$$
在样本分布 $p(x)$ 接近 $f(x)$ 时收敛更快的特性。
当 $p(x)=f(x)/\int f(x)\mathrm{d}x$ 时，一个样本就可以得到结果，方差为 0。这要求我们需要知道 $\int f(x)\mathrm{d}x$，而这正是我们使用蒙特卡洛方法所要计算的结果。蒙特卡洛方法和能够直接积分 $f(x)$ 是矛盾的。然而如果 $p(x)$ 与 $f(x)$ 相似的话，方差会减少。
如果采用了一个糟糕的采样分布，重要性采样会增加方差，比均匀采样还差。
幸运的是，找到一个好的采样分布不会太难。许多情况下，积分是多个函数的成绩，我们很难构建一个相似与整个乘积的 PDF，但是找到与一个乘数相似的分布也是有用的。
### 13.10.1 多重重要性采样 Multiple Impotance Sampling
当计算 $\int f(x)g(x)\mathrm{d}x$ 时，我们有分别针对 $f(x)$ 和 $g(x)$ 的两个重要性采样策略，我们该用哪个？（假设我们没法构建一个 PDF 使其正比于 $f(x)g(x)$）
考虑计算
$$
L _ { \mathrm { o } } \left( \mathrm { p } , \omega _ { \mathrm { o } } \right) = \int _ { \delta ^ { 2 } } f \left( \mathrm { p } , \omega _ { \mathrm { o } } , \omega _ { \mathrm { i } } \right) L _ { \mathrm { d } } \left( \mathrm { p } , \omega _ { \mathrm { i } } \right) \left| \cos \theta _ { \mathrm { i } } \right| \mathrm { d } \omega _ { \mathrm { i } }
$$
如果我们只基于 $f$ 或 $L_d$ 进行重要性采样，其中一个经常会表现糟糕。
不幸的是，从每个分布中抽取一些样本并对两个估计值求平均值的明显解决方案并不好。因为方差在这种情况下是可加性的，这种方法没有帮助——一旦方差进入一个估计值，我们不能通过将它加入另一个估计值来消除它，即使它本身方差很小。
多重重要性采样 MIS 解决了这个问题，只用了一个简单且易于实现的技术。基本思想是，当计算一个积分时，我们从多个分布中采样，期望至少一种一个能匹配上积分式的形状，即使我们并不知道具体是哪一个。MIS提供了一种方法来对来自每种技术的样本进行加权，这可以消除由于被积函数值和采样密度之间的不匹配而导致的大的方差尖峰。
如果两个采样分布 $p_f$ 和 $p_g$ 用于估计 $\int f(x)g(x)\mathrm{d}x$，则新的蒙特卡洛估计式为
$$
\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 _ { f } \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 _ { g } \left( Y _ { j } \right) }
$$
其中 $n_f$ 是按 $p_f$ 的采样个数，$n_g$ 是按 $p_g$ 的采样个数，并且 $w_f$ 和 $w_g$ 是特殊的权重函数，使得这个估计式的期望为 $\int f(x)g(x)\mathrm{d}x$。
权重函数对每个样本考虑了所有的采样方法，而不仅仅是该样本生成时所用的采样方法。
一个好的权重函数是平衡启发式 balance heuristic
$$
w _ { s } ( x ) = \frac { n _ { s } p _ { s } ( x ) } { \sum _ { i } n _ { i } p _ { i } ( x ) }
$$
`pbrt` 只需要两函数的权重函数，不需要更泛化的情况
实践中，幂启发式 power heuristic 能更好的减少方差
$$
w _ { s } ( x ) = \frac { \left( n _ { s } p _ { s } ( x ) \right) ^ { \beta } } { \sum _ { i } \left( n _ { i } p _ { i } ( x ) \right) ^ { \beta } }
$$
$\beta=2$ 是一个不错的选择