《学一辈子光线追踪》二点五一维MC积分（二）

118 阅读 0 评论 78 点赞

我是靠谱客的博主标致胡萝卜，这篇文章主要介绍《学一辈子光线追踪》二点五一维MC积分（二），现在分享给大家，希望可以做个参考。

蒙特卡洛光线追踪技术系列见蒙特卡洛光线追踪技术

我们现在可以对以前的积分进行采样

I = integral(x^2, 0, 2)

我们需要解释x的pdf的不均匀性。如果我们在某个位置的样本太多，我们应该降低其权重。

比如我们之前的例子，分两段，要计算从0-2的平均值，从1-2产生了150个数据，从0-1产生了50个数据，则最后的结果就是需要用1-2的150个数据的和除以150，再加上0-1的和除以50，而不是直接用所有的和除以200

pdf 是一个很好的度量方法，可以用来衡量在不同位置抽样的多少。所以权重函数应该与 1/pdf 成正比。实际上就是1/pdf：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "time.h"

double myRandom() {
	return rand() / (RAND_MAX + 1.0);
}
inline float pdf(float x) {
	return 0.5*x;
}
int main() {
	int N = 1000000;
	double sum = 0.0;
	for (int i = 0;i < N;i++) {
		double x = sqrt(4 * myRandom());
		if (x == 0)printf("errorn");
		sum += x*x / pdf(x);
	}
	printf("I = %lf n",sum/N);
	system("pause");
}

结果打印了一堆error，表示有的x是0，所以还是得修改一下原程序：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int main() {
	int N = 1000000;
	double sum = 0.0;
	for (int i = 0;i < N;i++) {
		double x = sqrt(4 * myRandom());
		if (x == 0) { 
			printf("errorn"); 
			
		}
		else {
			sum += x*x / pdf(x);
		}
	}
	printf("I = %lf n",sum/N);
	system("pause");
}

打印结果为：

说明估计是很准的。

由于我们在被积函数较大的地方采样更多，我们就可能得到更少的噪声，从而更快地收敛。这就是为什么使用精心选择的非均匀pdf通常被称为重要性抽样。

如果我们对相同的代码取相同的样本，使 pdf = ½ 在[0,2]范围内，我们可以使用机器得到x=2*drand48()，代码是：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
inline float pdf(float x) {
	return 0.5;
}
int main() {
	int N = 1000000;
	float sum = 0.0;
	for (int i = 0;i < N;i++) {
		float x = 2 * myRandom();
		if (x == 0) {
			printf("errorn");

		}
		else {
			sum += x*x / pdf(x);
		}
	}
	printf("I = %lf n", sum / N);
	system("pause");
}

请注意，我们不再需要2*sum/N中的2了-这是由pdf处理的，当您除以它时，pdf是0.5。你会注意到重要性抽样有点帮助，但不是很大。我们可以使pdf完全遵循被积函数：

p(x) = (3/8)x^2

我们得到相应的：

P(x) = (⅛)x^3

以及

P^-1(x) = pow(8x,⅓)

只有当我们已经知道答案（通过解析积分p得到P）时，这种完美的重要性抽样才是可能的，但是这是确保代码工作的一个很好的练习。

我们只需要一个样品：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
inline float pdf(float x) {
	return 3*x*x/8;
}
int main() {
	int N = 1;
	float sum = 0.0;
	for (int i = 0;i < N;i++) {
		float x = pow(8*myRandom(),1./3.);
		if (x == 0) {
			printf("errorn");

		}
		else {
			sum += x*x / pdf(x);
		}
	}
	printf("I = %lf n", sum / N);
	system("pause");
}