对随机变量的采样 | Blurred code

对随机变量的采样

2020/05/23

Updated:2020/05/23

Categories: CG

首先,先上Rendering Equation:

\[ L_o(p,\omega_o) = L_e(p,\omega_o) + \int_{\Omega +}L_i(p,\omega _i)f_r(p,\omega_i,\omega_o)(n \cdot \omega_i) \text{d}\omega_i \]

其中包含有含有递归定义的L项和对半球面的积分项,这个式子是没有解析解的,数值方法可以使用蒙特卡罗积分,而蒙特卡罗积分就涉及到对随机变量的采样.

逆采样方法

如果\(X\)是连续型随机变量,其概率密度函数\(p(x)\)已知, 则在定义域内可以积分得到其累积分布函数(CDF)为:

\[ P(x) = \int_0^x p(x')\text{d}x' \]

计算其反函数,令一随机变量\(u \in [0,1)\),则有

\[ P(P^{-1}(x)) = u \]

解可得:

\[ P^{-1}(x) = g(u) \]

其中\(P^{-1}(x)\)即为服从\(P(x)\)分布的样本, 而这个样本可以从均匀分布的\(u\)中导出.

一维随机变量的采样

设概率密度\(p(x) = 3x^2, x\in(0,1)\), 则可以得到其分布函数为\(P(x) = x^3, x\in(0,1)\). 令

\[ P(P^{-1}(x)) = u, u\in[0,1)\\ \rArr (P^{-1}(x))^3 = u, u\in[0,1]\\ \rArr X = P^{-1}(x) = u^{\frac{1}{3}}, u\in[0,1] \]

则可以通过对均匀分布的u采样(取(0,1)之间的随机数)来得到X的样本,而X的样本服从\(P(x)\)分布.

不同分布的转换

一维分布的转换

假设我们已经有了一个随机变量\(X\)且知道它的概率密度函数(PDF)\(p_x(x)\),那么如果我们再知道\(y=f(x)\)的映射,那么我们可以推断出\(Y\)变量的概率密度以及分布.

已知\(Y_i=f(X_i)\),则可以得出

\[ P(Y\le f(x)) = P(X \le x) \\ \rArr CDF: F(y) = F(f(x)) = F(X)\\ \text{两边对x微分求x的概率密度} \\ \rArr p(y)\frac{\text{d}y}{\text{d}x} = p(x) \text{移动微分到右边去,因为概率密度恒>0,所以加上绝对值} \\ \rArr p(y)= |\frac{\text{d}y}{\text{d}x}|^{-1} p(x) \]

多维分布的转换

同一维的算法是一样的,只是导数变成了偏微分,并且成为了矩阵形式,这个矩阵被称为Jacobian Matrix.

已知\(\bold{Y}=T(\bold{X}),\bold{X}=\{x_1,x_2,...,x_n\}\)

\[ p(y)= \frac{p(x)}{|J_T(x)|} \]

其中\(|J_T(x)|\)是Jaccobian矩阵的行列式的绝对值. 而Jaccobian矩阵为

\[ J_T(x) = \begin{pmatrix} \frac{\partial T_1}{\partial X_1} & \cdots & \frac{\partial T_1}{\partial X_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial T_n}{\partial X_1} & \cdots & \frac{\partial T_n}{\partial X_n} \end{pmatrix} \]

举个例子,已知三维坐标系下,球面的坐标为

\[ x = r\sin\theta\cos\phi\\ y = r\sin\theta\sin\phi\\ z = r\cos\theta \]

\[ J_T(r,\theta, \phi) = \begin{pmatrix} \sin\theta\cos\phi& r\cos\theta\cos\phi & -r\sin\theta\sin\phi\\ \sin\theta\sin\phi& r\cos\theta\sin\phi & r\sin\theta\cos\phi\\ \cos\theta& -r\sin\theta & 0 \end{pmatrix} \]

该行列式绝对值为\(|r^2\sin\theta|\),所以

\[ p(x,y,z) = \frac{p(r,\theta,\phi)}{|r^2\sin\theta|} \]

以上.

以下为球面几何的一些东西

假设在单位球体内,r=1,维度从\((r,\theta,\phi)\)降到二维\((\theta,\phi)\),又已知球面的立体角(solid Angle)微分为

\[ \text{d}\omega = \sin\theta\text{d}\theta\text{d}\phi \]

对球面做积分,积分的结果是点在球面上的分布函数,分别用两种不同形式的表示

\[ \int_{\Omega}p(\theta,\phi)d\theta d\phi= \int_{\Omega}p(\omega)d\omega \\ \rArr p(\theta,\phi)d\theta d\phi = p(\omega)d\omega\\ \rArr p(\theta,\phi) = \sin\theta p(\omega) \]

二维随机变量的采样

上面的例子已经介绍了球面上(x,y,z)和\((r,\theta,\phi)\)以及立体角\(\omega\)之间的概率密度函数之间的转换关系.假设\(r=1\)可以帮助我们把问题降到2维来, 问题转换成了,已知一个联合密度函数\(p(x,y)\),如何从中采样X,Y值,使X,Y样本服从该联合密度分布.

继续发挥基本的概率论知识,如果一个联合密度函数中,x,y变量是独立分布的,那么

\[ p(x,y) = p(x)p(y) \]

分别独立对X,Y采样即可. 然而这种情况很少见,不具备普遍性.

从联合密度函数中,我们可以算出边缘密度函数,即

\[ p(x) = \int p(x,y) \text{d}y \]

这样我们就得到了p(x)的单独分布. 似乎可以用同样的方法去求得p(y)的分布,然而概率论的知识告诉我们,以这样的方法求得的p(x)p(y),在x,y不独立的时候,并不等于p(x,y).

我们真正想要的结果是 \(p(y|x)\),y的条件分布.

\[ p(y|x) = \frac{p(x,y)}{p(x)} \]

以上.

继续对球面的采样,假设点在球面上均匀采样,则其立体角的密度函数应该等于一个常数c,对其进行积分

\[ \int_{\Omega}p(\omega)d\omega = \int_{\Omega}c d\omega = 1 \\ \rArr c = \frac{1}{4\pi} \\ \rArr p(\omega) = \frac{1}{4\pi} \\ \rArr p(\theta,\phi) = \frac{\sin\theta}{4\pi} \]

\(p(\theta,\phi)\)做积分

\[ p(\theta) = \int_0^{2\pi} p(\theta,\phi)d\phi = \frac{\sin\theta}{2} \\ p(\phi | \theta) = \frac{p(\theta,\phi)}{p(\theta)} = \frac{1}{2\pi} \]

求分布函数,

\[ P(\theta) = \int_0^\theta \frac{\sin\theta}{2}d\theta = \frac{1 - \cos\theta}{2} \\ P(\phi |\theta) = \frac{\phi}{2\pi} \]

逆采样,

\[ \phi = 2\pi u_1 \\ \theta = \arccos (1 -2u_2) \]

代入下式,即可采样得到(x,y,z)的值

\[ x = r\sin\theta\cos\phi\\ y = r\sin\theta\sin\phi\\ z = r\cos\theta \]