simulate zero-truncated Poisson distribution

(Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> writes:

On 01-May-05 Glazko, Galina wrote:

Dear All

I would like to know whether it is possible with R to 
generate random numbers from zero-truncated Poisson 
distribution.

Thanks in advance,
Galina

.....

(B) This has a deeper theoretical base.

Suppose the mean of the original Poisson distribution (i.e.
before the 0's are cut out) is T (notation chosen for intuitive
convenience).

The number of events in a Poisson process of rate 1 in the
interval (0,T) has a Poisson distribution with mean T.
The time to the first event of a Poisson process of rate 1
has the negative exponential distribution with density exp(-t).

Conditional on the first event lying in (0,T), the time
to it has the conditional distribution with density

 exp(-t)/(1 - exp(-T)) (0 <= t <= T)

and the PDF (cumulative distribution) is

 F(t) = (1 - exp(-t))/(1 - exp(-T))

If t is randomly sampled from this distribution, then
U = F(t) has a uniform distribution on (0,1). So, if
you sample U from runif(1), and then

 t = -log(1 - U*(1 - exp(-T)))

you will have a random variable which is the time to
the first event, conditional on it being in (0,T).

Next, the number of Poisson-process events in (t,T),
conditional on t, simply has a Poisson distribution
with mean (T-t).

So sample from rpois(1,(T-t)), and add 1 (corresponding to
the "first" event whose time you sampled as above) to this
value.

The result is a single value from a zero-truncated Poisson
distribution with (pre-truncation) mean T.

Something like the following code will do the job vectorially:

 n<-1000       # desired size of sample
 T<-3.5        # pre-truncation mean of Poisson
 U<-runif(n)   # the uniform sample
 t = -log(1 - U*(1 - exp(-T))) # the "first" event-times
 T1<-(T - t)   # the set of (T-t)

 X <- rpois(n,T1)+1 # the final truncated Poisson sample

The expected value of your truncated distribution is of course
related to the mean of the pre-truncated Poisson by

 E(X) = T/(1 - exp(-T))

There must be an easier way... Anything wrong with

 rtpois <- function(N, lambda)
     qpois(runif(N, dpois(0, lambda), 1), lambda)

 rtpois(100,5)

?