UA STAT675 统计计算I 随机数生成7 Envelope Accept-Reject Algorithm
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
算法分析
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)}α,βmin?nsup?n!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=βα?logu1?1?u1???if 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(λ)
總結
以上是生活随笔為你收集整理的UA STAT675 统计计算I 随机数生成7 Envelope Accept-Reject Algorithm的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: UA MATH567 高维统计专题1 稀
- 下一篇: UA STAT675 统计计算I 随机数