【光线追踪系列十四】蒙特卡洛积分与重要性采样

最新推荐文章于 2024-01-05 16:37:47 发布

沉默的舞台剧

最新推荐文章于 2024-01-05 16:37:47 发布

阅读量2.1k

点赞数 1

分类专栏：光线追踪文章标签：光线追踪

本文链接：https://blog.csdn.net/qq_35312463/article/details/116754314

版权

光线追踪专栏收录该内容

18 篇文章 18 订阅

订阅专栏

本文主要参照 Ray Tracing: The Rest of Your Life，其中只是主要精炼光追相关理论，具体实现可参照原文。
在这里插入图片描述
什么是蒙特卡洛积分？简而言之就是，在求积分时，如果找不到被积函数的原函数，那么利用经典积分方法是得不到积分结果的，但是蒙特卡洛积分方法告诉我们，利用一个随机变量对被积函数进行采样，并将采样值进行一定的处理，那么当采样数量很高时，得到的结果可以很好的近似原积分的结果。这样一来，我们就不用去求原函数的形式，就能求得积分的近似结果。

一、概率论基础

1.1 概率质量函数（PMF）

如果随机变量X只可能取有限个或至多可列个值，则称X为离散型随机变量。例如：
在这里插入图片描述
概率质量函数（probability mass function，简写为PMF），是离散随机变量在各特定取值上的概率。

PMF是对离散随机变量定义的，本身代表该值的概率。见下图：
在这里插入图片描述

1.2 概率密度函数（PDF）

连续型随机变量是指如果随机变量X的所有可能取值不可以逐个列举出来，而是取数轴上某一区间内的任一点的随机变量。

概率密度函数（probability distribution function，简写为PDF) 定义了一个连续的概率分布，例如正态分布函数。

PDF是对连续随机变量定义的，本身不是概率，只有对连续随机变量的PDF在某区间内进行积分后才是概率。它描述可能性的变化情况，如正态分布密度函数，在中间出现的情况最大，两端出现的情况较小。见下图：
在这里插入图片描述
PDF可以求一个区间内的概率：

PDF最重要的特点就是：整个概率区间的概率之和必定为1，即：
在这里插入图片描述

1.3 累计分布函数（CDF）

CDF是累计分布函数，描述发生某事件概率。

CDF将PDF的当前样本及其之前的所有样本都累加了起来。因为PDF是一个连续函数，所以CDF也是一个连续的函数。

下图以一个标准正态分布函数为例解释PDF与CDF之间关系：
在这里插入图片描述
而其CDF函数分布为：

不难发现其二者之间的关系为：

CDF是PDF的积分，积分区间为 −∞ 到 x
PDF是CDF的导数

1.4 连续型随机变量的期望与方差

在离散型随机变量中，我们通过分布列，求得加权的均值，即获得了离散型随机变量的期望。

而在连续性随机变量的场景下，假设一连续型随机变量X的样本空间为D，其概率密度分布函数为p(x)，则其数学期望为：
在这里插入图片描述
若另一连续随机变量Y满足Y=f(X)，则Y的数学期望E[Y]可由下式给出：

方差则是随机变量到期望的距离平方的期望，公式表达如下：

二、蒙特卡洛积分

2.1 简介

假设被积函数f(x)图像如下，积分区间为[a,b]，我们要求的积分为：
在这里插入图片描述
图示为：

(b−a)f(x)几何上可表示为一个宽为(b−a)，高为f(x)的矩形面积：

在这里插入图片描述
下图中，假设我们采样4次，则相当于计算了4个矩形的面积，并求平均值：

定义：
在这里插入图片描述

变形得：
在这里插入图片描述
其中， ⟨F^N⟩ 为样本数量为N时的F的平均值，也就是积分结果的近似值，以上公式叫做Basic Monte Carlo Estimator（基本蒙特卡洛估算器）。

