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

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:

(see

on a

probability,

interval,

more

seed

gives:

have

the

[mailto:owner-nmusers

On

Count

produces

--

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