From: Ted Harding (Ted.Harding@nessie.mcc.ac.uk)
Date: Thu 13 May 2004 - 07:55:47 EST
On 12-May-04 Bret Collier wrote:
> 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)
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.
Remail@example.com mailing list
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
This archive was generated by hypermail 2.1.3 : Mon 31 May 2004 - 23:05:09 EST