比率-的-制服分布中的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)) 

大量的代码可以被优化，但我认为这是不够好喜欢这个这个目的。我希望这有帮助。