# Generate values from a zero inflated negative binomial fit

I am using the pscl package to fit a simple zeroinf(...,dist = "negbin"). I would like to simulate the probability distribution from this fit.

Using another post on cross validated (Simulate from a zero-inflated poisson distribution) I see the following for the poisson case, but I am not sure what to do for the negative binomial case

library("VGAM")
library("pscl")
object <- zeroinf(...)
p <- predict(object, ..., type = "zero")
lambda <- predict(object, ..., type = "count")
rzipois(n, lambda = lambda, pstr0 = p)


In general, it is simple to generate any zero-inflated pseudo-random variable so long as you can already generate from the corresponding non-zero-inflated version. Suppose you want to generate $$n$$ independent values from the zero-inflated version of a random varible $$Y \sim F$$ with zero-inflation probability $$\theta$$. You can generate the pseudo-random variables $$X_1,...,X_n \sim \text{IID ZI}_\theta F$$ from the corresponding zero-inflated distribution as follows:

• Generate $$Y_1,...,Y_n \sim \text{IID } F$$.
• Generate $$U_1,...,U_n \sim \text{IID Bern}(1-\theta)$$.
• Set the values $$X_1,...,X_n$$ by $$X_i = U_i \cdot Y_i$$.

You can program this procedure in R using the underlying function rpois used for generating pseudo-random variables from the Poisson distribution.

rzipois <- function (n, lambda, zi.prob) {
Y <- rpois(n, lambda)
U <- sample(c(0, 1), size = n, prob = c(zi.prob, 1-zi.prob), replace = TRUE)
U*Y }


You can use the manual approach from Simulate from a zero-inflated poisson distribution and replace the rpois() call inside the ifelse() by an rnbinom() call.

We have also included a somewhat streamlined version of that code as rzinbinom() in the countreg package on R-Forge.