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 }

