+--------------------------------------------------------------------+
| |
| BIVARIATE FUNCTION |
| |
+--------------------------------------------------------------------+
This is a fully-worked out example of the use of the $ABBREVIATED
FUNCTION option to declare a user-defined function. In this example,
the function BIVARIATE is used to compute the integral of the bivari-
ate normal distribution of two correlated data points, DV1 and DV2,
with correlation RHO.
This is bivariate.ctl from the NONMEM 7 distribution medium. This
file, the data file bivariate.csv, and the function's source code
bivariate.f90 can be found in the examples directory.
$PROB BIVARIATE EXAMPLE
; THESE DECLARATIONS ALLOW ANY FUNCTION TO HAVE
; ALTERNATIVE DIMENSIONS FOR THEIR ARRAYS
; BUT, USER DEFINED DIMENSIONS ARE PASSED AS THE
; LAST ARGUMENT TO FUNC, SUCH AS:
; BV=BIVARIATE(VBI(1),FNC001_1(1,1),FNC001_2(1,1,1),5)
$ABBR FUNCTION BIVARIATE(VBI,5)
$INPUT SIM ID DOSE DV TIME
$DATA bivariate.csv IGNORE=
$SUBROUTINES OTHER=bivariate.f90
$PRED
B1=THETA(1)
B2=THETA(2)
B3=THETA(3)
K =LOG(2)/EXP(THETA(4))
ED50=EXP(THETA(5))
U =(1-EXP(-K*TIME ))
MU_1=B1+B3*DOSE/(DOSE+ED50)
MU_2=B2
MXB=MU_1+ETA(1)
MXU=MU_2+ETA(2)
MX =MXB + MXU*U ;***Current model prediction***;
PHIMX=PHI(MX)
IF(NEWIND.NE.2) THEN
TIMEP=0
MXP=0
DVP=0
PHIMXP=0.5
ENDIF
RHOB=(2/(1+EXP(-THETA(8)))-1)
IF(RHOB>0.0) RHO=RHOB**(TIME-TIMEP)
IF(RHOB==0.0) RHO=0.0
IF(RHOB<0.0) RHO=-(-RHOB)**(TIME-TIMEP)
PC =(1-PHIMX) *(1-DV ) + PHIMX*DV
IF(PC.LE.0.0) EXIT
V=SQRT(1+OMEGA(1,1)+OMEGA(2,2)*U**2)
POPP = (B1+B2*U +B3*DOSE/(DOSE+ED50))/V ;*Population mean prediction*;
IF (TIME.EQ.1) THEN
JP=PC
PCP=1.0
ELSE
;***Pass information to bivariate normal***;
VBI(1)=RHO
VBI(2)=MX
VBI(3)=MXP
VBI(4)=1 ;*0 = Upper tail as in Drezner & Wesolowsky; 1 = Bottom tail*;
VBI(5)=1 ;*0 = 3 pt approximation; 1 = 5 point approximation*;
BV=BIVARIATE(VBI)
JP=((DV-1)*(DVP-1)+(DV-1)*(1-2*DVP)*PHIMXP+(DVP-1) &
*(1-2*DV)*PHIMX+(1-2*DV)*(1-2*DVP)*BV)
ENDIF
IF(JP.LE.0.0) EXIT
LOGL=LOG(JP/PCP)
Y = -2*LOGL
MXP=MX
PCP=PC
DVP=DV
TIMEP=TIME
PHIMXP=PHIMX
$THETA
-1.7 ; 1 B1
1.2 ; 2 B2
2.9 ; 3 B3
1.4 ; 4 LOG(B4)
1.2 ; 5 LOG(B5)
(0.0 FIXED) ; 6 LOG SQRT VAR(ETA1)
(0.0 FIXED) ; 7 LOG SQRT VAR(ETA2)
2.2 ; 8 RHO parameter
$OMEGA DIAGONAL(2)
0.8 ; V1
0.8 ; V2
;$EST METHOD=IMP LAPLACE -2LL PRINT=1 NITER=300 ISAMPLE=300
;SIGL=6 CTYPE=3 NOHABORT
$EST MAX=0 PRINT=1 METHOD=1 LAPLACE -2LL SIGL=10 NOHABORT
$COV COMPRESS MATRIX=R PRINT=E UNCONDITIONAL
Here is a fragment of the data file bivariate.csv:
SIM,ID,DOSE,DV,TIME
1,1,0,0,1
1,1,0,0,2
1,1,0,0,4
1,1,0,0,8
1,1,0,0,16
1,1,0,0,24
1,1,0,0,36
This is the file bivariate.f90. (The extension .f90 in the file name
is not actually required because the contents of the file is copied
to FSUBS and compiled as part of FSUBS. It is not compiled indepen-
dently.)
! Integral of Bivariate normal distribution of two correlated data points,
! DV1, DV2, correlation RHO.
! BV1, BV2 and GAUSS FROM Drezner and Wesolowsky,
! J. Statistical Computat. Simul. 35, pp. 101-107, 1990.
! BV2 is more accurate for extreme values of rho
! Partial derivatives worked out by Robert Bauer.
!
! How to use in NONMEM 7.4
! VBI is a suggested vector name. Could be any name.
! But the function name must match what is in bivariate.f90
! $ABBR FUNCTION BIVARIATE(VBI,5)
! $SUBROUTINES ... OTHER=bivariate.f90
! $PK
! ...
! DV1 AND DV2 SHOULD BE TWO SEPARATE DATA ITEMS ON EACH RECORD.
! VBI(1)=RHO
! VBI(2)=DV1 ! H
! VBI(3)=DV2 !K
! VBI(4)=DTYPE !=0 FOR H TO INF, K TO INF, OR 1 FOR -INF TO H, -INF TO K
! VBI(5)=BVTYPE ! BVTYPE=0 USES SIMPLE INTEGRATOR, BVTYPE=1 USES MORE
! ACCURATE ONE FOR LARGE VALUES OF RHO
! F_FLAG=1
! IT IS SAFER TO AVOID BV BEING EXACTLY 0, JUST MAKE SURE
! IT IS REALLY SMALL, TO AVOID LOG() ERRORS
! BV=BIVARIATE(VBI)+1.0D-30
! IF YOU USE -2LL FORMAT (DEFAULT):
! Y=-2.0D+00*LOG(BV)
FUNCTION BIVARIATE(X,X1,X2,NDIM)
! NDIM SHOULD EQUAL 5
! IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT NONE
INTEGER NDIM
REAL*8 X(NDIM),X1(NDIM),X2(NDIM,NDIM)
REAL*8 R,H,K,ITYPE,PI,PI2,SPI2,R1,RSQR,HRK,KRH,PHRK,PKRH,XH,XK
REAL*8 XHRK,XKRH,XC,BV,BV2,BIVARIATE,GAUSS,BV1
EXTERNAL GAUSS,BV1,BV2
INTEGER BVTYPE
! unload the arguments
R=X(1)
H=X(2)
K=X(3)
ITYPE=X(4)
BVTYPE=X(5)
! ITYPE=0 INTGRATE H TO INF, K TO INF
! ITYPE=1 INTGRATE -INF TO H -INF TO K
! BVTYPE=0 USES THE SIMPLER BV1
! BVTYPE=1 USES THE MORE ACCURATE BV2
! SET UP SOME USEFUL CALCUATIONS FOR FIRST
! AND SECOND DERIVATIVES OF THE ARGUMENTS
! DERIVATIVES NOT NEEDED IF YOU USE IMP, SAEM, BAYES, OR IF
! YOU USE LAPLACE OR ITS WITH OPTMAP=1 ETADER=3 (NONMEM 7.3), OR
! DERIVATIVES NEEDED IF YOU USE STANDARD LAPLACE
PI=3.141592653589793238D+00
PI2=2.0D+00*pi
spi2=sqrt(pi2)
R1=1.0D+00-R*R
RSQR=SQRT(ABS(R1))
IF(RSQR.LE.0.0D+00) RSQR=1.0D-100
IF(R1.LE.0.0D+00) R1=1.0D-100
HRK=(H-R*K)/RSQR
KRH=(K-R*H)/RSQR
PHRK=1.0D+00-GAUSS(HRK)
PKRH=1.0D+00-GAUSS(KRH)
XH=EXP(-H*H/2.0D+00)/SPI2
XK=EXP(-K*K/2.0D+00)/SPI2
XHRK=EXP(-HRK*HRK/2.0D+00)/SPI2/RSQR
XKRH=EXP(-KRH*KRH/2.0D+00)/SPI2/RSQR
XC=EXP(-(H*H-2.0D+00*R*H*K+K*K)/2.0D+00/R1)/RSQR/PI2
! parital F WRT RHO, from Drezner and Wesolowsky, equation (4)
X1(1)=XC
! parital F WRT H
X1(2)=XH*PKRH+(ITYPE-1.0D+00)*XH
! parital F WRT K
X1(3)=XK*PHRK+(ITYPE-1.0D+00)*XK
! 2ND parital F WRT RHO,RHO
X2(1,1)=R/R1*XC+XC*(H*K*(1.0D+00-R)*(1.0D+00-R)-R*(H-K)*(H-K))/R1/R1
! 2ND parital F WRT RHO,H
X2(1,2)=XC*(R*K-H)/R1
! 2ND parital F WRT RHO,K
X2(1,3)=XC*(R*H-K)/R1
X2(2,1)=X2(1,2)
X2(3,1)=X2(1,3)
! 2ND parital F WRT H,K
X2(2,3)=XC
X2(3,2)=X2(2,3)
! 2ND parital F WRT H,H
X2(2,2)=-X1(2)*H-XC*R
! 2ND parital F WRT K,K
X2(3,3)=-X1(3)*K-XC*R
! write(50,*) 'A ',x1(1),x1(2),x1(3)
! write(50,*) 'B ',x2(1,1),x2(1,2),x2(1,3),x2(2,1), &
! x2(2,2),x2(2,3),x2(3,1),x2(3,2),x2(3,3)
! bV2 is more accurate for exterme rho values.
IF(BVTYPE==0) THEN
BV=BV1(H,K,R)
ELSE
BV=BV2(H,K,R)
ENDIF
IF(ITYPE==1.0D+00) BV=BV-GAUSS(H)-GAUSS(K)+1.0D+00
BIVARIATE=BV
RETURN
END
FUNCTION BV1(H1,H2,R)
IMPLICIT NONE
INTEGER I
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 H1,H2,R,H12,BV,BV1,RR,RR2,H3,GAUSS
EXTERNAL GAUSS
REAL*8 X(5),W(5),DVAL
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H12 =(H1*H1+ H2*H2)/2.0D+00
H3 = H1*H2
BV = 0.0D+00
DO 1 I= 1,5
RR = R*X(I)
RR2= 1.0D+00-RR*RR
IF(RR2.LE.0.0D+00) RR2=1.0D-100
DVAL=(RR*H3 - H12)/RR2
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
BV = BV + W(I)*EXP(DVAL)/SQRT(RR2)
1 CONTINUE
BV = BV*R + GAUSS(H1)*GAUSS(H2)
BV1=BV
RETURN
END
FUNCTION BV2(H1,HK,R)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 X(5),W(5)
REAL*8 H1,HK,R,H2,H12,BV,BV2,R2,R3,H3,H7,H6,H5,AA,AB,R1,RR,GAUSS,RR2
REAL*8 DVAL,DVAL2
EXTERNAL GAUSS
INTEGER I
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H2 = HK
H12 = (H1*H1+ H2*H2)/2.0D+00
BV = 0.0D+00
IF(ABS(R).LT.0.7D+00)GO TO 4
R2= 1.0D+00-R*R
R3 = SQRT(R2)
IF(R.LT.0.0D+00)H2 = -H2
H3=H1*H2
H7 = EXP( -H3/2.0D+00)
IF(R2.EQ.0.0D+00) GO TO 3
H6 = ABS(H1 - H2)
H5 = H6*H6/2.0D+00
H6 = H6/R3
AA = 0.5D+00 - H3/8.0D+00
AB = 3.0D+00 - 2.0D+00*AA*H5
BV = .13298076D+00*H6*AB*GAUSS(H6) &
-EXP(-H5/R2)*(AB + AA*R2)*.053051647D+00
DO 2 I= 1,5
R1= R3*X(I)
RR = R1*R1
R2 = SQRT(DABS(1.0D+00 - RR))
IF(R2.EQ.0.0D+00) R2=1.0D-100
IF(RR.EQ.0.0D+00) RR=1.0D-100
DVAL=-H5/RR
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
DVAL2=-H3/(1.0D+00 + R2)
IF(DVAL2.GT.300.0D+00) DVAL2=300.0D+00
BV = BV - W(I)*EXP(DVAL)*(EXP(DVAL2)/R2/H7 - 1.0D+00 - AA*RR)
2 CONTINUE
3 IF(R.GT.0.0D+00)BV = BV*R3*H7 + GAUSS(MAX(H1,H2))
IF(R.LT.0.0D+00)BV = MAX(0.0D+00,GAUSS(H1) - GAUSS(H2)) - BV*R3*H7
BV2=BV
RETURN
4 H3=H1*H2
DO 1 I= 1,5
R1= R*X(I)
RR2= 1.0D+00-R1*R1
IF(RR2.EQ.0.0D+00) RR2=1.0D-100
1 BV = BV + W(I)*EXP((R1*H3 -H12)/RR2)/SQRT(RR2)
BV2 =GAUSS(H1)*GAUSS(H2) + R*BV
RETURN
END
FUNCTION GAUSS(Z)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Z,GAUSS,X,G
REAL*8 A(4)
INTEGER I
DATA A/ -.72657601,.71070688,-.142248368,.127414796/
X= 1.0D+00/(1.0D+00 + .23164189D+00*ABS(Z))
G = .53070271D+00
DO 1 I= 1,4
1 G=G*X+A(I)
GAUSS = G*X*EXP( - Z*Z/2.0D+00)
IF(Z.LT.0.0D+00)GAUSS = 1.0D+00 -GAUSS
RETURN
END
! USE PROGRAM HEADER FOR STAND-ALONE EXECUTABLE TESTING
! PROGRAM BIVTEST
SUBROUTINE BIVTEST
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(9),A1(9),A2(9,9)
REAL*8 AH(9),AH1(9),AH2(9,9)
REAL*8 AJ(9),AJ1(9),AJ2(9,9)
REAL*8 B(9),B1(9),B2(9,9)
REAL*8 H,K
HDEL=1.0D-04
1 CONTINUE
WRITE(*,*) 'ENTER RHO,H,K,INTEGRAND-TYPE,BVTYPE'
READ(*,*) RHO,H,K,DTYPE,BVTYPE
A(1)=RHO
A(2)=H
A(3)=K
A(4)=DTYPE
A(5)=BVTYPE
BV=FUNCA(A,A1,A2)
WRITE(*,*) 'VALUE ',BV
WRITE(*,*) GAUSS(-1.96D+00),GAUSS(1.96D+00)
AH(1:5)=A(1:5)
AH(1)=RHO-HDEL
AJ(1:5)=A(1:5)
AJ(1)=RHO+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(1)=(BVJ-BVH)/2.0D+00/HDEL
B2(1,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(1,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(1,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG RHO',A1(1),A2(1,1),A2(1,2),A2(1,3)
WRITE(*,*) 'NUMER RHO',B1(1),B2(1,1),B2(1,2),B2(1,3)
AH(1:5)=A(1:5)
AH(2)=H-HDEL
AJ(1:5)=A(1:5)
AJ(2)=H+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(2)=(BVJ-BVH)/2.0D+00/HDEL
B2(2,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(2,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(2,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG H',A1(2),A2(2,1),A2(2,2),A2(2,3)
WRITE(*,*) 'NUMER H',B1(2),B2(2,1),B2(2,2),B2(2,3)
AH(1:5)=A(1:5)
AH(3)=K-HDEL
AJ(1:5)=A(1:5)
AJ(3)=K+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(3)=(BVJ-BVH)/2.0D+00/HDEL
B2(3,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(3,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(3,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG K',A1(3),A2(3,1),A2(3,2),A2(3,3)
WRITE(*,*) 'NUMER K',B1(3),B2(3,1),B2(3,2),B2(3,3)
GO TO 1
END
REFERENCES: Guide Introduction_7
Go to main index.
Created by nmhelp2html v. 1.0 written by Niclas Jonsson (Modified by AJB 5/2006,11/2007,10/2012)