递归生成指数随机变量

debugcn 发表于 Dev

用户名

我是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)

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

错误：评估嵌套太深：无限递归/选项（表达式=）？
包装过程中的错误：评估嵌套太深：无限递归/ options（expressions =）？

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

csgillespie

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

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次迭代。

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。