蒙特卡洛積分——指定pdf非均勻采樣

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

均勻采樣

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

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

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

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

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

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

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

  如何理解呢，\(b-a\) 還是寬，剩余部分，則是高的平均。

  當N趨於無窮時，根據大數定律，樣本的平均值，會無限趨近於期望值。即矩形的高無限趨近於其期望。

  但期望，等於真實值嗎，下面可以證明：

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

普通蒙特卡洛積分公式

\[\int_{a}^{b}f(x)dx=(b-a)\frac{1}{N}\sum_{i=1}^{N}f(X_i)\quad\quad X_i\sim{U(a,b)} \]

非均勻采樣

蒙特卡洛積分公式

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

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

  蒙特卡洛的積分形式為：

\[\int_{a}^{b}f(x)dx=\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}\quad\quad X_i\sim{p(x)} \]

  前面均勻采樣的公式，即是將\(p(x)=\frac{1}{b-a}\)代入。

幾何理解

  如何理解上述公式呢？

  我們要求\(f(x)\)在（a,b）區間內的積分,可以將整個面積划分成一個個寬度為\(△x\)，高度為\(f(x_{0})\)的小面積區域。

  假設我們已知該面積區域占總面積的比率為\(p(x_{0})\)，即

\[p(x)=\frac{S_{\Delta x}}{S_{all}} \]

  則總面積，即積分的值為

\[\int_{a}^{b}f(x)dx=\frac{S_{\Delta x}}{p(x)} \]

  其中\(p(x)\)被稱為概率密度函數。

  當然，實際某一塊面積的\(p(x)\)難以保證完全的精確，因此我們可以用多塊小面積，估計出的總面積取平均值，以達到更精確的結果，即有：

\[F_{N}=\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_{i})}{p(X_{i})} \]

  \(F_{N}\)即為\(\int_{a}^{b}f(x)dx\)的蒙特卡羅估計值。

  按照上述幾何理解，我們很容易發現，在\(p(x)\)的取值趨近與\(f(x)\)時，我們的采樣越容易接近真實值，更容易用較少的樣本收斂。

  因此，雖然理論上可以用任何pdf計算蒙特卡羅積分，但pdf的選取是有最優解的(\(cf(x)\))，這便是重要性采樣（Importance Sampling）。

代碼實現

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

#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

參考資料：https://zhuanlan.zhihu.com/p/61611088