由于随机数是等概率均匀产生的，上述方程的pdf为
在这里插入图片描述
上述公式仅在随机数是等概率均匀生成的情况下才生效。

如果将蒙特卡洛估算器扩展到任意的PDF，一个更加通用的⟨F^N⟩ 可以写成：

在这里插入图片描述
除以PDF，是为了抵消掉不均匀带来的影响。例如某个区域会有很多样本，那就除以该位置的PDF，以削弱其权重。

综上，使用蒙特卡洛算法求解积分的步骤如下：

找到在区间[a,b]内被积函数f(x)；
在区间[a,b]内挑选一个合适的非零的概率密度函数p；
生成一堆满足上述概率密度函数p分布的随机数r，从而得到一堆f®/p®，求这堆f®/p®的平均值。

2.2 推导

根据以上叙述，假设这里我们要计算一个一维积分式
在这里插入图片描述
根据经典的方法，我们需要求得f(x)的原函数F(x)，才能解出这个积分结果，但如果f(x)的原函数形式复杂，或者根本求不出来，总之在不知道F(x)的具体形式的情况下，如果我们还想计算这个积分，怎么办？这时候我们就需要借助蒙特卡洛积分(Monte Carlo Integration)方法。蒙特卡洛积分方法告诉我们，为求得积分结果，可以构造：
在这里插入图片描述
其中的每一个Xi(i=1,2,3,…,N)为[a,b]之间的均匀连续随机变量，所有的Xi组成一个随机变量集合。下面证明，FN的数学期望即为我们要求的结果A。

在这里插入图片描述
以上证明过程表明，若我们根据公式来构造一个新的随机变量FN，则FN的期望就是积分∫f(x)dx的结果，随着N的增加，FN就越逼近理论上A的值，即FN是A的一个无偏估计。这样我们实际上就避开了求f(x)的原函数F(x)的过程。

进一步对上述过程分析，我们发现这里的X被规定为与原积分区间相同的均匀分布随机变量。那么对于与原积分区间相同，但却不是均匀分布的一般随机变量，上述结论是否仍成立？结论是，如果这里的随机变量X的概率密度分布函数已知，那么上述结论还是成立的。假设其概率密度分布函数为p(x)，证明如下：

首先，这次我们构造：
在这里插入图片描述
再构造随机变量

在这里插入图片描述
首次构造的式子中所有量都是已知的。则FN的数学期望E[FN]即为

在这里插入图片描述
与之前的要求不同，这里仅要求随机变量X的概率密度分布函数p(x)已知且在X的样本空间内p(x)≠0。综合上述叙述，我们可以得到蒙特卡洛积分方法如下：

至此，蒙特卡洛积分推到结束。

三、重要性采样

3.1 方差

现在我们分析一下随着样本数量N的增加，估计值FN的方差σ2[FN]的变化情况，以便得出蒙特卡洛积分方法的收敛速度特性。在FN的数学期望E[FN]公式的基础上，我们继续计算σ2[FN]如下：
在这里插入图片描述
故，可得：

在这里插入图片描述

3.1 结论

上式告诉我们，估计值的不稳定来源于随机变量Y的取值不稳定。换句话说，如果f(Xi)/p(Xi)因不同Xi的取值变化地越剧烈，就会造成Y的方差较大，则会造成估计值的收敛速度越慢。这启示我们，若p(x)的形状越接近f(x)，则有益于最终结果的收敛。

上述思想即为“重要性采样”方法，即对积分值有重要贡献(f(x)较大)的被积函数区间，我们以较大概率生成处于这个区间附近的随机变量，用于快速逼近理论值。这也就是为什么我们要引入任意分布随机变量的蒙特卡洛积分方法，而不满足于利用均匀分布随机变量来求蒙特卡洛积分的原因。

