比率-的-制服分佈中的R

Question

我有一個運動，其中i的有如下創建一個算法：

比制服的是基於這樣的事實，與密度爲f的隨機變量X（X ），我們可以通過計算X = U/V爲一對（U生成所需的密度X，V）均勻地分佈在集

Af = {(u,v):0 < v ≤ f(u/v)} 

隨機

點能夠被均勻地在由房顫拒絕從MIN-採樣imal邊界矩形，即包含Af的最小可能矩形。

1. 生成隨機數： 它是由（U-，U +）×（0，V +），其中

   v+ = max f(x), x 
   
   u− = minx f(x), x 
   
   u+ = maxx f(x) 
   

   然後比率-的-制服的方法包括以下簡單步驟給出U統一在（u-，u +）中。

2. 在（0，v +）中統一生成隨機數V.

3. 設置X←U/V。

4. 如果V 2≤F（X）接受並返回X.

5. 否則再試一次。

我迄今爲止代碼：

x <- cnorm(1, mean = 0, sd=1) 

myrnorm <- function(pdf){ 
    ## call rou() n times 
    pdf <- function(x) {exp(-x^2/2)} 
    } 
rou <- function(u, v) { 
    uplus <- 1 
    vplus <- 1 
    n <- 100 
    u <- runif(n, min=0, max=uplus) 
    v <- runif(n, min=0, max=vplus) 
    xi <- v/u 
    while(v < sqrt(xi)) { 
    if(v^2 <= xi) 
     return(xi) 
    } 
} 


myx <- myrnorm(1000) 
hist(myx) 

但我真的不知道該怎麼走下去。我不喜歡這個練習。我會很感激任何建議。

Answer 1

繼此link和示例代碼的第8頁例1，我想出了這個解決方案：

ratioU <- function(nvals) 
{ 
    h_x = function(x) exp(-x) 
    # u- is b-, u+ is b+ and v+ is a in the example: 
    uminus = 0 
    uplus = 2/exp(1) 
    vplus = 1 
    X.vals <- NULL 
    i <- 0 
    repeat { 
    i <- i+1 
    u <- runif(1,0,vplus) 
    v <- runif(1,uminus,uplus) 
    X <- u/v 
    if(v^2 <= h_x(X)) { 
     tmp <- X 
    } 
    else { 
     next 
    } 
    X.vals <- c(X.vals,tmp) 
    if(length(X.vals) >= nvals) break 
    } 
    answer <- X.vals 
    answer 
} 

sol = ratioU(1000) 
par(mfrow=c(1,2)) 
hist(sol,breaks=50, main= "using ratioU",freq=F) 
hist(rexp(1000),breaks = 50, main="using rexp from R",freq=F) 
par(mfrow=c(1,1)) 

par(mfrow=c(1,2)) 
plot(density(sol)) 
plot(density(rexp(1000))) 
par(mfrow=c(1,1)) 

大量的代碼可以被優化，但我認爲這是不夠好喜歡這個這個目的。我希望這有幫助。