UA STAT675 统计计算I 随机数生成6 Accept-Reject Algorithm

随机模拟基本定理(Fundamental Theorem of Simulation)
- 根据随机模拟基本定理设计一元随机变量的随机数生成器
- 随机模拟基本定理的推论

上一讲我们介绍了生成随机数的general transformation method，那是以U(0,1)U(0,1)U(0,1)的随机数为基础，通过变换获得其他分布的随机数的方法，当我们知道各种分布之间的变换规则，或者知道分布函数并能比较容易地求出它的反函数时，这种方法就是最直观最简单的；但是当我们想进行抽样的总体分布比较复杂时，我们就需要设计一些其他的方法了。这一讲我们介绍第一类采样的算法：accept-reject methods。

随机模拟基本定理(Fundamental Theorem of Simulation)

Target density为fff，则从X∼fX \sim fX∼f中采样等价于从
(X,U)∼U{(x,u):0<u<f(x)}(X,U) \sim U\{(x,u):0<u<f(x)\}(X,U)∼U{(x,u):0<u<f(x)}

中采样。

证明
这个定理的证明非常简单，因为
f(x)=∫0f(x)duf(x) = \int_0^{f(x)}duf(x)=∫0f(x)du

所以f(x)f(x)f(x)是二元随机变量(X,U)∼U{(x,u):0<u<f(x)}(X,U) \sim U\{(x,u):0<u<f(x)\}(X,U)∼U{(x,u):0<u<f(x)}中XXX的边缘分布，因此对二元随机变量(X,U)(X,U)(X,U)采样得到的XXX的样本服从fff。

但是我们并不需要UUU的样本，所以称UUU是一个auxiliary variable。

根据随机模拟基本定理设计一元随机变量的随机数生成器

假设Target density是一元函数，满足

f(x)≤Mf(x)\le Mf(x)≤M（密度有界）
{x∈R:f(x)>0}⊂[a,b]\{x \in \mathbb{R}:f(x)>0\} \subset [a,b]{x∈R:f(x)>0}⊂[a,b]（支撑集有界）

则
P(a≤X<x)=∫axf(y)dy=∫ax∫0f(y)dudy=∫ax∫0f(y)dudy∫ab∫0f(y)dudy=P(Y≤x∣U<f(Y))P(a \le X< x) = \int_a^x f(y)dy = \int_a^x \int_0^{f(y)}dudy \\ = \frac{\int_a^x \int_0^{f(y)}dudy}{\int_a^b \int_0^{f(y)}dudy}=P(Y \le x|U<f(Y))P(a≤X<x)=∫axf(y)dy=∫ax∫0f(y)dudy=∫ab∫0f(y)dudy∫ax∫0f(y)dudy=P(Y≤x∣U<f(Y))

其中U∼U(0,M)U \sim U(0,M)U∼U(0,M), Y∼U(a,b)Y \sim U(a,b)Y∼U(a,b)，这个推导给了我们一种设计arget density的随机数生成器的思路：

Algorithm 1

Step 1: Generate y∼U(a,b)y \sim U(a,b)y∼U(a,b)
Step 2: Generate u∼U(0,M)u \sim U(0,M)u∼U(0,M)
Step 3: If u<f(y)u<f(y)u<f(y), accept yyy as a random number of fff; otherwise, repeat Step 1-Step 3

算法分析

算法适用条件：根据上面的推导，这个算法适用于值域与支撑集都有界的密度；用更直白的话讲，就是适用于在xxx轴和在yyy轴上都有界的分布；
算法几何解释：设想我们画出了fff的图像，并且找了[a,b]×[0,M][a,b] \times [0,M][a,b]×[0,M]这个矩形把它包围起来，fff的图像把这个矩形分成了上下两部分，接下来我们从(Y,U)(Y,U)(Y,U)中采样，得到的样本(y,u)(y,u)(y,u)其实是矩形中的点，yyy代表横坐标，uuu代表纵坐标，如果这个点位于矩形的下半部分，就认为yyy是fff的样本；
算法的效率：假设我们想要nnn个fff的样本，则我们平均至少需要生成nM(b−a)nM(b-a)nM(b−a)个随机数（因为[a,b][a,b][a,b]上fff围成的面积最大为1，矩形围成的面积为M(b−a)M(b-a)M(b−a)），这说明这个算法的效率取决于Target density的性质，如果Target density厚尾或者存在比较大的峰值，这个算法的效率就会非常低；
随机数的独立性分析：因为上面的算法中，每一步生成随机数与其他步骤都是可以互相独立的，所以最后得到的随机数可以有较强的独立性

随机模拟基本定理的推论

