《Fundamentals of Computer Graphics》第五版第十三章采样图形学中许多应用都需要对某种

积分（Integration）

粗略地讲，测度（measure）就是一个，以一种和长度、面积、体积一致的方式，将某个集合的子集映射到 $\mathbb{R}^{+}$ 的函数，比如，二维实平面 $\mathbb{R}^{2}$ 上的面积测度 $A$ 。测度的定义域是一个集合的幂集（power set）——所有子集构成的集合，比如，面积测度 $A$ 的定义域是幂集 $\mathcal{P}(\mathbb{R}^{2})$ ，因此面积测度可记为：

A: \mathcal{P}(\mathbb{R}^{2}) \rightarrow \mathbb{R}^{+}

平面上任一点的面积测度均为零，这种测度为零的集合称为零测集（zero measure sets）。

一般地，测度 $\mu: \mathcal{P}(\mathbb{S})\rightarrow\mathbb{R}^{+}$ 应满足如下条件：

$\mu(\emptyset)=0$
可数个两两不相交的集合的并满足可加性：
$\displaystyle\mu\left(\bigcup_{k=1}^{\infty}S_{k}\right)=\sum_{k=1}^{\infty}\mu(S_{k})$
对于两个可能相交的集合 $A$ 、 $B$ ：
$\mu(A\bigcup B)=\mu(A)+\mu(B)-\mu(A\bigcap B)$

通常使用积分（integration）来计算测度：

A(S) = \int_{x\in S}dA

给定集合 $\mathbb{S}$ 上的测度，总可以通过非负权重函数 $w:\mathbb{S}\rightarrow\mathbb{R}^{+}$ 来生成新测度：

\nu(S) = \int_{x\in S}w(x)d\mu

函数 $f$ 在集合 $S$ 上相对于测度 $\mu$ 的平均值：

\text{average}(f) \equiv \frac{\displaystyle\int_{x\in S}f(x)d\mu}{\displaystyle\int_{x\in S}d\mu}

这在积分几何——一个研究几何实体上的测度的领域——中经常出现。

对于二维平面上的直线，自然测度（natural measure）应该与坐标系无关，即旋转和平移不会改变自然测度。一条直线可以表示为对偶空间（dual space）中的一个点。比如，在 $(\phi, b)$ 空间中，自然测度为：

d\mu = \cos\phi\ d\phi\ db

其中， $b$ 是截距， $\phi\in[-\pi/2,\pi/2)$ 是直线与 x 轴的夹角。

自然测度也可以变换到其它参数空间中，比如， $(m,b)$ 空间：

d\mu = \frac{dm\ db}{(1 + m^{2})^{\frac{3}{2}}} = \frac{du\ dv}{(u^{2} + v^{2})^{\frac{3}{2}}} = \frac{ab\ da\ db}{(a^{2} + b^{2})^{\frac{3}{2}}} = dp\ d\theta

其中， $m$ 是直线的斜率； $ux+vy+1=0$ ； $a$ 、 $b$ 是 x 轴、y 轴截距； $x\cos\theta+y\sin\theta-p=0$ 。这里，简正坐标（normal coordinates） $p$ 、 $\theta$ 分别表示原点到直线的距离、法向量相对于 x 轴的偏角。

3D 直线可以表示为 $(x,y,\theta,\phi)$ ，其中 $(x,y)$ 是直线与 $xy$ 平面的交点， $(\theta,\phi)$ 是直线方向的球坐标，约定 $\theta\in[0,\pi/2]$ 。也就是说，3D 直线对应 4D 空间中的一个点。同样地，自然测度应该与 $(x,y)$ 无关，且正比于横截面，因此自然测度为：

d\mu = \cos\theta\ dx\ dy\ \sin\theta\ d\theta\ d\phi

原文中应该是少了横截面产生的 $\cos\theta$ 因子。

也可以使用直线与平面 $z=0$ 、 $z=1$ 的交点 $(u,v)$ 、 $(s,t)$ 所构成的四元组 $(u,v,s,t)$ 来描述直线，对于平行于 z 轴的直线束， $d\mu\approx du\ dv\ ds\ dt$ 。

给定一个半径为 $r$ 的球，对于那些与球相交的直线，可以使用两个交点的球坐标 $(\theta_{1},\phi_{1},\theta_{2},\phi_{2})$ 来表示。易知，垂直于直线方向的面元：

