From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>

Date: Mon 02 May 2005 - 05:35:22 EST

Date: Mon 02 May 2005 - 05:35:22 EST

(Ted Harding) <Ted.Harding@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)

?

-- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ______________________________________________ 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.htmlReceived on Mon May 02 05:41:57 2005

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