Search code examples
rexponential-distribution

Exponential distribution in R


I want to simulate some data from an exp(1) distribution but they have to be > 0.5 .so i used a while loop ,but it does not seem to work as i would like to .Thanks in advance for your responses !

x1<-c()

w<-rexp(1) 

while (length(x1) < 100) {

  if (w > 0.5) {

    x1<- w }

  else {

    w<-rexp(1)

  }

}

Solution

  • 1) The code in the question has these problems:

    • we need a new random variable on each iteration but it only generates new random variables if the if condition is FALSE

    • x1 is repeatedly overwritten rather than extended

    • although while could be used repeat seems better since having the test at the end is a better fit than the test at the beginning

    We can fix this up like this:

    x1 <- c()
    repeat {
      w <- rexp(1)
      if (w > 0.5) {
        x1 <- c(x1, w)
        if (length(x1) == 100) break
      }
    }
    

    1a) A variation would be the following. Note that an if whose condition is FALSE evaluates to NULL if there is no else leg so if the condition is FALSE on the line marked ## then nothing is concatenated to x1.

    x1 <- c()
    repeat {
      w <- rexp(1)
      x1 <- c(x1, if (w > 0.5) w)  ##
      if (length(x1) == 100) break
    }
    

    2) Alternately, this generates 200 exponential random variables keeping only those greater than 0.5. If fewer than 100 are generated then repeat. At the end it takes the first 100 from the last batch generated. We have chosen 200 to be sufficiently large that on most runs only one iteration of the loop will be needed.

    repeat {
      r <- rexp(200)
      r <- r[r > 0.5]
      if (length(r) >= 100) break
    }
    r <- head(r, 100)
    

    Alternative (2) is actually faster than (1) or (1a) because it is more highly vectorized. This is despite it throwing away more exponential random variables than the other solutions.