蒙特卡洛積分——非均勻取樣

mirroooor發表於2021-04-17

  為什麼需要蒙特卡洛法積分呢？數學上，積分的解析解，往往需要求出被積分函式的原函式，這對於計算機是相當困難的，因此有了求積分的數值方法。

均勻取樣

  假設我們現在要求\(x^2\)在\([0,2]\)上的積分

效果圖

  如何計算這塊面積呢，不妨將其看成“矩形”進行計算，矩形的寬為2，高為\(x^2\)在\([0,2]\)上的均值。

\[S=2*average(x^{2}) \]

  我們取越多的點來估算均值，獲得的結果也越精確。

  如何嚴格證明我們估算正確性呢？

  下面是我們要估算的真實值，即\(x^2\)在\([0,2]\)上的積分

效果圖

  蒙特卡洛積分表述為以下形式：

效果圖

  則有：

效果圖

  由於估計量的數學期望等於被估計引數的真實值，所以我們的估計是無偏的。

非均勻取樣

  如果不按均勻取樣，而是按\(p(x)\)的概率密度進行取樣，同樣也可以達到效果。

  只是此時，f(x)還需要除以p(x)，相當於出現概率更大的點，計算時賦予的權重就低一點。

  蒙特卡洛的積分形式為：

效果圖

  下面我們來證明其正確性

效果圖

  因此非均勻取樣的估計也是無偏的。

蒙特卡洛積分

  下面是一個蒙特卡洛積分的簡單示例

#include <iostream>
using namespace std;

//被積函式
double test_func(float x){
	return x * x;
}
inline double random() {
	return rand() / (RAND_MAX + 1.0);
}

inline double pdf(double x){
	return 3 * x*x / 8;
}

int main()
{
	int N = 10000;
	double sum = 0;
	for (int i = 0; i < N; i++){
		float x = pow(8 * random(), 1.0 / 3.0); //這裡取CDF的反函式，詳見上一篇文章
		sum += test_func(x) / pdf(x); // 本篇所述，1/pdf(x)的權重
	}
	cout << "I =" << sum / N << endl;
	return 0;
}

//輸出：I =2.66667

相關文章