利用蒙特卡洛方法，我们可以得到任意一个积分的结果，但是代价就是我们得不到其理论值，我们得到的只是一个对理论值的估计，估计值与理论值之间的误差可以通过增加样本数来减小，但收敛速率仅为O(√N)，也就是说，若想将误差降为现在的一半，我们需要再多计算4倍的计算量才可以达到。即便如此，原始的蒙特卡洛积分方法也不失为是一种经典有效的方法。

四、用PDF求解蒙特卡洛积分demo

求以下函数积分：

分析：上述要求本质上就是求0到2区间内的函数曲线下方和x轴所围成的面积。

4.1 均匀pdf求解

使用均匀采样的pdf，可得到公式：
在这里插入图片描述
带入具体公式即为：

写成代码如下：

_forceinline double random_double()
{
	int iRand = rand();
	return (iRand %= 10000) / 10000.0;
}

void pdf1()
{
	srand(time(NULL));
	int N = 1000000;
	float sum = 0;
	for (size_t i = 0; i < 10; i++)
	{
		sum = 0;
		for (int i = 0; i < N; i++) {
			float x = 2.0*random_double();
			sum += x * x;
		}
		std::cout << "I =" << 2 * sum / N << "\n";
	}
}

上述代码中，我们通过暴力增加样本数量的方式，使得平均值尽量接近真实答案。但是如果要结果误差尽量小，就需要很多样本，才能够模拟出正确的结果。接下来，问题变成了能否用最少的样本，模拟出最接近真实值的结果。

4.2 非均匀pdf求解

为 f(x)= x^2选择一个PDF，例如取 p ( r ) = r / 2，则可以用CDF反函数的方式来生成满足PDF概率分布的随机数。
根据PDF为 p ( r ) = r / 2，CDF为：
在这里插入图片描述
如果x是均匀生成的[0,1]区间内的随机数，求CDF的反函数，就可得到符合PDF的概率分布的随机数生成函数：

接下来套用蒙特卡洛积分公式：

写成代码如下：

//pdf:r/2
inline float pdf1(float x) {
	return 0.5*x;
}

void pdf2() {
	srand(time(NULL));
	int N = 100;
	float sum = 0;
	for (size_t i = 0; i < 10; i++)
	{
		sum = 0;
		for (int i = 0; i < N; i++) {
			float x = sqrt(4 * random_double());
			sum += x * x / pdf1(x);
		}
		std::cout << "I =" << sum / N << "\n";
	}
}

上述pdf似乎并不能帮我们快速求解积分。接下来我们人为定义一个更加合适的PDF：
在这里插入图片描述
对应有CDF为：

CDF的反函数为：

在这里插入图片描述
这是一个堪称完美的pdf，因为我们已经预先算出了题目中的积分值是8/3，从而使得一次蒙特卡洛积分采样就能给出结果：

//pdf: 3*x*x/8
inline float pdf2(float x) {
	return 3 * x*x / 8;
}

void pdf3() {

	srand(time(NULL));
	int N = 1;
	float sum = 0;
	for (size_t i = 0; i < 10; i++)
	{
		sum = 0;
		for (int i = 0; i < N; i++) {
			float x = pow(8 * random_double(), 1 / 3);
			sum += x * x / pdf2(x);
		}
		std::cout << "I =" << sum / N << "\n";
	}
}

上边三个pdf运行效率及结果如下：
在这里插入图片描述

由此可见选择一个合适的非均匀pdf，使得能够在被积函数较大的地方采样更多的样本，这样能够减少噪声、快速收敛，这也就叫做重要性采样。

沉默的舞台剧

关注

1
点赞
踩
22

收藏

觉得还不错? 一键收藏
2
评论
【光线追踪系列十四】蒙特卡洛积分与重要性采样

本文主要参照 Ray Tracing: The Rest of Your Life，其中只是主要精炼光追相关理论，具体实现可参照原文。困难积分的近似求解方式，我们开始正式求解渲染方程
复制链接

扫一扫

专栏目录