\begin{aligned} dS_{1\perp} &= \vec{n}_{1}\cdot\frac{\vec{n}_{1} - \vec{n}_{2}}{|\vec{n}_{1} - \vec{n}_{2}|}dS_{1} = \frac{1 - \vec{n}_{1}\cdot\vec{n}_{2}}{|\vec{n}_{1} - \vec{n}_{2}|}dS_{1} \\ dS_{2\perp} &= \vec{n}_{2}\cdot\frac{\vec{n}_{2} - \vec{n}_{1}}{|\vec{n}_{2} - \vec{n}_{1}|}dS_{2} = \frac{1 - \vec{n}_{1}\cdot\vec{n}_{2}}{|\vec{n}_{1} - \vec{n}_{2}|}dS_{2} \end{aligned}

$\vec{n}_{1}$ 、 $\vec{n}_{2}$ 是球面上面元 $dS_{1}$ 、 $dS_{2}$ 的单位方向向量。易知，自然测度应正比于：

$dS_{1\perp}$
$dS_{2\perp}$ 相对于 $dS_{1\perp}$ 上任一点所张开的立体角

\begin{aligned} d\mu &\propto dS_{1\perp}\frac{dS_{2\perp}}{r^{2}(\vec{n}_{1} - \vec{n}_{2})^{2}} \\ &=\frac{(1 - \vec{n}_{1}\cdot\vec{n}_{2})^{2}}{r^{2}(\vec{n}_{1} - \vec{n}_{2})^{4}}dS_{1}dS_{2} \\ &=\frac{(1 - \vec{n}_{1}\cdot\vec{n}_{2})^{2}}{r^{2}(2 - 2\vec{n}_{1}\cdot\vec{n}_{2})^{2}}dS_{1}dS_{2} \\ &\propto dS_{1}dS_{2} \end{aligned}

也就是说，在所有与给定球相交的直线构成的子集上，自然测度可以表示为：

d\mu = \sin\theta_{1}\ d\theta_{1}\ d\phi_{1}\ \sin\theta_{2}\ d\theta_{2}\ d\phi_{2}

上式表明，球面上均匀分布的两点的连线服从均匀分布。

连续型概率

连续型随机变量

连续型随机变量（continuous random variable） $x$ 的行为完全由它的概率密度函数（probability density function，简称 PDF） $p(x)$ 来描述：

\begin{gathered} p(x) \geqslant 0 \\ \int_{-\infty}^{+\infty}p(x)dx = 1 \end{gathered}

$x$ 处于区间 $[a,b]$ 上的概率可以表述为：

\text{Probability}(x\in [a, b]) = \int_{a}^{b}p(x)dx

规范随机变量（canonical random variable） $\xi$ 满足如下分布：

