递归生成指数随机变量

Question

我是 R 中递归的新手。我正在尝试生成满足 R 中特定条件的指数随机变量。这是一个简单的示例，展示了我尝试生成两个独立的指数随机变量的尝试。

代码

#### Function to generate two independent exponential random variable that meet two criteria
gen.exp<-function(q1,q2){
  a=rexp(1,q1)  # generate exponential random variable
  b=rexp(1,q2)  # generate exponential random variable
  if((a>10 & a<10.2) & (a+b>15)){ # criteria the random variables must meet
   return(c(a,b))
  }else{
   return(gen.exp(q1,q2)) #if above criteria is not met, repeat the process again
  } 

}

示例：q1=.25，q2=2

     gen.exp(.25,2)

当我运行上面的代码时，我得到以下错误：

1. 错误：求值嵌套太深：无限递归/选项（表达式=）？
2. 结束时出错：求值嵌套太深：无限递归/选项（表达式=）？

我已经尝试修改options(expressions=10000)，以允许R增加迭代次数。这似乎对我的情况没有帮助（也许我没有正确使用该选项）。我理解生成具有上述严格标准的连续分布可能是问题所在。话虽如此，有没有办法避免错误？或者至少在发生错误时重复递归？递归在这里过度杀戮吗？是否有更简单/更好的方法来生成所需的随机变量？ 我感谢任何见解。

Answer 1

您正在使用具有两个条件的简单拒绝采样器

• a > 10 & a < 10.2
• a+ b > 15

但是，两者都匹配的机会很低，即非常慢。但是，由于您对指数随机数感兴趣，我们可以避免模拟我们会拒绝的数字。

要生成指数随机数，我们使用公式

-rate * log(U)

其中U 是一个U(0,1) 随机数。因此，要从大于 10 的指数分布中生成值（例如），我们只需这样做

-log(U(0, exp(-10*rate))/rate

或在 R 代码中

-log(runif(1, 0, exp(-10*rate)))/rate

我们可以对上限使用类似的技巧。

使用上面的@Roland 函数，得到

gen.exp = function(q1, q2, maxiter = 1e3){
  i = 0
  repeat {
    i = i + 1
    upper = exp(-10*q1)
    lower = exp(-10.2*q1)
    a = -log(runif(1, lower, upper))/q1
    b = -log(runif(1, 0, exp(-4.8*q2)))/q2

    if((a>10 & a<10.2) & (a+b>15)) {message(i); return(c(a,b)) 
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries."))
  }
}

注意，我还打印了它需要多少次迭代。对于您的参数，我需要大约 2 次迭代。

Answer 2

不要为此使用递归：

gen.exp<-function(q1, q2, maxiter = 1e3){
  i <- 0
  repeat {
    i <- i + 1
    a=rexp(1,q1)  # generate exponential random variable
    b=rexp(1,q2)  # generate exponential random variable
    if((a>10 & a<10.2) & (a+b>15)) return(c(a,b)) # criteria the random variables must meet
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries."))
  }
}

set.seed(42)
gen.exp(.25, 2)
#Error in gen.exp(0.25, 2) : Conditions not fulfilled after 1000 tries.

b 大于 4.8 的概率为：

pexp(4.8, 2, lower.tail = FALSE)
#[1] 6.772874e-05

让我们尝试更多的迭代：

gen.exp(.25, 2, maxiter = 1e7)
#[1] 10.08664  5.55414

当然这个RNG太慢了，几乎没用。最好一次生产更大批量的a 和b。