From: Mats Karlsson <*mats.karlsson*>

Date: Fri, 7 Aug 2009 10:38:29 +0200

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 - 04:38:29 EDT

Date: Fri, 7 Aug 2009 10:38:29 +0200

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:

--

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 - 04:38:29 EDT