q(\xi) = \left\{ \begin{aligned} &1 && \text{if}\ 0 \leqslant\xi<1 \\ &0 && \text{otherwise} \end{aligned} \right.

对于随机变量 $x$ ，函数值 $f(x)$ 也是随机变量，它的期望值（expected value）定义如下：

E(f(x)) = \int_{-\infty}^{+\infty}f(x)p(x)dx

随机变量的期望具有如下性质：

E(f(x) + g(y)) = E(f(x)) + E(g(y))

不管随机变量 $x$ 、 $y$ 是否关联，上式均成立。

只要在随机变量的样本空间上定义一个测度 $\mu$ ，就可以很容易推广到多维随机变量的情况：

\begin{gathered} \text{Probability}(x\in S_{i}) = \int_{S_{i}}p(x)d\mu \\ E(f(x)) = \int_{S}f(x)p(x)d\mu \end{gathered}

一维随机变量的方差（variance）定义如下：

\begin{aligned} V(x) &\equiv E([x - E(x)]^{2}) \\ &=E(x^{2}) - [E(x)]^{2} \end{aligned}

一组相互独立的随机变量之和的方差等于随机变量方差之和。标准差（standard deviation） $\sigma=\sqrt{V}$ 。

样本均值

独立同分布（independent identically distributed，简称为 iid）随机变量 $x_{i}$ 的平均值——样本均值——可用于估计期望值：

E(x) \approx \frac{1}{N}\sum_{i=1}^{N}x_{i}

随着 $N$ 的增加，样本均值的方差 $\displaystyle V\left(\frac{1}{N}\sum_{i=1}^{N}x_{i}\right)=\frac{V(x_{i})}{N}$ 会逐渐减小。而大数定律（law of large numbers）则是使用样本均值估计期望值的基石：

\lim_{N\rightarrow\infty}\text{Probability}\left(\left|\frac{1}{N}\sum_{i=1}^{N}x_{i} - E(x)\right| < \epsilon\right) = 1

$\epsilon$ 为任意正实数。

蒙特卡罗积分（Monte Carlo Integration）

函数积分可以看作随机变量的期望值，进而用样本均值来近似：

\int_{x\in S}g(x)d\mu \approx \frac{1}{N}\sum_{i=1}^{N}\frac{g(x_{i})}{p(x_{i})}

其中，随机变量 $x_{i}$ 服从 $p(x)$ 分布。由于上式右侧的方差为 $\displaystyle\frac{1}{N}V\left(\frac{g(x)}{p(x)}\right)$ ，因此为了得到较好的结果，可以从以下两点入手：

尽可能增加样本数目 $N$ ；
让 $g(x)/p(x)$ 的方差尽可能小，即函数 $g$ 、 $p$ 形状近似，这种选取 $p$ 的方式称为重要性采样（importance sampling）。

样本均值的方差也暗含了蒙特卡罗积分（Monte Carlo integration）的一个基本问题：收益递减（diminishing return）。由于标准差 $\sigma\propto 1/\sqrt{N}$ ，因此，为了降低一半误差，需将样本数增加为原来的四倍。

另一种降低方差的方法是将积分区域 $S$ 划分为多个小区域 $S_{i}$ ，然后对 $S_{i}$ 分别计算积分值，这也称为分层采样（stratified sampling）。通常，每个 $S_{i}$ 中只以分布 $p_{i}$ 采一个样本，这种估计的方差为：

V\left(\sum_{i=1}^{N}\frac{g(x_{i})}{p_{i}(x_{i})}\right) = \sum_{i=1}^{N}V\left(\frac{g(x_{i})}{p_{i}(x_{i})}\right)

如果等概率分层，即 $\displaystyle\int_{S_{i}}p(x)d\mu=1/N$ ，可以证明，分层采样的方差一定不超过未分层采样。对于等概率分层：

p_{i}(x) = \left\{ \begin{aligned} &Np(x) && x\in S_{i} \\ &0 && x\notin S_{i} \end{aligned} \right.

由柯西不等式可知：

N\sum_{i=1}^{N}\left(\int_{S_{i}}g(x)dx\right)^{2} \geqslant \left(\int_{S}g(x)dx\right)^{2}

因此：

\begin{aligned} \sum_{i=1}^{N}V\left(\frac{g(x_{i})}{p_{i}(x_{i})}\right) &= \frac{1}{N}\sum_{i=1}^{N}\int_{S_{i}}\left(\frac{g(x)}{p(x)}\right)^{2}p(x)dx - \sum_{i=1}^{N}\left(\int_{S_{i}}g(x)dx\right)^{2} \\ &\leqslant \frac{1}{N}\left[\int_{S}\left(\frac{g(x)}{p(x)}\right)^{2}p(x)dx - \left(\int_{S}g(x)dx\right)^{2}\right] \\ &= \frac{1}{N}V\left(\frac{g(x)}{p(x)}\right) \end{aligned}

由于平均上每个分区的方差通常会小于整个区域的方差，因此对于大多数函数，分层采样都会改善其采样效率。分层采样通常远优于重要性采样。

用准随机点（quasi-random point）代替随机点是一种常用的方法，以此所做的蒙特卡罗积分称为准蒙特卡罗积分（quasi–Monte Carlo integration）。准随机点是确定性的，但分布均匀，一个区域上的点数正比于该区域的测度，比如晶格上的规则样本点。准随机点可以改善性能，但需注意不要引入混淆。

随机采样

生成任意分布的随机点或伪随机点有多种方法：反函数法（function inversion）、舍选法（rejection）、Metropolis。

高维随机变量可以使用边缘概率分布逐个生成每个分量。

反函数法（Function Inversion）

若 $\xi$ 是 $[0,1]$ 上均匀分布的随机变量，则可以通过单调增函数 $x(t)$ ，将 $\xi$ 变换为区间 $[x_{\text{min}},x_{\text{max}}]$ 上服从分布 $f(x)$ 的随机变量 $\alpha$ 。注意到下式成立：

\begin{aligned} & &\text{Probability}(\xi < t) &= \text{Probability}(\alpha < x) \\ \Rightarrow& & \int_{0}^{t}dt' &= \int_{x_{\text{min}}}^{x}f(x')dx' \\ \Rightarrow& & t &= P(x) \end{aligned}

因此，随机变量 $\alpha$ 可由 $\xi$ 表示为：

\alpha = P^{-1}(\xi)

其中， $P$ 是累积概率分布函数（cumulative probability distribution function）。

对于 $[x_{\text{min}},x_{\text{max}}]\times[y_{\text{min}},y_{\text{max}}]$ 上满足分布 $f(x,y)$ 的二维随机变量 $(\alpha_{x},\alpha_{y})$ ，分布函数：

\text{Probability}(\alpha_{x} < x, \alpha_{y} < y) = F(x, y) = \int_{y_{\text{min}}}^{y}dy\int_{x_{\text{min}}}^{x}dxf(x', y')

$[0,1]$ 上均匀分布的随机变量 $\xi_{1}$ 、 $\xi_{2}$ 与随机变量 $\alpha_{x}$ 、 $\alpha_{y}$ 之间，由下式所确定的函数 $x(t_{1},t_{2})$ 、 $y(t_{1},t_{2})$ 联系起来：

\begin{aligned} t_{1} &= F(x, y_{\text{max}}) \\ t_{2} &= \frac{F'_{x}(x, y)}{F'_{x}(x, y_{\text{max}})} \end{aligned}

其中， $F(x, y_{\text{max}})$ 是边缘分布（marginal distribution）， $F'_{x}$ 是分布函数 $F$ 对 $x$ 的一阶偏导。注意到：

\begin{gathered} \frac{\partial(t_{1}, t_{2})}{\partial(x, y)} = F''_{xy}(x, y) = f(x, y) \\ \Rightarrow\quad dt_{1}dt_{2} = f(x, y)dxdy \end{gathered}

可见由此生成的随机变量 $\alpha_{x}$ 、 $\alpha_{y}$ 满足分布 $f(x,y)$ 。

通常在球面上采样时，需要将角度 $(\theta,\phi)$ 转为方向向量：

\vec{a} = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)

