Re: Re: count simulations
Mats,
A uniform random deviate from the interval 0-1 is by convention allowed to include 0 but not include 1. So in theory R could be 0 and LOG(0) would be invalid, however 1-R if would be OK if R was 0.
I don't understand what triggered the error message that Elodie shows below. It implies that the language processor (is it NM-TRAN or FORTRAN?) 'knows' that R could be 0 and therefore does not allow LOG(R). But surely both NM-TRAN and FORTRAN allow the LOG of any variable and it is only at run-time that a variable <=0 would throw an exception?
In practice the probability of having R=0 is approaching 0 so I cannot see that there is a real risk of LOG(R) being a problem. On the other hand there is certainly a performance overhead in computing 1-R everytime when R would do just as well.
Nick
Mats Karlsson wrote:
> Hi Nick,
>
> Elodie came up with one reason why it may have been 1-R rather than R (see
> below). I rationalize the code as simulating one event after another on a
> time interval standardized by lambda. You sample a survival probability,
> translate that into a time, check if it is beyond the standardized interval,
> if not increase N and add the event time to the elapsed time in the
> interval.
> SUR = EXP(-LAMB*T)
> R = EXP(-LAMB*T)
> LOG(R) = -LAMB*T
> T = -LOG(R)/LAMB
>
> Here is elodies mail
> ;------------------
> Really in theory having (R) or (1-R) cannot give simulated datasets more
> different to each other than if they were simulated with 2 different seed
> number.
>
> But there was a reason in practice in nonmem because having only (R) gives:
> FSUBS.f: In subroutine `pred':
>
> FSUBS.f:34: B00006=DLOG(R) ^
>
> Reference to intrinsic `DLOG' at (^) invalid -- one or more arguments have
> incorrect type
> No nonmem execution.
>
> Is it that nonmem can give random number 0 and not 1? I can see only the
> numerical issue of log(0).
>
> ;------------------
>
> Best regards,
> Mats
>
> Mats Karlsson, PhD
> Professor of Pharmacometrics
> Dept of Pharmaceutical Biosciences
> Uppsala University
> Box 591
> 751 24 Uppsala Sweden
> phone: +46 18 4714105
> fax: +46 18 471 4003
>
Quoted reply history
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Nick Holford
> Sent: Thursday, August 06, 2009 9:48 PM
> To: nmusers
> Subject: [NMusers] Re: count simulations
>
> Mats,
>
> Thanks for pointing out that R and 1-R are equivalent when R is a uniform 0-1 random deviate.
>
> There is an NM-TRAN example using 1-R in this paper:
>
> Frame B, Miller R, Lalonde RL. Evaluation of Mixture Modeling with Count Data using NONMEM. Journal of Pharmacokinetics and Pharmacodynamics. 2003;30(3):167-83.
>
> I have to admit to having cut and pasted this example and used it to show others how to simulate count data so it may have propogated that way too.
>
> Do you know of a clear explanation of why this simple algorithm produces Poisson distribution samples?
>
> Nick
>
> Mats Karlsson wrote:
>
> > Dear both,
> >
> > You have both simulated count data using the code below (or very similar). My question is why do you use LOG(1-R) rather than the simpler LOG(R)? If you've done it because you inherited the code, where did you get the code.
> >
> > *IF (ICALL.EQ.4) THEN*
> >
> > * T=0*
> >
> > * N=0*
> >
> > * DO WHILE (T.LT.1)*
> >
> > * CALL RANDOM (2,R)*
> >
> > * T=T-LOG(1-R)/LAMB*
> >
> > * IF (T.LT.1) N=N+1*
> >
> > * END DO*
> >
> > * DV=N*
> >
> > *ENDIF*
> >
> > Best regards,
> >
> > Mats
> >
> > Mats Karlsson, PhD
> >
> > Professor of Pharmacometrics
> >
> > Dept of Pharmaceutical Biosciences
> >
> > Uppsala University
> >
> > Box 591
> >
> > 751 24 Uppsala Sweden
> >
> > phone: +46 18 4714105
> >
> > fax: +46 18 471 4003
--
Nick Holford, Professor Clinical Pharmacology
Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
[email protected] tel:+64(9)923-6730 fax:+64(9)373-7090
mobile: +64 21 46 23 53
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford