RE: Re: count simulations
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