simulate zero-truncated Poisson distribution
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
Maybe someone has written an efficient function for this.
If not, you anyway have at least two options with basic R.
(A) "Rejection sampling": OK if the probability of 0 is small.
In that case, suppose you want n zero-truncated Poisson values.
Sample n values from rpois. Reject any which are 0, leaving
you with say r to find. Sample r from rpois, and continue
in the same way until you have all n.
This could be done by something on the following lines:
n<-1000; T<-3.5
Y<-rpois(n,T); Y0<-Y[Y>0]; r<-(n - length(Y0))
while(r>0){
Y<-rpois(r,T); Y0<-c(Y0,Y[Y>0]); r<-(n - length(Y0))
}
and Y is then the required sample of n=1000 from a Poisson
distribution with mean T=3.5, after rejecting all zeros.
(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))
Hoping this helps,
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 01-May-05 Time: 17:20:09
------------------------------ XFMail ------------------------------