# Re: Re: count simulations

From: Nick Holford <n.holford>
Date: Fri, 07 Aug 2009 23:43:04 +1200

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
> 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 - 07:43:04 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.