蒙特卡洛积分入门

最近脑子都快生锈了，复习一下基本方法。

数值求积规则

给定某定义域 $Ω$ （比如 $s$ 维立方体 $[0, 1]^{s}$ ）和其上的函数 $f : Ω \to R$ ，令测度函数为 $d μ (x) = d x^{1} \dots d x^{s}$ ，则积分：

I = \int_{Ω} f (x) d μ (x)

常被用以下规则估值：

\hat{I} = i = 1 \sum N w_{i} f (x_{i})

如何选择 $w_{i}$ 和 $x_{i}$ ？这可谓一个深不见底的大坑，常见的、针对一维情形的梯形法则、辛普森规则等大多具有形如 $O (n^{- r}) (r \in N^{+})$ 的收敛率。在高维情形下，这一规则可以被自然地推广为：

\hat{I} = i_{1} = 1 \sum n i_{2} = 1 \sum n \dots i_{s} = 1 \sum n w_{i_{1}} w_{i_{2}} \dots w_{i_{s}} f (x_{i_{1}}, x_{i_{2}}, \dots, x_{i_{s}})

推广版本的规则依然保持着 $O (n^{- r})$ 的收敛率，但这要求 $N = n^{s}$ 个采样点。换言之，若保持采样数量 $N$ 不变，则收敛率其实降低到了 $O (N^{- r / s})$ ，维数一上去，效率就撑不住了。

当然，不是所有高维求积规则都使用这种张量积的形式，但Bakhvalov定理已经说了，在座的各位都是垃圾：

Bakhvalow’s Theorem. 任给 $s$ 维求积规则，存在具有 $r$ 阶有界连续导数的函数 $f$ ， $f$ 在该规则下的误差是 $O (N^{- r / s})$ 。

估计量

给定积分：

I = \int_{Ω} f (x) d μ (x)

令估计量 $F_{N}$ 为：

F_{N} = \frac{1}{N} i = 1 \sum N \frac{f ( X _{i} )}{p ( X _{i} )}

其中 $p$ 是用于采样 $X_{i}$ 的概率分布函数。 $F_{N}$ 的方差为

V [F_{N}] = V [\frac{1}{N} i = 1 \sum N \frac{f ( X _{i} )}{p ( X _{i} )}] = \frac{1}{N ^{2}} i = 1 \sum N V [\frac{f ( X _{i} )}{p ( X _{i} )}] = \frac{1}{N} V [\frac{f ( X )}{p ( X )}]

可见 $F_{N}$ 的收敛速度是 $O (N)$ 。即使 $V [f (X) / p (X)]$ 是无穷，只要 $E [f (X) / p (X)]$ 存在，根据强大数定律 $F_{N}$ 也一定会收敛（只不过要慢一些）。

随机变量采样

已知分布函数 $P$ ，令 $U$ 是 $[0, 1]$ 上的均匀随机变量， $X = P^{- 1} (U)$ ，则 $X$ 符合分布 $P$ 。
已知概率密度函数 $p$ ，若存在概率密度函数 $q$ 和常量 $M$ 使得 $p (x) \leq Mq (x)$ ，则可先按 $q$ 采样得到一个 $X_{i}$ ，然后在 $[0, 1]$ 内均匀采样得到一个 $U_{i}$ ，若 $U_{i} \leq p (X_{i}) / (Mq (X_{i}))$ 则令采样结果为 $X_{i}$ ，否则重来一遍。这样得到的采样点符合 $p$ 。
Metropolis方法，这里不详述。

设 $F_{N}$ 是某个量 $Q$ 的一个估计量，则 $F_{N} - Q$ 称为误差，误差的期望值称为偏差：

β [F_{N}] = E [F_{N} - Q]

偏差为0的估计量被称为是无偏的（unbiased）。另一方面，一个估计量是一致的（consistent）当且仅当误差随 $N$ 的增大以概率1收敛到0，即：

P r {N \to \infty lim F_{N} = Q} = 1

定义均方误差 $MSE [F] = E [(F - Q)^{2}]$ ，则：

MSE [F] = E [(F - Q)^{2}] = (E [F^{2}] - 2 E^{2} [F] + E^{2} [F]) + (E^{2} [F] - 2 E [F] Q + Q^{2}) = E [(F^{2} - 2 FE [F] + E^{2} [F])] + E [(E^{2} [F] - 2 E [F] Q + Q^{2})] = E [(F - E [F])^{2}] + (E [F] - Q)^{2} = V [F] + β^{2} [F]

对无偏估计量而言，其MSE就等于方差。于是，对某个无偏估计量 $Y$ ，设 $Y_{1}, Y_{2}, \dots, Y_{N}$ 是其 $N$ 个独立采样值，则：

\hat{V} [F_{N}] = \frac{1}{N - 1} ⎩ ⎨ ⎧ (\frac{1}{N} i = 1 \sum N Y_{i}^{2}) - (\frac{1}{N} i = 1 \sum N Y_{i})^{2} ⎭ ⎬ ⎫

是无偏估计量

F_{N} = \frac{1}{N} i = 1 \sum N Y_{i}

的方差的无偏估计量。

一些常见的减小方差的方法：

用条件期望降维

若能计算出：

f (X) = \int f (x, y) d y p (x) = \int p (x, y) d y

