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.
1596 lines
56 KiB
Fortran
1596 lines
56 KiB
Fortran
MODULE TRIVARIATEVAR
|
|
! Global variables used in calculation of TRIVARIATE
|
|
! normal and TRIVARIATE student T probabilties
|
|
INTEGER :: NU
|
|
DOUBLE PRECISION :: H1, H2, H3, R23, RUA, RUB, AR, RUC
|
|
END MODULE TRIVARIATEVAR
|
|
!
|
|
! FIMOD contains functions for calculating 1D, 2D and 3D Normal and student T probabilites
|
|
! and 1D expectations
|
|
MODULE FIMOD
|
|
! USE ERFCOREMOD
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
PUBLIC :: NORMPRB, FI, FIINV, MVNLIMITS, MVNLMS, BVU,BVNMVN
|
|
PUBLIC :: GAUSINT, GAUSINT2, EXLMS, EXINV
|
|
PUBLIC :: STUDNT, BVTL, TVTL, TVNMVN
|
|
|
|
INTERFACE NORMPRB
|
|
MODULE PROCEDURE NORMPRB
|
|
END INTERFACE
|
|
|
|
INTERFACE FI
|
|
MODULE PROCEDURE FI
|
|
END INTERFACE
|
|
|
|
INTERFACE FI2
|
|
MODULE PROCEDURE FI2
|
|
END INTERFACE
|
|
|
|
INTERFACE FIINV
|
|
MODULE PROCEDURE FIINV
|
|
END INTERFACE
|
|
|
|
INTERFACE MVNLIMITS
|
|
MODULE PROCEDURE MVNLIMITS
|
|
END INTERFACE
|
|
|
|
INTERFACE MVNLMS
|
|
MODULE PROCEDURE MVNLMS
|
|
END INTERFACE
|
|
|
|
INTERFACE BVU
|
|
MODULE PROCEDURE BVU
|
|
END INTERFACE
|
|
|
|
INTERFACE BVNMVN
|
|
MODULE PROCEDURE BVNMVN
|
|
END INTERFACE
|
|
|
|
INTERFACE STUDNT
|
|
MODULE PROCEDURE STUDNT
|
|
END INTERFACE
|
|
|
|
INTERFACE BVTL
|
|
MODULE PROCEDURE BVTL
|
|
END INTERFACE
|
|
|
|
INTERFACE TVTL
|
|
MODULE PROCEDURE TVTL
|
|
END INTERFACE
|
|
|
|
INTERFACE GAUSINT
|
|
MODULE PROCEDURE GAUSINT
|
|
END INTERFACE
|
|
|
|
INTERFACE GAUSINT2
|
|
MODULE PROCEDURE GAUSINT2
|
|
END INTERFACE
|
|
|
|
INTERFACE EXLMS
|
|
MODULE PROCEDURE EXLMS
|
|
END INTERFACE
|
|
|
|
INTERFACE EXINV
|
|
MODULE PROCEDURE EXINV
|
|
END INTERFACE
|
|
|
|
CONTAINS
|
|
FUNCTION FIINV(P) RESULT (VAL)
|
|
IMPLICIT NONE
|
|
*
|
|
* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
|
|
*
|
|
* Produces the normal deviate Z corresponding to a given lower
|
|
* tail area of P.
|
|
* Absolute error less than 1e-13
|
|
* Relative error less than 1e-15 for abs(VAL)>0.1
|
|
*
|
|
* The hash sums below are the sums of the mantissas of the
|
|
* coefficients. They are included for use in checking
|
|
* transcription.
|
|
*
|
|
DOUBLE PRECISION, INTENT(in) :: P
|
|
DOUBLE PRECISION :: VAL
|
|
!local variables
|
|
DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2, ONE, ZERO, HALF,
|
|
& A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
|
|
& C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
|
|
& E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
|
|
& Q, R
|
|
PARAMETER ( SPLIT1 = 0.425D0, SPLIT2 = 5.D0,
|
|
& CONST1 = 0.180625D0, CONST2 = 1.6D0,
|
|
& ONE = 1.D0, ZERO = 0.D0, HALF = 0.5D0 )
|
|
*
|
|
* Coefficients for P close to 0.5
|
|
*
|
|
PARAMETER (
|
|
* A0 = 3.38713 28727 96366 6080D0,
|
|
* A1 = 1.33141 66789 17843 7745D+2,
|
|
* A2 = 1.97159 09503 06551 4427D+3,
|
|
* A3 = 1.37316 93765 50946 1125D+4,
|
|
* A4 = 4.59219 53931 54987 1457D+4,
|
|
* A5 = 6.72657 70927 00870 0853D+4,
|
|
* A6 = 3.34305 75583 58812 8105D+4,
|
|
* A7 = 2.50908 09287 30122 6727D+3,
|
|
* B1 = 4.23133 30701 60091 1252D+1,
|
|
* B2 = 6.87187 00749 20579 0830D+2,
|
|
* B3 = 5.39419 60214 24751 1077D+3,
|
|
* B4 = 2.12137 94301 58659 5867D+4,
|
|
* B5 = 3.93078 95800 09271 0610D+4,
|
|
* B6 = 2.87290 85735 72194 2674D+4,
|
|
* B7 = 5.22649 52788 52854 5610D+3 )
|
|
* HASH SUM AB 55.88319 28806 14901 4439
|
|
*
|
|
* Coefficients for P not close to 0, 0.5 or 1.
|
|
*
|
|
PARAMETER (
|
|
* C0 = 1.42343 71107 49683 57734D0,
|
|
* C1 = 4.63033 78461 56545 29590D0,
|
|
* C2 = 5.76949 72214 60691 40550D0,
|
|
* C3 = 3.64784 83247 63204 60504D0,
|
|
* C4 = 1.27045 82524 52368 38258D0,
|
|
* C5 = 2.41780 72517 74506 11770D-1,
|
|
* C6 = 2.27238 44989 26918 45833D-2,
|
|
* C7 = 7.74545 01427 83414 07640D-4,
|
|
* D1 = 2.05319 16266 37758 82187D0,
|
|
* D2 = 1.67638 48301 83803 84940D0,
|
|
* D3 = 6.89767 33498 51000 04550D-1,
|
|
* D4 = 1.48103 97642 74800 74590D-1,
|
|
* D5 = 1.51986 66563 61645 71966D-2,
|
|
* D6 = 5.47593 80849 95344 94600D-4,
|
|
* D7 = 1.05075 00716 44416 84324D-9 )
|
|
* HASH SUM CD 49.33206 50330 16102 89036
|
|
*
|
|
* Coefficients for P near 0 or 1.
|
|
*
|
|
PARAMETER (
|
|
* E0 = 6.65790 46435 01103 77720D0,
|
|
* E1 = 5.46378 49111 64114 36990D0,
|
|
* E2 = 1.78482 65399 17291 33580D0,
|
|
* E3 = 2.96560 57182 85048 91230D-1,
|
|
* E4 = 2.65321 89526 57612 30930D-2,
|
|
* E5 = 1.24266 09473 88078 43860D-3,
|
|
* E6 = 2.71155 55687 43487 57815D-5,
|
|
* E7 = 2.01033 43992 92288 13265D-7,
|
|
* F1 = 5.99832 20655 58879 37690D-1,
|
|
* F2 = 1.36929 88092 27358 05310D-1,
|
|
* F3 = 1.48753 61290 85061 48525D-2,
|
|
* F4 = 7.86869 13114 56132 59100D-4,
|
|
* F5 = 1.84631 83175 10054 68180D-5,
|
|
* F6 = 1.42151 17583 16445 88870D-7,
|
|
* F7 = 2.04426 31033 89939 78564D-15 )
|
|
* HASH SUM EF 47.52583 31754 92896 71629
|
|
*
|
|
Q = ( P - HALF)
|
|
IF ( ABS(Q) .LE. SPLIT1 ) THEN ! Central range.
|
|
R = CONST1 - Q*Q
|
|
VAL = Q*( ( ( ((((A7*R + A6)*R + A5)*R + A4)*R + A3)
|
|
* *R + A2 )*R + A1 )*R + A0 )
|
|
* /( ( ( ((((B7*R + B6)*R + B5)*R + B4)*R + B3)
|
|
* *R + B2 )*R + B1 )*R + ONE)
|
|
ELSE ! near the endpoints
|
|
R = MIN( P, ONE - P )
|
|
IF (R .GT.ZERO) THEN ! ( 2.d0*R .GT. CFxCutOff) THEN ! R .GT.0.d0
|
|
R = SQRT( -LOG(R) )
|
|
IF ( R .LE. SPLIT2 ) THEN
|
|
R = R - CONST2
|
|
VAL = ( ( ( ((((C7*R + C6)*R + C5)*R + C4)*R + C3)
|
|
* *R + C2 )*R + C1 )*R + C0 )
|
|
* /( ( ( ((((D7*R + D6)*R + D5)*R + D4)*R + D3)
|
|
* *R + D2 )*R + D1 )*R + ONE )
|
|
ELSE
|
|
R = R - SPLIT2
|
|
VAL = ( ( ( ((((E7*R + E6)*R + E5)*R + E4)*R + E3)
|
|
* *R + E2 )*R + E1 )*R + E0 )
|
|
* /( ( ( ((((F7*R + F6)*R + F5)*R + F4)*R + F3)
|
|
* *R + F2 )*R + F1 )*R + ONE )
|
|
END IF
|
|
ELSE
|
|
VAL = 37.0d0 !9.D0 !XMAX 9.d0
|
|
END IF
|
|
IF ( Q < ZERO ) VAL = - VAL
|
|
END IF
|
|
RETURN
|
|
END FUNCTION FIINV
|
|
! *********************************
|
|
SUBROUTINE NORMPRB(Z, P, Q)
|
|
! USE ERFCOREMOD
|
|
! USE GLOBALDATA, ONLY : XMAX
|
|
! Normal distribution probabilities accurate to 18 digits between
|
|
! -XMAX and XMAX
|
|
!
|
|
! Z = no. of standard deviations from the mean.
|
|
! P, Q = probabilities to the left & right of Z. P + Q = 1.
|
|
!
|
|
! by pab 23.03.2003
|
|
!
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(IN) :: Z
|
|
DOUBLE PRECISION, INTENT(OUT) :: P
|
|
DOUBLE PRECISION, INTENT(OUT), OPTIONAL :: Q
|
|
!Local variables
|
|
DOUBLE PRECISION :: PP, ZABS
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0
|
|
DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0
|
|
DOUBLE PRECISION, PARAMETER :: SQ2M1 = 0.70710678118655D0 ! 1/SQRT(2)
|
|
DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0
|
|
DOUBLE PRECISION, PARAMETER :: XMAX = 37D0
|
|
ZABS = ABS(Z)
|
|
!
|
|
! |Z| > 37 (or XMAX)
|
|
!
|
|
IF ( ZABS .GT. XMAX ) THEN
|
|
IF (Z > ZERO) THEN
|
|
P = ONE
|
|
IF (PRESENT(Q)) Q = ZERO
|
|
ELSE
|
|
P = ZERO
|
|
IF (PRESENT(Q)) Q = ONE
|
|
END IF
|
|
ELSE
|
|
!
|
|
! |Z| <= 37
|
|
!
|
|
PP = DERFC(ZABS*SQ2M1)*HALF
|
|
|
|
IF (Z < ZERO) THEN
|
|
P = PP
|
|
IF (PRESENT(Q)) Q = ONE - PP
|
|
ELSE
|
|
P = ONE - PP
|
|
IF (PRESENT(Q)) Q = PP
|
|
END IF
|
|
END IF
|
|
|
|
RETURN
|
|
END SUBROUTINE NORMPRB
|
|
FUNCTION FI( Z ) RESULT (VALUE)
|
|
! USE ERFCOREMOD
|
|
! USE GLOBALDATA, ONLY : XMAX
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: Z
|
|
DOUBLE PRECISION :: VALUE
|
|
! Local variables
|
|
DOUBLE PRECISION :: ZABS
|
|
DOUBLE PRECISION, PARAMETER:: SQ2M1 = 0.70710678118655D0 ! 1/SQRT(2)
|
|
DOUBLE PRECISION, PARAMETER:: HALF = 0.5D0
|
|
DOUBLE PRECISION, PARAMETER:: XMAX = 37.D0
|
|
ZABS = ABS(Z)
|
|
*
|
|
* |Z| > 37 (or XMAX)
|
|
*
|
|
IF ( ZABS .GT. XMAX ) THEN
|
|
IF (Z < 0.0D0) THEN
|
|
VALUE = 0.0D0
|
|
ELSE
|
|
VALUE = 1.0D0
|
|
ENDIF
|
|
ELSE
|
|
VALUE = DERFC(-Z*SQ2M1)*HALF
|
|
ENDIF
|
|
RETURN
|
|
END FUNCTION FI
|
|
|
|
FUNCTION FI2( Z ) RESULT (VALUE)
|
|
! USE GLOBALDATA, ONLY : XMAX
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: Z
|
|
DOUBLE PRECISION :: VALUE
|
|
*
|
|
* Normal distribution probabilities accurate to 1.e-15.
|
|
* relative error less than 1e-8;
|
|
* Z = no. of standard deviations from the mean.
|
|
*
|
|
* Based upon algorithm 5666 for the error function, from:
|
|
* Hart, J.F. et al, 'Computer Approximations', Wiley 1968
|
|
*
|
|
* Programmer: Alan Miller
|
|
*
|
|
* Latest revision - 30 March 1986
|
|
*
|
|
DOUBLE PRECISION :: P0, P1, P2, P3, P4, P5, P6,
|
|
* Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,XMAX,
|
|
* P, EXPNTL, CUTOFF, ROOTPI, ZABS, Z2
|
|
PARAMETER(
|
|
* P0 = 220.20 68679 12376 1D0,
|
|
* P1 = 221.21 35961 69931 1D0,
|
|
* P2 = 112.07 92914 97870 9D0,
|
|
* P3 = 33.912 86607 83830 0D0,
|
|
* P4 = 6.3739 62203 53165 0D0,
|
|
* P5 = 0.70038 30644 43688 1D0,
|
|
* P6 = 0.035262 49659 98910 9D0 )
|
|
PARAMETER(
|
|
* Q0 = 440.41 37358 24752 2D0,
|
|
* Q1 = 793.82 65125 19948 4D0,
|
|
* Q2 = 637.33 36333 78831 1D0,
|
|
* Q3 = 296.56 42487 79673 7D0,
|
|
* Q4 = 86.780 73220 29460 8D0,
|
|
* Q5 = 16.064 17757 92069 5D0,
|
|
* Q6 = 1.7556 67163 18264 2D0,
|
|
* Q7 = 0.088388 34764 83184 4D0 )
|
|
PARAMETER( ROOTPI = 2.5066 28274 63100 1D0 )
|
|
PARAMETER( CUTOFF = 7.0710 67811 86547 5D0 )
|
|
PARAMETER( XMAX = 37.D0 )
|
|
*
|
|
ZABS = ABS(Z)
|
|
*
|
|
* |Z| > 37 (or XMAX)
|
|
*
|
|
IF ( ZABS .GT. XMAX ) THEN
|
|
P = 0.d0
|
|
ELSE
|
|
*
|
|
* |Z| <= 37
|
|
*
|
|
Z2 = ZABS * ZABS
|
|
EXPNTL = EXP( -Z2 * 0.5D0 )
|
|
*
|
|
* |Z| < CUTOFF = 10/SQRT(2)
|
|
*
|
|
IF ( ZABS < CUTOFF ) THEN
|
|
P = EXPNTL*( (((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
|
|
* + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
|
|
* + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
|
|
* + Q0 )
|
|
*
|
|
* |Z| >= CUTOFF.
|
|
*
|
|
ELSE
|
|
P = EXPNTL/( ZABS + 1.d0/( ZABS + 2.d0/( ZABS + 3.d0/( ZABS
|
|
* + 4.d0/( ZABS + 0.65D0 ) ) ) ) )/ROOTPI
|
|
END IF
|
|
END IF
|
|
IF ( Z .GT. 0.d0 ) P = 1.d0 - P
|
|
VALUE = P
|
|
RETURN
|
|
END FUNCTION FI2
|
|
|
|
SUBROUTINE MVNLIMITS( A, B, INFIN, AP, PRB, AQ)
|
|
! RETURN probabilities for being between A and B
|
|
! WHERE
|
|
! AP = FI(A), AQ = 1 - FI(A)
|
|
! BP = FI(B), BQ = 1 - FI(B)
|
|
! PRB = BP-AP IF BP+AP<1
|
|
! = AQ-BQ OTHERWISE
|
|
!
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: A, B
|
|
DOUBLE PRECISION, INTENT(out) :: AP
|
|
DOUBLE PRECISION, INTENT(out),OPTIONAL :: PRB,AQ
|
|
INTEGER,INTENT(in) :: INFIN
|
|
! LOCAL VARIABLES
|
|
DOUBLE PRECISION :: BP,AQQ, BQQ
|
|
DOUBLE PRECISION, PARAMETER :: ONE=1.D0, ZERO = 0.D0
|
|
|
|
SELECT CASE (infin)
|
|
CASE (:-1)
|
|
AP = ZERO
|
|
! BP = ONE
|
|
IF (PRESENT(PRB)) PRB = ONE
|
|
IF (PRESENT(AQ)) AQ = ONE
|
|
! IF (PRESENT(BQ)) BQ = ZERO
|
|
CASE (0)
|
|
AP = ZERO
|
|
CALL NORMPRB(B,BP) !,BQQ)
|
|
IF (PRESENT(PRB)) PRB = BP
|
|
IF (PRESENT(AQ)) AQ = ONE
|
|
! IF (PRESENT(BQ)) BQ = BQQ
|
|
CASE (1)
|
|
! BP = ONE
|
|
CALL NORMPRB(A,AP,AQQ)
|
|
IF (PRESENT(PRB)) PRB = AQQ
|
|
IF (PRESENT(AQ)) AQ = AQQ
|
|
! IF (PRESENT(BQ)) BQ = ZERO
|
|
CASE (2:)
|
|
CALL NORMPRB(A,AP,AQQ)
|
|
CALL NORMPRB(B,BP,BQQ)
|
|
IF (PRESENT(PRB)) THEN
|
|
IF (AP+BP < ONE) THEN
|
|
PRB = BP - AP
|
|
ELSE
|
|
PRB = AQQ - BQQ
|
|
END IF
|
|
ENDIF
|
|
IF (PRESENT(AQ)) AQ = AQQ
|
|
! IF (PRESENT(BQ)) BQ = BQQ
|
|
END SELECT
|
|
RETURN
|
|
END SUBROUTINE MVNLIMITS
|
|
|
|
|
|
SUBROUTINE MVNLMS( A, B, INFIN, LOWER, UPPER )
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: A, B
|
|
DOUBLE PRECISION, INTENT(out) :: LOWER, UPPER
|
|
INTEGER,INTENT(in) :: INFIN
|
|
|
|
LOWER = 0.0D0
|
|
UPPER = 1.0D0
|
|
IF ( INFIN < 0 ) RETURN
|
|
IF ( INFIN .NE. 0 ) LOWER = FI(A)
|
|
IF ( INFIN .NE. 1 ) UPPER = FI(B)
|
|
RETURN
|
|
END SUBROUTINE MVNLMS
|
|
|
|
FUNCTION TVNMVN(A, B, INFIN, R, EPSI ) RESULT (VAL)
|
|
IMPLICIT NONE
|
|
*
|
|
* A function for computing trivariate normal probabilities.
|
|
*
|
|
* Parameters
|
|
*
|
|
* A REAL, array of lower integration limits.
|
|
* B REAL, array of upper integration limits.
|
|
! R REAL, array of correlation coefficents
|
|
! R = [r12 r13 r23]
|
|
! EPSI = REAL tolerance
|
|
* INFIN INTEGER, array of integration limits flags:
|
|
* if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
|
|
* if INFIN(I) = 1, Ith limits are [A(I), infinity);
|
|
* if INFIN(I) = 2, Ith limits are [A(I), B(I)].
|
|
DOUBLE PRECISION, DIMENSION(:), INTENT (IN) :: A, B , R
|
|
DOUBLE PRECISION, INTENT (IN) :: EPSI
|
|
INTEGER, DIMENSION(:), INTENT (IN) :: INFIN
|
|
DOUBLE PRECISION :: VAL
|
|
|
|
IF ( INFIN(1) .EQ. 2 ) THEN
|
|
IF ( INFIN(2) .EQ. 2 ) THEN
|
|
IF (INFIN(3) .EQ. 2 ) THEN !OK
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI )
|
|
& + TVNL( A(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& + TVNL( A(1), B(2), A(3),R(1),R(2),R(3),EPSI )
|
|
$ + TVNL( B(1), A(2), A(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( A(1), A(2), A(3),R(1),R(2),R(3),EPSI )
|
|
ELSE IF (INFIN(3) .EQ. 1 ) THEN ! B(3) = inf ok
|
|
VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
& - TVNL( A(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
& - TVNL( B(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
& + TVNL( A(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
ELSE IF (INFIN(3) .EQ. 0 ) THEN !OK A(3) = -inf
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& + TVNL( A(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
ELSE ! INFIN(1:2)=2
|
|
VAL = BVNMVN( A ,B, INFIN, R(1) )
|
|
ENDIF
|
|
ELSE IF (INFIN(2) .EQ. 1 ) THEN ! B(2) = inf
|
|
IF (INFIN(3) .EQ. 2 ) THEN
|
|
VAL = TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
|
|
& - TVNL( A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
|
|
& - TVNL( B(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI )
|
|
& + TVNL( A(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI )
|
|
ELSE IF (INFIN(3) .EQ. 1 ) THEN
|
|
VAL = TVNL( B(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI )
|
|
$ - TVNL( A(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI )
|
|
ELSE IF (INFIN(3) .EQ. 0 ) THEN
|
|
VAL = TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
|
|
$ - TVNL( A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
|
|
ELSE
|
|
VAL = BVNMVN( A ,B, INFIN, R(1) )
|
|
ENDIF
|
|
ELSE IF (INFIN(2) .EQ. 0 ) THEN
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2:) ! % % A(2)=-INF
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
$ - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
$ - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI )
|
|
$ + TVNL( A(1), B(2), A(3),R(1),R(2),R(3),EPSI )
|
|
CASE (1) !% % A(2)=-INF B(3) = INF
|
|
VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
|
|
$ - TVNL(A(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
|
|
CASE (0) ! % % A(2)=-INF A(3) = -INF
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
$ - TVNL( A(1), B(2), B(3),R(1),R(2),R(3),EPSI)
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
ELSE
|
|
VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
|
|
ENDIF
|
|
ELSE IF ( INFIN(1) .EQ. 1 ) THEN
|
|
SELECT CASE (INFIN(2))
|
|
CASE (2)
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2) !% B(1) = INF %OK
|
|
VAL = TVNL(-A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
|
|
$ - TVNL(-A(1), B(2), A(3),-R(1),-R(2),R(3),EPSI )
|
|
$ - TVNL(-A(1), A(2), B(3),-R(1),-R(2),R(3),EPSI )
|
|
$ + TVNL(-A(1), A(2), A(3),-R(1),-R(2),R(3),EPSI )
|
|
CASE (1) ! % B(1) = INF B(3) = INF %OK
|
|
VAL = TVNL(-A(1), B(2), -A(3),-R(1),R(2),-R(3),EPSI )
|
|
$ - TVNL(-A(1), A(2), -A(3),-R(1),R(2),-R(3),EPSI)
|
|
CASE (0) ! % B(1) = INF A(3) = -INF %OK
|
|
VAL = TVNL(-A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
|
|
$ - TVNL(-A(1), A(2), B(3),-R(1),-R(2),R(3),EPSI)
|
|
CASE (-1)
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE (1) !%B(2) = INF
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2) ! % B(1) = INF B(2) = INF % OK
|
|
VAL = TVNL( -A(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI )
|
|
& - TVNL( -A(1), -A(2),A(3),-R(1),R(2),-R(3),EPSI )
|
|
CASE (1) ! % B(1:3) = INF %OK
|
|
VAL = TVNL( -A(1), -A(2), -A(3),R(1),R(2),R(3),EPSI)
|
|
CASE (0) ! % B(1:2) = INF A(3) = -INF %OK
|
|
VAL = TVNL( -A(1), -A(2), B(3),R(1),-R(2),-R(3),EPSI )
|
|
CASE (:-1)
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE (0) ! A(2) = -INF
|
|
SELECT CASE ( INFIN(3))
|
|
CASE (2) ! B(1) = INF , A(2) = -INF %OK
|
|
VAL = TVNL( -A(1), B(2), B(3),-R(1),R(2),-R(3),EPSI )
|
|
& - TVNL( -A(1), B(2),A(3),-R(1),R(2),-R(3),EPSI )
|
|
CASE (1) ! B(1) = INF , A(2) = -INF B(3) = INF % OK
|
|
VAL = TVNL( -A(1), B(2), -A(3),-R(1),-R(2),R(3),EPSI)
|
|
CASE (0) !% B(1) = INF , A(2:3) = -INF
|
|
VAL = TVNL( -A(1), B(2), B(3),-R(1),-R(2),R(3),EPSI )
|
|
CASE (:-1)
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
|
|
END SELECT
|
|
ELSE IF ( INFIN(1) .EQ. 0 ) THEN
|
|
SELECT CASE (INFIN(2))
|
|
CASE (2)
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2:) ! A(1) = -INF %OK
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI)
|
|
& - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& + TVNL( B(1), A(2), A(3),R(1),R(2),R(3),EPSI )
|
|
CASE (1) ! % A(1) = -INF , B(3) = INF %OK
|
|
VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
$ - TVNL( B(1), A(2), -A(3),R(1),-R(2),-R(3),EPSI )
|
|
CASE (0) ! A(1) = -INF , A(3) = -INF %OK
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
& - TVNL( B(1), A(2), B(3),R(1),R(2),R(3),EPSI )
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE (1) ! B(2) = INF
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2:) ! A(1) = -INF B(2) = INF %OK
|
|
VAL = TVNL( B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
|
|
$ - TVNL( B(1), -A(2), A(3),-R(1),R(2),-R(3),EPSI)
|
|
CASE (1) ! A(1) = -INF B(2) = INF B(3) = INF %OK
|
|
VAL = TVNL( B(1), -A(2), -A(3),-R(1),-R(2),R(3),EPSI)
|
|
CASE (0) ! % A(1) = -INF B(2) = INF A(3) = -INF %OK
|
|
VAL = TVNL(B(1), -A(2), B(3),-R(1),R(2),-R(3),EPSI)
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE (0) ! A(2) = -INF
|
|
SELECT CASE (INFIN(3))
|
|
CASE (2:) ! % A(1:2) = -INF
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI)
|
|
$ - TVNL( B(1), B(2), A(3),R(1),R(2),R(3),EPSI)
|
|
CASE (1) ! A(1:2) = -INF B(3) = INF
|
|
VAL = TVNL( B(1), B(2), -A(3),R(1),-R(2),-R(3),EPSI)
|
|
CASE (0) ! % A(1:3) = -INF
|
|
VAL = TVNL( B(1), B(2), B(3),R(1),R(2),R(3),EPSI )
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A,B,INFIN,R(1))
|
|
END SELECT
|
|
CASE DEFAULT
|
|
VAL = BVNMVN(A(1:3:2),B(1:3:2),INFIN(1:3:2),R(2))
|
|
END SELECT
|
|
ELSE
|
|
VAL = BVNMVN(A(2:3),B(2:3),INFIN(2:3),R(3))
|
|
END IF
|
|
CONTAINS
|
|
DOUBLE PRECISION FUNCTION TVNL(H1,H2,H3, R12,R13,R23, EPSI )
|
|
!Returns Trivariate Normal CDF
|
|
DOUBLE PRECISION, INTENT(IN) :: R12,R13,R23
|
|
DOUBLE PRECISION, INTENT(IN) :: H1,H2,H3, EPSI
|
|
! Locals
|
|
INTEGER, PARAMETER :: NU = 0
|
|
DOUBLE PRECISION,DIMENSION(3) :: H,R
|
|
H(:) = (/ H1, H2, H3 /)
|
|
R(:) = (/ R12, R13, R23 /)
|
|
TVNL = TVTL(NU,H,R,EPSI)
|
|
END FUNCTION TVNL
|
|
END FUNCTION TVNMVN
|
|
FUNCTION BVNMVN( LOWER, UPPER, INFIN, CORREL ) RESULT (VAL)
|
|
IMPLICIT NONE
|
|
*
|
|
* A function for computing bivariate normal probabilities.
|
|
*
|
|
* Parameters
|
|
*
|
|
* LOWER REAL, array of lower integration limits.
|
|
* UPPER REAL, array of upper integration limits.
|
|
* INFIN INTEGER, array of integration limits flags:
|
|
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
|
|
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
|
|
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
|
|
* CORREL REAL, correlation coefficient.
|
|
*
|
|
|
|
DOUBLE PRECISION, DIMENSION(:), INTENT (IN) :: LOWER, UPPER
|
|
DOUBLE PRECISION, INTENT (IN) :: CORREL
|
|
INTEGER, DIMENSION(:), INTENT (IN) :: INFIN
|
|
DOUBLE PRECISION :: VAL
|
|
DOUBLE PRECISION :: E
|
|
SELECT CASE (INFIN(1))
|
|
CASE (2:)
|
|
SELECT CASE ( INFIN(2) )
|
|
CASE (2:)
|
|
VAL = BVU ( LOWER(1), LOWER(2), CORREL )
|
|
& - BVU ( UPPER(1), LOWER(2), CORREL )
|
|
& - BVU ( LOWER(1), UPPER(2), CORREL )
|
|
& + BVU ( UPPER(1), UPPER(2), CORREL )
|
|
|
|
CASE (1)
|
|
VAL = BVU ( LOWER(1), LOWER(2), CORREL )
|
|
& - BVU ( UPPER(1), LOWER(2), CORREL )
|
|
CASE (0)
|
|
VAL = BVU ( -UPPER(1), -UPPER(2), CORREL )
|
|
& - BVU ( -LOWER(1), -UPPER(2), CORREL )
|
|
CASE DEFAULT
|
|
CALL MVNLIMITS(LOWER(1),UPPER(1),INFIN(1),E,VAL)
|
|
END SELECT
|
|
CASE (1)
|
|
SELECT CASE ( INFIN(2))
|
|
CASE ( 2: )
|
|
VAL = BVU ( LOWER(1), LOWER(2), CORREL )
|
|
& - BVU ( LOWER(1), UPPER(2), CORREL )
|
|
CASE (1)
|
|
VAL = BVU ( LOWER(1), LOWER(2), CORREL )
|
|
CASE (0)
|
|
VAL = BVU ( LOWER(1), -UPPER(2), -CORREL )
|
|
CASE DEFAULT
|
|
CALL MVNLIMITS(LOWER(2),UPPER(2),INFIN(2),E,VAL)
|
|
END SELECT
|
|
CASE (0)
|
|
SELECT CASE ( INFIN(2))
|
|
CASE ( 2: )
|
|
VAL = BVU ( -UPPER(1), -UPPER(2), CORREL )
|
|
& - BVU ( -UPPER(1), -LOWER(2), CORREL )
|
|
CASE ( 1 )
|
|
VAL = BVU ( -UPPER(1), LOWER(2), -CORREL )
|
|
CASE (0)
|
|
VAL = BVU ( -UPPER(1), -UPPER(2), CORREL )
|
|
CASE DEFAULT
|
|
CALL MVNLIMITS(LOWER(1),UPPER(1),INFIN(1),E,VAL)
|
|
END SELECT
|
|
CASE DEFAULT !ELSE !INFIN(1)<0
|
|
CALL MVNLIMITS(LOWER(2),UPPER(2),INFIN(2),E,VAL)
|
|
END SELECT
|
|
END FUNCTION BVNMVN
|
|
FUNCTION BVU( SH, SK, R ) RESULT (VAL)
|
|
! USE GLOBALDATA, ONLY: XMAX
|
|
IMPLICIT NONE
|
|
*
|
|
! A function for computing bivariate normal probabilities.
|
|
!
|
|
! Yihong Ge
|
|
! Department of Computer Science and Electrical Engineering
|
|
! Washington State University
|
|
! Pullman, WA 99164-2752
|
|
! and
|
|
! Alan Genz
|
|
! Department of Mathematics
|
|
! Washington State University
|
|
! Pullman, WA 99164-3113
|
|
! Email : alangenz@wsu.edu
|
|
!
|
|
! This function is based on the method described by
|
|
! Drezner, Z and G.O. Wesolowsky, (1989),
|
|
! On the computation of the bivariate normal integral,
|
|
! Journal of Statist. Comput. Simul. 35, pp. 101-107,
|
|
! with major modifications for double precision, and for |R| close to 1.
|
|
!
|
|
! BVU - calculate the probability that X > SH and Y > SK.
|
|
! (to accuracy of 1e-16?)
|
|
!
|
|
! Parameters
|
|
!
|
|
! SH REAL, lower integration limit
|
|
! SK REAL, lower integration limit
|
|
! R REAL, correlation coefficient
|
|
!
|
|
! LG INTEGER, number of Gauss Rule Points and Weights
|
|
!
|
|
! Revised pab added check on XMAX
|
|
DOUBLE PRECISION, INTENT(IN) :: SH, SK, R
|
|
DOUBLE PRECISION :: VAL
|
|
! Local variables
|
|
DOUBLE PRECISION :: ZERO,ONE,FOUR
|
|
DOUBLE PRECISION :: SQTWOPI ,TWOPI1,FOURPI1
|
|
DOUBLE PRECISION :: HALF,ONETHIRD,ONEEIGHT,ONESIXTEEN
|
|
DOUBLE PRECISION :: TWELVE, EXPMIN, XMAX
|
|
INTEGER :: I, LG, NG
|
|
PARAMETER ( ZERO = 0.D0,ONE=1.0D0,HALF=0.5D0)
|
|
PARAMETER (FOUR = 4.0D0, TWELVE = 12.0D0)
|
|
PARAMETER (EXPMIN = -100.0D0)
|
|
PARAMETER (ONESIXTEEN = 0.0625D0) !1/16
|
|
PARAMETER (ONEEIGHT = 0.125D0 ) !1/8
|
|
PARAMETER (ONETHIRD = 0.3333333333333333333333D0)
|
|
! PARAMETER (TWOPI = 6.283185307179586D0 )
|
|
PARAMETER (TWOPI1 = 0.15915494309190D0 ) !1/(2*pi)
|
|
PARAMETER (FOURPI1 = 0.0795774715459476D0 ) !/1/(4*pi)
|
|
PARAMETER (SQTWOPI = 2.50662827463100D0) ! SQRT(2*pi)
|
|
PARAMETER (XMAX = 8.3D0)
|
|
DOUBLE PRECISION, DIMENSION(10,3) :: X, W
|
|
DOUBLE PRECISION :: AS, A, B, C, D, RS, XS
|
|
DOUBLE PRECISION :: SN, ASR, H, K, BS, HS, HK
|
|
! Gauss Legendre Points and Weights, N = 6
|
|
DATA ( W(I,1), X(I,1), I = 1,3) /
|
|
* 0.1713244923791705D+00,-0.9324695142031522D+00,
|
|
* 0.3607615730481384D+00,-0.6612093864662647D+00,
|
|
* 0.4679139345726904D+00,-0.2386191860831970D+00/
|
|
! Gauss Legendre Points and Weights, N = 12
|
|
DATA ( W(I,2), X(I,2), I = 1,6) /
|
|
* 0.4717533638651177D-01,-0.9815606342467191D+00,
|
|
* 0.1069393259953183D+00,-0.9041172563704750D+00,
|
|
* 0.1600783285433464D+00,-0.7699026741943050D+00,
|
|
* 0.2031674267230659D+00,-0.5873179542866171D+00,
|
|
* 0.2334925365383547D+00,-0.3678314989981802D+00,
|
|
* 0.2491470458134029D+00,-0.1252334085114692D+00/
|
|
! Gauss Legendre Points and Weights, N = 20
|
|
DATA ( W(I,3), X(I,3), I = 1,10) /
|
|
* 0.1761400713915212D-01,-0.9931285991850949D+00,
|
|
* 0.4060142980038694D-01,-0.9639719272779138D+00,
|
|
* 0.6267204833410906D-01,-0.9122344282513259D+00,
|
|
* 0.8327674157670475D-01,-0.8391169718222188D+00,
|
|
* 0.1019301198172404D+00,-0.7463319064601508D+00,
|
|
* 0.1181945319615184D+00,-0.6360536807265150D+00,
|
|
* 0.1316886384491766D+00,-0.5108670019508271D+00,
|
|
* 0.1420961093183821D+00,-0.3737060887154196D+00,
|
|
* 0.1491729864726037D+00,-0.2277858511416451D+00,
|
|
* 0.1527533871307259D+00,-0.7652652113349733D-01/
|
|
SAVE W, X
|
|
VAL = ZERO
|
|
HK = MIN(SH,SK)
|
|
IF ( HK < -XMAX) THEN ! pab 24.05.2003
|
|
VAL = FI(-MAX(SH,SK))
|
|
RETURN
|
|
ELSE IF ( XMAX < MAX(SH,SK)) THEN
|
|
RETURN
|
|
ENDIF
|
|
IF ( ABS(R) < 0.3D0 ) THEN
|
|
NG = 1
|
|
LG = 3
|
|
ELSE IF ( ABS(R) < 0.75D0 ) THEN
|
|
NG = 2
|
|
LG = 6
|
|
ELSE
|
|
NG = 3
|
|
LG = 10
|
|
ENDIF
|
|
H = SH
|
|
K = SK
|
|
HK = H*K
|
|
|
|
IF ( ABS(R) < 0.925D0 ) THEN
|
|
IF (ABS(R) .GT. ZERO ) THEN
|
|
HS = ( H*H + K*K )*HALF
|
|
ASR = ASIN(R)
|
|
DO I = 1, LG
|
|
SN = SIN(ASR*(ONE + X(I,NG))*HALF)
|
|
VAL = VAL + W(I,NG)*EXP( ( SN*HK - HS )/( ONE - SN*SN ) )
|
|
SN = SIN(ASR*(ONE - X(I,NG))*HALF)
|
|
VAL = VAL + W(I,NG)*EXP( ( SN*HK - HS )/( ONE - SN*SN ) )
|
|
END DO
|
|
VAL = VAL*ASR*FOURPI1
|
|
ENDIF
|
|
VAL = VAL + FI(-H)*FI(-K)
|
|
ELSE
|
|
IF ( R < ZERO ) THEN
|
|
K = -K
|
|
HK = -HK
|
|
ENDIF
|
|
IF ( ABS(R) < ONE ) THEN
|
|
AS = ( ONE - R )*( ONE + R )
|
|
A = SQRT(AS)
|
|
B = ABS( H - K ) !**2
|
|
BS = B * B
|
|
C = ( FOUR - HK ) * ONEEIGHT !/8D0
|
|
D = ( TWELVE - HK ) * ONESIXTEEN !/16D0
|
|
ASR = -(BS/AS + HK)*HALF
|
|
IF (ASR.GT.EXPMIN) THEN
|
|
VAL = A*EXP( ASR ) *
|
|
& ( ONE - C*(BS - AS)*(ONE - D*BS*0.2D0)*ONETHIRD +
|
|
& C*D*AS*AS*0.2D0 )
|
|
ENDIF
|
|
IF ( HK .GT. EXPMIN ) THEN
|
|
VAL = VAL - EXP(-HK*HALF)*SQTWOPI*FI(-B/A)*B
|
|
+ *( ONE - C*BS*( ONE - D*BS*0.2D0 )*ONETHIRD )
|
|
ENDIF
|
|
A = A * HALF
|
|
DO I = 1, LG
|
|
XS = ( A * (ONE + X(I,NG)) ) !**2
|
|
XS = XS * XS
|
|
RS = SQRT( ONE - XS )
|
|
ASR = -(BS / XS + HK) * HALF
|
|
IF (ASR.GT.EXPMIN) THEN
|
|
VAL = VAL + A*W(I,NG)*EXP( ASR )
|
|
& * ( EXP( - HALF*HK*( ONE - RS )/( ONE + RS ) )/RS
|
|
$ -( ONE + C*XS*( ONE + D*XS ) ) )
|
|
ENDIF
|
|
XS = ( A * (ONE - X(I,NG)) ) !**2
|
|
XS = XS * XS
|
|
RS = SQRT( ONE - XS )
|
|
ASR = -(BS / XS + HK) * HALF
|
|
IF (ASR.GT.EXPMIN) THEN
|
|
VAL = VAL + A*W(I,NG)*EXP( ASR )
|
|
& *( EXP( - HALF*HK*( ONE - RS )/( ONE + RS ) )/RS-
|
|
$ ( ONE + C*XS*( ONE + D*XS ) ) )
|
|
ENDIF
|
|
END DO
|
|
VAL = -VAL*TWOPI1
|
|
ENDIF
|
|
IF ( R .GT. ZERO ) THEN
|
|
VAL = VAL + FI( -MAX( H, K ) )
|
|
ELSE
|
|
VAL = -VAL
|
|
IF ( H < K ) VAL = VAL + FI(K)-FI(H)
|
|
ENDIF
|
|
ENDIF
|
|
RETURN
|
|
END FUNCTION BVU
|
|
DOUBLE PRECISION FUNCTION STUDNT( NU, T )
|
|
IMPLICIT NONE
|
|
!
|
|
! Student t Distribution Function
|
|
!
|
|
! T
|
|
! STUDNT = C I ( 1 + y*y/NU )**( -(NU+1)/2 ) dy
|
|
! NU -INF
|
|
!
|
|
INTEGER, INTENT(IN) :: NU
|
|
DOUBLE PRECISION, INTENT(IN) :: T
|
|
! Locals
|
|
INTEGER :: J
|
|
DOUBLE PRECISION :: ZRO, ONE
|
|
PARAMETER ( ZRO = 0.0D0, ONE = 1.0D0 )
|
|
DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979D0
|
|
DOUBLE PRECISION :: CSSTHE, SNTHE, POLYN, TT, TS, RN
|
|
|
|
IF ( NU < 1 ) THEN
|
|
STUDNT = FI( T )
|
|
ELSE IF ( NU .EQ. 1 ) THEN
|
|
STUDNT = ( ONE + 2.0D0*ATAN(T)/PI )*0.5D0
|
|
ELSE IF ( NU .EQ. 2 ) THEN
|
|
STUDNT = ( ONE + T/SQRT( 2.0D0 + T*T ))*0.5D0
|
|
ELSE
|
|
RN = NU ! convert to double
|
|
TT = T * T
|
|
CSSTHE = ONE/( ONE + TT/RN )
|
|
POLYN = 1
|
|
DO J = NU-2, 2, -2
|
|
POLYN = ONE + ( J - 1 )*CSSTHE*POLYN/J
|
|
END DO
|
|
|
|
IF ( MOD( NU, 2 ) .EQ. 1 ) THEN
|
|
TS = T/SQRT(RN)
|
|
STUDNT = ( ONE + 2.0D0*( ATAN(TS) +
|
|
& TS*CSSTHE*POLYN )/PI )*0.5D0
|
|
ELSE
|
|
SNTHE = T/SQRT( RN + TT )
|
|
STUDNT = ( ONE + SNTHE*POLYN )*0.5D0
|
|
END IF
|
|
STUDNT = MAX( ZRO, MIN( STUDNT, ONE ) )
|
|
ENDIF
|
|
END FUNCTION STUDNT
|
|
DOUBLE PRECISION FUNCTION BVTL( NU, DH, DK, R )
|
|
IMPLICIT NONE
|
|
!*
|
|
!* A function for computing bivariate t probabilities.
|
|
!*
|
|
!* Alan Genz
|
|
!* Department of Mathematics
|
|
!* Washington State University
|
|
!* Pullman, WA 99164-3113
|
|
!* Email : alangenz@wsu.edu
|
|
!*
|
|
!* This function is based on the method described by
|
|
!* Dunnett, C.W. and M. Sobel, (1954),
|
|
!* A bivariate generalization of Student's t-distribution
|
|
!* with tables for certain special cases,
|
|
!* Biometrika 41, pp. 153-169.
|
|
!*
|
|
!* BVTL - calculate the probability that X < DH and Y < DK.
|
|
!*
|
|
!* parameters
|
|
!*
|
|
!* NU number of degrees of freedom (NOTE: NU = 0 gives bivariate normal prb)
|
|
!* DH 1st lower integration limit
|
|
!* DK 2nd lower integration limit
|
|
!* R correlation coefficient
|
|
!*
|
|
INTEGER, INTENT(IN) ::NU
|
|
DOUBLE PRECISION, INTENT(IN) :: DH, DK, R
|
|
! Locals
|
|
INTEGER :: J, HS, KS
|
|
DOUBLE PRECISION :: ORS, HRK, KRH, BVT
|
|
DOUBLE PRECISION :: DH2, DK2, SNU ,DNU, DHDK
|
|
!, BVND, STUDNT
|
|
DOUBLE PRECISION :: GMPH, GMPK, XNKH, XNHK, QHRK, HKN, HPK, HKRN
|
|
DOUBLE PRECISION :: BTNCKH, BTNCHK, BTPDKH, BTPDHK
|
|
DOUBLE PRECISION :: ZERO, ONE, EPS, PI,TPI
|
|
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, EPS = 1.0D-15 )
|
|
PARAMETER (PI = 3.14159265358979D0, TPI = 6.28318530717959D0)
|
|
IF ( NU < 1 ) THEN
|
|
BVTL = BVU( -DH, -DK, R )
|
|
ELSE IF ( ONE - R .LE. EPS .OR. 1.0D+16<MAX(ABS(DH),ABS(DK))) THEN
|
|
BVTL = STUDNT( NU, MIN( DH, DK ) )
|
|
ELSE IF ( R + ONE .LE. EPS ) THEN
|
|
IF ( DH .GT. -DK ) THEN
|
|
BVTL = STUDNT( NU, DH ) - STUDNT( NU, -DK )
|
|
ELSE
|
|
BVTL = ZERO
|
|
END IF
|
|
ELSE
|
|
!PI = ACOS(-ONE)
|
|
!TPI = 2*PI
|
|
DNU = NU ! convert to double
|
|
SNU = SQRT(DNU)
|
|
ORS = ONE - R * R
|
|
HRK = DH - R * DK
|
|
KRH = DK - R * DH
|
|
DK2 = DK * DK
|
|
DH2 = DH * DH
|
|
IF ( ABS(HRK) + ORS .GT. ZERO ) THEN
|
|
XNHK = HRK**2/( HRK**2 + ORS*( DNU + DK2 ) )
|
|
XNKH = KRH**2/( KRH**2 + ORS*( DNU + DH2 ) )
|
|
ELSE
|
|
XNHK = ZERO
|
|
XNKH = ZERO
|
|
END IF
|
|
HS = INT(DSIGN( ONE, HRK)) !DH - R*DK )
|
|
KS = INT(DSIGN( ONE, KRH)) !DK - R*DH )
|
|
IF ( MOD( NU, 2 ) .EQ. 0 ) THEN
|
|
BVT = ATAN2( SQRT(ORS), -R )/TPI
|
|
GMPH = DH/SQRT( 16.0D0*( DNU + DH2 ) )
|
|
GMPK = DK/SQRT( 16.0D0*( DNU + DK2 ) )
|
|
BTNCKH = 2*ATAN2( SQRT( XNKH ), SQRT( ONE - XNKH ) )/PI
|
|
BTPDKH = 2*SQRT( XNKH*( ONE - XNKH ) )/PI
|
|
BTNCHK = 2*ATAN2( SQRT( XNHK ), SQRT( ONE - XNHK ) )/PI
|
|
BTPDHK = 2*SQRT( XNHK*( ONE - XNHK ) )/PI
|
|
DO J = 1, NU/2
|
|
BVT = BVT + GMPH*( ONE + KS*BTNCKH )
|
|
BVT = BVT + GMPK*( ONE + HS*BTNCHK )
|
|
BTNCKH = BTNCKH + BTPDKH
|
|
BTPDKH = 2*J*BTPDKH*( ONE - XNKH )/( 2*J + 1 )
|
|
BTNCHK = BTNCHK + BTPDHK
|
|
BTPDHK = 2*J*BTPDHK*( ONE - XNHK )/( 2*J + 1 )
|
|
GMPH = GMPH*( 2*J - 1 )/( 2*J*( ONE + DH2/DNU ) )
|
|
GMPK = GMPK*( 2*J - 1 )/( 2*J*( ONE + DK2/DNU ) )
|
|
END DO
|
|
ELSE ! NU is ODD
|
|
DHDK = DH*DK
|
|
QHRK = SQRT( DH2 + DK2 - 2.0D0*R*DHDK + DNU*ORS )
|
|
HKRN = DHDK + R*DNU
|
|
HKN = DHDK - DNU
|
|
HPK = DH + DK
|
|
BVT = ATAN2( -SNU*( HKN*QHRK + HPK*HKRN ),
|
|
& HKN*HKRN-DNU*HPK*QHRK )/TPI
|
|
IF ( BVT < -EPS ) BVT = BVT + ONE
|
|
GMPH = DH/( TPI*SNU*( ONE + DH2/DNU ) )
|
|
GMPK = DK/( TPI*SNU*( ONE + DK2/DNU ) )
|
|
BTNCKH = SQRT( XNKH )
|
|
BTPDKH = BTNCKH
|
|
BTNCHK = SQRT( XNHK )
|
|
BTPDHK = BTNCHK
|
|
DO J = 1, ( NU - 1 )/2
|
|
BVT = BVT + GMPH*( ONE + KS*BTNCKH )
|
|
BVT = BVT + GMPK*( ONE + HS*BTNCHK )
|
|
BTPDKH = ( 2*J - 1 )*BTPDKH*( ONE - XNKH )/( 2*J )
|
|
BTNCKH = BTNCKH + BTPDKH
|
|
BTPDHK = ( 2*J - 1 )*BTPDHK*( ONE - XNHK )/( 2*J )
|
|
BTNCHK = BTNCHK + BTPDHK
|
|
GMPH = 2*J*GMPH/( ( 2*J + 1 )*( ONE + DH2/DNU ) )
|
|
GMPK = 2*J*GMPK/( ( 2*J + 1 )*( ONE + DK2/DNU ) )
|
|
END DO
|
|
END IF
|
|
BVTL = BVT
|
|
END IF
|
|
END FUNCTION BVTL
|
|
|
|
DOUBLE PRECISION FUNCTION TVTL( NU1, H, R, EPSI )
|
|
USE TRIVARIATEVAR !Block to transfer variables to TVTMFN
|
|
IMPLICIT NONE
|
|
!
|
|
! A function for computing trivariate normal and t-probabilities.
|
|
! This function uses algorithms developed from the ideas
|
|
! described in the papers:
|
|
! R.L. Plackett, Biometrika 41(1954), pp. 351-360.
|
|
! Z. Drezner, Math. Comp. 62(1994), pp. 289-294.
|
|
! with adaptive integration from (0,0,1) to (0,0,r23) to R.
|
|
!
|
|
! Calculate the probability that X(I) < H(I), for I = 1,2,3
|
|
! NU INTEGER degrees of freedom; use NU = 0 for normal cases.
|
|
! H REAL array of uppoer limits for probability distribution
|
|
! R REAL array of three correlation coefficients, R should
|
|
! contain the lower left portion of the correlation matrix r.
|
|
! R should contains the values r21, r31, r23 in that order.
|
|
! EPSI REAL required absolute accuracy; maximum accuracy for most
|
|
! computations is approximately 1D-14
|
|
!
|
|
! The software is based on work described in the paper
|
|
! "Numerical Computation of Rectangular Bivariate and Trivariate
|
|
! Normal and t Probabilities", by the code author:
|
|
!
|
|
! Alan Genz
|
|
! Department of Mathematics
|
|
! Washington State University
|
|
! Pullman, WA 99164-3113
|
|
! Email : alangenz@wsu.edu
|
|
!
|
|
! EXTERNAL TVTMFN
|
|
INTEGER, INTENT(IN) :: NU1
|
|
DOUBLE PRECISION,DIMENSION(:), INTENT(IN) :: H, R
|
|
DOUBLE PRECISION, INTENT(IN) :: EPSI
|
|
!Locals
|
|
DOUBLE PRECISION :: R12, R13, TVT
|
|
DOUBLE PRECISION :: ONE, ZERO, EPS, PT
|
|
! DOUBLE PRECISION RUA, RUB, AR, RUC,
|
|
! BVTL, PHID, ADONET
|
|
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
|
|
PARAMETER ( PT = 1.57079632679489661923132169163975D0 ) !pi/2
|
|
|
|
! COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NU
|
|
EPS = MAX( 1.0D-13, EPSI )
|
|
! PT = ASIN(ONE)
|
|
NU = NU1
|
|
H1 = H(1)
|
|
H2 = H(2)
|
|
H3 = H(3)
|
|
R12 = R(1)
|
|
R13 = R(2)
|
|
R23 = R(3)
|
|
!
|
|
! Sort R's and check for special cases
|
|
!
|
|
IF ( ABS(R12) .GT. ABS(R13) ) THEN
|
|
H2 = H3
|
|
H3 = H(2)
|
|
R12 = R13
|
|
R13 = R(1)
|
|
END IF
|
|
IF ( ABS(R13) .GT. ABS(R23) ) THEN
|
|
H1 = H2
|
|
H2 = H(1)
|
|
R23 = R13
|
|
R13 = R(3)
|
|
END IF
|
|
TVT = 0
|
|
IF ( ABS(H1) + ABS(H2) + ABS(H3) < EPS ) THEN
|
|
TVT = (ONE + (ASIN(R12) + ASIN(R13) + ASIN(R23))/PT )*0.125D0
|
|
ELSE IF ( NU < 1 .AND. ABS(R12) + ABS(R13) < EPS ) THEN
|
|
TVT = FI(H1)*BVTL( NU, H2, H3, R23 )
|
|
ELSE IF ( NU < 1 .AND. ABS(R13) + ABS(R23) < EPS ) THEN
|
|
TVT = FI(H3)*BVTL( NU, H1, H2, R12 )
|
|
ELSE IF ( NU < 1 .AND. ABS(R12) + ABS(R23) < EPS ) THEN
|
|
TVT = FI(H2)*BVTL( NU, H1, H3, R13 )
|
|
ELSE IF ( ONE - R23 < EPS ) THEN
|
|
TVT = BVTL( NU, H1, MIN( H2, H3 ), R12 )
|
|
ELSE IF ( R23 + ONE < EPS ) THEN
|
|
IF ( H2 .GT. -H3 )
|
|
& TVT = BVTL( NU, H1, H2, R12 ) - BVTL( NU, H1, -H3, R12 )
|
|
ELSE
|
|
!
|
|
! Compute singular TVT value
|
|
!
|
|
IF ( NU < 1 ) THEN
|
|
TVT = BVTL( NU, H2, H3, R23 )*FI(H1)
|
|
ELSE IF ( R23 .GE. ZERO ) THEN
|
|
TVT = BVTL( NU, H1, MIN( H2, H3 ), ZERO )
|
|
ELSE IF ( H2 .GT. -H3 ) THEN
|
|
TVT = BVTL( NU, H1, H2, ZERO ) - BVTL( NU, H1, -H3, ZERO )
|
|
END IF
|
|
!
|
|
! Use numerical integration to compute probability
|
|
!
|
|
!
|
|
RUA = ASIN( R12 )
|
|
RUB = ASIN( R13 )
|
|
AR = ASIN( R23)
|
|
RUC = SIGN( PT, AR ) - AR
|
|
TVT = TVT + ADONET(TVTMFN, ZERO, ONE, EPS )/( 4.D0*PT )
|
|
END IF
|
|
TVTL = MAX( ZERO, MIN( TVT, ONE ) )
|
|
END FUNCTION TVTL
|
|
! CONTAINS
|
|
!
|
|
DOUBLE PRECISION FUNCTION TVTMFN( X )
|
|
USE TRIVARIATEVAR
|
|
IMPLICIT NONE
|
|
!
|
|
! Computes Plackett formula integrands
|
|
!
|
|
DOUBLE PRECISION, INTENT(IN) :: X
|
|
! Locals
|
|
DOUBLE PRECISION R12, RR2, R13, RR3, R, RR, ZRO !, PNTGND
|
|
! Parameters transfeered from TRIVARIATEVAR
|
|
! INTEGER :: NU
|
|
! DOUBLE PRECISION :: H1, H2, H3, R23, RUA, RUB, AR, RUC
|
|
PARAMETER ( ZRO = 0.0D0 )
|
|
! COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NU
|
|
TVTMFN = 0.0D0
|
|
CALL SINCS( RUA*X, R12, RR2 )
|
|
CALL SINCS( RUB*X, R13, RR3 )
|
|
IF ( ABS(RUA) .GT. ZRO )
|
|
& TVTMFN = TVTMFN + RUA*PNTGND( NU, H1,H2,H3, R13,R23,R12,RR2 )
|
|
IF ( ABS(RUB) .GT. ZRO )
|
|
& TVTMFN = TVTMFN + RUB*PNTGND( NU, H1,H3,H2, R12,R23,R13,RR3 )
|
|
IF ( NU .GT. 0 ) THEN
|
|
CALL SINCS( AR + RUC*X, R, RR )
|
|
TVTMFN = TVTMFN - RUC*PNTGND( NU, H2, H3, H1, ZRO, ZRO, R, RR )
|
|
END IF
|
|
END FUNCTION TVTMFN
|
|
!
|
|
SUBROUTINE SINCS( X, SX, CS )
|
|
!
|
|
! Computes SIN(X), COS(X)^2, with series approx. for |X| near PI/2
|
|
!
|
|
DOUBLE PRECISION, INTENT(IN) :: X
|
|
DOUBLE PRECISION, INTENT(OUT) :: SX, CS
|
|
!Locals
|
|
DOUBLE PRECISION :: EE, PT, KS, KC, ONE, SMALL, HALF, ONETHIRD
|
|
PARAMETER (ONE = 1.0D0, SMALL = 5.0D-5, HALF = 0.5D0 )
|
|
PARAMETER ( PT = 1.57079632679489661923132169163975D0 )
|
|
PARAMETER ( KS = 0.0833333333333333333333333333333D0) !1/12
|
|
PARAMETER ( KC = 0.1333333333333333333333333333333D0) !2/15
|
|
PARAMETER ( ONETHIRD = 0.33333333333333333333333333333333D0) !1/3
|
|
EE = ( PT - ABS(X) )
|
|
EE = EE * EE
|
|
IF ( EE < SMALL ) THEN
|
|
SX = SIGN( ONE - EE*( ONE - EE*KS )*HALF, X )
|
|
CS = EE *( ONE - EE*( ONE - EE*KC )*ONETHIRD)
|
|
ELSE
|
|
SX = SIN(X)
|
|
CS = ONE - SX*SX
|
|
END IF
|
|
END SUBROUTINE SINCS
|
|
!
|
|
DOUBLE PRECISION FUNCTION PNTGND( NU, BA, BB, BC, RA, RB, R, RR )
|
|
IMPLICIT NONE
|
|
!
|
|
! Computes Plackett formula integrand
|
|
!
|
|
INTEGER, INTENT(IN) :: NU
|
|
DOUBLE PRECISION, INTENT(IN) :: BA, BB, BC, RA, RB, R, RR
|
|
! Locals
|
|
DOUBLE PRECISION :: DT, FT, BT,RAB2, BARB, rNU!, PHID, STUDNT
|
|
PNTGND = 0.0D0
|
|
FT = ( RA - RB )
|
|
RAB2 = FT*FT
|
|
DT = RR*( RR - RAB2 - 2.0D0*RA*RB*( 1.D0 - R ) )
|
|
IF ( DT .GT. 0.0D0 ) THEN
|
|
BT = ( BC*RR + BA*( R*RB - RA ) + BB*( R*RA -RB ) )/SQRT(DT)
|
|
BARB = ( BA - R*BB )
|
|
FT = ( BARB * BARB ) / RR + BB * BB
|
|
IF ( NU < 1 ) THEN
|
|
IF ( BT .GT. -10.0D0 .AND. FT < 100.0D0 ) THEN
|
|
PNTGND = EXP( -FT * 0.5D0 ) * FI( BT )
|
|
! PNTGND = EXP( -FT*0.5D0)
|
|
! IF ( BT < 10.0D0 ) PNTGND = PNTGND * FI(BT)
|
|
END IF
|
|
ELSE
|
|
rNU = NU
|
|
FT = SQRT( 1.0D0 + FT/rNU )
|
|
PNTGND = STUDNT( NU, BT/FT )/FT**NU
|
|
END IF
|
|
END IF
|
|
END FUNCTION PNTGND
|
|
!
|
|
DOUBLE PRECISION FUNCTION ADONET( F, A, B, TOL )
|
|
IMPLICIT NONE
|
|
!
|
|
! One Dimensional Globally Adaptive Integration Function
|
|
!
|
|
! EXTERNAL F
|
|
DOUBLE PRECISION, INTENT(IN) :: A, B, TOL
|
|
INTEGER :: NL, I, IM, IP
|
|
PARAMETER ( NL = 100 )
|
|
DOUBLE PRECISION, DIMENSION(NL) :: EI, AI, BI, FI
|
|
DOUBLE PRECISION :: FIN, ERR !, KRNRDT
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0
|
|
DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0
|
|
DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0
|
|
DOUBLE PRECISION, PARAMETER :: FOUR = 4.0D0
|
|
INTERFACE
|
|
FUNCTION F(Z) RESULT (VAL)
|
|
DOUBLE PRECISION, INTENT(IN) :: Z
|
|
DOUBLE PRECISION :: VAL
|
|
END FUNCTION F
|
|
END INTERFACE
|
|
! COMMON /ABLK/ ERR, IM
|
|
AI(1) = A
|
|
BI(1) = B
|
|
ERR = ONE
|
|
IP = 1
|
|
IM = 1
|
|
DO WHILE ( (ERR .GT. TOL) .AND. (IM < NL) )
|
|
IM = IM + 1
|
|
BI(IM) = BI(IP)
|
|
AI(IM) = ( AI(IP) + BI(IP) ) * HALF
|
|
BI(IP) = AI(IM)
|
|
FI(IP) = KRNRDT( AI(IP), BI(IP), F, EI(IP) )
|
|
FI(IM) = KRNRDT( AI(IM), BI(IM), F, EI(IM) )
|
|
ERR = ZERO
|
|
FIN = ZERO
|
|
DO I = 1, IM
|
|
IF ( EI(I) .GT. EI(IP) ) IP = I
|
|
FIN = FIN + FI(I)
|
|
ERR = ERR + EI(I)*EI(I)
|
|
END DO
|
|
ERR = FOUR * SQRT( ERR )
|
|
END DO
|
|
ADONET = FIN
|
|
END FUNCTION ADONET
|
|
!
|
|
DOUBLE PRECISION FUNCTION KRNRDT( A, B, F, ERR )
|
|
!
|
|
! Kronrod Rule
|
|
!
|
|
DOUBLE PRECISION, intent(in) :: A, B
|
|
DOUBLE PRECISION, intent(out) :: ERR
|
|
DOUBLE PRECISION T, CEN, FC, WID, RESG, RESK
|
|
!
|
|
! The abscissae and weights are given for the interval (-1,1);
|
|
! only positive abscissae and corresponding weights are given.
|
|
!
|
|
! XGK - abscissae of the 2N+1-point Kronrod rule:
|
|
! XGK(2), XGK(4), ... N-point Gauss rule abscissae;
|
|
! XGK(1), XGK(3), ... optimally added abscissae.
|
|
! WGK - weights of the 2N+1-point Kronrod rule.
|
|
! WG - weights of the N-point Gauss rule.
|
|
!
|
|
INTEGER :: J, N
|
|
PARAMETER ( N = 11 )
|
|
DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0
|
|
DOUBLE PRECISION WG(0:(N+1)/2), WGK(0:N), XGK(0:N)
|
|
SAVE WG, WGK, XGK
|
|
INTERFACE
|
|
FUNCTION F(Z) RESULT (VAL)
|
|
DOUBLE PRECISION, INTENT(IN) :: Z
|
|
DOUBLE PRECISION :: VAL
|
|
END FUNCTION F
|
|
END INTERFACE
|
|
DATA WG( 0)/ 0.2729250867779007D+00/
|
|
DATA WG( 1)/ 0.5566856711617449D-01/
|
|
DATA WG( 2)/ 0.1255803694649048D+00/
|
|
DATA WG( 3)/ 0.1862902109277352D+00/
|
|
DATA WG( 4)/ 0.2331937645919914D+00/
|
|
DATA WG( 5)/ 0.2628045445102478D+00/
|
|
!
|
|
DATA XGK( 0)/ 0.0000000000000000D+00/
|
|
DATA XGK( 1)/ 0.9963696138895427D+00/
|
|
DATA XGK( 2)/ 0.9782286581460570D+00/
|
|
DATA XGK( 3)/ 0.9416771085780681D+00/
|
|
DATA XGK( 4)/ 0.8870625997680953D+00/
|
|
DATA XGK( 5)/ 0.8160574566562211D+00/
|
|
DATA XGK( 6)/ 0.7301520055740492D+00/
|
|
DATA XGK( 7)/ 0.6305995201619651D+00/
|
|
DATA XGK( 8)/ 0.5190961292068118D+00/
|
|
DATA XGK( 9)/ 0.3979441409523776D+00/
|
|
DATA XGK(10)/ 0.2695431559523450D+00/
|
|
DATA XGK(11)/ 0.1361130007993617D+00/
|
|
!
|
|
DATA WGK( 0)/ 0.1365777947111183D+00/
|
|
DATA WGK( 1)/ 0.9765441045961290D-02/
|
|
DATA WGK( 2)/ 0.2715655468210443D-01/
|
|
DATA WGK( 3)/ 0.4582937856442671D-01/
|
|
DATA WGK( 4)/ 0.6309742475037484D-01/
|
|
DATA WGK( 5)/ 0.7866457193222764D-01/
|
|
DATA WGK( 6)/ 0.9295309859690074D-01/
|
|
DATA WGK( 7)/ 0.1058720744813894D+00/
|
|
DATA WGK( 8)/ 0.1167395024610472D+00/
|
|
DATA WGK( 9)/ 0.1251587991003195D+00/
|
|
DATA WGK(10)/ 0.1312806842298057D+00/
|
|
DATA WGK(11)/ 0.1351935727998845D+00/
|
|
!
|
|
! Major variables
|
|
!
|
|
! CEN - mid point of the interval
|
|
! WID - half-length of the interval
|
|
! RESG - result of the N-point Gauss formula
|
|
! RESK - result of the 2N+1-point Kronrod formula
|
|
!
|
|
! Compute the 2N+1-point Kronrod approximation to
|
|
! the integral, and estimate the absolute error.
|
|
!
|
|
WID = ( B - A ) * HALF
|
|
CEN = ( B + A ) * HALF
|
|
FC = F(CEN)
|
|
RESG = FC * WG(0)
|
|
RESK = FC * WGK(0)
|
|
DO J = 1, N
|
|
T = WID * XGK(J)
|
|
FC = F( CEN - T ) + F( CEN + T )
|
|
RESK = RESK + WGK(J) * FC
|
|
IF( MOD( J, 2 ) .EQ. 0 ) RESG = RESG + WG(J/2) * FC
|
|
END DO
|
|
KRNRDT = WID * RESK
|
|
ERR = ABS( WID * ( RESK - RESG ) )
|
|
END FUNCTION KRNRDT
|
|
! END FUNCTION TVTL
|
|
|
|
FUNCTION GAUSINT (X1, X2, A, B, C, D) RESULT (value)
|
|
! USE GLOBALDATA,ONLY: xCutOff
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: X1,X2,A,B,C,D
|
|
DOUBLE PRECISION :: value
|
|
!Local variables
|
|
DOUBLE PRECISION :: Y1,Y2,Y3
|
|
DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
|
|
DOUBLE PRECISION, PARAMETER :: XMAX = 37.d0
|
|
! Let X be standardized Gaussian variable,
|
|
! i.e., X=N(0,1). The function calculate the
|
|
! following integral E[I(X1<X<X2)(A+BX)(C+DX)
|
|
! where I(X1<X<X2) is an indicator function of
|
|
! the set {X1<X<X2}.
|
|
IF (X1.GE.X2) THEN
|
|
value = 0.d0
|
|
RETURN
|
|
ENDIF
|
|
IF (ABS (X1) .GT.XMAX) THEN
|
|
Y1 = 0.d0
|
|
ELSE
|
|
Y1 = (A * D+B * C + X1 * B * D) * EXP ( - 0.5d0 * X1 * X1)
|
|
ENDIF
|
|
IF (ABS (X2) .GT.XMAX) THEN
|
|
Y2 = 0.d0
|
|
ELSE
|
|
Y2 = (A * D+B * C + X2 * B * D) * EXP ( - 0.5d0 * X2 * X2)
|
|
ENDIF
|
|
Y3 = (A * C + B * D) * (FI (X2) - FI (X1) )
|
|
value = Y3 + SQTWOPI1 * (Y1 - Y2)
|
|
RETURN
|
|
END FUNCTION GAUSINT
|
|
|
|
|
|
FUNCTION GAUSINT2 (X1, X2, A, B) RESULT (value)
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: X1,X2,A,B
|
|
DOUBLE PRECISION :: value
|
|
! Local variables
|
|
DOUBLE PRECISION :: X0,Y0,Y1,Y2
|
|
DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
|
|
! Let X be standardized Gaussian variable,
|
|
! i.e., X=N(0,1). The function calculate the
|
|
! following integral E[I(X1<X<X2)ABS(A+BX)]
|
|
! where I(X1<X<X2) is an indicator function of
|
|
! the set {X1<X<X2}.
|
|
IF (X1.GE.X2) THEN
|
|
value = 0.d0
|
|
RETURN
|
|
ENDIF
|
|
IF (ABS(B).EQ.0.d0) THEN
|
|
value = ABS(A)*(FI(X2)-FI(X1))
|
|
RETURN
|
|
ENDIF
|
|
|
|
Y1 = -A*FI(X1)+SQTWOPI1*B*EXP(-0.5d0*X1*X1)
|
|
Y2 = A*FI(X2)-SQTWOPI1*B*EXP(-0.5d0*X2*X2)
|
|
IF ((B*X1 < -A).AND.(-A < B*X2))THEN
|
|
X0 = -A/B
|
|
Y0 = 2.d0*(A*FI(X0)-SQTWOPI1*B*EXP(-0.5d0*X0*X0))
|
|
value=ABS(Y2-Y1-Y0)
|
|
ELSE
|
|
value=ABS(Y1+Y2)
|
|
ENDIF
|
|
RETURN
|
|
END FUNCTION GAUSINT2
|
|
|
|
SUBROUTINE EXLMS(A, X1, X2, INFIN, LOWER, UPPER, Ca,Pa)
|
|
IMPLICIT NONE
|
|
DOUBLE PRECISION, INTENT(in) :: A, X1, X2
|
|
DOUBLE PRECISION, INTENT(out) :: LOWER, UPPER,Ca,Pa
|
|
INTEGER,INTENT(in) :: INFIN
|
|
! Local variables
|
|
DOUBLE PRECISION :: P1
|
|
! DOUBLE PRECISION, PARAMETER :: AMAX = 5D15
|
|
|
|
! Let X be standardized Gaussian variable,
|
|
! i.e., X=N(0,1). The function calculate the
|
|
! following integral E[I(X1<X<X2)ABS(A+X)]/C(A)
|
|
! where I(X1<X<X2) is an indicator function of
|
|
! the set {X1<X<X2} and C(A) is a normalization factor
|
|
! i.e. E(I(-inf<X<inf)*ABS(A+X)) = C(A)
|
|
!
|
|
! Pa = probability at the inflection point of the CDF
|
|
! Ca = C(A) normalization parameter
|
|
! A = location parameter of the inflection point
|
|
! X1,X2 = integration limits < infinity
|
|
LOWER = 0.d0
|
|
UPPER = 1.d0
|
|
|
|
P1 = EXFUN(-A,A)
|
|
Ca = A+2D0*P1
|
|
Pa = P1/Ca
|
|
IF ( INFIN < 0 ) RETURN
|
|
IF ( INFIN .NE. 0 ) THEN
|
|
IF (X1.LE.-A) THEN
|
|
LOWER = EXFUN(X1,A)/Ca
|
|
ELSE
|
|
LOWER = 1D0-EXFUN(-X1,-A)/Ca
|
|
ENDIF
|
|
ENDIF
|
|
IF ( INFIN .NE. 1 ) THEN
|
|
IF (X2.LE.-A) THEN
|
|
UPPER = EXFUN(X2,A)/Ca
|
|
ELSE
|
|
UPPER = 1D0-EXFUN(-X2,-A)/Ca
|
|
ENDIF
|
|
ENDIF
|
|
|
|
RETURN
|
|
CONTAINS
|
|
FUNCTION EXFUN(X,A) RESULT (P)
|
|
DOUBLE PRECISION, INTENT(in) :: X,A
|
|
DOUBLE PRECISION :: P
|
|
DOUBLE PRECISION, PARAMETER :: SQTWOPI1=3.9894228040143d-1 !=1/sqrt(2*pi)
|
|
P = EXP(-X*X*0.5D0)*SQTWOPI1-A*FI(X)
|
|
END FUNCTION EXFUN
|
|
END SUBROUTINE EXLMS
|
|
|
|
|
|
FUNCTION EXINV(P,A,Ca,Pa) RESULT (VAL)
|
|
!EXINV calculates the inverse of the CDF of abs(x+A)*exp(-x^2/2)/SQRT(2*pi)/C(A)
|
|
! where C(A) is a normalization parameter depending on A
|
|
!
|
|
! CALL: val = exinv(P,A,Ca,Pa)
|
|
!
|
|
! val = quantiles
|
|
! p = probabilites
|
|
! A = location parameter
|
|
! Ca = normalization parameter
|
|
! Pa = excdf(-A,A) probability at the inflection point of the CDF
|
|
double precision, intent(in) :: P,A,Ca,Pa
|
|
double precision :: val
|
|
! local variables
|
|
double precision, parameter :: amax = 5.D15
|
|
double precision, parameter :: epsl = 5.D-15
|
|
double precision, parameter :: xmax = 8
|
|
double precision :: P1,Xk,Ak,Zk,SGN
|
|
|
|
! if (P<0.D0.OR.P.GT.1.D0) PRINT *,'warning P<0 or P>1'
|
|
|
|
! The inverse cdf of 0 is -inf, and the inverse cdf of 1 is inf.
|
|
if (P.LE.EPSL.OR.P+EPSL.GE.1.D0) THEN
|
|
VAL = SIGN(xmax,P-0.5D0)
|
|
return
|
|
endif
|
|
Ak = ABS(A)
|
|
if (EPSL < Ak .AND. Ak < amax) THEN
|
|
IF (ABS(p-Pa).LE.EPSL) THEN
|
|
VAL = SIGN(MIN(Ak,xmax),-A)
|
|
RETURN
|
|
ENDIF
|
|
IF (Ak < 1D-2) THEN ! starting guess always less than 0.2 from the true value
|
|
IF (P.GE.0.5D0) THEN
|
|
xk = SQRT(-2D0*log(2D0*(1D0-P)))
|
|
ELSE
|
|
xk = -SQRT(-2D0*log(2D0*P))
|
|
ENDIF
|
|
ELSE
|
|
xk = FIINV(P) ! starting guess always less than 0.8 from the true value
|
|
! Modify starting guess if possible in order to speed up Newtons method
|
|
IF (1D-3.LE.P.AND. P.LE.0.99D0.AND.
|
|
& 3.5.LE.Ak.AND.Ak.LE.1D3 ) THEN
|
|
SGN = SIGN(1.d0,-A)
|
|
Zk = xk*SGN
|
|
xk = SGN*(Zk+((1D0/(64.9495D0*Ak-178.3191D0)-0.02D0/Ak)*
|
|
& Zk+1D0/(-0.99679234298211D0*Ak-0.07195350071872D0))/
|
|
& (Zk/(-1.48430620263825D0*Ak-0.33340759016175D0)+1D0))
|
|
ELSEIF ((P < 1D-3.AND.A.LE.-3.5D0).OR.
|
|
& (3.5D0.LE.A.AND.P.GT.0.99D0)) THEN
|
|
SGN = SIGN(1.d0,-A)
|
|
Zk = xk*SGN
|
|
P1 = -2.00126182192701D0*Ak-2.57306603933111D0
|
|
xk = SGN*Zk*(1D0+
|
|
& P1/((-0.99179258785909D0*Ak-0.21359746002397D0)*
|
|
& (Zk+P1)))
|
|
ENDIF
|
|
ENDIF
|
|
! Check if the starting guess is on the correct side of the inflection point
|
|
IF (xk.LE.-A .AND. P.GT.Pa) xk = 1.D-2-A
|
|
IF (xk.GE.-A .AND. P < Pa) xk = -1.D-2-A
|
|
|
|
|
|
IF (P < Pa) THEN
|
|
VAL = funca(xk,A,P*Ca)
|
|
ELSE ! exploit the symmetry of the CDF
|
|
VAL = -funca(-xk,-A,(1.D0-P)*Ca)
|
|
ENDIF
|
|
|
|
ELSEIF (ABS(A).LE.EPSL) THEN
|
|
IF (P>=0.5D0) THEN
|
|
VAL = SQRT(-2D0*log(2D0*(1.D0-P)))
|
|
ELSE
|
|
VAL = -SQRT(-2D0*log(2D0*P))
|
|
ENDIF
|
|
ELSE ! ABS(A) > AMAX
|
|
VAL = FIINV(P)
|
|
ENDIF
|
|
!CALL EXLMS(A,0.d0,VAL,0,ak,P1,zk,sgn)
|
|
!If (ABS(p-P1).GT.0.0001) PRINT *,'excdf(x,a)-p',p-P1
|
|
RETURN
|
|
|
|
CONTAINS
|
|
|
|
function funca(xk0,ak,CaP) RESULT (xk)
|
|
double precision, intent(in) :: xk0,ak,CaP ! =Ca*P
|
|
DOUBLE PRECISION :: xk
|
|
!Local variables
|
|
INTEGER, PARAMETER :: ixmax = 25
|
|
double precision, parameter :: crit = 7.1D-08 ! = sqrt(1e-15)
|
|
double precision, parameter :: SQTWOPI1 = 0.39894228040143D0 !=1/SQRT(2*pi)
|
|
double precision, parameter :: SQTWOPI = 2.50662827463100D0 !=SQRT(2*pi)
|
|
INTEGER :: IX
|
|
DOUBLE PRECISION :: H,H1,tmp0,tmp1,XNEW
|
|
! Newton's Method or Fixed point iteration to find the inverse of the EXCDF.
|
|
! Assumption: xk0 < -ak and xk < -ak
|
|
! Permit no more than IXMAX iterations.
|
|
IX = 0
|
|
H = 1.D0
|
|
xk = xk0 ! starting guess for the iteration
|
|
|
|
|
|
! Break out of the iteration loop for the following:
|
|
! 1) The last update is very small (compared to x).
|
|
! 2) The last update is very small (compared to sqrt(eps)=crit).
|
|
! 3) There are more than 15 iterations. This should NEVER happen.
|
|
IF (.TRUE..OR.ABS(ak) < 1.D-2) THEN
|
|
! Newton's method
|
|
!~~~~~~~~~~~~~~~~~
|
|
DO WHILE( ABS(H).GT.MIN(crit*ABS(xk),crit).AND.IX < IXMAX)
|
|
|
|
IX = IX+1
|
|
!print *,'Iteration ',IX
|
|
|
|
tmp0 = FI(xk)
|
|
tmp1 = EXP(-xk*xk*0.5D0)*SQTWOPI1 ! =normpdf(x)
|
|
H1 = (tmp1-ak*tmp0-CaP)/(ABS(xk+ak)*tmp1)
|
|
H = DSIGN(MIN(ABS(H1),0.7D0/DBLE(IX)),H1) ! Only allow smaller and smaller steps
|
|
|
|
xnew = xk - H
|
|
! Make sure that the current guess is less than -a.
|
|
! When Newton's Method suggests steps that lead to -a guesses
|
|
! take a step 9/10ths of the way to -a:
|
|
IF (xnew.GT.-ak-crit) THEN
|
|
xnew = (xk - 9.D0*ak)*1D-1
|
|
H = xnew - xk
|
|
ENDIF
|
|
xk = xnew
|
|
END DO
|
|
ELSE ! FIXED POINT iteration
|
|
!~~~~~~~~~~~~~~~~~~~~~~~
|
|
DO WHILE (ABS(H).GT.MIN(crit*ABS(xk),crit).AND.IX < IXMAX)
|
|
IX = IX+1
|
|
tmp0 = SQTWOPI1*EXP(-xk*xk*0.5D0)/FI(xk)
|
|
tmp1 = -2.D0*LOG(SQTWOPI*CaP*tmp0/(tmp0-ak))
|
|
SGN = sign(1.D0,tmp1)
|
|
xnew = -SQRT(SGN*tmp1)*SGN
|
|
! Make sure that the current guess is less than -a.
|
|
! When this method suggests steps that lead to -a guesses
|
|
! take a step 9/10ths of the way to -a:
|
|
IF (xnew.GT.-ak-crit) xnew = (xk - 9.D0*ak)*1.D-1
|
|
|
|
H = xnew - xk
|
|
xk = xnew
|
|
END DO
|
|
ENDIF
|
|
|
|
!print *,'EXINV total number of iterations ',IX
|
|
if (IX.GE.IXMAX) THEN
|
|
! print *, 'Warning: EXINV did not converge. Cap=',Cap
|
|
! print *, 'The last step was: ', h, ' value=,',xk,' ak=',ak
|
|
endif
|
|
return
|
|
END FUNCTION FUNCA
|
|
END FUNCTION EXINV
|
|
END MODULE FIMOD
|