Simulating zero-truncated Poisson

I have been working on simulating data from Bob Dorazio’s version of the Royle N-mixture model for repeated count data. This version of the model assumes that the abundance at site is generated from a Hurdle model such that and , . This can be interpreted as specifying an occupancy rate of with abundance given by a zero-truncated Poisson with mean at the occupied sites.

I was originally generating the zero-truncated Poisson values using a while loop to simulate from a Poisson with mean until a non-zero value was generated. This is inefficient and runs into problems if is sufficiently small that the probability of a non-zero value is below computer precision.

To remedy this I found this nice solution from Peter Daalgard in response to a question from Spencer Graves on the R-help list. This solution essentially uses the inverse probability transform but truncates the CDF so that the values generated are guaranteed to be greater than 1. The same could be done for any distribution (that has a distribution function implemented in R). My code below implements this trick but adds a tolerance that avoids problems if the value of lambda is too small.

rtpois <- function(n, lambda,tol=1e-10){
## Simulate from zero-truncated Poisson

## Initialize output
x <- rep(NA,n)

## Identify lambda values below tolerance
low <- which(lambda < tol)
nlow <- length(low)

if(nlow > 0){
x[low] <- 1

if(nlow < n)
x[-low] <- qpois(runif(n-nlow, dpois(0, lambda[-low]), 1), lambda[-low])
}
else
x <- qpois(runif(n-nlow, dpois(0, lambda), 1), lambda)

return(x)
}