rejection：接受拒绝采样（Acceptance-Rejection Sampling）

2024-04-20 09:00:57 0 0

我们所说的抽样，其实是指从一个概率分布中生成观察值（observations）的方法。而这个分布通常是由其概率密度函数（PDF）来表示的。而且，即使在已知PDF的情况下，让计算机自动生成观测值也不是一件容易的事情。从本质上来说，计算机只能实现对均匀分布（Uniform distribution）的采样。那如何实现计算机很好的采样数据样本呢？今天我们一起来看看实现方法。

在采样问题上我们可能会面对这些问题：

计算机只能实现对均匀分布的采样，但我们仍然可以在此基础上对更为复杂的分布进行采样，那具体该如何操作呢？
随机分布的某些数字特征可能需要通过积分的形式来求解，但是某些积分可能没有（或者很难求得）解析解，彼时我们该如何处理呢？
在贝叶斯推断中，后验概率的分布是正⽐于先验和似然函数之积的，但是先验和似然函数的乘积形式可能相对复杂，我们又该如何对这种形式复杂的分布进行采样呢？

针对这些问题衍生出一系列求解的方法。

一、接受拒绝采样（Acceptance-Rejection Sampling）

在数学中，拒绝抽样是用来从分布产生观测值的基本技术。它也被称为接受拒绝方法或“接受 - 拒绝算法”，是一种蒙特卡罗方法，这种方法与Metropolis-Hastings算法也有一定关系。

1. 简单认识

下图是一个随机变量的密度函数曲线，试问如何获得这个随机变量的样本呢？

利用这个曲线的形状来抽取样本,用一个矩形将这个密度曲线套起来，把密度曲线框在一个矩形里，如下：

然后，向这个矩形里随机投点。随机投点意味着在矩形这块区域内，这些点是满足均匀分布的。投了大概10000个点,如下面这个样子：

显然，有的点落在了密度曲线下侧，有的点落在了密度曲线的上侧。上侧的点用绿色来表示，下侧的点用蓝色来表示，如下图：

只保留密度曲线下侧的点，即蓝色的点：

到这里，提一个问题，在密度曲线以下的这块区域里，这些点满足什么分布？均匀分布！这是拒绝采样最关键的部分，搞个矩形、向矩形里投点等等，所做的一切都是为了获得一个密度曲线所围成区域的均匀分布。只要能获得这样一个在密度曲线下满足均匀分布的样本，我们就可以获得与该密度曲线相匹配的随机变量的采样样本。方法是，只需把每个蓝点的横坐标提取出来，这些横坐标所构成的样本就是我们的目标样本。下图左侧，是按照以上方法获得的一个样本的直方图以及核密度估计，下图右侧，是开始的密度曲线。

可见，采样样本的核密度估计与目标密度曲线基本一致，可以肯定这个样本就是目标样本。

最开始时候用到了一个矩形，这个矩形就是一个满足均匀分布的建议分布，建议分布只是获得目标密度函数曲线下均匀分布样本的一个辅助工具。采用均匀分布作为建议分布有时效率很低，为什么这么说？从上例就可以看出，均匀分布的好多点（那些绿点）都被剔除了，造成了一种浪费。可以选择一些其他曲线来把密度曲线框起来，效率会提高一点，如下图：

数曲线为h(x), 对应于下图中的蓝线，建议分布密度曲线为g(x),我们把g(x)乘上一个常数因子c,然后用cg(x)这条曲线将目标密度曲线框起来。

我们假定满足g(x)的随机变量易于采样，那么拒绝采样的步骤如下：

从g(x)采到一个样本数据，记为 x ⋆ x^{\star} x⋆，我们把它作为一个建议

要不要接受这个建议，作为满足h(x)分布的一个样本数据呢？我们定义一个接受概率:

也就是说，我们以 α \alpha α的概率接受 x ⋆ x^{\star} x⋆作为h(x)分布的一个样本数据。实际操作中，我们是取一个 U ( 0 , 1 ) U(0, 1) U(0,1)的随机数 μ \mu μ,如果 μ < α \mu<\alpha μ<α,就接受 x ⋆ x^{\star} x⋆作为h(x)的一个样本数据，否则，把它舍弃掉，回到1步继续循环。最终可以获得一个样本。

文章开头是一下子抽取10000个点，到后来怎么成了一个个抽了呢？其实它们是对应的，把蓝点去掉的过程就相当于你做是否拒绝判断的过程。
如果有密度曲线下的均匀分布样本，就可以得到与密度曲线相匹配的分布的一个样本。
如果建议分布的形状和目标分布越接近，采样的效率就越高。

2. Acceptance-Rejection Sampling过程

3. Acceptance-Rejection Sampling的直观解释

4. Acceptance-Rejection Sampling有效性证明(待)

5.python实现

2. 生成代码如下：

import randomimport mathimport matplotlib.pyplot as pltimport seaborn as snsimport numpy as np%matplotlib inlinesns.set_style('darkgrid')plt.rcParams['figure.figsize'] = (12, 8)def AceeptReject(split_val): global c global power while True: x = random.uniform(0, 1) y = random.uniform(0, 1) if y*c <= math.pow(x - split_val, power): return xpower = 4t = 0.4 sum_ = (math.pow(1-t, power + 1) - math.pow(-t, power + 1)) / (power + 1) #求积分x = np.linspace(0, 1, 100)#常数值cc = 0.6**4/sum_cc = [c for xi in x]plt.plot(x, cc, '--',label='c*f(x)')#目标概率密度函数的值f(x)y = [math.pow(xi - t, power)/sum_ for xi in x]plt.plot(x, y,label='f(x)')#采样10000个点samples = []for i in range(10000): samples.append(AceeptReject(t))plt.hist(samples, bins=50, normed=True,label='sampling')plt.legend()plt.show()