0%

蒙特卡洛方法

《全局光照技术》中的蒙特卡洛方法总结

概率论基础

概率密度函数与累计分布函数

随机变量$X$每个$x$值关联一个概率,概率密度函数表示随机变量所有可能值组成的概率分布函数$p(x)$。

累计分布函数定义为所有随机变量值中小于或等于$y$的随机变量概率的积分:
$$
P(y)=\int^y_{-\infty}p(x)dx
$$

期望与方差

离散形式

$$
E[X]=\sum^n_{i=1}p_ix_i
$$

$$
\sigma^2=\sum_i(x_i-E[x])^2p_i
$$

连续形式

$$
E[X]=\int^\infty_{-\infty}xp(x)dx
$$

$$
\sigma^2=\int(x-E[X])^2p(x)dx
$$

大数定律

当独立同分布的随机变量抽样和除以抽样数时,可得到对该随机变量期望的一个估计:
$$
E[x]\approx\bar{X}=\frac{1}{N}\sum^N_{i=1}x_i
$$
当$N\to\infty$时,随机变量统计平均值接近期望:
$$
P{E[X]=\lim_{N\to\infty}\frac{1}{N}\sum^N_{i=1}x_i}=1
$$

蒙特卡洛方法

在渲染方程中的被积函数往往比较复杂,不能简单求出数值解。蒙特卡洛积分的目的是对积分进行近似计算。

由大数定律可知:
$$
E[f(x)]=\int_{x\in S}f(x)p(x)d\mu\approx\frac{1}{N}\sum^N_{i=1}f(x_i)
$$
令$g(x)=f(x)p(x)$,上述期望近似变为积分$I$的估值$I_N$:
$$
\int_{x\in S}g(x)d\mu\approx\frac{1}{N}\sum^N_{i=1}\frac{g(x_i)}{p(x_i)}
$$
该估值的期望为:
$$
E[I_N]=E[\frac{1}{N}\sum^N_{i=1}\frac{g(x_i)}{p(x_i)}]\
=\frac{1}{N}\sum^N_{i=1}E[\frac{g(x_i)}{p(x_i)}]\
\frac{1}{N}N\int\frac{g(x)}{p(x)}p(x)dx\
=\int g(x)dx\
=I
$$
由此可得该估计可以用作$I$的近似,该估计的方差为
$$
\sigma^2=\frac{1}{N}\int(\frac{g(x)}{p(x)}-I)^2p(x)dx
$$
可见方差随N的增加线性降低,又由于误差与标准差成正比,因此误差随着$\sqrt{N}$线性降低。理论上蒙特卡洛方法没有对$p(x)$进行要求,但由上式可知$p(x)$分布越接近$g(x)$误差会更小。

综上,蒙特卡洛方法归纳为以下步骤:

  • 对满足某一概率分布的随机数进行抽样
  • 使用抽样值计算$\frac{g(x_i)}{p(x_i)}$,称为样本贡献值
  • 对所有的抽样结果求平均值得到积分的估计

抽样方法

抽样过程可以描述为从$p(x)$对应的随机变量$X$中产生一系列随机数$X_1,X_2,…$使得对于任意$\Omega$属于定义空间$\Omega_0$满足:
$$
P{X_k\in\Omega}=\int_\Omega p(x)dx \le1
$$
实际上不能直接从$p(x)$进行抽样,计算机进行抽样过程时首先求取某些基础随机数序列。通常使用均匀分布的随机数,由这个序列作为抽样所需的基础随机数对于概率分布函数确定的随机变量,有逆变换算法和取舍算法。

逆变换算法

算法步骤:

  1. 首先计算$p(x)$的累计分布函数:$P(x)=\int^x_0p(x’)dx’$
  2. 计算累计分布函数的反函数:$P^{-1}(x)$
  3. 从$[0,1]$上的均匀分布产生一个随机数$\xi$
  4. 最后求出满足$p(x)$分布的随机数:$X_i=P^{-1}(\xi)$

取舍算法

考虑一个限制在$[a,b]$上的分布函数$f(x)$,其最大值为$c$,其他区间上为0。

算法步骤:

  1. 产生一个$[a,b]$上的均匀分布随机数:$X\sim U(a,b)$。
  2. 产生一个$[0,c]$上独立于$X$的均匀分布随机数:$Y\sim U(0,c)$。
  3. 如果$Y\le f(X)$则接受,接受的随机数为$Z=X$,否则拒绝返回第一步。

基于马尔可夫链的梅特罗波利算法

此方法是对于随机变量维度高或是概率分布不完全已知的概率分布时使用的方法。

详见另一篇博客专门讲解,《梅特罗波利斯算法》

方差缩减

以下为图形学中常用的方差缩减方法。

复合重要性抽样

复合重要性抽样是从多个不同的分布中抽样,然后对这些不同的抽样结果进行加权组合。其估计表达式为:
$$
I_N=\sum^n_{i=1}\frac{1}{n_j}\sum^{n_i}_{j=1}\omega_i(X_{x,j})\frac{X_{i,j}}{p_i(X_{i,j})}
$$
它将整个空间按一定的形式划分,然后对每个划分的区域使用不同的抽样技术$\sum^n_{i=1}$就是不同抽样技术的和;$\sum^{n_i}_{j=1}$是每种抽样技术在所对应空间中对$I_N$的估计,每种技术只贡献部分值故乘上加权系数$\omega_i(X_{x,j})$;$\omega_i(x)$满足$\sum^n_{i=1}\omega_i(x)=1$。

一种加权系数函数–平衡启发式:
$$
\hat{\omega_i}=\frac{c_ip_i(x)}{\sum_jc_jp_j(x)}
$$
其中$c_i$是每个抽样分布$p_i$对应的抽样数量比例,$c_i=\frac{n_i}{N}$。

这种加权系数函数对应的积分估计为:
$$
I_N=\frac{1}{N}\sum^n_{i=1}\sum^{n_j}{j=1}\frac{f(X{i,j})}{\bar{P}(X_{i,j})}\
\bar{P}(x)=\sum^n_{i=1}c_ip_i(x)
$$

分层抽样

分层抽样将被积函数定义域$\Omega$划分为多个彼此不相交的子定义域$\Omega_1,\Omega_2,…$,满足$\bigcup^n_{i=1}\Omega_i=\Omega$。每个阶层内运用不同的抽样分布进行抽样,每个阶层的蒙特卡洛估计为:

$$
F_i=\frac{1}{n_i}\sum^{n_i}{j=1}\frac{f(X{i,j})}{p_i(X_{i,j})}
$$

因此,整个被积函数的估计为:
$$
I_N=\sum^n_{i=1}v_iF_i
$$
其中$i$为每个阶层占整个定义域空间的体积比例。