recursively generate exponential random variables

Question

I am new to recursion in R. I am trying to generate exponential random variables that meet a certain condition in R. Here is a simple example showing my attempt to generate two independent exponential random variables.

CODE

#### 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
  } 

}

EXAMPLE: q1=.25, q2=2

     gen.exp(.25,2)

When I run the above code, I get the following errors:

1. Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
2. Error during wrapup: evaluation nested too deeply: infinite recursion / options(expressions=)?

I have already tried modifying options(expressions=10000) and in order to allow R to increase the number of iterations. That has does not seem to help with my case (maybe I am not using the option correctly). I understand generating continuous distributions with stringent criteria as above maybe the problem. Having said that, is there anyway to avoid the errors? Or at least repeat the recursion whenever an error occurs? Is recursion an over kill here? Are there simpler/better ways of generating the desired random variables? I appreciate any insights.

Answer 1

You are using a simple rejection sampler with two conditions

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

However, the chance of matching both of them is low, i.e. very slow. However, since you are interested in exponential random numbers, we can avoid simulating numbers we would reject.

To generate an exponential random number, we use the formula

-rate * log(U)

where U is a U(0,1) random number. So to generate values from the exponential distribution larger than 10 (say), we just do

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

or in R code

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

We can use a similar trick for upper bounds.

Using @Roland's function from above, this gives

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."))
  }
}

Notice, I've also printed out how many iterations it took. For your parameters, I need around 2 iterations.

Answer 2

Don't use recursion for this:

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.

The probability for b to be larger than 4.8 is:

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

Let's try with more iterations:

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

Of course this RNG is so slow that it is almost useless. It would be better to produce larger batches of a and b at once.