机器学习小组知识点30：接受-拒绝采样（Acceptance-Rejection Sampling）

蒙特卡洛数值积分

如果我们要求f(x)的积分，如

∫abf(x)dx

而f(x)的形式比较复杂积分不好求，则可以通过数值解法来求近似的结果。
常用的方法是蒙特卡洛积分：

∫abf(x)q(x)q(x)dx

这样把q(x)看做是x在区间内的概率分布，而把前面的分数部门看做一个函数，然后在q(x)下抽取n个样本，当n足够大时，可以用采用均值来近似：

1n∑if(xi)q(xi)

因此只要q(x)比较容易采到数据样本就行了。随机模拟方法的核心就是如何对一个概率分布得到样本，即抽样（sampling）。

均匀分布，Box-Muller 变换

在计算机中生成[0,1]之间的伪随机数序列，就可以看成是一种均匀分布。而随机数生成方法有很多，最简单的如：

xn+1=(axn+c) mod m

当然计算机产生的随机数都是伪随机数。

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于Uniform(0,1) 的样本生成。例如正态分布可以通过著名的Box−Muller变换得到

[Box-Muller 变换]

如果随机变量 U1,U2 独立且U1,U2∼Uniform[0,1]，

Z0=−2lnU1−−−−−−√cos(2πU2)

Z1=−2lnU1−−−−−−√sin(2πU2)

则 Z0,Z1独立且服从标准正态分布。

其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来，大家可以参考统计计算的书，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。

Monte Carlo principle

Monte Carlo 抽样计算随即变量的期望值是接下来内容的重点：X表示随即变量，服从概率分布 p(x), 那么要计算 f(x) 的期望，只需要我们不停从 p(x)中抽样xi，然后对这些f（xi）取平均即可近似f(x)的期望。

EN(f)=1N∑i=1Nf(xi)

接受-拒绝抽样

很多实际问题中，p(x)是很难直接采样的的，因此，我们需要求助其他的手段来采样。既然p(x) 太复杂在程序中没法直接采样，那么我设定一个程序可抽样的分布 q(x)比如高斯分布，然后按照一定的方法拒绝某些样本，达到接近p(x)分布的目的，其中q(x)叫做 proposal distribution，但是要满足p(x)≤Mq(x)，其中M<∞。

然后重复如下步骤，直到获得m个被接受的采样点：

• 从 q(x) 中获得一个随机采样点 xi
• 对于 x_i 计算接受概率（acceptance probability）
  
  α=p(xi)Mq(xi)
• 从Uniform(0,1)中随机生成一个值，用u 表示
• 如果 α≥u，则接受 xi作为一个来自p(x)的采样值，否则就拒绝 xi并回到第一步。

你当然可以采用严密的数学推导来证明Reject Sampling的可行性。但它的原理从直观上来解释也是相当容易理解的。你可以想象一下在上图的例子中，从哪些位置抽出的点会比较容易被接受。显然，Mq(x)和p(x)所示之函数更加接近的地方接受概率较高，也即是更容易被接受，所以在这样的地方采到的点就会比较多，而在接受概率较低（即两个函数差距较大）的地方采到的点就会比较少，这也就保证了这个方法的有效性。
但是该算法也面临着几个限制。一方面并不是较为容易找到一个M去限制住p(x)q(x)。如果M太大，那么接收概率α为

Pα=P(u<p(x)Mq(x))=1M

将会比较小，这将在高维情形中变得不实际。

参考文献：
1. http://www.cnblogs.com/xbinworld/p/4266146.html
2. http://blog.csdn.net/baimafujinji/article/details/51407703
3. https://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/
4. An Introduction to MCMC for Machine Learning，2003