蒙特卡洛积分——指定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