RE: [R] simulate zero-truncated Poisson distribution

From: Ted Harding <Ted.Harding_at_nessie.mcc.ac.uk>
Date: Mon 02 May 2005 - 02:24:39 EST


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@nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861
Date: 01-May-05                                       Time: 17:20:09
------------------------------ XFMail ------------------------------

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon May 02 02:51:54 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:31:31 EST