则可作以下替换：

F = \frac{f ( X , Y )}{p ( X , Y )} \Rightarrow F^{'} = \frac{f ( X )}{p ( X )}

这是因为：

E_{Y} [\frac{f ( X , Y )}{p ( X , Y )}] = \int \frac{f ( X , y )}{p ( X , y )} p (y ∣ X) d y = \int \frac{f ( X , y )}{p ( X , y )} \frac{p ( X , y )}{\int p ( X , y ^{'} ) d y ^{'}} d y = \frac{f ( X )}{p ( X )}

即 $F^{'} (X) = E_{Y} [F (X, Y)]$ 。又注意到对随机变量 $G$ ：

V [G] = E [G^{2}] - E^{2} [G] = E_{X} [E_{Y} [G^{2}]] - (E_{X} [E_{Y} [G]])^{2} = E_{X} [E_{Y} [G^{2}] - E_{Y}^{2} [G]] + (E_{X} [E_{Y}^{2} [G]] - (E_{X} [E_{Y} [G]])^{2}) = E_{X} [V_{Y} [G]] + V_{X} [E_{Y} [G]]

于是 $V [F] - V [F^{'}] = E_{X} [V_{Y} [F]]$ ，可见这一变换消去了 $Y$ 带来的方差。

重要性采样

如果采样所使用的概率密度函数和被积函数只差一个常数因子，那么方差甚至可以是0。然而要做到这一点需要事先知道积分结果，因此不具有实践意义。尽管如此，概率密度函数和被积函数形态还是越相似越好。一种常见的方法是随便拿什么拟合一下被积函数（比方说，分段线性），然后用归一化的拟合函数作为采样的概率密度。

Control Variates

这个翻译过来就是控制变量？总觉得怪怪的。若能找到一个形态与被积函数 $f$ 相似但可以积出解析结果的函数 $g$ ，则可以把原积分改写为：

I = \int_{Ω} f (x) d μ (x) = \int_{Ω} g (x) d μ (x) + \int_{Ω} (f (x) - g (x)) d μ (x)

估计量也就变成了：

F = \int_{Ω} g (x) d μ (x) + \frac{1}{N} i = 1 \sum N \frac{f ( X _{i} ) - g ( X _{i} )}{p ( X _{i} )}

这或许能有效地降低方差：

V [\frac{f ( X _{i} ) - g ( X _{i} )}{p ( X _{i} )}] \leq V [\frac{f ( X _{i} )}{p ( X _{i} )}]

事实上这个方法和重要性采样作用的途径相仿，因此在有了一个之后再用另一个往往效果不佳。Control variates可能在 $f$ 严格非负的情况下带来负样本值，还有一些别的毛病，因此在图形学中没什么存在感。

分块采样

将积分域分成多块，每块分别估值，最后合起来。这东西不适合被积函数有奇异点或变化剧烈的情形，因而也在图形学中没什么存在感。

Quasi-Monte Carlo

想办法把采样点均匀地散布开来，这就是耳熟能详的QMC方法。

首先定义一个衡量采样方案好坏的标准：差异度（Discrepancy）。令：

P B^{*} = {x_{1}, \dots, x_{N}} \subset [0, 1]^{s} = {[0, u_{1}] \times \dots [0, u_{s}] ∣ 0 \leq u_{i} \leq 1}

$P$ 是一堆采样点， $B$ 则是一堆有个角位于原点的轴对齐盒子。在（不存在的）理想情况下，对每个体积为 $λ (B)$ 的 $B \in B$ ， $B$ 中都恰好包含着 $λ (B) N$ 个点。于是，可以把这个理想和现实的差距定义为差异度：

D_{N}^{*} (P) = B \in B^{*} sup \frac{# { P \cap B }}{N} - λ (B)

低差异度序列：称点列 $x_{1}, x_{2}, \dots$ 是一个低差异度序列，当且仅当对其任意前缀 $P = {x_{1}, \dots, x_{N}}$ 均有：

D_{N}^{*} (P) = O (\frac{lo g ^{s} N}{N})

可以证明使用低差异度序列的蒙特卡罗积分的误差为 $O (lo g^{s} N / N)$ 。尽管这看起来并不是很突出（ $lo g^{s} N$ 可能会很大，只有在 $N$ 极大的时候才有优势），但QMC在实践中表现非常优秀，这大概是因为图形学的被积函数对QMC而言性质良好。

Halton序列：设 $i$ 的 $b$ 进制展开形式为

i = k \geq 0 \sum d_{i, k} b^{k}

将 $ϕ_{b} (i)$ 定义为

ϕ_{b} (i) = k \geq 0 \sum d_{i, k} b^{- 1 - k}

称序列：

{x_{i}} = {(ϕ_{b_{1}} (i), ϕ_{b_{2}} (i), \dots, ϕ_{b_{s}} (i))}

为Halton序列，其中 $b_{1}, \dots, b_{s}$ 两两互质（最常用的就是前 $s$ 个素数）。特别地，若 $N$ 在一开始就是已知的，则点集：

{x_{i}} = {(i / N, ϕ_{2} (i), ϕ_{3} (i), \dots, ϕ_{p_{s - 1}} (i))}

被称做Hammersley点集，其差异度为 $O (lo g^{s - 1} N / N)$ 。

Automorphism