UA STAT675 統(tǒng)計計算I 隨機(jī)數(shù)生成7 Envelope Accept-Reject Algorithm

• • Squeeze Principle
  • Atkinson's Poisson Simulation

上一講我們推導(dǎo)了Accept-Reject Algorithm

---

Algorithm 3: Accept-Reject Algorithm

Step 1: Generate $y～g$, $u～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

---

它對所有的密度都適用，但前提是找到另一個密度作為工具密度，工具密度必須是目標(biāo)密度的強(qiáng)函數(shù)；要提高這個算法的效率，最好的做法是設(shè)計一個$g$，它比$f$稍微大一點(diǎn)點(diǎn)但又特別接近，使得$M\approx 1$，這一講我們就介紹一些常用的設(shè)計工具密度的方法。

Squeeze Principle

上一講我們分析了Accept-Reject Algorithm的效率，也就是要生成1個目標(biāo)密度的隨機(jī)數(shù)，平均需要$M$個均勻分布與工具密度的隨機(jī)數(shù)；一種可行的改進(jìn)是給目標(biāo)密度再找一個下界，進(jìn)一步降低采樣器的計算成本。

---

Algorithm 4: Envelope Accept-Reject Algorithm

Step 1: Generate $y～g_u$, $u～U(0,1)$
Step 2: If $u\le \frac{g_l(y)}{M{g}_u(y)}$, accept $y$ as a random number of $f$; otherwise, go to Step 3;
Step 3: if $u\le \frac{f(y)}{M{g}_u(y)}$, accept $y$ as a random number of $f$; otherwise, repeat Step 1-Step 3

---

算法分析

算法適用條件：構(gòu)造$g_l$與$g_u$滿足$$g_l\le f\le g_u$$其中$g_u$是工具密度(instrumental density)、$f$是目標(biāo)密度(target density)、$g_l$是輔助函數(shù)

算法幾何解釋：$g_l,M{g}_u$構(gòu)成了一個帶狀區(qū)域，$g_l$為下界，$M{g}_u$為上界，我們在$M{g}_u$下面的區(qū)域內(nèi)均勻采樣，一部分落在帶狀區(qū)域內(nèi)，另一部分落在帶狀區(qū)域下方；直接接受帶狀區(qū)域下方的點(diǎn)，把帶狀區(qū)域內(nèi)的點(diǎn)作為候選；目標(biāo)密度$f$把帶狀區(qū)域分割為上下兩部分，落在下部分的候選點(diǎn)就是$f$的隨機(jī)樣本

算法的效率：假設(shè)我們生成了一個$g_u$與$U(0,1)$的樣本，它被直接接受的概率是$$\frac{{\int }_{\mathbb{R}}{g}_l(x)dx}{{\int }_{\mathbb{R}}M{g}_u(x)dx}=\frac{{\int }_{\mathbb{R}}{g}_l(x)dx}{M}$$成為候選點(diǎn)并被接受為隨機(jī)樣本的概率是$$\frac{{\int }_{\mathbb{R}}(f?g_l)(x)dx}{{\int }_{\mathbb{R}}M{g}_u(x)dx}=\frac{1?{\int }_{\mathbb{R}}{g}_l(x)dx}{M}$$于是整體的接受率為$$\frac{{\int }_{\mathbb{R}}{g}_l(x)dx}{M}+\frac{1?{\int }_{\mathbb{R}}{g}_l(x)dx}{M}=\frac{1}{M}$$因此Envelope Accept-Reject Algorithm與Accept-Reject Algorithm的接受率一樣，在效率上沒有改進(jìn)，但是直接接受的樣本，也就是$\frac{{\int }_{\mathbb{R}}{g}_l(x)dx}{M}$這么多的樣本可以用$g_l(x)$的計算代替$f(x)$的計算（這個原則叫Squeeze Principle），因此只要$g_l$計算復(fù)雜度足夠低，并且與$f$足夠接近，Envelope Accept-Reject Algorithm的計算效率仍然可以得到很大提高

隨機(jī)數(shù)的獨(dú)立性分析：因?yàn)樯厦娴乃惴ㄖ?#xff0c;每一步生成隨機(jī)數(shù)與其他步驟都是可以互相獨(dú)立的，所以最后得到的隨機(jī)數(shù)可以有較強(qiáng)的獨(dú)立性

Squeeze Principle給我們提供了另一種改進(jìn)采樣器的思路，從Algorithm 1的簡單rejection采樣到Algorithm 3 Accept-Reject Algorithm，我們通過提高接受率來提高采樣器的效率；從Accept-Reject Algorithm到Envelope Accept-Reject Algorithm，我們通過降低計算量來提高采樣器的效率；提高接受率、降低計算量是兩種改進(jìn)采樣器的主流思路。

Atkinson’s Poisson Simulation

早期從Poisson分布中采樣有兩種主流方法，第一種方法是反函數(shù)法；第二種方法是使用Poisson過程法，這兩種方法都比較低效。Atkinson (1979)基于Accept-Reject Algorithm提出了一種從Poisson分布中采樣的方法。

考慮服從Logistics分布的隨機(jī)變量$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 }}}$$

定義$N=?x+0.5?$，考慮left-truncated Logistics分布，限制$x\in [?0.5,+\infty )$，則
$$P(N=n)=?\begin{array}{l}\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{array}$$

Atkinson的想法是使用這個分布作為從Poisson分布（記參數(shù)為$\lambda$）中采樣的工具密度，根據(jù)Accept-Reject Algorithm的算法效率，$\alpha ,\beta$的參數(shù)選擇最好使得
$$\frac{{\lambda }^n{e}^{?\lambda }}{n!P(N=n)}$$

的上確界越小越好，也就是要求解
$$\underset{\alpha ,\beta }{\min }\underset{n}{\sup }\frac{{\lambda }^n{e}^{?\lambda }}{n!P(N=n)}$$

由此Atkinson推導(dǎo)出了Poisson分布的采樣器：

---

Algorithm 5: Atkinson’s Poisson Sampler

Step 0: Define $\beta =\frac{\pi }{\sqrt{3\lambda }},\alpha =\lambda \beta ,c=0.767?3.36/\lambda ,k=\log c?\lambda ?\log \beta$;
Step 1: Generate $u_1～U(0,1)$ and calculate $$x=\frac{\alpha ?\log \frac{1?u_1}{u_1}}{\beta }$$if $x>?0.5$, go to Step 2; otherwise, repeat Step 1;
Step 2: Define $N=?x+0.5?$ and generate $u_2～U(0,1)$
Step 3: if $\alpha ?\beta x+\log \frac{u_2}{(1+e^{\alpha ?\beta x}{)}^2}\le k+N\log \lambda ?\log N!$, accept $N$ as a random number from $Poisson(\lambda )$

---

總結(jié)

以上是生活随笔為你收集整理的UA STAT675 统计计算I 随机数生成7 Envelope Accept-Reject Algorithm的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。