Generating Random Zero-truncated Negative Binomial Values using Rejection Sampling

Question

I am interested in generating zero-truncated negative binomial random variables using some sort of rejection sampling. My first thought was to simply draw from a negative binomial distribution, and then get rid of any zeroes that may arise. I'd then repeat said process until I reached $N$, the desired number of random values. So something along the lines of

rztnegbin = function(n,k,mu){
  size = k
  fill = rep(NA,n)
  for(i in 1:n){
    ztrv = 0
    while(ztrv==0){
      ztrv = rnbinom(1,size=k,mu=mu)
    }
    fill[i] = ztrv
  }
  fill
}

all(rztnegbin(n=100,k=1,mu=5)>0)

Just wanted to make sure that this was appropriate. I'm sure there are more efficient ways out there.

Answer 1

The basic approach* will produce zero-truncated negative binomials.

*(I'm not checking your code, but commenting on the description of the algorithm)

So as far as correctness goes, the approach is completely appropriate.

If $p_0=P(X=0)$ is not large, it may be reasonably efficient. If that number gets much above 50% you'll probably want to consider some other approach, but it depends how many such variables you need (if you only need one thousand or even ten thousand, it probably doesn't matter if it's only 10% efficient - but if you'll need billions over many uses, you'll definitely want to look at other methods).

If speed is your main issue (e.g. you only need a moderate number, but for some reason you need them to arrive fast), you probably should consider coding something up in C or C++ and calling it, though there may already be a function for this somewhere.

You might consider having two different methods, one for when the probability of a zero is not large and another for when it is.

Answer 2

Here is an algorithm that avoids the for() loop, but still allows different means or probabilities:

Y <-rnbinom(n = n, mu = mu, size =Theta )

    while( sum(Y==0) >0 ){ # rejection sampling
    Y.new <- rnbinom(n = sum(Y==0) , mu = mu[which(Y==0)], size =Theta )
    Y[which(Y==0)] <- Y.new
    }
Y