Re: Re: count simulations

> -----Original Message----- > From: [email protected] [mailto:[email protected]] On > Behalf Of Nick Holford > Sent: Friday, August 07, 2009 1:43 PM > To: nmusers > Subject: Re: [NMusers] 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 > > > > -----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