!
! ************************************************************************
!
! This program (Fortran 90) is used to evaluate the critical values for the testing
! procedure given in Huang, Wen, Cheung, and Kwong.
! "Non-inferiority studies with multiple reference treatments" (submitted to
! Statistical Methods in Medical Research)
!
! Input Variables:
!
! (1) alpha: familywise error rate (FWE) (real)
! (2) k: the number of new (experimantal) treatments (integer)
! (3) sss: the number of reference treatments (integer)
! (4) nE1, ..., nEk: the sample sizes of new treatments (integers)
! (5) nR1, ..., nRs: the sample sizes of reference treatments (integers)
! (6) np: the sample size ofthe placebo (integer)
!
! Output variable: required critical value (real)
!
! Modified on January 28, 2015
!
! Remark: Possible values of sss = 1, 2, 3, 4 (Note that in practice sss=4 is very unusual)
! For sss = 1, the obtained critical values are identical to those povided
! in the paper by Kwong KS, Cheung SH, Hayter A and Wen MJ (2012)
! Extension of three-arm non-inferiority studies to trials with multiple new treatments.
! Statistics in Medicine, 31, 2833-2843
! For sss = 4, it will take much longer to compute the critical value due to
! the need to compute very high dimensional integration. The exact time depends
! on the machine being used. Again, sss=4 is extremely unusual in practice. For sss = 3,
! the computational time, with an ordinary PC
! is around 10 minutes and for sss = 2, it is about 10 seconds.
! *************************************************************************
! Main program
PROGRAM main
IMPLICIT NONE
INTEGER i, sss, k, M, j
REAL alpha, v, nnk(80), nns(80), np, D, cp
REAL C0rho(80), C0rhos(80), C1rho(80), C1rhos(80)
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
REAL Time_start, Time_end
WRITE(*,*) 'Please input a familywise error rate (alpha):'
READ(*,*) alpha
WRITE(*,*) 'Please input a number for the group number of new treatments (k):'
READ(*,*) k
WRITE(*,*) 'Please input a number for the group number of reference treatments (s):'
READ(*,*) sss
DO i = 1, k
WRITE(*,*) "Please input the sample sizes for new treatment" ,i, ":"
READ(*,*) nnk(i)
END DO
DO i = 1, sss
WRITE(*,*) 'Please input the sample sizes for reference treatment',i, ':'
READ(*,*) nns(i)
END DO
WRITE(*,*) 'Please input the sample for the placebo group:'
READ(*,*) np
WRITE(*,*)
WRITE(*,*) 'Please wait...'
WRITE(*,*)
OPEN(unit = 11, file = '11cp.txt')
CALLCPU_Time(Time_start)
alpha = real(alpha)
k = int(k)
sss = int(sss)
nnk = real(nnk)
nns = real(nns)
np = real(np)
M = k + 1
v = sum(nnk) + sum(nns) + np - (k + sss + 1)
DO j = 1, sss
DO i = 1, k
RM(i, j) = sqrt(nnk(i)/(nnk(i) + nns(j)))
RsM(i, j) = sqrt(nns(j)/(nnk(i) + nns(j)))
END DO
GM(j) = sqrt(nns(j)/(nns(j) + np))
GsM(j) = sqrt(np/(nns(j) + np))
END DO
D = 1.96
CALL GGEN(D, RM, RsM, Gm, Gsm, alpha, M, sss, v)
cp = D
CALLCPU_Time(Time_end)
WRITE(*, '(a30, f10.4)') "critical point = ", cp
! WRITE(*, '(a30, f10.4)') "Using time (sec) =", Time_end - Time_start
WRITE(11, '(a30, f10.4)') "critical point = ", cp
! WRITE(11, '(a30, f10.4)') "Using time (sec) =", Time_end - Time_start
STOP
ENDPROGRAM
!
! ************************************************************************
!
! MODULE
!
! ************************************************************************
MODULE premain
INTEGER, PARAMETER :: N = 80
REAL, PARAMETER :: SQRT2PI = sqrt(2*3.141592654)
REAL*8 QX(N), QW(N)
DATA QX /-1.000000000000000000000000000000D0, -0.998838663101196289062500000000D0, &
-0.996108710765838623046875000000D0, -0.991822957992553710937500000000D0, &
-0.985988378524780273437500000000D0, -0.978614449501037597656250000000D0, &
-0.969712495803833007812500000000D0, -0.959296405315399169921875000000D0, &
-0.947382450103759765625000000000D0, -0.933989346027374267578125000000D0, &
-0.919137895107269287109375000000D0, -0.902851343154907226562500000000D0, &
-0.885155081748962402343750000000D0, -0.866076767444610595703125000000D0, &
-0.845646142959594726562500000000D0, -0.823895156383514404296875000000D0, &
-0.800857722759246826171875000000D0, -0.776569962501525878906250000000D0, &
-0.751069545745849609375000000000D0, -0.724396467208862304687500000000D0, &
-0.696592390537261962890625000000D0, -0.667700529098510742187500000000D0, &
-0.637766242027282714843750000000D0, -0.606836140155792236328125000000D0, &
-0.574958503246307373046875000000D0, -0.542183101177215576171875000000D0, &
-0.508561253547668457031250000000D0, -0.474145233631134033203125000000D0, &
-0.438989043235778808593750000000D0, -0.403147280216217041015625000000D0, &
-0.366676121950149536132812500000D0, -0.329632431268692016601562500000D0, &
-0.292074084281921386718750000000D0, -0.254059702157974243164062500000D0, &
-0.215648576617240905761718750000D0, -0.176900759339332580566406250000D0, &
-0.137876763939857482910156250000D0, -0.098637439310550689697265625000D0, &
-0.059244122356176376342773437500D0, -0.019758338108658790588378906250D0, &
0.019758308306336402893066406250D0, 0.059244148433208465576171875000D0, &
0.098637476563453674316406250000D0, 0.137876763939857482910156250000D0, &
0.176900774240493774414062500000D0, 0.215648591518402099609375000000D0, &
0.254059672355651855468750000000D0, 0.292074084281921386718750000000D0, &
0.329632490873336791992187500000D0, 0.366676181554794311523437500000D0, &
0.403147339820861816406250000000D0, 0.438989043235778808593750000000D0, &
0.474145323038101196289062500000D0, 0.508561253547668457031250000000D0, &
0.542183220386505126953125000000D0, 0.574958503246307373046875000000D0, &
0.606836199760437011718750000000D0, 0.637766242027282714843750000000D0, &
0.667700469493865966796875000000D0, 0.696592330932617187500000000000D0, &
0.724396467208862304687500000000D0, 0.751069605350494384765625000000D0, &
0.776569962501525878906250000000D0, 0.800857841968536376953125000000D0, &
0.823895215988159179687500000000D0, 0.845646202564239501953125000000D0, &
0.866076767444610595703125000000D0, 0.885155081748962402343750000000D0, &
0.902851343154907226562500000000D0, 0.919137954711914062500000000000D0, &
0.933989346027374267578125000000D0, 0.947382569313049316406250000000D0, &
0.959296464920043945312500000000D0, 0.969712615013122558593750000000D0, &
0.978614509105682373046875000000D0, 0.985988497734069824218750000000D0, &
0.991823017597198486328125000000D0, 0.996108710765838623046875000000D0, &
0.998838663101196289062500000000D0, 1.000000000000000000000000000000D0/
DATA QW /0.000316442543407902121543884277D0, 0.001950567937456071376800537109D0, &
0.003508417401462793350219726562D0, 0.005061226896941661834716796875D0, &
0.006605997681617736816406250000D0, 0.008140315301716327667236328125D0, &
0.009661571122705936431884765625D0, 0.011167781427502632141113281250D0, &
0.012656762264668941497802734375D0, 0.014125877059996128082275390625D0, &
0.015572983771562576293945312500D0, 0.016996013000607490539550781250D0, &
0.018392089754343032836914062500D0, 0.019759664312005043029785156250D0, &
0.021096242591738700866699218750D0, 0.022400053218007087707519531250D0, &
0.023668784648180007934570312500D0, 0.024900555610656738281250000000D0, &
0.026093553751707077026367187500D0, 0.027245631441473960876464843750D0, &
0.028355251997709274291992187500D0, 0.029420763254165649414062500000D0, &
0.030440155416727066040039062500D0, 0.031412057578563690185546875000D0, &
0.032334897667169570922851562500D0, 0.033207248896360397338867187500D0, &
0.034027762711048126220703125000D0, 0.034795157611370086669921875000D0, &
0.035508289933204650878906250000D0, 0.036165662109851837158203125000D0, &
0.036767188459634780883789062500D0, 0.037310685962438583374023437500D0, &
0.037796217948198318481445312500D0, 0.038222733885049819946289062500D0, &
0.038589447736740112304687500000D0, 0.038896244019269943237304687500D0, &
0.039141800254583358764648437500D0, 0.039326533675193786621093750000D0, &
0.039449851959943771362304687500D0, 0.039511464536190032958984375000D0, &
0.039511471986770629882812500000D0, 0.039449658244848251342773437500D0, &
0.039326600730419158935546875000D0, 0.039141751825809478759765625000D0, &
0.038896277546882629394531250000D0, 0.038589425384998321533203125000D0, &
0.038222678005695343017578125000D0, 0.037796292454004287719726562500D0, &
0.037310808897018432617187500000D0, 0.036766897886991500854492187500D0, &
0.036166030913591384887695312500D0, 0.035508446395397186279296875000D0, &
0.034795094281435012817382812500D0, 0.034027997404336929321289062500D0, &
0.033207222819328308105468750000D0, 0.032334662973880767822265625000D0, &
0.031412266194820404052734375000D0, 0.030440336093306541442871093750D0, &
0.029420448467135429382324218750D0, 0.028355462476611137390136718750D0, &
0.027245484292507171630859375000D0, 0.026093257591128349304199218750D0, &
0.024901069700717926025390625000D0, 0.023668944835662841796875000000D0, &
0.022399920970201492309570312500D0, 0.021095737814903259277343750000D0, &
0.019760118797421455383300781250D0, 0.018391948193311691284179687500D0, &
0.016995908692479133605957031250D0, 0.015573025681078433990478515625D0, &
0.014125997200608253479003906250D0, 0.012656771577894687652587890625D0, &
0.011167856864631175994873046875D0, 0.009661483578383922576904296875D0, &
0.008140278980135917663574218750D0, 0.006605920847505331039428710938D0, &
0.005061406176537275314331054688D0, 0.003508394351229071617126464844D0, &
0.001950456295162439346313476562D0, 0.000316502759233117103576660156/
END MODULE
!
! ************************************************************************
!
! This subroutine is to find the root such that the Equation (3) is
! satisfied.
!
! ************************************************************************
! GGEN(D, RM, RsM, Gm, Gsm, ALPHA, M, sss, DF)
! ======
SUBROUTINE GGEN(D, RM, RsM, Gm, Gsm, ALPHA, M, sss, DF)
IMPLICIT NONE
INTEGER M, sss
REAL D, ALPHA, DF, PR, EEPS
REAL VALX0, VALX1, DD, CV, XABS, GINT1
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
EXTERNAL GINT1
PR = 1.0 - ALPHA
EEPS = 0.001
VALX0 = GINT1(D, RM, RsM, Gm, Gsm, M, sss, DF) - PR
DD = D - 0.1
VALX1 = GINT1(DD, RM, RsM, Gm, Gsm, M, sss, DF) - PR
90 CV = DD - VALX1*(DD - D)/(VALX1 - VALX0)
VALX0 = VALX1
D = DD
VALX1 = GINT1(CV, RM, RsM, Gm, Gsm, M, sss, DF) - PR
DD = CV
XABS = ABS(DD - D)
IF (XABS .GT. EEPS) THEN
GOTO 90
END IF
D = DD
RETURN
END SUBROUTINE
!
! ************************************************************************
!
! Computing the probability of the left hand side of Equation (3)
!
! ************************************************************************
! GINT1(D1, RM, RsM, Gm, Gsm, NC, sss, DF)
! ======
! integerate chi square distribution w (that is X1)
! D1 is critical point (that is D)
! NC is k + 1 (that is M)
! DF
FUNCTION GINT1(D1, RM, RsM, Gm, Gsm, NC, sss, DF)
IMPLICIT NONE
INTEGER, PARAMETER :: NLEGQ = 16
INTEGER, PARAMETER :: IHALFQ = 8
REAL, PARAMETER :: EPS = 1.0
REAL, PARAMETER :: EPS1 = -30.0
REAL, PARAMETER :: EPS2 = 0.0000000000001
REAL, PARAMETER :: DHAF = 100.0
REAL, PARAMETER :: DQUAR = 800.0
REAL, PARAMETER :: DEIGH = 5000.0
REAL, PARAMETER :: ULEN1 = 1.0
REAL, PARAMETER :: ULEN2 = 0.5
REAL, PARAMETER :: ULEN3 = 0.25
REAL, PARAMETER :: ULEN4 = 0.125
INTEGER NC, i, j, jj, sss
REAL GINT1, GINT, QLGAMA, DF, D2(8), D1, F2, F2LF, F21, FF4
REAL ULEN, QTSUM, RQTSUM, TWA1, T1, XXX1, XX2
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
REAL*8 R2, XLEGQ(IHALFQ), ALEGQ(IHALFQ), qlgamma
EXTERNAL GINT, qlgamma
DATA R2 /0.693147180559945309417232121458D0/
DATA XLEGQ / 0.98940093499164993259615417345D0, 0.944575023073232576077988415535D0, &
0.865631202387831743880467897712D0, 0.755404408355003033895101194847D0, &
0.617876244402643748446671764049D0, 0.458016777657227386342419442984D0, &
0.281603550779258913230460501460D0, 0.095012509837637440185319335425D0 /
DATA ALEGQ /0.0271524594117540948517805724560D0, 0.0622535239386478928628438369944D0, &
0.0951585116824927848099251076022D0, 0.124628971255533872052476282192D0, &
0.149595988816576732081501730547D0, 0.169156519395002538189312079030D0, &
0.182603415044923588866763667969D0, 0.189450610455068496285396723208D0 /
IF (DF .GT. 2000.0) THEN
GINT1 = GINT(1.0, D1, RM, RsM, Gm, Gsm, NC, sss)
GOTO 900
END IF
F2 = DF*0.5
FF4 = DF*0.25
F2LF = F2*ALOG(DF) - DF*R2 - qlgamma(F2)
F21 = F2 - 1.0
IF (DF .LE. DHAF) THEN
ULEN = ULEN1
ELSE IF (Df .GT. DHAF .AND. DF .LE. DQUAR) THEN
ULEN = ULEN2
ELSE IF (DF .GT. DQUAR .AND. DF .LE. DEIGH) THEN
ULEN = ULEN3
ELSE
ULEN = ULEN4
END IF
F2LF = F2LF + ALOG(ULEN)
GINT1 = 0.0D0
DO i = 1, 50
QTSUM = 0.0D0
TWA1 = (2.0*REAL(i) - 1.0)*ULEN
DO jj = 1, NLEGQ
IF (IHALFQ .LT. jj) THEN
j = jj - IHALFQ
T1 = F2LF + F21*DLOG(TWA1 + XLEGQ(j)*ULEN) - (XLEGQ(j)*ULEN + TWA1)*FF4
ELSE
j = jj
T1 = F2LF + F21*DLOG(TWA1 - XLEGQ(j)*ULEN) + (XLEGQ(j)*ULEN - TWA1)*FF4
END IF
IF (T1 .GE. EPS1) THEN
IF (IHALFQ .LT. jj) THEN
XXX1 = SQRT((XLEGQ(j)*ULEN + TWA1)*0.5)
ELSE
XXX1 = SQRT((-XLEGQ(j)*ULEN + TWA1)*0.5)
END IF
XX2 = GINT(XXX1, D1, RM, RsM, Gm, Gsm, NC, sss)
RQTSUM = XX2*ALEGQ(j)*EXP(T1)
QTSUM = RQTSUM + QTSUM
END IF
END DO
IF ((REAL(i)*ULEN .GE. 1.0) .AND. (QTSUM .LE. EPS2)) GOTO 900
GINT1 = GINT1 + QTSUM
END DO
900RETURN
END FUNCTION
!
! GINT(X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
FUNCTION GINT(X1, D, RM, RsM, Gm, Gsm, M, sss)
IMPLICIT NONE
INTEGER sss, M, j
REAL GINT, X1, D, XX(80)
REAL GINTs1, GINTs2, GINTs3, GINTs4
REAL XX2, XX3, XX4
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
EXTERNAL GINTs1, GINTs2, GINTs3, GINTs4
IF (sss == 1) THEN
GINT = GINTs1(XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
ELSE IF (sss == 2) THEN
GINT = GINTs2(XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
ELSE IF (sss == 3) THEN
GINT = GINTs3(XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
ELSE
GINT = GINTs4(X1, D, RM, RsM, Gm, Gsm, M, sss)
END IF
RETURN
END FUNCTION
! GINTs4(X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
! integrate for normal y1, ... ys
FUNCTION GINTs4(X1, D, RM, RsM, Gm, Gsm, M, sss)
USE premain
IMPLICIT NONE
INTEGER i, M, sss, js
REAL GINTs3, SUM, C, A, X1, D, GF, GINTs4
REAL XX(80)
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
REAL XX4
EXTERNAL GINTs3
C = 8.0
A = -8.0
SUM = 0.0
DO i = 1, N
XX4 = ((C -A)*QX(i) + C + A)*0.5
SUM = SUM + GINTs3(XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)*QW(i)
END DO
GINTs4 = SUM*(C - A)*0.5
RETURN
END FUNCTION
! GINTs3(XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
! integrate for normal y1, ... ys
FUNCTION GINTs3(XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
USE premain
IMPLICIT NONE
INTEGER i, M, sss, js
REAL GINTs3, SUM, C, A, X1, D, GF, GINTs2
REAL XX3, XX4
REAL XX(80)
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
EXTERNAL GINTs2
C = 8.0
A = -8.0
SUM = 0.0
DO i = 1, N
XX3 = ((C -A)*QX(i) + C + A)*0.5
SUM = SUM + GINTs2(XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)*QW(i)
END DO
GINTs3 = SUM*(C - A)*0.5
RETURN
END FUNCTION
! GINTs2(XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
! integrate for normal y1, ... ys
FUNCTION GINTs2(XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
USE premain
IMPLICIT NONE
INTEGER i, M, sss, js
REAL GINTs1, SUM, C, A, X1, D, GF, GINTs2
REAL XX2, XX3, XX4
REAL XX(80)
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
EXTERNAL GINTs1
C = 8.0
A = -8.0
SUM = 0.0
DO i = 1, N
XX2 = ((C -A)*QX(i) + C + A)*0.5
SUM = SUM + GINTs1(XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)*QW(i)
END DO
GINTs2 = SUM*(C - A)*0.5
RETURN
END FUNCTION
! GINTs1(XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
! integrate for normal y1, ... ys
FUNCTION GINTs1(XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
USE premain
IMPLICIT NONE
INTEGER i, M, sss, js
REAL GINTs1, SUM, C, A, X1, D, GF
REAL XX1, XX2, XX3, XX4
REAL XX(80)
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
EXTERNAL GF
C = 8.0
A = -8.0
SUM = 0.0
DO i = 1, N
XX1 = ((C -A)*QX(i) + C + A)*0.5
SUM = SUM + GF(XX1, XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)*QW(i)
END DO
GINTs1 = SUM*(C - A)*0.5
RETURN
END FUNCTION
! GF(XX1, XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
! ======
! prod_i=1^k cdfNormal()*[1 - cdfNormal()]*pdf(y1)*...*pdf(ys)
! XX for normal that is y1, ..., ys
! X1 for chi sqruare that is w
! D for critical point that is c
! nnk for (nE1, nE2, ..., nEk)
! nns for (nR1, nR2, ..., nRs)
! np for nP
! M is k + 1
! sss is the dim for XX, that is the s
! DTT is w*c
FUNCTION GF(XX1, XX2, XX3, XX4, X1, D, RM, RsM, Gm, Gsm, M, sss)
IMPLICIT NONE
REAL GF, normalcdf, normalpdf
REAL XX1, XX2, XX3, XX4, XX(80)
REAL X1, D, DTT, F1, A1, A2, C1, C2
REAL C0rho, C0rhos, C0rho1, c0rhos1, C1rho, C1rhos, C1rho1, C1rhos1
REAL RM(80, 80), RsM(80, 80), Gm(80), Gsm(80)
INTEGER i, M, j, sss
EXTERNAL normalcdf, normalpdf
XX(1) = XX1
XX(2) = XX2
XX(3) = XX3
XX(4) = XX4
DTT = D*X1
F1 = 1.0
DO j = 1, (M - 1) ! k
C0rho = RM(j, 1)
C0rhos = RsM(j, 1)
C1 = (DTT + XX(1)*C0rho)/C0rhos
IF (sss == 1) THEN
GOTO 1001
ELSE
DO i = 2, sss ! s
C0rho1 = RM(j, i)
C0rhos1 = RsM(j, i)
C1 = min(C1, (DTT + XX(i)*C0rho1)/C0rhos1)
END DO
END IF
1001 A1 = normalcdf(C1)
F1 = F1*A1
END DO
C1rho = Gm(1)
C1rhos = Gsm(1)
C2 = (XX(1)*C1rhos - DTT)/C1rho
IF (sss == 1) THEN
GOTO 1002
ELSE
DO i = 2, sss
C1rho1 = Gm(i)
C1rhos1 = Gsm(i)
C2 = max(C2, (XX(i)*C1rhos1 - DTT)/C1rho1)
END DO
END IF
1002 A2 = normalcdf(-C2)
GF = F1*A2
DO i = 1, sss
GF = GF*normalpdf(XX(i))
END DO
GF = GF
RETURN
END FUNCTION
! nature log gamma function
! qlgamma(x)
! ======
FUNCTION qlgamma(x)
IMPLICIT NONE
REAL x
REAL*8 xx, a, alog, frac, y
REAL*8 qlgamma
REAL*8 sqrtpi, logsqrt2pi, d012, d030,d035, d013,d071
PARAMETER (xx = 21.0)
DATA sqrtpi /1.772453850905516027298167D0/
DATA logsqrt2pi /0.9189385332046727417803297D0/
DATA d012 /0.833333333333333333333333333333D-1/
DATA d030 /30.0D0/
DATA d035 /3.5D0/
DATA d013 /1.333333333333333333333333333333D0/
DATA d071 /0.707142857142857142857142857143D0/
frac = x - int(x)
IF ((x .LT. xx) .AND. (frac .EQ. 0.0)) THEN
qlgamma = 1.0
y = x - 1.0
DO a = 1.0, y
qlgamma = qlgamma*a
END DO
qlgamma = dlog(qlgamma)
ELSE IF ((x .LT. xx) .AND. (frac .EQ. 0.5)) THEN
qlgamma = 1.0
y = 2.0*x - 1.0
DO a = 1.0, y, 2.0
qlgamma = qlgamma*a*0.5
END DO
qlgamma = dlog(qlgamma*sqrtpi)
ELSE
y = 1.0/(x*x)
qlgamma = d012*(1.0 - y/d030*(1.0 - y/d035*(1.0 - y/d013*(1.0 - y/d071))))/x
qlgamma = qlgamma + logsqrt2pi + (x - 0.5)*alog(x) - x
END IF
RETURN
END FUNCTION
! cdf of normal distribtuion form -8 to x
! ======
FUNCTION normalcdf(x)
USE premain
IMPLICIT NONE
INTEGER i
REAL normalcdf, x
normalcdf = 0.0
DO i = 1, N
normalcdf = normalcdf + QW(i)*0.5*(x + 8)*exp(-0.125*((x + 8)*QX(i) + x - 8)**2)/SQRT2PI
END DO
normalcdf = normalcdf
RETURN
END FUNCTION
!! pdf of normal distribtuion
!! ======
FUNCTION normalpdf(x)
USE premain
IMPLICIT NONE
REAL normalpdf, x
normalpdf = exp((-1)*0.5*x*x)/SQRT2PI
RETURN
END FUNCTION