NONMEM Users Network Archive

Hosted by Cognigen

RE: Re: count simulations

From: Elodie Plan <elodie.plan>
Date: Fri, 7 Aug 2009 14:41:34 +0200

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: owner-nmusers
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: owner-nmusers
[mailto:owner-nmusers
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
n.holford
mobile: +64 21 46 23 53
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford

Received on Fri Aug 07 2009 - 08:41:34 EDT

The NONMEM Users Network is maintained by ICON plc. Requests to subscribe to the network should be sent to: nmusers-request@iconplc.com.

Once subscribed, you may contribute to the discussion by emailing: nmusers@globomaxnm.com.