Generate random numbers in R satisfying constraints

2
votes

I need help with a code to generate random numbers according to constraints. Specifically, I am trying to simulate random numbers ALFA and BETA from, respectively, a Normal and a Gamma distribution such that ALFA - BETA < 1.

Here is what I have written but it does not work at all.

set.seed(42) 
n <- 0 
repeat {
  n <- n + 1
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) 
  alfa[n] <- a
  beta[n] <- b
  l = length(alfa)
  if (l == 10000) break
}
2
rrandomconstraints
You should use curly braces {} after your first if statement to enclose the dependent lines. By default, only the first line after the if will depend on the if. – Gregor Thomas

2 Answers

3
votes

Due to vectorization, it will be faster to generate the numbers "all at once" rather than in a loop:

set.seed(42)
N = 1e5
a = rnorm(N, 10, 2)
b = rgamma(N, 8, 1)
d = a - b
alfa = a[d < 1]
beta = b[d < 1]
length(alfa)
# [1] 36436

This generated 100,000 candidates, 36,436 of which met your criteria. If you want to generate n samples, try setting N = 4 * n and you'll probably generate more than enough, keep the first n.

Your loop has 2 problems: (a) you need curly braces to enclose multiple lines after an if statement. (b) you are using n as an attempt counter, but it should be a success counter. As written, your loop will only stop if the 10000th attempt is a success. Move n <- n + 1 inside the if statement to fix:

set.seed(42) 
n <- 0 
alfa = numeric(0)
beta = numeric(0)
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) {
    n <- n + 1
    alfa[n] <- a
    beta[n] <- b
    l = length(alfa)
    if (l == 500) break
  }
}

But the first way is better... due to "growing" alfa and beta in the loop, and generating numbers one at a time, this method takes longer to generate 500 numbers than the code above takes to generate 30,000.

0
votes

As commented by @Gregor Thomas, the failure of your attempt is due to the missing of curly braces to enclose the if statement. If you would like to skip {} for if control, maybe you can try the code below

set.seed(42) 
r <- list()
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) r[[length(r)+1]] <- cbind(alfa = a, beta = b)
  if (length(r) == 100000) break
}
r <- do.call(rbind,r)

such that

> head(r)
          alfa      beta
[1,]  9.787751 12.210648
[2,]  9.810682 14.046190
[3,]  9.874572 11.499204
[4,]  6.473674  8.812951
[5,]  8.720010  8.799160
[6,] 11.409675 10.602608