1

We have a big for loop in R for simulating various data where for some iterations the data generate in such a way that a quantity comes 0 inside the loop, which is not desirable and we should skip that step of data generation. But at the same time we also need to increase the number of iterations by one step because of such skip, otherwise we will have fewer observations than required.

For example, while running the following code, we get z=0 in iteration 1, 8 and 9. 

rm(list=ls())
n <- 10
z <- NULL
for(i in 1:n){
  set.seed(i)
  a <- rbinom(1,1,0.5)
  b <- rbinom(1,1,0.5)
  z[i] <- a+b
}
z
[1] 0 1 1 1 1 2 1 0 0 1

We desire to skip these steps so that we do not have any z=0 but we also want a vector z of length 10. It may be done in many ways. But what I particularly want to see is how we can stop the iteration and skip the current step when z=0 is encountered and go to the next step, ultimately obtaining 10 observations for z.

Blain Waan
  • 393
  • 2
  • 4
  • 17

1 Answers1

4

Normally we do this via a while loop, as the number of iterations required is unknown beforehand.

n <- 10L
z <- integer(n)
m <- 1L; i <- 0L
while (m <= n) {
  set.seed(i)
  z_i <- sum(rbinom(2L, 1, 0.5))
  if (z_i > 0L)  {z[m] <- z_i; m <- m + 1L}
  i <- i + 1L
  }

Output:

z
# [1] 1 1 1 1 1 2 1 1 1 1

i
# [1] 14

So we sample 14 times, 4 of which are 0 and the rest 10 are retained.

More efficient vectorized method

set.seed(0)
n <- 10L
z <- rbinom(n, 1, 0.5) + rbinom(n, 1, 0.5)
m <- length(z <- z[z > 0L])  ## filtered samples
p <- m / n  ## estimated success probability
k <- round(1.5 * (n - m) / p)   ## further number of samples to ensure successful (n - m) non-zero samples
z_more <- rbinom(k, 1, 0.5) + rbinom(k, 1, 0.5)
z <- c(z, z_more[which(z_more > 0)[seq_len(n - m)]])

Some probability theory of geometric distribution has been used here. Initially we sample n samples, m of which are retained. So the estimated probability of success in accepting samples is p <- m/n. According to theory of Geometric distribution, on average, we need at least 1/p samples to observe a success. Therefore, we should at least sample (n-m)/p more times to expect (n-m) success. The 1.5 is just an inflation factor. By sampling 1.5 times more samples we hopefully can ensure (n-m) success.

According to Law of large numbers, the estimate of p is more precise when n is large. Therefore, this approach is stable for large n.

If you feel that 1.5 is not large enough, use 2 or 3. But my feeling is that it is sufficient.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • No. They should use vectorization, create more values than needed, and subset to the desired number of non-zero values. – Roland Aug 16 '16 at 18:54
  • 2
    @Roland It's not obvious that the OP's true code can be easily vectorized. Also, the approach you suggest may still fail to reach `n` successes (non-zero results). – Frank Aug 16 '16 at 19:08
  • 2
    @Frank As the update to the answer shows you can get a good estimate of the number of extra values needed. If you don't get the necessary number of "good" values, just run again. You will still be faster than the loop approach for any simulation of decent size. – Roland Aug 16 '16 at 19:43
  • @Roland Sure, I read the OP's "For example..." as just that (an example and not their actual problem), but I guess it's open to interpretation. Obvs, the OP should just `n = 10; nd = 2; sample(nd, n, prob=dbinom(1:nd, nd, .5), replace=TRUE)` if this is their true problem... – Frank Aug 16 '16 at 19:50