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

Squeeze Principle
Atkinson's Poisson Simulation

上一讲我们推导了Accept-Reject Algorithm

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

它对所有的密度都适用，但前提是找到另一个密度作为工具密度，工具密度必须是目标密度的强函数；要提高这个算法的效率，最好的做法是设计一个ggg，它比fff稍微大一点点但又特别接近，使得M≈1M \approx 1M≈1，这一讲我们就介绍一些常用的设计工具密度的方法。

Squeeze Principle

上一讲我们分析了Accept-Reject Algorithm的效率，也就是要生成1个目标密度的随机数，平均需要MMM个均匀分布与工具密度的随机数；一种可行的改进是给目标密度再找一个下界，进一步降低采样器的计算成本。

Algorithm 4: Envelope Accept-Reject Algorithm

Step 1: Generate y∼guy \sim g_uy∼gu, u∼U(0,1)u \sim U(0,1)u∼U(0,1)
Step 2: If u≤gl(y)Mgu(y)u\le \frac{g_l(y)}{Mg_u(y)}u≤Mgu(y)gl(y), accept yyy as a random number of fff; otherwise, go to Step 3;
Step 3: if u≤f(y)Mgu(y)u \le \frac{f(y)}{Mg_u(y)}u≤Mgu(y)f(y), accept yyy as a random number of fff; otherwise, repeat Step 1-Step 3

算法分析

算法适用条件：构造glg_lgl与gug_ugu满足gl≤f≤gug_l \le f \le g_ugl≤f≤gu其中gug_ugu是工具密度(instrumental density)、fff是目标密度(target density)、glg_lgl是辅助函数
算法几何解释：gl,Mgug_l,Mg_ugl,Mgu构成了一个带状区域，glg_lgl为下界，MguMg_uMgu为上界，我们在MguMg_uMgu下面的区域内均匀采样，一部分落在带状区域内，另一部分落在带状区域下方；直接接受带状区域下方的点，把带状区域内的点作为候选；目标密度fff把带状区域分割为上下两部分，落在下部分的候选点就是fff的随机样本
算法的效率：假设我们生成了一个gug_ugu与U(0,1)U(0,1)U(0,1)的样本，它被直接接受的概率是∫Rgl(x)dx∫RMgu(x)dx=∫Rgl(x)dxM\frac{\int_{\mathbb{R}}g_l(x)dx}{\int_{\mathbb{R}}Mg_u(x)dx}=\frac{\int_{\mathbb{R}}g_l(x)dx}{M}∫RMgu(x)dx∫Rgl(x)dx=M∫Rgl(x)dx成为候选点并被接受为随机样本的概率是∫R(f−gl)(x)dx∫RMgu(x)dx=1−∫Rgl(x)dxM\frac{\int_{\mathbb{R}}(f-g_l)(x)dx}{\int_{\mathbb{R}}Mg_u(x)dx}=\frac{1-\int_{\mathbb{R}}g_l(x)dx}{M}∫RMgu(x)dx∫R(f−gl)(x)dx=M1−∫Rgl(x)dx于是整体的接受率为∫Rgl(x)dxM+1−∫Rgl(x)dxM=1M\frac{\int_{\mathbb{R}}g_l(x)dx}{M}+\frac{1-\int_{\mathbb{R}}g_l(x)dx}{M}=\frac{1}{M}M∫Rgl(x)dx+M1−∫Rgl(x)dx=M1因此Envelope Accept-Reject Algorithm与Accept-Reject Algorithm的接受率一样，在效率上没有改进，但是直接接受的样本，也就是∫Rgl(x)dxM\frac{\int_{\mathbb{R}}g_l(x)dx}{M}M∫Rgl(x)dx这么多的样本可以用gl(x)g_l(x)gl(x)的计算代替f(x)f(x)f(x)的计算（这个原则叫Squeeze Principle），因此只要glg_lgl计算复杂度足够低，并且与fff足够接近，Envelope Accept-Reject Algorithm的计算效率仍然可以得到很大提高
随机数的独立性分析：因为上面的算法中，每一步生成随机数与其他步骤都是可以互相独立的，所以最后得到的随机数可以有较强的独立性

Squeeze Principle给我们提供了另一种改进采样器的思路，从Algorithm 1的简单rejection采样到Algorithm 3 Accept-Reject Algorithm，我们通过提高接受率来提高采样器的效率；从Accept-Reject Algorithm到Envelope Accept-Reject Algorithm，我们通过降低计算量来提高采样器的效率；提高接受率、降低计算量是两种改进采样器的主流思路。

