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～f$中采樣等價于從
$$(X,U)～U\{(x,u):0<u<f(x)\}$$

中采樣。

證明
這個定理的證明非常簡單，因為
$$f(x)={\int }_0^{f(x)}du$$

所以$f(x)$是二元隨機變量$(X,U)～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\}?[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)$，這個推導給了我們一種設計arget density的隨機數生成器的思路：

---

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

---

算法分析

算法適用條件：根據上面的推導，這個算法適用于值域與支撐集都有界的密度；用更直白的話講，就是適用于在$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)$，$?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}$$

這個式子說明，在$U<f(Y)$的條件下，$X$的分布與$Y$的分布是相同的，于是此時的$Y$的隨機數服從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

---

這樣關于均勻分布的處理就比較標準化了，定義
$$\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的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。