You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1100 lines
37 KiB
Fortran
1100 lines
37 KiB
Fortran
C
|
|
C f2py -m mvnprd -h mvnprd.pyf mvnprd.f only: mvnprd
|
|
C edit mvnprd.pyf with input and output and then
|
|
C f2py mvnprd.pyf mvnprd.f -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
|
|
C
|
|
C f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m mvnprd -c mvnprd.f
|
|
|
|
C Altarnative: compile mvnprd and link to it through mvnprd_interface.f
|
|
C
|
|
C gfortran -fPIC -c mvnprd.f
|
|
C f2py -m mvnprdmod -c mvnprd.o mvnprd_interface.f --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
|
|
C
|
|
C df -c mvnprd.f
|
|
C f2py -m mvnprdmod -c mvnprd.obj mvnprd_interface.f --fcompiler=compaqv --compiler=mingw32 -lmsvcr71
|
|
|
|
! This is a MEX-file for MATLAB.
|
|
! and contains a mex-interface to Charles W. Dunnett's programs
|
|
! ,MVNPRD and MVSTUD subroutines for computing multivariate normal
|
|
! or student T probabilities with product correlation structure. The
|
|
! file should compile without errors on (Fortran77) standard Fortran
|
|
! compilers.
|
|
*
|
|
* The mex-interface was written by
|
|
* Per Andreas Brodtkorb
|
|
* Norwegian Defence Research Establishment
|
|
* P.O. Box 115m
|
|
* N-3191 Horten
|
|
* Norway
|
|
* Email: Per.Brodtkorb@ffi.no
|
|
*
|
|
* Charles Dunnett
|
|
C Dept. of Mathematics and Statistics
|
|
C McMaster University
|
|
C Hamilton, Ontario L8S 4K1
|
|
C Canada
|
|
C E-mail: dunnett@mcmaster.ca
|
|
C Tel.: (905) 525-9140 (Ext. 27104)
|
|
*
|
|
* MVNPRDMEX Computes multivariate normal or student T probability
|
|
* with product correlation structure.
|
|
|
|
*
|
|
* CALL [value,bound,inform] = mvnprdmex(RHO,A,B,D,NDF,abseps,IERC,HNC)
|
|
*
|
|
* RHO REAL, array of coefficients defining the correlation
|
|
* coefficient by:
|
|
* correlation(I,J) = RHO(I)*RHO(J) for J/=I
|
|
* where
|
|
* 1 < RHO(I) < 1
|
|
* A REAL, array of lower integration limits.
|
|
* B REAL, array of upper integration limits.
|
|
* NOTE: any values greater the 37, are considered as
|
|
* infinite values.
|
|
* D Real array of means
|
|
* NDF Degrees of freedom, NDF<=0 gives normal probabilities
|
|
* ABSEPS REAL absolute error tolerance.
|
|
* IERC INTEGER 1 if strict error control based on fourth
|
|
* derivative
|
|
* 0 if intuitive error control based on halving the
|
|
* intervals
|
|
* HINC REAL start interval width of simpson rule
|
|
*
|
|
* OUTPUT:
|
|
* VALUE REAL estimated value for the integral
|
|
* BOUND REAL bound on the error of the approximation
|
|
* INFORM INTEGER, termination status parameter:
|
|
* 0, if normal completion with ERROR < EPS;
|
|
* 1, if N > 100 or N < 1.
|
|
* 2, IF any abs(rho)>=1
|
|
* 4, if ANY(B(I)<=A(i))
|
|
* 5, if number of terms computed exceeds maximum number of
|
|
* evaluation points
|
|
* 6, if fault accurs in normal subroutines
|
|
* 7, if subintervals are too narrow or too many
|
|
* 8, if bounds exceeds abseps
|
|
*
|
|
*
|
|
* MVNPRDMEX calculates multivariate normal or student T probability
|
|
* with product correlation structure for rectangular regions.
|
|
* The accuracy is up to around single precision, i.e., about 1e-7.
|
|
*
|
|
* This file was successfully compiled for matlab 5.3
|
|
* using Compaq Visual Fortran 6.1, and Windows 2000.
|
|
* The example here uses Fortran77 source.
|
|
* First, you will need to modify your mexopts.bat file.
|
|
* To find it, issue the command prefdir(1) from the Matlab command line,
|
|
* the directory it answers with will contain your mexopts.bat file.
|
|
* Open it for editing. The first section will look like:
|
|
*
|
|
*rem ********************************************************************
|
|
*rem General parameters
|
|
*rem ********************************************************************
|
|
*set MATLAB=%MATLAB%
|
|
*set DF_ROOT=C:\Program Files\Microsoft Visual Studio
|
|
*set VCDir=%DF_ROOT%\VC98
|
|
*set MSDevDir=%DF_ROOT%\Common\msdev98
|
|
*set DFDir=%DF_ROOT%\DF98
|
|
*set PATH=%MSDevDir%\bin;%DFDir%\BIN;%VCDir%\BIN;%PATH%
|
|
*set INCLUDE=%DFDir%\INCLUDE;%DFDir%\IMSL\INCLUDE;%INCLUDE%
|
|
*set LIB=%DFDir%\LIB;%VCDir%\LIB
|
|
*
|
|
* then you are ready to compile this file at the matlab prompt using the
|
|
* following command:
|
|
* mex -O mvnprdmex.f
|
|
|
|
|
|
|
|
C The rest of this file contains:
|
|
C 1. A Readme file provided by Charles Dunnett, the author of AS 251.
|
|
C 2. The published algorithm AS 251 together with the two other AS algorithms
|
|
C which it calls (AS 66 and AS 241).
|
|
C 3. A driver program (MVTIN) for either multivariate normal or t.
|
|
C 4. MVSTUD for calculating multivariate t probabilities.
|
|
C ***************************************************************************
|
|
|
|
C Date: Mon, 10 Apr 1995 16:49:10 +0059 (EDT)
|
|
C From: "Charles W. Dunnett" <dunnett@mcmail.cis.mcmaster.ca>
|
|
C Subject: Readme for AS251 (extended version incl. multivariate t)
|
|
|
|
C MVTIN is a driver program for computing multivariate normal or t
|
|
C probability integrals over arbitrary rectangular regions. The
|
|
C correlation structure is assumed to be of product form, rho_ij =
|
|
C b_i x b_j, where -1 < b_i < +1.
|
|
|
|
C It requires the following:-
|
|
|
|
C 1. MVNPRD (published as algorithm AS 251 in Applied
|
|
C Statistics (1989), 38: 564-579; see also the correction
|
|
C note in Applied Statistics (1993), 42: 709),
|
|
|
|
C 2. ALNORM and PPND7 (published as algorithms AS 66 and
|
|
C AS 241, respectively, in Applied Statistics, and
|
|
|
|
C 3. MVSTUD ( which Studentizes MVNPRD).
|
|
|
|
|
|
|
|
SUBROUTINE MVNPRD(A, B, BPD, EPS, N, INF, IERC, HINC, PROB, BOUND,
|
|
* IFAULT)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 251.1 APPL.STATIST. (1989), VOL.38, NO.3
|
|
C
|
|
C FOR A MULTIVARIATE NORMAL VECTOR WITH CORRELATION STRUCTURE
|
|
C DEFINED BY RHO(I,J) = BPD(I) * BPD(J), COMPUTES THE PROBABILITY
|
|
C THAT THE VECTOR FALLS IN A RECTANGLE IN N-SPACE WITH ERROR
|
|
C LESS THAN EPS.
|
|
C
|
|
INTEGER NN
|
|
PARAMETER (NN = 100)
|
|
DOUBLE PRECISION A(*), B(*), BPD(*), ESTT(22), FV(5), FD(5),
|
|
& F1T(22), F2T(22), F3T(22), G1T(22), G3T(22), PSUM(22), H(NN)
|
|
$ , HL(NN),BB(NN)
|
|
INTEGER INF(*), INFT(NN), LDIR(22)
|
|
DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, SIX, PT1, PT24,
|
|
* SMALL, DXMIN, SQRT2, PROB, ERRL, BI, START,
|
|
* Z, HINC, ADDN, EPS2, EPS1, EPS, ZU, Z2, Z3, Z4, Z5, ZZ,
|
|
* ERFAC, EL, EL1, BOUND, PART0, PART2, PART3, FUNC0, FUNC2,
|
|
* FUNCN, WT, CONTRB, DLG, DX, DA, ESTL, ESTR, TSUM, EXCESS, ERROR,
|
|
* PROB1, SAFE, ONEP5,X2880
|
|
INTEGER N, IERC, IFAULT, I, NTM, NMAX, LVL, NR, NDIM
|
|
DOUBLE PRECISION ALNORM, PPND7
|
|
EXTERNAL ALNORM, PPND7
|
|
DATA ZERO, HALF, ONE, TWO, FOUR, SIX /0.0, 0.5, 1.0, 2.0,
|
|
* 4.0, 6.0/
|
|
DATA PT1, PT24, ONEP5, X2880 /0.1, 0.24, 1.5, 2880.0/
|
|
DATA SMALL, DXMIN, SQRT2 /1.0E-10, 0.0000001, 1.41421356237310/
|
|
C
|
|
C CHECK FOR INPUT VALUES OUT OF RANGE.
|
|
C
|
|
PROB = ZERO
|
|
BOUND = ZERO
|
|
IFAULT = 1
|
|
IF (N .LT. 1 .OR. N .GT. NN) RETURN
|
|
DO 10 I = 1, N
|
|
BI = ABS(BPD(I))
|
|
IFAULT = 2
|
|
IF (BI .GE. ONE) RETURN
|
|
IFAULT = 3
|
|
IF (INF(I) .LT. 0 .OR. INF(I) .GT. 2) RETURN
|
|
IFAULT = 4
|
|
IF (INF(I) .EQ. 2 .AND. A(I) .LE. B(I)) RETURN
|
|
10 CONTINUE
|
|
IFAULT = 0
|
|
PROB = ONE
|
|
C
|
|
C CHECK WHETHER ANY BPD(I) = 0.
|
|
C
|
|
NDIM = 0
|
|
DO 20 I = 1, N
|
|
IF (BPD(I) .NE. ZERO) THEN
|
|
NDIM = NDIM + 1
|
|
H(NDIM) = A(I)
|
|
HL(NDIM) = B(I)
|
|
BB(NDIM) = BPD(I)
|
|
INFT(NDIM) = INF(I)
|
|
ELSE
|
|
C
|
|
C IF ANY BPD(I) = 0, THE CONTRIBUTION TO PROB FOR THAT
|
|
C VARIABLE IS COMPUTED FROM A UNIVARIATE NORMAL.
|
|
C
|
|
IF (INF(I) .LT. 1) THEN
|
|
PROB = PROB * (ONE - ALNORM(B(I), .FALSE.))
|
|
ELSE IF (INF(I) .EQ. 1) THEN
|
|
PROB = PROB * ALNORM(A(I), .FALSE.)
|
|
ELSE
|
|
PROB = PROB * (ALNORM(A(I), .FALSE.) -
|
|
* ALNORM(B(I), .FALSE.))
|
|
END IF
|
|
IF (PROB .LE. SMALL) PROB = ZERO
|
|
END IF
|
|
20 CONTINUE
|
|
IF (NDIM .EQ. 0 .OR. PROB .EQ. ZERO) RETURN
|
|
C
|
|
C IF NOT ALL BPD(I) = 0, PROB IS COMPUTED BY SIMPSON'S RULE.
|
|
C BUT FIRST, INITIALIZE THE VARIABLES.
|
|
C
|
|
Z = ZERO
|
|
IF (HINC .LE. ZERO) HINC = PT24
|
|
ADDN = -ONE
|
|
DO 30 I = 1, NDIM
|
|
IF (INFT(I) .EQ. 2 .OR.
|
|
* (INFT(I) .NE. INFT(1) .AND. BB(I) * BB(1) .GT. ZERO) .OR.
|
|
* (INFT(I) .EQ. INFT(1) .AND. BB(I) * BB(1) .LT. ZERO))
|
|
* ADDN = ZERO
|
|
30 CONTINUE
|
|
C
|
|
C THE VALUE OF ADDN IS TO BE ADDED TO THE PRODUCT EXPRESSIONS IN
|
|
C THE INTEGRAND TO INSURE THAT THE LIMITING VALUE IS ZERO.
|
|
C
|
|
PROB1 = ZERO
|
|
NTM = 0
|
|
NMAX = 400
|
|
IF (IERC .EQ. 0) NMAX = NMAX * 2
|
|
CALL PFUNC (Z, H, HL, BB, NDIM, INFT, ADDN, SAFE, FUNC0, NTM,
|
|
* IERC, PART0)
|
|
EPS2 = EPS * PT1 * HALF
|
|
C
|
|
C SET UPPER BOUND ON Z AND APPORTION EPS.
|
|
C
|
|
ZU = -PPND7(EPS2, IFAULT) / SQRT2
|
|
IF (IFAULT .NE. 0) THEN
|
|
IFAULT = 6
|
|
RETURN
|
|
END IF
|
|
!NR = IFIX(ZU / HINC) + 1
|
|
NR = NINT(ZU / HINC) + 1
|
|
ERFAC = ONE
|
|
IF (IERC .NE. 0) ERFAC = X2880 / HINC ** 5
|
|
EL = (EPS - EPS2) / FLOAT(NR) * ERFAC
|
|
EL1 = EL
|
|
C
|
|
C START COMPUTATIONS FOR THE INTERVAL (Z, Z + HINC).
|
|
C
|
|
40 ERROR = ZERO
|
|
LVL = 0
|
|
FV(1) = PART0
|
|
FD(1) = SAFE
|
|
START = Z
|
|
DA = HINC
|
|
Z3 = START + HALF * DA
|
|
CALL PFUNC(Z3, H, HL, BB, NDIM, INFT, ADDN, FD(3), FUNCN, NTM,
|
|
* IERC, FV(3))
|
|
Z5 = START + DA
|
|
CALL PFUNC(Z5, H, HL, BB, NDIM, INFT, ADDN, FD(5), FUNC2, NTM,
|
|
* IERC, FV(5))
|
|
PART2 = FV(5)
|
|
SAFE = FD(5)
|
|
WT = DA / SIX
|
|
CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5))
|
|
DLG = ZERO
|
|
IF (IERC .NE. 0) THEN
|
|
CALL WMAX(FD(1), FD(3), FD(5), DLG)
|
|
IF (DLG .LE. EL) GO TO 90
|
|
DX = DA
|
|
GO TO 60
|
|
END IF
|
|
LVL = 1
|
|
LDIR(LVL) = 2
|
|
PSUM(LVL) = ZERO
|
|
C
|
|
C BISECT INTERVAL. IF IERC = 1, COMPUTE ESTIMATE ON LEFT
|
|
C HALF; IF IERC = 0, ON BOTH HALVES.
|
|
C
|
|
50 DX = HALF * DA
|
|
WT = DX / SIX
|
|
Z2 = START + HALF * DX
|
|
CALL PFUNC(Z2, H, HL, BB, NDIM, INFT, ADDN, FD(2), FUNCN, NTM,
|
|
* IERC,FV(2))
|
|
ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3))
|
|
IF (IERC .EQ. 0) THEN
|
|
Z4 = START + ONEP5 * DX
|
|
CALL PFUNC(Z4, H, HL, BB, NDIM, INFT, ADDN, FD(4), FUNCN,
|
|
* NTM, IERC, FV(4))
|
|
ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5))
|
|
TSUM = ESTL + ESTR
|
|
DLG = ABS(CONTRB - TSUM)
|
|
EPS1 = EL / TWO ** (LVL - 1)
|
|
ERRL = DLG
|
|
ELSE
|
|
FV(3) = FV(2)
|
|
FD(3) = FD(2)
|
|
CALL WMAX(FD(1), FD(3), FD(5), DLG)
|
|
ERRL = DLG / TWO ** (5 * LVL)
|
|
TSUM = ESTL
|
|
EPS1 = EL * (TWO ** LVL) ** 4
|
|
END IF
|
|
C
|
|
C STOP SUBDIVIDING INTERVAL WHEN ACCURACY IS SUFFICIENT,
|
|
C OR IF INTERVAL TOO NARROW OR SUBDIVIDED TOO OFTEN.
|
|
C
|
|
IF (DLG .LE. EPS1 .OR. DLG .LT. SMALL) GO TO 70
|
|
IF (IFAULT .EQ. 0 .AND. NTM .GE. NMAX) IFAULT = 5
|
|
IF (ABS(DX) .LE. DXMIN .OR. LVL .GT. 21) IFAULT = 7
|
|
IF (IFAULT .NE. 0) GO TO 70
|
|
C
|
|
C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF AND APPLY
|
|
C SIMPSON'S RULE TO LEFT HALF.
|
|
C
|
|
60 LVL = LVL + 1
|
|
LDIR(LVL) = 1
|
|
F1T(LVL) = FV(3)
|
|
F3T(LVL) = FV(5)
|
|
DA = DX
|
|
FV(5) = FV(3)
|
|
IF (IERC .EQ. 0) THEN
|
|
F2T(LVL) = FV(4)
|
|
ESTT(LVL) = ESTR
|
|
CONTRB = ESTL
|
|
FV(3) = FV(2)
|
|
ELSE
|
|
G1T(LVL) = FD(3)
|
|
G3T(LVL) = FD(5)
|
|
FD(5) = FD(3)
|
|
END IF
|
|
GO TO 50
|
|
C
|
|
C ACCEPT APPROXIMATE VALUE FOR INTERVAL.
|
|
C RESTORE SAVED INFORMATION TO PROCESS
|
|
C RIGHT HALF INTERVAL.
|
|
C
|
|
70 ERROR = ERROR + ERRL
|
|
80 IF (LDIR(LVL) .EQ. 1) THEN
|
|
PSUM(LVL) = TSUM
|
|
LDIR(LVL) = 2
|
|
IF (IERC .EQ. 0) DX = DX * TWO
|
|
START = START + DX
|
|
DA = HINC / TWO ** (LVL - 1)
|
|
FV(1) = F1T(LVL)
|
|
IF (IERC .EQ. 0) THEN
|
|
FV(3) = F2T(LVL)
|
|
CONTRB = ESTT(LVL)
|
|
ELSE
|
|
FV(3) = F3T(LVL)
|
|
FD(1) = G1T(LVL)
|
|
FD(5) = G3T(LVL)
|
|
END IF
|
|
FV(5) = F3T(LVL)
|
|
GO TO 50
|
|
END IF
|
|
TSUM = TSUM + PSUM(LVL)
|
|
LVL = LVL - 1
|
|
IF (LVL .GT. 0) GO TO 80
|
|
CONTRB = TSUM
|
|
LVL = 1
|
|
DLG = ERROR
|
|
90 PROB1 = PROB1 + CONTRB
|
|
BOUND = BOUND + DLG
|
|
EXCESS = EL - DLG
|
|
EL = EL1
|
|
IF (EXCESS .GT. ZERO) EL = EL1 + EXCESS
|
|
IF ((FUNC0 .GT. ZERO .AND. FUNC2 .LE. FUNC0) .OR.
|
|
* (FUNC0 .LT. ZERO .AND. FUNC2 .GE. FUNC0)) THEN
|
|
ZZ = -SQRT2 * Z5
|
|
PART3 = ABS(FUNC2) * ALNORM(ZZ, .FALSE.) + BOUND / ERFAC
|
|
IF (PART3 .LE. EPS .OR. NTM .GE. NMAX .OR. Z5 .GE. ZU) GOTO 100
|
|
END IF
|
|
Z = Z5
|
|
PART0 = PART2
|
|
FUNC0 = FUNC2
|
|
IF (Z .LT. ZU .AND. NTM .LT. NMAX) GO TO 40
|
|
100 PROB = (PROB1 - ADDN * HALF) * PROB
|
|
BOUND = PART3
|
|
IF (NTM .GE. NMAX .AND. IFAULT .EQ. 0) IFAULT = 5
|
|
IF (BOUND .GT. EPS .AND. IFAULT .EQ. 0) IFAULT = 8
|
|
RETURN
|
|
END
|
|
SUBROUTINE PFUNC(Z, A, B, BPD, N, INF, ADDN, DERIV, FUNCN, NTM,
|
|
* IERC, RESULT)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 251.2 APPL.STATIST. (1989), VOL.38, NO.3
|
|
C
|
|
C
|
|
C COMPUTE FUNCTION IN INTEGRAND AND ITS 4TH DERIVATIVE.
|
|
C
|
|
INTEGER NN
|
|
PARAMETER (NN = 100)
|
|
DOUBLE PRECISION A(*), B(*), BPD(*), FOU(NN), FOU1(4, NN), TMP(4),
|
|
& GOU(NN), GOU1(4, NN), FF(4), GF(4), TERM(4), GERM(4)
|
|
INTEGER INF(*)
|
|
DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE,
|
|
& SIXTN, SMALL, Z, U, U1, U2, BI, HI, HLI, BP, ADDN, DERIV,
|
|
$ FUNCN,RESULT, RSLT1, RSLT2, DEN, SQRT2, SQRTPI, PHI, PHI1,
|
|
$ PHI2,PHI3, PHI4, FRM, GRM
|
|
INTEGER N, NTM, IERC, INFI, I, J, K, M, L, IK
|
|
DOUBLE PRECISION ALNORM
|
|
EXTERNAL ALNORM
|
|
DATA ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE, SIXTN,
|
|
* SMALL /0.0, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 12.0, 16.0, 0.1E-12/
|
|
DATA SQRT2, SQRTPI /1.41421356237310, 1.77245385090552/
|
|
DERIV = ZERO
|
|
NTM = NTM + 1
|
|
RSLT1 = ONE
|
|
RSLT2 = ONE
|
|
BI = ONE
|
|
HI = A(1) + ONE
|
|
HLI = B(1) + ONE
|
|
INFI = -1
|
|
DO 60 I = 1, N
|
|
IF (BPD(I) .EQ. BI .AND. A(I) .EQ. HI .AND. B(I) .EQ. HLI .AND.
|
|
* INF(I) .EQ. INFI) THEN
|
|
FOU(I) = FOU(I - 1)
|
|
GOU(I) = GOU(I - 1)
|
|
DO 10 IK = 1, 4
|
|
FOU1(IK, I) = FOU1(IK, I - 1)
|
|
GOU1(IK, I) = GOU1(IK, I - 1)
|
|
10 CONTINUE
|
|
ELSE
|
|
BI = BPD(I)
|
|
HI = A(I)
|
|
HLI = B(I)
|
|
INFI = INF(I)
|
|
IF (BI .EQ. ZERO) THEN
|
|
IF (INFI .LT. 1) THEN
|
|
FOU(I) = ONE - ALNORM(HLI, .FALSE.)
|
|
ELSE IF (INFI .EQ. 1) THEN
|
|
FOU(I) = ALNORM(HI, .FALSE.)
|
|
ELSE
|
|
FOU(I) = ALNORM(HI, .FALSE.) - ALNORM(HLI, .FALSE.)
|
|
END IF
|
|
GOU(I) = FOU(I)
|
|
DO 20 IK = 1, 4
|
|
FOU1(IK, I) = ZERO
|
|
GOU1(IK, I) = ZERO
|
|
20 CONTINUE
|
|
ELSE
|
|
DEN = SQRT(ONE - BI * BI)
|
|
BP = BI * SQRT2 / DEN
|
|
IF (INFI .LT. 1) THEN
|
|
U = -HLI / DEN + Z * BP
|
|
FOU(I) = ALNORM(U, .FALSE.)
|
|
CALL ASSIGN (U, BP, FOU1(1, I))
|
|
BP = -BP
|
|
U = -HLI / DEN + Z * BP
|
|
GOU(I) = ALNORM(U, .FALSE.)
|
|
CALL ASSIGN (U, BP, GOU1(1, I))
|
|
ELSE IF (INFI .EQ. 1) THEN
|
|
U = HI / DEN + Z * BP
|
|
GOU(I) = ALNORM(U, .FALSE.)
|
|
CALL ASSIGN (U, BP, GOU1(1, I))
|
|
BP = -BP
|
|
U = HI / DEN + Z * BP
|
|
FOU(I) = ALNORM(U, .FALSE.)
|
|
CALL ASSIGN (U, BP, FOU1(1, I))
|
|
ELSE
|
|
U2 = -HLI / DEN + Z * BP
|
|
CALL ASSIGN (U2, BP, FOU1(1, I))
|
|
BP = -BP
|
|
U1 = HI / DEN + Z * BP
|
|
CALL ASSIGN (U1, BP, TMP(1))
|
|
FOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.) - ONE
|
|
DO 30 IK = 1, 4
|
|
FOU1(IK, I) = FOU1(IK, I) + TMP(IK)
|
|
30 CONTINUE
|
|
IF (-HLI .EQ. HI) THEN
|
|
GOU(I) = FOU(I)
|
|
DO 40 IK = 1, 4
|
|
GOU1(IK, I) = FOU1(IK, I)
|
|
40 CONTINUE
|
|
ELSE
|
|
U2 = -HLI / DEN + Z * BP
|
|
CALL ASSIGN (U2, BP, GOU1(1, I))
|
|
BP = -BP
|
|
U1 = HI / DEN + Z * BP
|
|
GOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.)-ONE
|
|
CALL ASSIGN (U1, BP, TMP(1))
|
|
DO 50 IK = 1, 4
|
|
GOU1(IK, I) = GOU1(IK, I) + TMP(IK)
|
|
50 CONTINUE
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
RSLT1 = RSLT1 * FOU(I)
|
|
RSLT2 = RSLT2 * GOU(I)
|
|
IF (RSLT1 .LE. SMALL) RSLT1 = ZERO
|
|
IF (RSLT2 .LE. SMALL) RSLT2 = ZERO
|
|
60 CONTINUE
|
|
FUNCN = RSLT1 + RSLT2 + ADDN
|
|
RESULT = FUNCN * EXP(-Z * Z) / SQRTPI
|
|
C
|
|
C IF 4TH DERIVATIVE IS NOT WANTED, STOP HERE.
|
|
C OTHERWISE, PROCEED TO COMPUTE 4TH DERIVATIVE.
|
|
C
|
|
IF (IERC .EQ. 0) RETURN
|
|
DO 70 IK = 1, 4
|
|
FF(IK) = ZERO
|
|
GF(IK) = ZERO
|
|
70 CONTINUE
|
|
DO 100 I = 1, N
|
|
FRM = ONE
|
|
GRM = ONE
|
|
DO 80 J = 1, N
|
|
IF (J .EQ. 1) GO TO 80
|
|
FRM = FRM * FOU(J)
|
|
GRM = GRM * GOU(J)
|
|
IF (FRM .LE. SMALL) FRM = ZERO
|
|
IF (GRM .LE. SMALL) GRM = ZERO
|
|
80 CONTINUE
|
|
DO 90 IK = 1, 4
|
|
FF(IK) = FF(IK) + FRM * FOU1(IK, I)
|
|
GF(IK) = GF(IK) + GRM * GOU1(IK, I)
|
|
90 CONTINUE
|
|
100 CONTINUE
|
|
IF (N .LE. 2) GO TO 230
|
|
DO 130 I = 1, N
|
|
DO 120 J = I + 1, N
|
|
TERM(2) = FOU1(1, I) * FOU1(1, J)
|
|
GERM(2) = GOU1(1, I) * GOU1(1, J)
|
|
TERM(3) = FOU1(2, I) * FOU1(1, J)
|
|
GERM(3) = GOU1(2, I) * GOU1(1, J)
|
|
TERM(4) = FOU1(3, I) * FOU1(1, J)
|
|
GERM(4) = GOU1(3, I) * GOU1(1, J)
|
|
TERM(1) = FOU1(2, I) * FOU1(2, J)
|
|
GERM(1) = GOU1(2, I) * GOU1(2, J)
|
|
DO 110 K = 1, N
|
|
IF (K .EQ. I .OR. K .EQ. J) GO TO 110
|
|
CALL TOOSML (1, TERM, FOU(K))
|
|
CALL TOOSML (1, GERM, GOU(K))
|
|
110 CONTINUE
|
|
FF(2) = FF(2) + TWO * TERM(2)
|
|
FF(3) = FF(3) + TWO * TERM(3) * THREE
|
|
FF(4) = FF(4) + TWO * (TERM(4) * FOUR + TERM(1) * THREE)
|
|
GF(2) = GF(2) + TWO * GERM(2)
|
|
GF(3) = GF(3) + TWO * GERM(3) * THREE
|
|
GF(4) = GF(4) + TWO * (GERM(4) * FOUR + GERM(1) * THREE)
|
|
120 CONTINUE
|
|
130 CONTINUE
|
|
DO 170 I = 1, N
|
|
DO 160 J = I + 1, N
|
|
DO 150 K = J + 1, N
|
|
TERM(3) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K)
|
|
TERM(4) = FOU1(2, I) * FOU1(1, J) * FOU1(1, K)
|
|
GERM(3) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K)
|
|
GERM(4) = GOU1(2, I) * GOU1(1, J) * GOU1(1, K)
|
|
IF (N .GT. 3) THEN
|
|
DO 140 M = 1, N
|
|
IF (M .EQ. I .OR. M .EQ. J .OR. M .EQ. K) GO TO 140
|
|
CALL TOOSML (3, TERM, FOU(M))
|
|
CALL TOOSML (3, GERM, GOU(M))
|
|
140 CONTINUE
|
|
END IF
|
|
FF(3) = FF(3) + SIX * TERM(3)
|
|
FF(4) = FF(4) + SIX * TERM(4) * SIX
|
|
GF(3) = GF(3) + SIX * GERM(3)
|
|
GF(4) = GF(4) + SIX * GERM(4) * SIX
|
|
150 CONTINUE
|
|
160 CONTINUE
|
|
170 CONTINUE
|
|
IF (N .LE. 3) GO TO 230
|
|
DO 220 I = 1, N
|
|
DO 210 J = I + 1, N
|
|
DO 200 K = J + 1, N
|
|
DO 190 M = K + 1, N
|
|
TERM(4) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K) * FOU1(1, M)
|
|
GERM(4) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K) * GOU1(1, M)
|
|
IF (N .GT. 4) THEN
|
|
DO 180 L = 1, N
|
|
IF (L .EQ. I .OR. L .EQ. J .OR. L .EQ. K .OR. L .EQ. M)GOTO 180
|
|
CALL TOOSML (4, TERM, FOU(L))
|
|
CALL TOOSML (4, GERM, GOU(L))
|
|
180 CONTINUE
|
|
END IF
|
|
FF(4) = FF(4) + FOUR * SIX * TERM(4)
|
|
GF(4) = GF(4) + FOUR * SIX * GERM(4)
|
|
190 CONTINUE
|
|
200 CONTINUE
|
|
210 CONTINUE
|
|
220 CONTINUE
|
|
C
|
|
230 CONTINUE
|
|
PHI = EXP(-Z * Z) / SQRTPI
|
|
PHI1 = -TWO * Z * PHI
|
|
PHI2 = (FOUR * Z ** 2 - TWO) * PHI
|
|
PHI3 = (-EIGHT * Z ** 3 + TWELVE * Z) * PHI
|
|
PHI4 = (SIXTN * Z ** 2 * (Z ** 2 - THREE) + TWELVE) * PHI
|
|
DERIV = PHI * (FF(4) + GF(4)) + FOUR * PHI1 * (FF(3) + GF(3))
|
|
* + SIX * PHI2 * (FF(2) + GF(2)) + FOUR * PHI3 * (FF(1) + GF(1))
|
|
* + PHI4 * FUNCN
|
|
RETURN
|
|
END
|
|
SUBROUTINE ASSIGN (U, BP, FF)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 251.3 APPL.STATIST. (1989), VOL.38, NO.3
|
|
C
|
|
C
|
|
C COMPUTE DERIVATIVES OF NORMAL CDF'S.
|
|
C
|
|
DOUBLE PRECISION FF(4)
|
|
DOUBLE PRECISION U, U2, BP, HALF, ONE, THREE, SQ2PI, T1, T2, T3
|
|
$ ,ZERO, UMAX, SMALL
|
|
INTEGER I
|
|
DATA HALF, ONE, THREE, SQ2PI /0.5, 1.0, 3.0, 2.50662827463100/
|
|
DATA ZERO, UMAX, SMALL /0.0, 8.0, 0.1E-07/
|
|
IF (ABS(U) .GT. UMAX) THEN
|
|
DO 10 I = 1, 4
|
|
FF(I) = ZERO
|
|
10 CONTINUE
|
|
ELSE
|
|
U2 = U * U
|
|
T1 = BP * EXP(-HALF * U2) / SQ2PI
|
|
T2 = BP * T1
|
|
T3 = BP * T2
|
|
FF(1) = T1
|
|
FF(2) = -U * T2
|
|
FF(3) = (U2 - ONE) * T3
|
|
FF(4) = (THREE - U2) * U * BP * T3
|
|
DO 20 I = 1, 4
|
|
IF(ABS(FF(I)) .LT. SMALL) FF(I) = ZERO
|
|
20 CONTINUE
|
|
END IF
|
|
RETURN
|
|
END
|
|
SUBROUTINE WMAX(W1, W2, W3, DLG)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 251.4 APPL.STATIST. (1989), VOL.38, NO.3
|
|
C
|
|
C
|
|
C LARGEST ABSOLUTE VALUE OF QUADRATIC FUNCTION FITTED
|
|
C TO THREE POINTS.
|
|
C
|
|
DOUBLE PRECISION W1, W2, W3, DLG, QUAD, QLIM, QMIN, ONE, TWO, B2C
|
|
DATA ONE, TWO, QMIN /1.0, 2.0, 0.00001/
|
|
DLG = MAX( ABS(W1), ABS(W3) )
|
|
QUAD = W1 - W2 * TWO + W3
|
|
QLIM = MAX( ABS(W1 - W3) / TWO , QMIN)
|
|
IF (ABS(QUAD) .LE. QLIM) RETURN
|
|
B2C = (W1 - W3) / QUAD / TWO
|
|
IF (ABS(B2C) .GE. ONE) RETURN
|
|
DLG = MAX( DLG, ABS(W2 - B2C * QUAD * B2C / TWO) )
|
|
RETURN
|
|
END
|
|
SUBROUTINE TOOSML (N, FF, F)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 251.5 APPL.STATIST. (1989), VOL.38, NO.3
|
|
C
|
|
C
|
|
C MULTIPLY FF(I) BY F FOR I = N TO 4. SET TO ZERO IF TOO SMALL.
|
|
C
|
|
DOUBLE PRECISION FF(4), F, ZERO, SMALL
|
|
INTEGER N, I
|
|
DATA ZERO, SMALL /0.0, 0.1E-12/
|
|
DO 10 I = N, 4
|
|
FF(I) = FF(I) * F
|
|
IF (ABS(FF(I)) .LE. SMALL) FF(I) = ZERO
|
|
10 CONTINUE
|
|
RETURN
|
|
END
|
|
DOUBLE PRECISION FUNCTION ALNORM(X, UPPER)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, P.424
|
|
C
|
|
C EVALUATES THE TAIL AREA OF THE STANDARDIZED NORMAL CURVE
|
|
C FROM X TO INFINITY IF UPPER IS .TRUE. OR
|
|
C FROM MINUS INFINITY TO X IF UPPER IS .FALSE.
|
|
C
|
|
DOUBLE PRECISION LTONE, UTZERO, ZERO, HALF, ONE, CON, A1, A2, A3,
|
|
$ A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7, B8, B9,
|
|
$ B10, B11, B12, X, Y, Z, ZEXP
|
|
LOGICAL UPPER, UP
|
|
C
|
|
C LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR COMPUTER
|
|
C (SEE INTRODUCTORY TEXT)
|
|
C
|
|
DATA LTONE, UTZERO /7.0, 18.66/
|
|
DATA ZERO, HALF, ONE, CON /0.0, 0.5, 1.0, 1.28/
|
|
DATA A1, A2, A3,
|
|
$ A4, A5, A6,
|
|
$ A7
|
|
$ /0.398942280444, 0.399903438504, 5.75885480458,
|
|
$ 29.8213557808, 2.62433121679, 48.6959930692,
|
|
$ 5.92885724438/
|
|
DATA B1, B2, B3,
|
|
$ B4, B5, B6,
|
|
$ B7, B8, B9,
|
|
$ B10, B11, B12
|
|
$ /0.398942280385, 3.8052E-8, 1.00000615302,
|
|
$ 3.98064794E-4, 1.98615381364, 0.151679116635,
|
|
$ 5.29330324926, 4.8385912808, 15.1508972451,
|
|
$ 0.742380924027, 30.789933034, 3.99019417011/
|
|
C
|
|
ZEXP(Z) = EXP(Z)
|
|
C
|
|
UP = UPPER
|
|
Z = X
|
|
IF (Z .GE. ZERO) GOTO 10
|
|
UP = .NOT. UP
|
|
Z = -Z
|
|
10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20
|
|
ALNORM = ZERO
|
|
GOTO 40
|
|
20 Y = HALF * Z * Z
|
|
IF (Z .GT. CON) GOTO 30
|
|
C
|
|
ALNORM = HALF - Z * (A1 - A2 * Y / (Y + A3 - A4 / (Y + A5 +
|
|
$ A6 / (Y + A7))))
|
|
GOTO 40
|
|
C
|
|
30 ALNORM = B1 * ZEXP(-Y) / (Z - B2 + B3 / (Z + B4 + B5 / (Z -
|
|
$ B6 + B7 / (Z + B8 - B9 / (Z + B10 + B11 / (Z + B12))))))
|
|
C
|
|
40 IF (.NOT. UP) ALNORM = ONE - ALNORM
|
|
RETURN
|
|
END
|
|
DOUBLE PRECISION FUNCTION PPND7 (P, IFAULT)
|
|
implicit none
|
|
C
|
|
C ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
|
|
C
|
|
C PRODUCES THE NORMAL DEVIATE Z CORRESPONDING TO A GIVEN LOWER
|
|
C TAIL AREA OF P; Z IS ACCURATE TO ABOUT 1 PART IN 10**7.
|
|
C
|
|
C THE HASH SUMS BELOW ARE THE SUMS OF THE MANTISSAS OF THE
|
|
C COEFFICIENTS. THEY ARE INCLUDED FOR USE IN CHECKING
|
|
C TRANSCRIPTION.
|
|
C
|
|
INTEGER IFAULT
|
|
DOUBLE PRECISION ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, CONST2,
|
|
* A0, A1, A2, A3, B1, B2, B3, C0, C1, C2, C3, D1, D2,
|
|
* E0, E1, E2, E3, F1, F2, P, Q, R
|
|
PARAMETER (ZERO = 0.0E0, ONE = 1.0E0, HALF = 0.5E0,
|
|
* SPLIT1 = 0.425E0, SPLIT2 = 5.0E0,
|
|
* CONST1 = 0.180625E0, CONST2 = 1.6E0)
|
|
C
|
|
C COEFFICIENTS FOR P CLOSE TO 1/2
|
|
PARAMETER (A0 = 3.3871327179D0,
|
|
* A1 = 5.0434271938D1,
|
|
* A2 = 1.5929113202D2,
|
|
* A3 = 5.9109374720D1,
|
|
* B1 = 1.7895169469D1,
|
|
* B2 = 7.8757757664D1,
|
|
* B3 = 6.7187563600D1)
|
|
C HASH SUM AB 32.3184577772
|
|
C
|
|
C COEFFICIENTS FOR P NEITHER CLOSE TO 1/2 NOR 0 OR 1
|
|
PARAMETER (C0 = 1.4234372777D0,
|
|
* C1 = 2.7568153900D0,
|
|
* C2 = 1.3067284816D0,
|
|
* C3 = 1.7023821103D-1,
|
|
* D1 = 7.3700164250D-1,
|
|
* D2 = 1.2021132975D-1)
|
|
C HASH SUM CD 15.7614929821
|
|
C
|
|
C COEFFICIENTS FOR P NEAR 0 OR 1
|
|
PARAMETER (E0 = 6.6579051150E0,
|
|
* E1 = 3.0812263860E0,
|
|
* E2 = 4.2868294337E-1,
|
|
* E3 = 1.7337203997E-2,
|
|
* F1 = 2.4197894225E-1,
|
|
* F2 = 1.2258202635E-2)
|
|
C HASH SUM EF 19.4052910204
|
|
C
|
|
IFAULT = 0
|
|
Q = P - HALF
|
|
IF (ABS(Q) .LE. SPLIT1) THEN
|
|
R = CONST1 - Q * Q
|
|
PPND7 = Q * (((A3 * R + A2) * R + A1) * R + A0) /
|
|
* (((B3 * R + B2) * R + B1) * R + ONE)
|
|
RETURN
|
|
ELSE
|
|
IF (Q .LT. 0) THEN
|
|
R = P
|
|
ELSE
|
|
R = ONE - P
|
|
ENDIF
|
|
IF (R .LE. ZERO) THEN
|
|
IFAULT = 1
|
|
PPND7 = ZERO
|
|
RETURN
|
|
ENDIF
|
|
R = SQRT(-LOG(R))
|
|
IF (R .LE. SPLIT2) THEN
|
|
R = R - CONST2
|
|
PPND7 = (((C3 * R + C2) * R + C1) * R + C0) /
|
|
* ((D2 * R + D1) * R + ONE)
|
|
ELSE
|
|
R = R - SPLIT2
|
|
PPND7 = (((E3 * R + E2) * R + E1) * R + E0) /
|
|
* ((F2 * R + F1) * R + ONE)
|
|
ENDIF
|
|
IF (Q .LT. 0) PPND7 = -PPND7
|
|
RETURN
|
|
ENDIF
|
|
END
|
|
SUBROUTINE MVSTUD(NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB,
|
|
* BND,IFLT)
|
|
implicit none
|
|
C
|
|
C COMPUTE MULTIVARIATE STUDENT INTEGRAL,
|
|
C USING MVNPRD (DUNNETT, APPL. STAT., 1989)
|
|
C IF RHO(I,J) = BPD(I)*BPD(J).
|
|
C
|
|
C IF RHO(I,J) HAS GENERAL STRUCTURE, USE
|
|
C MULNOR (SCHERVISH, APPL. STAT., 1984) AND REPLACE
|
|
C CALL MVNPRD(A,B,BPD,EPS,N,INF,IERC,HNC,PROB,BND,IFLT)
|
|
C BY CALL MULNOR(A,B,SIG,EPS,N,INF,PROB,BND,IFLT).
|
|
C
|
|
C AUTHOR: C.W. DUNNETT, MCMASTER UNVERSITY
|
|
C
|
|
C BASED ON ADAPTIVE SIMPSON'S RULE ALGORITHM
|
|
C DESCRIBED IN SHAMPINE & ALLEN: "NUMERICAL
|
|
C COMPUTING", (1974), PAGE 240.
|
|
C
|
|
C PARAMETERS ARE SAME AS IN ALGORITHM AS 251
|
|
C IN APPL. STAT. (1989), VOL. 38: 564-579
|
|
C WITH THE FOLLOWING ADDITIONS:
|
|
C NDF INTEGER INPUT DEGREES OF FREEDOM
|
|
C D REAL ARRAY INPUT NON-CENTRALITY VECTOR
|
|
C (PUT NDF = 0 FOR INFINITE D.F.)
|
|
C
|
|
DOUBLE PRECISION :: HNC,PROB,BND
|
|
INTEGER :: NN, MAXDF,I,IERC,NDF,N,IFLT
|
|
PARAMETER (NN=100, MAXDF = 150)
|
|
integer :: INF(*)
|
|
DOUBLE PRECISION :: A(*),B(*),BPD(*),D(*),F(3),
|
|
& AA(NN),BB(NN)
|
|
DOUBLE PRECISION :: ERB2, ERRB, AX,BX,XX
|
|
DOUBLE PRECISION,SAVE :: ZERO,HALF,TWO,THREE,FOUR
|
|
INTEGER :: NF
|
|
!DIMENSION A(*),B(*),BPD(*),INF(*),D(*),F(3),AA(NN),BB(NN)
|
|
DATA ZERO,HALF,TWO,THREE,FOUR / 0.0, 0.5, 2.0, 3.0, 4.0 /
|
|
!external float
|
|
DO 10 I = 1, N
|
|
AA(I) = A(I) - D(I)
|
|
BB(I) = B(I) - D(I)
|
|
10 CONTINUE
|
|
IF (NDF .LE. 0) THEN
|
|
CALL MVNPRD(AA,BB,BPD,ERRB,N,INF,IERC,HNC,PROB,BND,IFLT)
|
|
RETURN
|
|
ENDIF
|
|
BND = ZERO
|
|
IFLT = 0
|
|
|
|
ERB2 = ERRB
|
|
C
|
|
C CHECK IF D.F. EXCEED MAXDF; IF YES, THEN PROB
|
|
C IS COMPUTED BY QUADRATIC INTERPOLATION ON 1./D.F.
|
|
C
|
|
IF (NDF .LE. MAXDF) GO TO 20
|
|
CALL MVNPRD(AA,BB,BPD,ERB2,N,INF,IERC,HNC,F(1),BND,IFLT)
|
|
NF = MAXDF / 2
|
|
CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(3),BND,IFLT)
|
|
NF = NF * 2
|
|
CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(2),BND,IFLT)
|
|
XX = DBLE(NF) / DBLE(NDF)
|
|
AX = F(3) - F(2)*TWO + F(1)
|
|
BX = F(2)*FOUR - F(3) - F(1)*THREE
|
|
PROB = F(1) + XX * (AX * XX + BX) * HALF
|
|
RETURN
|
|
20 CALL SIMPSN (NDF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,PROB,BND,IFLT)
|
|
RETURN
|
|
END
|
|
SUBROUTINE SIMPSN (NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB,
|
|
* BND,IFLT)
|
|
implicit none
|
|
C
|
|
C STUDENTIZES A MULTIVARIATE INTEGRAL USING SIMPSON'S RULE.
|
|
C
|
|
double precision :: A,B,BPD,D,
|
|
* FV,F1T,F2T,F3T,
|
|
* LDIR,PSUM,ESTT,ERRR,GV,G1T,G2T,
|
|
* G3T,GSUM
|
|
double precision :: PROB, BOUNDA, BOUNDG,sTART, DAX, ERB2,ERRB,
|
|
$ EPS1, F0,G0, HNC, ERROR,DA,Z3,WT, CONTRG,DX,Z2,Z4,ESTL,ESTR,
|
|
$ ESTGL,ESTGR,CONTRB,TSUM,SUMG,DLG,ERRL,EXCESS,BND
|
|
double precision :: ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN
|
|
INTEGER :: IFLAG, IER, NDF,N, INF, IERC,LVL,IFLT
|
|
DIMENSION A(*),B(*),BPD(*),INF(*),D(*),
|
|
* FV(5),F1T(30),F2T(30),F3T(30),
|
|
* LDIR(30),PSUM(30),ESTT(30),ERRR(30),GV(5),G1T(30),G2T(30),
|
|
* G3T(30),GSUM(30)
|
|
DATA ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN /0.0,0.5,1.0,1.5,
|
|
* 2.0,4.0,6.0,0.000004/
|
|
PROB = ZERO
|
|
BOUNDA = ZERO
|
|
BOUNDG = ZERO
|
|
IFLAG = 0
|
|
IER = 0
|
|
START = -ONE
|
|
DAX = ONE
|
|
ERB2 = ERRB * HALF
|
|
EPS1 = ERB2 * HALF
|
|
CALL FUN (ZERO,NDF,A,B,BPD,ERB2,N,INF,D,F0,G0,IERC,HNC,IER)
|
|
10 FV(1) = ZERO
|
|
GV(1) = ZERO
|
|
ERROR = ZERO
|
|
DA = DAX
|
|
LVL = 1
|
|
Z3 = START + HALF*DA
|
|
CALL FUN(Z3,NDF,A,B,BPD,ERB2,N,INF,D,FV(3),GV(3),IERC,HNC,IER)
|
|
FV(5) = F0
|
|
GV(5) = G0
|
|
WT = ABS(DA) / SIX
|
|
CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5))
|
|
CONTRG = WT * (GV(1) + FOUR * GV(3) + GV(5))
|
|
LDIR(LVL) = 2
|
|
PSUM(LVL) = ZERO
|
|
GSUM(LVL) = ZERO
|
|
C
|
|
C BISECT INTERVAL; COMPUTE ESTIMATES FOR EACH HALF.
|
|
C
|
|
20 DX = HALF * DA
|
|
WT = ABS(DX) / SIX
|
|
Z2 = START + HALF * DX
|
|
CALL FUN(Z2,NDF,A,B,BPD,ERB2,N,INF,D,FV(2),GV(2),IERC,HNC,IER)
|
|
Z4 = START + ONEP5 * DX
|
|
CALL FUN(Z4,NDF,A,B,BPD,ERB2,N,INF,D,FV(4),GV(4),IERC,HNC,IER)
|
|
ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3))
|
|
ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5))
|
|
ESTGL = WT * (GV(1) + FOUR * GV(2) + GV(3))
|
|
ESTGR = WT * (GV(3) + FOUR * GV(4) + GV(5))
|
|
TSUM = ESTL + ESTR
|
|
SUMG = ESTGL + ESTGR
|
|
DLG = ABS(CONTRB - TSUM)
|
|
ERRL = DLG
|
|
C
|
|
C STOP BISECTING WHEN ACCURACY SUFFICIENT, OR IF
|
|
C INTERVAL TOO NARROW OR BISECTED TOO OFTEN.
|
|
C
|
|
30 IF (DLG .LE. EPS1) GO TO 50
|
|
IF (ABS(DX) .LE. DXMIN .OR. LVL .GE. 30) GO TO 40
|
|
C
|
|
C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF
|
|
C AND APPLY SIMPSON'S RULE TO LEFT HALF.
|
|
C
|
|
LVL = LVL + 1
|
|
LDIR(LVL) = 1
|
|
F1T(LVL) = FV(3)
|
|
F2T(LVL) = FV(4)
|
|
F3T(LVL) = FV(5)
|
|
G1T(LVL) = GV(3)
|
|
G2T(LVL) = GV(4)
|
|
G3T(LVL) = GV(5)
|
|
DA = DX
|
|
FV(5) = FV(3)
|
|
FV(3) = FV(2)
|
|
GV(5) = GV(3)
|
|
GV(3) = GV(2)
|
|
ESTT(LVL) = ESTR
|
|
CONTRB = ESTL
|
|
CONTRG = ESTGL
|
|
EPS1 = EPS1 * HALF
|
|
ERRR(LVL) = EPS1
|
|
GO TO 20
|
|
C
|
|
C ACCEPT APPROXIMATE VALUE FOR INTERVAL.
|
|
C
|
|
40 IFLAG = 11
|
|
50 ERROR = ERROR + ERRL
|
|
60 IF (LDIR(LVL) .EQ. 1) GO TO 70
|
|
TSUM = TSUM + PSUM(LVL)
|
|
SUMG = SUMG + GSUM(LVL)
|
|
LVL = LVL - 1
|
|
IF (LVL .GT. 0) GO TO 60
|
|
CONTRB = TSUM
|
|
CONTRG = SUMG
|
|
LVL = 1
|
|
DLG = ERROR
|
|
GO TO 80
|
|
C
|
|
C RESTORE SAVED INFORMATION TO PROCESS RIGHT HALF.
|
|
C
|
|
70 PSUM(LVL) = TSUM
|
|
GSUM(LVL) = SUMG
|
|
LDIR(LVL) = 2
|
|
DA = DAX / TWO**(LVL-1)
|
|
START = START + DX * TWO
|
|
FV(1) = F1T(LVL)
|
|
FV(3) = F2T(LVL)
|
|
FV(5) = F3T(LVL)
|
|
GV(1) = G1T(LVL)
|
|
GV(3) = G2T(LVL)
|
|
GV(5) = G3T(LVL)
|
|
CONTRB = ESTT(LVL)
|
|
EXCESS = EPS1 - DLG
|
|
EPS1 = ERRR(LVL)
|
|
IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS
|
|
GO TO 20
|
|
80 PROB = PROB + CONTRB
|
|
BOUNDG = BOUNDG + CONTRG
|
|
BOUNDA = BOUNDA + DLG
|
|
IF (Z4 .LE. ZERO) GO TO 90
|
|
IF (IFLT .EQ. 0) IFLT = IER
|
|
IF (IFLT .EQ. 0) IFLT = IFLAG
|
|
BOUNDA = BOUNDA + BOUNDG
|
|
IF (BND .LT. BOUNDA) BND = BOUNDA
|
|
RETURN
|
|
90 EPS1 = ERB2 * HALF
|
|
EXCESS = EPS1 - BND
|
|
IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS
|
|
START = ONE
|
|
DAX = -ONE
|
|
GO TO 10
|
|
END
|
|
DOUBLE PRECISION FUNCTION SDIST(Y,N)
|
|
implicit none
|
|
C
|
|
C COMPUTE Y**(N/2 - 1) EXP(-Y) / GAMMA(N/2)
|
|
C
|
|
C (Revised: 1994-01-19)
|
|
C
|
|
DOUBLE PRECISION :: Y,XN,TEST, ZERO, HALF, ONE, X23,SQRTPI
|
|
INTEGER :: N,JJ,JK,JKP,J
|
|
DATA ZERO, HALF, ONE, X23 / 0.0, 0.5, 1.0, -23.0 /
|
|
DATA SQRTPI / 1.77245385090552 /
|
|
SDIST = ZERO
|
|
IF (Y .LE. ZERO) RETURN
|
|
JJ = N/2 - 1
|
|
JK = 2 * JJ - N + 2
|
|
JKP = JJ - JK
|
|
SDIST = ONE
|
|
IF (JK .LT. 0) SDIST = SDIST / SQRT(Y) / SQRTPI
|
|
IF (JKP .EQ. 0) GO TO 20
|
|
XN = DBLE(N) * HALF
|
|
TEST = LOG(Y) - Y / DBLE(JKP)
|
|
IF ( TEST .LT. X23 ) THEN
|
|
SDIST = ZERO
|
|
RETURN
|
|
ENDIF
|
|
SDIST = LOG ( SDIST )
|
|
DO 10 J = 1, JKP
|
|
XN = XN - ONE
|
|
SDIST = SDIST + TEST - LOG(XN)
|
|
10 CONTINUE
|
|
IF ( SDIST .LT. X23 ) THEN
|
|
SDIST = ZERO
|
|
ELSE
|
|
SDIST = EXP( SDIST )
|
|
ENDIF
|
|
RETURN
|
|
20 SDIST = SDIST * EXP(-Y)
|
|
RETURN
|
|
END
|
|
SUBROUTINE FUN (Z,NDF,H,HL,BPD,ERB2,N,INF,D,F0,G0,IERC
|
|
* ,HNC,IER)
|
|
implicit none
|
|
double precision :: ZERO, ONE, TWO, SMALL, Z, arg, term , f0, g0
|
|
$ ,df,ERB2,HNC,BND,PROB
|
|
INTEGER NN,NDF,N, I, IER,IERC,IFLT
|
|
PARAMETER (NN=100)
|
|
DOUBLE precision :: A,B,H,HL,BPD,D, SDIST
|
|
integer :: INF
|
|
DIMENSION A(NN),B(NN),H(*),HL(*),BPD(*),INF(*),D(*)
|
|
DATA ZERO, ONE, TWO, SMALL / 0.0, 1.0, 2.0, 1.0E-08 /
|
|
external SDIST
|
|
F0 = ZERO
|
|
G0 = ZERO
|
|
IF (Z .LE. -ONE .OR. Z .GE. ONE) RETURN
|
|
DF = DBLE(NDF)
|
|
ARG = (ONE + Z) / (ONE - Z)
|
|
TERM = ARG * DF * TWO / (ONE-Z)**2 * SDIST(DF/TWO*ARG*ARG,NDF)
|
|
IF (TERM .LE. SMALL) RETURN
|
|
DO 10 I = 1, N
|
|
A(I) = ARG * H(I) - D(I)
|
|
B(I) = ARG * HL(I) - D(I)
|
|
10 CONTINUE
|
|
CALL MVNPRD (A,B,BPD,ERB2,N,INF,IERC,HNC,PROB,BND,IFLT)
|
|
IF (IER .EQ. 0) IER = IFLT
|
|
G0 = TERM * BND
|
|
F0 = TERM * PROB
|
|
RETURN
|
|
END
|
|
|
|
C * * * * * * * * * * * * * * * * * * * * * * * * * * * *
|
|
C Charles Dunnett
|
|
C Dept. of Mathematics and Statistics
|
|
C McMaster University
|
|
C Hamilton, Ontario L8S 4K1
|
|
C Canada
|
|
C E-mail: dunnett@mcmaster.ca
|
|
C Tel.: (905) 525-9140 (Ext. 27104)
|
|
C * * * * * * * * * * * * * * * * * * * * * * * * * * * *
|
|
|
|
|