RE: [R] Setting max values for rpois

About this list Date view Thread view Subject view Author view Attachment view

From: Ted Harding (Ted.Harding@nessie.mcc.ac.uk)
Date: Thu 13 May 2004 - 07:55:47 EST


Message-id: <XFMail.040512224508.Ted.Harding@nessie.mcc.ac.uk>

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.

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Mon 31 May 2004 - 23:05:09 EST