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
count simulations
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
Quoted reply history
-----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
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
>
Quoted reply history
> -----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
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
----------------------------------------------------
Quoted reply history
-----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
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
> ----------------------------------------------------
>
Quoted reply history
> -----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
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