蒙特卡罗模拟中的采样方法

背景知识

随机模拟也可以叫做蒙特卡罗模拟(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis，在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法,并在最早的计算机上进行编程实现。

蒙特卡洛数值积分

如果我们要求𝑓(𝑥)的积分，可化成 $\int_a^b f(x)=\int_a^b \frac{f(x)}{p(x)}p(x)dx$ 记 $g(x)=\frac{f(x)}{p(x)}$ ， $p(x)$ 为x的概率密度函数，在[a,b]区间内取n个样本[x1,x2,x3....xn]，则
$\int_a^b f(x)=\int_a^b \frac{f(x)}{p(x)}p(x)dx=\frac{g(x_1)+g(x_2)+...+g(x_n)}{n}$
蒙特卡洛数值积分将积分的计算转化为期望的计算，进而可以用统计量来得到。
然而有一个重要的问题——怎么在给定的区间内采样？

采样就是给定一个概率分布，根据该分布从样本空间中生成样本。

比如你可以通过抛硬币近似完成对 $p=\frac12$ 的0-1分布的采样。再举个最简单的通过计算机程序对低维离散样本空间进行采样的例子，比如一维变量X取值的样本空间为{1,2,3}，且取1,2,3三个值的概率分别为 $\frac12,\frac14,\frac14$ 。那么这时候你要如何从这个分布中进行采样？我想大多数人自己都写过这样简单的程序，首先根据各离散取值的概率大小对[0,1]区间进行等比例划分，如划分为[0,0.5],[0,5,0.75],[0.75,1]这三个区间，再通过计算机产生[0,1]之间的伪随机数，根据伪随机数的落点即可完成一次采样。

那么问题来了，如果我们要对一个连续分布（即给定一个已知的概率密度函数 $p(x)$ ）进行采样，那么上述对[0,1]区间划分的方式显然就失效了（当然，最简单的方法，可能会有人考虑将概率密度函数均匀分段积分，然后继续采用之前的做法，不过这样永远只能得到近似的采样结果，永远不可能得到原始分布的采样结果，并且高维情况下积分的计算代价以及是否可积本身也是个问题）。所以，在给定概率密度函数的连续分布下要如何采样呢？

对于概率密度函数的采样，相信一些人可以很直观的想到可以利用累积分布函数（P(x<t)，即 $\int_{-\infty}^tp(x)dx$ ，即在[0,1]间随机生成一个数a，然后求使得P(x<t)=a成立的t，t即可以视作从该分部中得到的一个采样结果。

采样方法

（0）Box-Muller方法

用平均分布的一对随机变量生成高斯分布的随机变量,如果随机变量 U1,U2 独立且U1,U2∼Uniform[0,1]
$Z_0=\sqrt{-2lnU_1}cos(2\pi U_2), Z_1=\sqrt{-2lnU_1}cos(2\pi U_2)$
$Z_0,Z_1$ 满足高斯分布。

但是这是基于累积分布函数可积，且P(x<t)=a可解，即累积分布函数具有反函数的条件下的。假如累积分布函数没有反函数呢？

（1）拒绝采样

很多实际问题中，既然 p(x) 太复杂在程序中没法直接采样，那么我们可以设定一个程序可抽样的分布 q(x) 比如高斯分布，然后按照一定的方法拒绝某些样本，达到接近 p(x) 分布的目的，其中q(x)叫做 proposal distribution 。

拒绝采样

具体操作如下，设定一个方便抽样的函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。

x 轴方向：从 q(x) 分布抽样得到 a。
y 轴方向：从均匀分布（0, kq(a)) 中抽样得到 u。
如果刚好落到灰色区域： u > p(a), 拒绝，否则接受这次抽样。
重复以上过程

问题：在高维的情况下，Rejection Sampling 会出现两个问题，第一是合适的 q 分布比较难以找到，第二是很难确定一个合理的 k 值。这两个问题会导致拒绝率很高，无用计算增加。

（2）MCMC采样

从马尔可夫链到PageRank再到RWR
马尔可夫链蒙特卡洛(MCMC)采样

最后编辑于：2020.01.14 20:19:59