舍选法（Rejection）

舍选法（rejection）就是从一系列随机数中选出一部分符合条件的，舍去那些不符合条件的。要想生成满足分布 $p:[a,b]\mapsto\mathbb{R}$ 的随机数，假定 $p(x)<m$ ，可以先生成矩形 $[a,b]\times[0,m]$ 上均匀分布的随机数，然后选出 $y<p(x)$ 的 $x$ ：

\begin{aligned} &\text{done = false} \\ &\textbf{while}\ (\text{not done})\ \textbf{do} \\ &\quad\begin{aligned} &x = a + r()(b - a) \\ &y = r()m \\ &\textbf{if}\ (y < p(x))\ \textbf{then} \\ &\quad\text{done = true} \end{aligned} \end{aligned}

这里， $r$ 是标准随机数函数。使用舍选法生成球面上均匀分布的随机点：

\begin{aligned} &\text{done = false} \\ &\textbf{while}\ (\text{not done})\ \textbf{do} \\ &\quad\begin{aligned} &x = -1 + 2r() \\ &y = -1 + 2r() \\ &z = -1 + 2r() \\ &\textbf{if}\ ((l = \sqrt{x^{2} + y^{2} + z^{2}}) < 1)\ \textbf{then} \\ &\quad\text{done = true} \end{aligned} \\ &x = x / l \\ &y = y / l \\ &z = z / l \end{aligned}

舍选法易实现，但收敛较慢，因此主要用于调试或较困难的场景。

Metropolis

Metropolis 方法使用随机跃迁来生成指定分布的随机样本。假定当前样本点为 $x$ ，先按照跃迁概率密度 $t(x\rightarrow y)$ 生成一个候选点 $y$ ，然后以概率 $a(x\rightarrow y)$ 决定是否接受，依次类推，生成一个随机序列。通过约定细致平衡条件——任意两点间沿相反方向的流相等，让随机过程达到指定分布 $f(x)$ 的稳态：

\begin{aligned} & &\text{flow}(x\rightarrow y) &= \text{flow}(y\rightarrow x) \\ \Rightarrow& & f(x)t(x\rightarrow y)a(x\rightarrow y) &= f(y)t(y\rightarrow x)a(y\rightarrow x) \\ \Rightarrow& & \frac{a(y\rightarrow x)}{a(x\rightarrow y)} &= \frac{f(x)t(x\rightarrow y)}{f(y)t(y\rightarrow x)} \end{aligned}

这样就可以生成满足分布 $f(x)$ 的随机变量。实际操作中通常会把接受率 $a(x\rightarrow y)$ 、 $a(y\rightarrow x)$ 中较大者设为 $1$ 。Metropolis 采样的难点在于估计暂态过程的样本点数。

《Fundamentals of Computer Graphics》第五版 第十三章 采样