正如我们在算法分析中讨论的一样，基于随机模拟基本定理设计的算法效率取决于Target density的形状，如果Target density形状比较差，比如支撑集为R\mathbb{R}R或者有比较严重的concentration，上面的算法效率就会很差。不难发现上述算法局限在于我们总是在试图用一个矩形去包围一个面积固定但形状可以千奇百怪的区域，那么是否可以放弃矩形包围的思路，针对不同形状的区域设计不一样的包围方法呢？

随机模拟基本定理的推论
Target density f(x)f(x)f(x)
Instrumental density g(x)g(x)g(x)
假设f(x)≤Mg(x)f(x) \le Mg(x)f(x)≤Mg(x)，∃M≥1\exists M\ge 1∃M≥1，则从X∼fX \sim fX∼f中抽样可以用下面的算法：

Algorithm 2

Step 1: Generate y∼gy \sim gy∼g
Step 2: Generate u∼U(0,Mg(y))u \sim U(0,Mg(y))u∼U(0,Mg(y))
Step 3: If u<f(y)u<f(y)u<f(y), accept yyy as a random number of fff; otherwise, repeat Step 1-Step 3

证明
如果X∼fX \sim fX∼f，∀B∈B(R)\forall B \in \mathcal{B}(\mathbb{R})∀B∈B(R)，
P(X∈B)=∫Bf(y)dy=∫B∫0f(y)1Mg(y)dudy=∫B∫0f(y)1Mg(y)dudy∫R∫0f(y)1Mg(y)dudy=P(Y∈B∣U<f(Y))P(X \in B) = \int_{B} f(y)dy = \int_B\int_0^{f(y)}\frac{1}{Mg(y)}dudy\\ = \frac{\int_B \int_0^{f(y)}\frac{1}{Mg(y)}dudy}{\int_{\mathbb{R}} \int_0^{f(y)}\frac{1}{Mg(y)}dudy}=P(Y \in B|U<f(Y)) P(X∈B)=∫Bf(y)dy=∫B∫0f(y)Mg(y)1dudy=∫R∫0f(y)Mg(y)1dudy∫B∫0f(y)Mg(y)1dudy=P(Y∈B∣U<f(Y))

这个式子说明，在U<f(Y)U<f(Y)U<f(Y)的条件下，XXX的分布与YYY的分布是相同的，于是此时的YYY的随机数服从target density；

算法分析
首先，我们把算法2的第2、3步合并一下：

Algorithm 3: Accept-Reject Algorithm

Step 1: Generate y∼gy \sim gy∼g, u∼U(0,1)u \sim U(0,1)u∼U(0,1)
Step 2: If u<f(y)Mg(y)u<\frac{f(y)}{Mg(y)}u<Mg(y)f(y), accept yyy as a random number of fff; otherwise, repeat Step 1-Step 2

这样关于均匀分布的处理就比较标准化了，定义
α(y)=f(y)Mg(y)\alpha(y) = \frac{f(y)}{Mg(y)}α(y)=Mg(y)f(y)

称α\alphaα为acceptance rate；在fff与MgMgMg比较接近的区域，acceptance rate较高。

算法适用条件：Accept-Reject Algorithm对所有的密度都适用，但前提是找到另一个密度作为工具密度，工具密度必须是目标密度的强函数；
算法几何解释：与算法1不同，现在我们放松了支撑集有界的假设，改成了用Mg(x)Mg(x)Mg(x)来包围f(x)f(x)f(x)；
算法的效率：不难发现Accept-Reject Algorithm取决于f(x)≤Mg(x)f(x)\le Mg(x)f(x)≤Mg(x)这个不等式有多tight，也就是Mg(x)Mg(x)Mg(x)与f(x)f(x)f(x)的距离有多近，可以简单计算一下∫RMg(x)dx∫Rf(x)dx=M\frac{\int_{\mathbb{R}}Mg(x)dx}{\int_{\mathbb{R}}f(x)dx}=M∫Rf(x)dx∫RMg(x)dx=M所以要得到nnn个服从fff的随机数，平均需要MMM个均匀分布的随机变量，因此要提高这个算法的效率，最好的做法是设计一个ggg，它比fff稍微大一点点但又特别接近，使得M≈1M \approx 1M≈1，这种Accept-Reject sampler就会具有非常高的效率，一个非常好的例子是Horseshoe estimation的算法中的一个rejection sampler，参考James Johndrow, Paulo Orenstein, Anirban Bhattacharya; 21(73):1−61, 2020. appendix S1.
随机数的独立性分析：因为上面的算法中，每一步生成随机数与其他步骤都是可以互相独立的，所以最后得到的随机数可以有较强的独立性