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

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

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

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

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

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

中采样。

证明
这个定理的证明非常简单，因为
$$f(x)={\int }_0^{f(x)}du$$

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

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

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

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

$f(x)\le M$（密度有界）

$\{x\in \mathbb{R}:f(x)>0\}\subset [a,b]$（支撑集有界）

则
$$\begin{gathered} 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)) \end{gathered}$$

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

---

Algorithm 1

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

---

算法分析

算法适用条件：根据上面的推导，这个算法适用于值域与支撑集都有界的密度；用更直白的话讲，就是适用于在$x$轴和在$y$轴上都有界的分布；

算法几何解释：设想我们画出了$f$的图像，并且找了$[a,b]\times [0,M]$这个矩形把它包围起来，$f$的图像把这个矩形分成了上下两部分，接下来我们从$(Y,U)$中采样，得到的样本$(y,u)$其实是矩形中的点，$y$代表横坐标，$u$代表纵坐标，如果这个点位于矩形的下半部分，就认为$y$是$f$的样本；

算法的效率：假设我们想要$n$个$f$的样本，则我们平均至少需要生成$nM(b-a)$个随机数（因为$[a,b]$上$f$围成的面积最大为1，矩形围成的面积为$M(b-a)$），这说明这个算法的效率取决于Target density的性质，如果Target density厚尾或者存在比较大的峰值，这个算法的效率就会非常低；

随机数的独立性分析：因为上面的算法中，每一步生成随机数与其他步骤都是可以互相独立的，所以最后得到的随机数可以有较强的独立性

随机模拟基本定理的推论

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

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

---

Algorithm 2

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

---

证明
如果$X\sim f$，$\forall B\in \mathcal{B}(\mathbb{R})$，
$$\begin{gathered} 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)) \end{gathered}$$

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

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

---

Algorithm 3: Accept-Reject Algorithm

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

---

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

称$\alpha$为acceptance rate；在$f$与$Mg$比较接近的区域，acceptance rate较高。

算法适用条件：Accept-Reject Algorithm对所有的密度都适用，但前提是找到另一个密度作为工具密度，工具密度必须是目标密度的强函数；

算法几何解释：与算法1不同，现在我们放松了支撑集有界的假设，改成了用$Mg(x)$来包围$f(x)$；

算法的效率：不难发现Accept-Reject Algorithm取决于$f(x)\le Mg(x)$这个不等式有多tight，也就是$Mg(x)$与$f(x)$的距离有多近，可以简单计算一下$$\frac{{\int }_{\mathbb{R}}Mg(x)dx}{{\int }_{\mathbb{R}}f(x)dx}=M$$所以要得到$n$个服从$f$的随机数，平均需要$M$个均匀分布的随机变量，因此要提高这个算法的效率，最好的做法是设计一个$g$，它比$f$稍微大一点点但又特别接近，使得$M\approx 1$，这种Accept-Reject sampler就会具有非常高的效率，一个非常好的例子是Horseshoe estimation的算法中的一个rejection sampler，参考James Johndrow, Paulo Orenstein, Anirban Bhattacharya; 21(73):1−61, 2020. appendix S1.

随机数的独立性分析：因为上面的算法中，每一步生成随机数与其他步骤都是可以互相独立的，所以最后得到的随机数可以有较强的独立性

总结

以上是生活随笔为你收集整理的UA STAT675 统计计算I 随机数生成6 Accept-Reject Algorithm的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。