RE: Re: count simulations

From: Bob Leary Date: August 10, 2009 technical Source: mail-archive.com
I suspect that the reason that log(1-R) works but log(R) does not is the real*4 type problem, not the fact that R may take on the value 0 (I do not know of any pseudorandom uniform number generators where this can happen). The 1 may be treated as a double precision constant, in which case fortran casts the mixed quantity 1-R (real*8 - real*4) as a real*8 quantity, which agrees in type with what dlog is expecting. Robert H. Leary, PhD Fellow Pharsight - A Certara(tm) Company 5625 Dillard Dr., Suite 205 Cary, NC 27511 Phone/Voice Mail: (919) 852-4625, Fax: (919) 859-6871 Email: [email protected] > This email message (including any attachments) is for the sole use of the > intended recipient and may contain confidential and proprietary information. > Any disclosure or distribution to third parties that is not specifically > authorized by the sender is prohibited. If you are not the intended > recipient, please contact the sender by reply email and destroy all copies of > the original message.
Quoted reply history
-----Original Message----- From: [email protected] [mailto:[email protected]]on Behalf Of Nick Holford Sent: Monday, August 10, 2009 1:51 AM To: 'nmusers' Subject: Re: [NMusers] Re: count simulations Elodie, I've tried using LOG(R) in NONMEM 6 with the Intel V 10 and G77 compilers. T=0 NEVENT=0 DO WHILE (T.LT.1) CALL RANDOM (2,R) T=T-LOG(R)/COUNT IF (T.LT.1) NEVENT=NEVENT+1 ENDDO DV=NEVENT The Intel compiler gave this error: FSUBS.obj : error LNK2019: unresolved external symbol _DLOG referenced in function _ERROR count_sim_R.exe : fatal error LNK1120: 1 unresolved externals The G77 compiler gave this error: FSUBS.for: In subroutine `error': FSUBS.for:221: B00009=DLOG(R) ^ Reference to intrinsic `DLOG' at (^) invalid -- one or more arguments have incorrect type With NONMEM 7 the Intel V 10 and Gfortran compilers had no problems. This seems to be a problem caused by the code generated by NONMEM 6. R is declared as a REAL by NONMEM 6 but the DLOG function wants a DOUBLE and the compilers won't play the game. NONMEM 7 declares REAL(KIND=DPSIZE) :: R which seems to let DLOG work OK. In any event its safer (to avoid LOG(0)) to use LOG(1-R) even in NONMEM 7. Nick Elodie Plan wrote: > Dear Nick, > > I agree, it must be "remembered" somewhere that log(R) is forbidden, because > the simple following trick works (but won't once every 2 billion times, when > 0 is picked): > IF (ICALL.EQ.4) THEN > T=0 > N=0 > DO WHILE (T.LT.1) > CALL RANDOM (2,R) > RA=R > T=T-LOG(RA)/LAMB > IF (T.LT.1) N=N+1 > END DO > DV=N > ENDIF > > Best regards, > Elodie > > Elodie L. Plan, PharmD, MSc, > PhDstudent > ******************************************* > Department of Pharmaceutical Biosciences, > Faculty of Pharmacy, Uppsala University > Box 591 - 751 24 Uppsala - SWEDEN > Office +46 18 4714385 - Fax +46 18 4714003 > ---------------------------------------------------- > > > -----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
Aug 06, 2009 Nick Holford Re: count simulations
Aug 07, 2009 Mats Karlsson RE: Re: count simulations
Aug 07, 2009 Nick Holford Re: Re: count simulations
Aug 07, 2009 Elodie Plan RE: Re: count simulations
Aug 10, 2009 Nick Holford Re: Re: count simulations
Aug 10, 2009 Bob Leary RE: Re: count simulations