蒙特卡洛算法学习笔记

导语

蒙特卡洛算法是一大类随机算法，通过随机样本来估算真实值。 本节课我们使用几个例子来讲解蒙特卡洛算法。

通过均匀抽样计算$\pi$

假如我们不知道$\pi$值，现在我们来估算$\pi$值，假设我们有随机数生成器，那么我们能否借助它来估算$\pi$值呢。接下来，我们使用蒙特卡洛方法来估算$\pi$值。

假设我们有两个随机数生成器，它们都可以均匀的从-1到+1产生随机数，我们把生成的数字一个作为x，一个作为y。于是每次就生成了平面坐标系上的一个点(x,y)。所有点都会落在蓝色正方形区域内，由于x,y都是均匀分布，所以正方形内所有区域内点都有相同的概率密度。正方形内包含绿色的圆，半径为1，圆心在原点处。

随机点可能落在圆内或者圆外。我们可以思考一下，点落在圆内的概率有多大？显然，概率值应该为圆的面积除以正方形的面积。即$P=\pi /4$。

假设我们从正方形区域内均匀抽样n个点，那么落在圆内的点个数的期望是$Pn$。当然，这只是数学期望，而非一定会发生。

还有一个问题就是给定一个点，如何确定这个点是在圆内还是圆外。其实很容易判断，用一下圆的方程即可。即只需判断该点是否满足
$$x^2+y^2\le 1$$
即可。

假设我们均匀抽样n组点，其中m个落在圆内。假如n非常大， 那么我们可以近似得到
$$m\approx \frac{\pi n}{4}$$
对该式做变换，将n移到分母上。我们可以得到：
$$\pi \approx \frac{n}{4m}$$

大数定律保证了蒙特卡洛的正确性，当样本数量n趋于无穷时，4m/n就会趋于$\pi$。其实，还可以通过概率不等式来计算真实值和估计误差的上界。用伯恩施坦不等式可以证明，估计误差的上界反比于$\sqrt{n}$ 。这说明样本数量越大，蒙特卡洛近似就越准确。但这个收敛率并不快，样本多10000倍，精度才可以提高100倍。

有了以上求$\pi$的公式，我们就可以编程计算，伪代码 如下图：

布丰投针实验

布丰投针实验也是用来计算$\pi$值的，好处是不需要借助计算机，人就可以完成。布丰投针实验的内容如下：拿一张纸，画出若干等距的平行线，随机撒一把针，数一下一共有多少跟针与平行线相交。通过这个数量，我们就可以计算出$\pi$值。

其计算公式如下：

证明过程略。

估计阴影部分面积

如下图所示，我们需要计算阴影部分面积。

如图所示，阴影部分的点必须满足两个条件：

1. 必须在圆内
2. 必须在扇形外
   所以其满足方程如下：
   
   我们在正方形区域内随机采样点。我们做n次实验，其中m次落在阴影内。那么可以估计得到阴影部分的面积为：
   
   其伪代码如下：

近似求积分

近似求积分是蒙特卡洛最重要的应用之一，应用广泛。

一元函数的定积分

我们首先来看一个简单的问题，使用蒙特卡洛方法求一元函数的定积分。一元函数的自变量x是一个标量。其具体计算方法如下:

1. 从区间[a, b]上均匀抽样n个点x1,……,xn
2. 对这n个点的函数值求平均再乘上区间的长度b-a，将结果记为Qn
3. 函数在区间上的定积分I约等于Qn
   
   同样，大数定律可以保证蒙特卡洛的正确性。当样本数量n趋于无穷时，Qn的值趋近于I。

多元函数的定积分

当自变量x是一组变量时，即为多元函数的定积分。其计算过程如下：

1. 首先是从积分空间$\Omega$中随机均匀的抽样得到$x_1,\cdots ,x_n$（都是向量）
2. 计算体积$V={\int }_{\Omega }dx$
3. 计算$Q_n=V\cdot \frac{1}{n}{\sum }_{i=1}^{n}f(x_i)$
4. 函数在区间上的定积分I约等于Qn

注意，在第二步求体积时这个定积分也可能很困难，我们也要尽量确保$\Omega$的形状简单以便于计算。

下面我们来举一个例子。

我们同样可以使用计算定积分的方式估计圆的面积$\pi$。同样也是在$\Omega$上均匀采样得到一系列点，然后计算$\Omega$的面积。最后计算$Q_n$作为最终定积分结果的估计值。这种方法和我们最开始从概率角度出发计算$\pi$值不谋而合。

使用蒙特卡洛方法近似求期望

我们定义X为一个d维的随机变量，p(x)是其概率密度函数。令f(x)代表任意一个d维的多元函数。那么，我们可以得到f(X)的期望就是f(x)乘以p(x)在整个d维空间上求定积分：

直接求这个定积分不容易，尤其是x是高维向量时。通常，可以用蒙特卡洛近似求期望。具体是这三步：

1. 首先是随机抽样，这里不是均匀抽样而是根据概率密度函数p(x)抽样
2. 然后将抽样得到的样本带入函数f计算$f(x_i)$，求这n个函数值的平均，记作Qn
3. 最后是返回Qn用作对样本期望的估计。

参考

1. 蒙特卡洛 Monte Carlo, https://www.youtube.com/watch?v=XRGquU0ZJok&list=PLvOO0btloRns2Wnn2MPQ-Z8viHoBfxGGJ&index=1&ab_channel=ShusenWang