Atkinson’s Poisson Simulation

早期从Poisson分布中采样有两种主流方法，第一种方法是反函数法；第二种方法是使用Poisson过程法，这两种方法都比较低效。Atkinson (1979)基于Accept-Reject Algorithm提出了一种从Poisson分布中采样的方法。

考虑服从Logistics分布的随机变量XXX：
f(x)=1βe−x−αβ(1+e−x−αβ)2,F(x)=11+e−x−αβf(x) = \frac{1}{\beta}\frac{e^{-\frac{x-\alpha}{\beta}}}{(1+e^{-\frac{x-\alpha}{\beta}})^2},F(x)=\frac{1}{1+e^{-\frac{x-\alpha}{\beta}}}f(x)=β1(1+e−βx−α)2e−βx−α,F(x)=1+e−βx−α1

定义N=⌊x+0.5⌋N=\lfloor x+0.5 \rfloorN=⌊x+0.5⌋，考虑left-truncated Logistics分布，限制x∈[−0.5,+∞)x \in [-0.5,+\infty)x∈[−0.5,+∞)，则
P(N=n)={11+e−n+0.5−αβ−11+e−n−0.5−αβ,x>1/2(11+e−n+0.5−αβ−11+e−n−0.5−αβ)1+e−0.5+αβe−0.5+αβ,−1/2<x≤1/2P(N=n)=\begin{cases} \frac{1}{1+e^{-\frac{n+0.5-\alpha}{\beta}}}- \frac{1}{1+e^{-\frac{n-0.5-\alpha}{\beta}}},x>1/2 \\ (\frac{1}{1+e^{-\frac{n+0.5-\alpha}{\beta}}}- \frac{1}{1+e^{-\frac{n-0.5-\alpha}{\beta}}})\frac{1+e^{-\frac{0.5+\alpha}{\beta}}}{e^{-\frac{0.5+\alpha}{\beta}}},-1/2<x \le 1/2\end{cases}P(N=n)=⎩⎪⎨⎪⎧1+e−βn+0.5−α1−1+e−βn−0.5−α1,x>1/2(1+e−βn+0.5−α1−1+e−βn−0.5−α1)e−β0.5+α1+e−β0.5+α,−1/2<x≤1/2

Atkinson的想法是使用这个分布作为从Poisson分布（记参数为λ\lambdaλ）中采样的工具密度，根据Accept-Reject Algorithm的算法效率，α,β\alpha,\betaα,β的参数选择最好使得
λne−λn!P(N=n)\frac{\lambda^ne^{-\lambda}}{n!P(N=n)}n!P(N=n)λne−λ

的上确界越小越好，也就是要求解
min⁡α,βsup⁡nλne−λn!P(N=n)\min_{\alpha,\beta} \sup_{n}\frac{\lambda^ne^{-\lambda}}{n!P(N=n)}α,βminnsupn!P(N=n)λne−λ

由此Atkinson推导出了Poisson分布的采样器：

Algorithm 5: Atkinson’s Poisson Sampler

Step 0: Define β=π3λ,α=λβ,c=0.767−3.36/λ,k=log⁡c−λ−log⁡β\beta = \frac{\pi}{\sqrt{3\lambda}},\alpha=\lambda \beta,c=0.767-3.36/\lambda,k=\log c-\lambda - \log \betaβ=3λπ,α=λβ,c=0.767−3.36/λ,k=logc−λ−logβ;
Step 1: Generate u1∼U(0,1)u_1 \sim U(0,1)u1∼U(0,1) and calculate x=α−log⁡1−u1u1βx=\frac{\alpha-\log \frac{1-u_1}{u_1}}{\beta}x=βα−logu11−u1if x>−0.5x>-0.5x>−0.5, go to Step 2; otherwise, repeat Step 1;
Step 2: Define N=⌊x+0.5⌋N=\lfloor x+0.5 \rfloorN=⌊x+0.5⌋ and generate u2∼U(0,1)u_2 \sim U(0,1)u2∼U(0,1)
Step 3: if α−βx+log⁡u2(1+eα−βx)2≤k+Nlog⁡λ−log⁡N!\alpha-\beta x+\log \frac{u_2}{(1+e^{\alpha-\beta x})^2} \le k+ N\log \lambda - \log N!α−βx+log(1+eα−βx)2u2≤k+Nlogλ−logN!, accept NNN as a random number from Poisson(λ)Poisson(\lambda)Poisson(λ)