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) <[EMAIL PROTECTED]> 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)
?
______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
