[R] Setting max values for rpois

[(Ted Harding)]

On 12-May-04 Bret Collier wrote:
> R-users,
>          I am simulating a birth process for 4 classes of individuals
> with l[i] being the average No. fetuses per individual.  However, I
> need to bound the resulting values for each generated rpois to be <=3
> (no individual can have > 3 offspring).  I have not been able to
> figure out how to incorporate this into the below example.  Any
> suggestions on integrating  would be appreciated.
> 
> 
> recruit.f <- c(12, 12, 25, 51)  #No. females in each age class
> l <- c(.05, 1.22, 1.6, 1.8)  #mean No. fetuses for each age class
> x <- sapply(lapply(1:4, function(i) rpois(recruit.f[i], l[i])), sum)

It looks as though you are seeking to sample randomly from a Poisson
truncated at 3. This is in effect a multinomial with n=1, so one
approach could be (using your (51,1.8) case as illustration)

  prob<-dpois((0:3),1.8); prob<-prob/sum(prob);
  Nfetuses <- t((0:3))*rmultinom(51,1,prob)

which would give you 51 draws of 1 from the distribution
P(0)=0.1854599, P(1)=0.3338279, P(2)=0.3004451, P(3)=0.1802671
(though there may be a more direct way).

Anyway, whatever method you use to get the sample, you can wrap it
in a function to replace 'rpois' in your expression above.

Ted.