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

• • 隨機(jī)模擬基本定理(Fundamental Theorem of Simulation)
  • • 根據(jù)隨機(jī)模擬基本定理設(shè)計(jì)一元隨機(jī)變量的隨機(jī)數(shù)生成器
    • 隨機(jī)模擬基本定理的推論

上一講我們介紹了生成隨機(jī)數(shù)的general transformation method，那是以$U(0,1)$的隨機(jī)數(shù)為基礎(chǔ)，通過變換獲得其他分布的隨機(jī)數(shù)的方法，當(dāng)我們知道各種分布之間的變換規(guī)則，或者知道分布函數(shù)并能比較容易地求出它的反函數(shù)時(shí)，這種方法就是最直觀最簡(jiǎn)單的；但是當(dāng)我們想進(jìn)行抽樣的總體分布比較復(fù)雜時(shí)，我們就需要設(shè)計(jì)一些其他的方法了。這一講我們介紹第一類采樣的算法：accept-reject methods。

隨機(jī)模擬基本定理(Fundamental Theorem of Simulation)

Target density為$f$，則從$X～f$中采樣等價(jià)于從
$$(X,U)～U\{(x,u):0<u<f(x)\}$$

中采樣。

證明
這個(gè)定理的證明非常簡(jiǎn)單，因?yàn)?br /> $$f(x)={\int }_0^{f(x)}du$$

所以$f(x)$是二元隨機(jī)變量$(X,U)～U\{(x,u):0<u<f(x)\}$中$X$的邊緣分布，因此對(duì)二元隨機(jī)變量$(X,U)$采樣得到的$X$的樣本服從$f$。

但是我們并不需要$U$的樣本，所以稱$U$是一個(gè)auxiliary variable。

根據(jù)隨機(jī)模擬基本定理設(shè)計(jì)一元隨機(jī)變量的隨機(jī)數(shù)生成器

假設(shè)Target density是一元函數(shù)，滿足

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

$\{x\in \mathbb{R}:f(x)>0\}?[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～U(0,M)$, $Y～U(a,b)$，這個(gè)推導(dǎo)給了我們一種設(shè)計(jì)arget density的隨機(jī)數(shù)生成器的思路：

---

Algorithm 1

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

---

算法分析

算法適用條件：根據(jù)上面的推導(dǎo)，這個(gè)算法適用于值域與支撐集都有界的密度；用更直白的話講，就是適用于在$x$軸和在$y$軸上都有界的分布；

算法幾何解釋：設(shè)想我們畫出了$f$的圖像，并且找了$[a,b]\times [0,M]$這個(gè)矩形把它包圍起來，$f$的圖像把這個(gè)矩形分成了上下兩部分，接下來我們從$(Y,U)$中采樣，得到的樣本$(y,u)$其實(shí)是矩形中的點(diǎn)，$y$代表橫坐標(biāo)，$u$代表縱坐標(biāo)，如果這個(gè)點(diǎn)位于矩形的下半部分，就認(rèn)為$y$是$f$的樣本；

算法的效率：假設(shè)我們想要$n$個(gè)$f$的樣本，則我們平均至少需要生成$nM(b?a)$個(gè)隨機(jī)數(shù)（因?yàn)?span id="ozvdkddzhkzd" class="katex--inline">$[a,b]$上$f$圍成的面積最大為1，矩形圍成的面積為$M(b?a)$），這說明這個(gè)算法的效率取決于Target density的性質(zhì)，如果Target density厚尾或者存在比較大的峰值，這個(gè)算法的效率就會(huì)非常低；

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

隨機(jī)模擬基本定理的推論

正如我們?cè)谒惴ǚ治鲋杏懻摰囊粯?#xff0c;基于隨機(jī)模擬基本定理設(shè)計(jì)的算法效率取決于Target density的形狀，如果Target density形狀比較差，比如支撐集為$\mathbb{R}$或者有比較嚴(yán)重的concentration，上面的算法效率就會(huì)很差。不難發(fā)現(xiàn)上述算法局限在于我們總是在試圖用一個(gè)矩形去包圍一個(gè)面積固定但形狀可以千奇百怪的區(qū)域，那么是否可以放棄矩形包圍的思路，針對(duì)不同形狀的區(qū)域設(shè)計(jì)不一樣的包圍方法呢？

隨機(jī)模擬基本定理的推論
Target density $f(x)$
Instrumental density $g(x)$
假設(shè)$f(x)\le Mg(x)$，$?M\ge 1$，則從$X～f$中抽樣可以用下面的算法：

---

Algorithm 2

Step 1: Generate $y～g$
Step 2: Generate $u～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～f$，$?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}$$

這個(gè)式子說明，在$U<f(Y)$的條件下，$X$的分布與$Y$的分布是相同的，于是此時(shí)的$Y$的隨機(jī)數(shù)服從target density；

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

---

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

---

這樣關(guān)于均勻分布的處理就比較標(biāo)準(zhǔn)化了，定義
$$\alpha (y)=\frac{f(y)}{Mg(y)}$$

稱$\alpha$為acceptance rate；在$f$與$Mg$比較接近的區(qū)域，acceptance rate較高。

算法適用條件：Accept-Reject Algorithm對(duì)所有的密度都適用，但前提是找到另一個(gè)密度作為工具密度，工具密度必須是目標(biāo)密度的強(qiáng)函數(shù)；

算法幾何解釋：與算法1不同，現(xiàn)在我們放松了支撐集有界的假設(shè)，改成了用$Mg(x)$來包圍$f(x)$；

算法的效率：不難發(fā)現(xiàn)Accept-Reject Algorithm取決于$f(x)\le Mg(x)$這個(gè)不等式有多tight，也就是$Mg(x)$與$f(x)$的距離有多近，可以簡(jiǎn)單計(jì)算一下$$\frac{{\int }_{\mathbb{R}}Mg(x)dx}{{\int }_{\mathbb{R}}f(x)dx}=M$$所以要得到$n$個(gè)服從$f$的隨機(jī)數(shù)，平均需要$M$個(gè)均勻分布的隨機(jī)變量，因此要提高這個(gè)算法的效率，最好的做法是設(shè)計(jì)一個(gè)$g$，它比$f$稍微大一點(diǎn)點(diǎn)但又特別接近，使得$M\approx 1$，這種Accept-Reject sampler就會(huì)具有非常高的效率，一個(gè)非常好的例子是Horseshoe estimation的算法中的一個(gè)rejection sampler，參考James Johndrow, Paulo Orenstein, Anirban Bhattacharya; 21(73):1?61, 2020. appendix S1.

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

總結(jié)

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

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