[R] simulate zero-truncated Poisson distribution
Spencer Graves
spencer.graves at pdf.com
Sun May 1 22:39:20 CEST 2005
Hi, Peter:
Your solution is simple, elegant, very general yet sufficiently subtle
that I (and I suspect many others) might not have thought of it until
you mentioned it. It is worth remembering.
Thanks. spencer graves
Peter Dalgaard wrote:
> (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)
>
> ?
>
More information about the R-help
mailing list