蒙特卡洛积分——采样方法

蒙特卡洛积分

• 目的： 通过计算机进行采样近似求解复杂的积分
• 理论基础： 大数定律，当$n$足够大时，$X$的均值将收敛于它的期望（不严谨表述）
• 一般形式： $$\theta ={\int }_a^{b}f(x)dx={\int }_a^b\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum _{i=1}^n\frac{f(x_i)}{p(x_i)}$$ $p(x)$为$x$的概率密度，$x_i$为从$x$的概率分布$p(x)$中随机采样得到的第$i$个样本，然而，当$p(x)$是一个比较复杂的概率密度时，从$p(x)$中随机采样并不容易，那么我们现在的问题就转化为了"如何采样出$x$的分布$p(x)$对应的若干个样本"
• 重点强调一下，我们所说的采样指的是对$x$采样，而$p(x)$是$x$的概率密度，影响每一个$x$被采样到的概率。

概率分布采样

• 虽然对于复杂分布$p(x)$而言，采样是一件困难的事情，但在计算机中, 得到一个均匀分布$U(0,1)$上的随机数是很简单的, 因此，我们自然而然的会想要通过这个简单的均匀分布来解决我们的问题。

离散型随机变量情形

• 我们这里用$\theta$表示服从$U(0,1)$的随机变量. 现在我们就探讨如何通过对服从$U(0,1)$的随机变量$\theta$采样，来得到我们想要的服从其它离散分布的随机变量$X$的样本.
• 设有一个离散型随机变量$X\sim \left(\begin{array}{cccc}{x}_1& x_2& \cdots & x_n\\ p_1& p_2& \cdots & p_n\end{array}\right)$，$X$与$\theta$之间存在如下映射关系$X=G(\theta )$
• 由于$\theta$恰好服从$U(0,1)$，因此，我们只需要将$\theta$依概率映射到
  样本空间$X$中每个值即可，即$$X=G(\theta )=\{\begin{array}{cc}{x}_1& 0\le \theta <p_1\\ x_2& p_1\le \theta <p_1+p_2\\ \cdots & \cdots \\ x_n& {\sum }_1^{n-1}{p}_i\le \theta <1\end{array}$$
  其实可以将右边的不等式看作累积分布概率作用下的结果，比如第二行$p_1\le \theta <p_1+p_2$，其实它区间的长度就等于$p_2+p_1-p_1=p_2$，也就是$x_2$对应的概率$p_2$
• 这样就实现了在离散型随机变量的情形之下，找到了服从$U(0,1)$的随机变量$\theta$与$X$之间的映射关系$G$，也就可以通过对$\theta$采样的结果转化为对$X$的采样结果

连续型随机变量情形

• 在连续型随机变量的情形之下，同样的，我们也探讨如何通过一个服从$U(0,1)$的随机变量$\theta$来得到我们想要的服从其它概率分布的随机变量的样本.
• 设有一个连续型随机变量$X$，其概率密度为$f(x)$，累积分布函数为$F(x)$，$X$与$\theta$之间存在如下映射关系$X=G(\theta )$
• 由累积分布函数的定义可得 P { X < a } = F ( a ) P\{X<a\}=F(a) P{X<a}=F(a)将$X=G(\theta )$代入可得 P { G ( θ ) < a } = F ( a ) P\{G(\theta)<a\}=F(a) P{G(θ)<a}=F(a)
• 由累积分布函数的性质可得，$F(x)$为一个单调递增函数，因此，对不等式两边同时取$G$的反函数可得 P { θ < G − 1 ( a ) } = F ( a ) P\{\theta<G^{-1}(a)\}=F(a) P{θ<G−1(a)}=F(a)
• 又已知$\theta$服从$U(0,1)$，因此有 P { θ < b } = b , b ∈ [ 0 , 1 ] P\{\theta<b\}=b,b\in[0,1] P{θ<b}=b,b∈[0,1]
• 令$b=G^{-1}(a)$可得 P { θ < G − 1 ( a ) } = G − 1 ( a ) P\{\theta<G^{-1}(a)\}=G^{-1}(a) P{θ<G−1(a)}=G−1(a)
• 又因为刚才证得 P { θ < G − 1 ( a ) } = F ( a ) P\{\theta<G^{-1}(a)\}=F(a) P{θ<G−1(a)}=F(a)，所以有$$F(a)=G^{-1}(a)$$
• 因此$$G(a)=F^{-1}(a)$$
• 所以$$X=G(\theta )=F^{-1}(\theta )$$
• 根据以上结论，我们可以得出在连续型随机变量的情形下，$\theta$与随机变量$X$的映射关系就是$x$的累积分布函数的反函数

接收—拒绝采样

• 刚刚介绍了概率分布采样方法，对于连续型随机变量的情形而言，我们可以通过使用随机变量$X$的的累积分布函数的反函数，也就是如下式子$$X=G(\theta )=F^{-1}(\theta )$$来将$\theta$转化为$X$，但众所周知，累积分布函数是不好求的，尤其是对于复杂的概率密度$p(x)$而言，难以通过积分得到累积分布函数$F(\theta )$。
• 对于概率分布不是常见的分布的情况，一个可行的办法是牺牲采样的效率，采用接受-拒绝采样来得到该分布的样本。
• 现在假设我们想要得到从一个复杂的分布$p(x)$中进行采样，既然 $p(x)$太复杂在程序中没法直接采样，那么我设定一个便于程序采样的分布 $q(x)$，比如高斯分布，然后从$q(x)$中进行采样，并按照一定的方法拒绝某些不合理的样本，以达到接近$p(x)$分布的目的，其中$q(x)$被称为提议分布
  
• 采样步骤如下：
  • 我们设置一个常量$k$，使得$kq(x)$总是在$p(x)$的上方
  • 首先，我们使用概率分布采样方法对$q(x)$进行采样，得到一个样本$z_0$
  • 然后，从均匀分布$U(0,1)$中随机采样得到一个样本$u_0$
  • 如果$u_0\le \frac{p(z_0)}{kq(z_0)}$，那么就接受$z_0$这个样本，否则拒绝这次采样
  • 重复以上过程，得到$n$个接受的样本$z_0,z_1,...z_{n-1}$
• 直观上理解这个采样过程就是：
  • 从均匀分布$U(0,1)$中生成一个随机概率$u_0$
  • 由于我们需要近似的概率分布$p(x)$是图中的红色线条，对于我们采集到的样本点$z_0$而言，它有$\frac{p(z_0)}{kq(z_0)}$的概率可以被认为是一个合理的样本（$\frac{p(z_0)}{kq(z_0)}$代表的就是在$x=z_0$这条竖直的线上$p(z_0)$与$kq(z_0)$的比值，可以按几何概率的思想来理解）
  • 随机概率$u_0$如果没有超过$\frac{p(z_0)}{kq(z_0)}$这个阈值，则认为这次采样的结果是合理的，并接受$z_0$
• 总之，通过我们通过提议分布$q(x)$，最终采集到的样本为$z_0,z_1,...z_{n-1}$，我们可以通过这组样本，应用蒙特卡洛方法解决最开始的积分问题$$\theta ={\int }_a^{b}f(x)dx={\int }_a^b\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum _{i=1}^n\frac{f(z_i)}{p(z_i)}$$
• 但此采样方法的效率深受提议分布的影响，如果提议分布与原始分布相差巨大，那么采样结果被拒绝的概率就会非常高，导致采样的效率低下，同时，对于一些高维的复杂非常见分布$p(x_1,x_2,...,x_n)$，我们要找到一个合适的$q(x)$和$k$非常困难。

总结

• 当我们需要求解一个复杂积分时，我们可以对此积分稍作变换$${\int }_a^{b}f(x)dx={\int }_a^b\frac{f(x)}{p(x)}p(x)dx$$以便于应用蒙特卡洛方法，因为蒙特卡洛方法最重要的思想便是使用均值来逼近期望，这就要求我们抽取一系列的$x$的样本，然后令$${\int }_a^b\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum _{i=1}^n\frac{f(x_i)}{p(x_i)}$$ 以此来近似积分的结果.
• 在抽取$x$的样本的过程中，又出现了新的问题，对于常见的概率分布，我们可以通过概率分布采样方法进行采样，然而对于不常见的分布而言，我们很难积分得到它的累积分布函数，更不必说它的反函数，因此我们又引入了接收—拒绝采样方法,使用提议分布为中介进行采样.

References

• https://www.cnblogs.com/pinard/p/6625739.html
• https://zhuanlan.zhihu.com/p/191487550
• https://zhuanlan.zhihu.com/p/259280292