Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
3690 lines (3684 sloc) 150 KB
MODULE PRECISION_MODEL
IMPLICIT NONE
INTEGER, PARAMETER, PUBLIC :: STND = SELECTED_REAL_KIND(12, 60)
END MODULE PRECISION_MODEL
!
MODULE MVSTAT
!
!
! This package contains Fortran 90 software for the MVT distribution,
! plus supporting software. The main subroutine is MVDIST.
! The package is self-contained and should compile without errors on
! standard Fortran(90) compilers.
!
USE PRECISION_MODEL
IMPLICIT NONE
PUBLIC :: MVDIST, MVDSTD, MVCRIT, LWRUPR, LWUPRH, LWUPRT, COVSCL
! MVT, MVT+Deriv, Critical, Bounds, Bounds,
! MVDIST: for computing non-central multivariate t probabilities.
! MVDSTD: for computing a multivariate t probability and derivative.
! MVCRIT: for computing critical values. The subroutine tries to
! compute the critical value TALPHA satisfying P(TALPHA) = 1 - ALPHA.
! P(T) is a multivariate T probability integral, with covariance
! matrix COVRNC, for an integration region determined by
! TALPHA*LOWER < DELTA + CONSTR*X < TALPHA*UPPER .
! LWUPRH: for computing Second Order Bounds
! LWUPRH: for computing Hybrid Second-Third Order Bounds
! LWUPRT: for computing Third Order Bounds
!
PUBLIC :: MVKBRV, UNIFRM
! KOROBOV, Uniform
!
PRIVATE :: MVINIT, MVINTD, MVSORT, MVSWAP, MVFUNC, MVLIMS
PRIVATE :: MVKRSV, MVSPCL, MVVRNC, MVCHNT
!
INTEGER, PRIVATE, PARAMETER :: ML = 500, NL = 200
REAL(KIND=STND), PRIVATE :: SQTNU
REAL(KIND=STND), PRIVATE, DIMENSION(ML,NL) :: CNSTRN
REAL(KIND=STND), PRIVATE, DIMENSION(ML) :: A, B, D
INTEGER, PRIVATE, DIMENSION(ML) :: INF
INTEGER, PRIVATE, DIMENSION(ML) :: CONLNG
INTEGER, PRIVATE :: NP, MP, NUP
!
PUBLIC :: MVDNST, MVUVT, MVSTDT, MVSTNV
! T-Dens, 1-D Dist, T-Dist, Inv T
PUBLIC :: MVPHI, MVPHNV, MVCHNV, MVCHI
! Phi, Inv Phi, Inv Chi, Chi
PUBLIC :: MVBVT, MVBVTL, MVTVT, MVTVTL, MVBVU
! BVT, TVT, BVNU
PUBLIC :: MVSTDC, MVBVTC, MVTVTC
! Comp: 1-D, 2-D, 3-D
PUBLIC :: MVADON, MVKRND
! 1-D Adaptive, Kronrod Rule
PRIVATE :: TVTFNC, MVCHNC
INTEGER, PRIVATE :: NUF
REAL(KIND=STND), PRIVATE :: BF1, BF2, BF3, RF31, RF32
CONTAINS
SUBROUTINE MVDIST( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DELTA, &
MAXPTS, ABSEPS, RELEPS, ERROR, VALUE, NEVALS, INFORM )
!
! A subroutine for computing non-central multivariate t probabilities.
! This subroutine uses an algorithm (QRSVN) described in the paper
! "Methods for the Computation of Multivariate t-Probabilities",
! by Alan Genz and Frank Bretz
!
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, WA 99164-3113
! Email : AlanGenz@wsu.edu
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL NxN covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! LOWER < DELTA + CONSTR*X < UPPER.
! UPPER REAL, array of upper integration limits.
! INFIN INTEGER, array of integration limits flags:
! if INFIN(I) < 0, Ith limits are (-infinity, infinity);
! 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)].
! DELTA REAL array of non-centrality parameters.
! MAXPTS INTEGER, maximum number of function values allowed. This
! parameter can be used to limit the time. A sensible
! strategy is to start with MAXPTS = 1000*N, and then
! increase MAXPTS if ERROR is too large.
! ABSEPS REAL absolute error tolerance.
! RELEPS REAL relative error tolerance.
! ERROR REAL estimated absolute error, with 99% confidence level.
! VALUE REAL estimated value for the integral
! NEVALS Integer number of integrand values used for this call.
! INFORM INTEGER, termination status parameter:
! if INFORM = 0, normal completion with ERROR < EPS;
! if INFORM = 1, completion with ERROR > EPS and MAXPTS
! function vaules used; increase MAXPTS to
! decrease ERROR;
! if INFORM = 2, N > NL or N < 1.
! if INFORM = 3, covariance matrix not positive semidefinite.
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, NU, M, MAXPTS
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(:,:), INTENT(IN) :: COVRNC, CONSTR
INTEGER, DIMENSION(:), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: RELEPS, ABSEPS
INTEGER, INTENT(OUT) :: INFORM, NEVALS
REAL(KIND=STND), INTENT(OUT) :: ERROR, VALUE
!
!
!
REAL(KIND=STND), DIMENSION(1) :: V
!
NEVALS = 0
CALL MVINIT( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DELTA, &
VALUE, ERROR, INFORM )
IF ( INFORM == 0 .AND. NP > 0 ) THEN
CALL MVKBRV( NP - 1 + MIN(1,NUP), 0, MAXPTS, 1, MVFUNC, &
ABSEPS, RELEPS, ERROR, V, NEVALS, INFORM )
VALUE = V(1)
ENDIF
!
END SUBROUTINE MVDIST
!
FUNCTION MVFUNC( NF, W ) RESULT(VALUE)
!
! Integrand subroutine
!
!
! Global Variables
!
INTEGER, INTENT(IN) :: NF
REAL(KIND=STND), INTENT(IN), DIMENSION(:) :: W
REAL(KIND=STND), DIMENSION(NF) :: VALUE
!
! Local Variables
!
REAL(KIND=STND), DIMENSION(NP) :: Y
INTEGER :: I, J, L, LI, FINITA, FINITB
REAL(KIND=STND) :: R, S, AI, AT, BI, BT, DI, EI, DE
!
VALUE(1:NF) = 0
IF ( NUP > 0 ) THEN
CALL MVCHNT( NUP, W(NP+NF-1), R, VALUE(1) )
R = R/SQTNU
ELSE
VALUE(1) = 1
R = 1
END IF
LI = 0
DO I = 1, NP
FINITA = 0
FINITB = 0
DO L = LI + 1, LI + CONLNG(I)
S = D(L) + SUM( CNSTRN(L,1:I-1)*Y(1:I-1) )
IF ( INF(L) /= 0 ) THEN
AI = R*A(L) - S
IF ( FINITA > 0 ) AI = MAX( AI, AT )
FINITA = 1
AT = AI
END IF
IF ( INF(L) /= 1 ) THEN
BI = R*B(L) - S
IF ( FINITB > 0 ) BI = MIN( BI, BT )
FINITB = 1
BT = BI
END IF
END DO
CALL MVLIMS( 0, AI, BI, FINITA + FINITA + FINITB - 1, DI, EI, DE )
VALUE(1) = VALUE(1)*DE
IF ( DE == 0 ) EXIT
IF ( I < NP .OR. NF == 2 ) Y(I) = MVPHNV( DI + W(I)*DE )
LI = LI + CONLNG(I)
END DO
IF ( VALUE(1) > 0 .AND. NF == 2 ) VALUE(2) = VALUE(1)*( NP - SUM(Y*Y) )
!
END FUNCTION MVFUNC
!
!
SUBROUTINE MVCHNT( N, P, R, S )
!
! R
! P = 1 - K I exp(-t*t/2) t**(N-1) dt, for N >= 1.
! N 0
!
! with Normal approximation when N >= NR
!
INTEGER, INTENT(IN) :: N
REAL(KIND=STND), INTENT(IN) :: P
REAL(KIND=STND), INTENT(OUT) :: R
REAL(KIND=STND), INTENT(OUT) :: S ! Integrand scale factor for N >= 100
!
INTEGER, PARAMETER :: NR = 1000
INTEGER :: I
!
REAL(KIND=STND), PARAMETER :: LRP = -.22579135264472743235E0_STND
! LRP = LOG( SQRT( 2/PI ) )
REAL(KIND=STND), PARAMETER :: LRT = .91893853320467274177E0_STND
! LRT = LOG( SQRT( 2*PI ) )
REAL(KIND=STND) :: RI
REAL(KIND=STND), SAVE :: LKN
INTEGER, SAVE :: NO = 0
!
IF ( N < NR ) THEN
R = MVCHNV( N, P )
S = 1
ELSE
!
! Use Normal approximation for Chi
!
IF ( N /= NO ) THEN
!
! First call computes LKN = LOG(K )
! N
NO = N
LKN = 0
DO I = N - 2, 2, -2
RI = I
LKN = LKN - LOG( RI )
END DO
IF ( MODULO( N, 2 ) == 1 ) LKN = LKN + LRP
END IF
RI = N - 1
RI = SQRT(RI)
R = RI - MVPHNV(P)
S = EXP( LRT + LKN + RI*( RI*LOG(R) - R + RI/2 ) )
!
! S = K R**(N-1)*EXP( -( R*R - (R-SQRT(N-1))**2 )/2 )*SQRT(2*PI)
! N
END IF
!
END SUBROUTINE MVCHNT
!
SUBROUTINE MVINIT( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DELTA, &
VALUE, ERROR, INFORM )
!
! Initialization and computation of covariance Cholesky factor.
!
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, NU, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
INTEGER, DIMENSION(M), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(OUT) :: VALUE, ERROR
INTEGER, INTENT(OUT) :: INFORM
!
! Initialization and computation of covariance Cholesky factor.
!
!
IF ( 0 < N .AND. N <= NL ) THEN
NUP = MAX( NU, 0 )
CALL MVSORT( N, COVRNC, M, LOWER, CONSTR, UPPER, INFIN, DELTA, &
.TRUE., INFORM )
ELSE
INFORM = 2
END IF
CALL MVSPCL( NU, VALUE, ERROR, INFORM )
!
END SUBROUTINE MVINIT
!
SUBROUTINE MVSPCL( NU, VALUE, ERROR, INFORM )
!
! Special cases subroutine
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: NU, INFORM
REAL(KIND=STND), INTENT(OUT) :: VALUE, ERROR
!
! Local Variables
!
REAL(KIND=STND) :: DS, R22, R33
!
IF ( INFORM > 0 ) THEN
VALUE = 0
ERROR = 1
ELSE
SQTNU = SQRT( REAL( MAX( NUP, 0 ), STND ) )
!
! Special cases
!
IF ( NP == 0 ) THEN
ERROR = 0
VALUE = 1
ELSE IF ( NP < 4 ) THEN
DS = SUM ( ABS( D( 1:CONLNG(1) ) ) )
IF ( NP == 1 .AND. ( NUP < 1 .OR. DS == 0 ) ) THEN
!
! 1-d case for normal or central t
!
VALUE = MVUVT( NUP, A(1)-D(1), B(1)-D(1), INF(1) )
ERROR = 2E-16_STND
NP = 0
ELSE IF ( NP == 2 .AND. CONLNG(1) == 1 .AND. CONLNG(2) < 3 ) THEN
IF ( NUP < 1 .OR. SUM( ABS(D(1:1+CONLNG(2))) ) == 0 ) THEN
IF ( INF(1) /= 0 ) A(1) = A(1) - D(1)
IF ( INF(1) /= 1 ) B(1) = B(1) - D(1)
R22 = SQRT( 1 + CNSTRN(2,1)**2 )
IF ( INF(2) /= 0 ) A(2) = ( A(2) - D(2) )/R22
IF ( INF(2) /= 1 ) B(2) = ( B(2) - D(2) )/R22
IF ( CONLNG(2) == 1 ) THEN
!
! 2-d case for normal or central t
!
CNSTRN(2,1) = CNSTRN(2,1)/R22
VALUE = MVBVT( NUP, A, B, INF, CNSTRN(2,1))
ERROR = 1E-15_STND
ELSE
!
! 3-d singular case for normal or central t
!
R33 = SQRT( 1 + CNSTRN(3,1)**2 )
IF ( INF(3) /= 0 ) A(3) = ( A(3) - D(3) )/R33
IF ( INF(3) /= 1 ) B(3) = ( B(3) - D(3) )/R33
CNSTRN(3,2) = CNSTRN(3,2) + CNSTRN(2,1)*CNSTRN(3,1)
CNSTRN(2,1) = CNSTRN(2,1)/R22
CNSTRN(3,1) = CNSTRN(3,1)/R33
CNSTRN(3,2) = CNSTRN(3,2)/( R33*R22 )
VALUE = MVTVT( NUP, A, B, INF, &
(/CNSTRN(2,1),CNSTRN(3,1),CNSTRN(3,2)/) )
ERROR = 1E-14_STND
END IF
NP = 0
END IF
ELSE IF ( NP == 3 .AND. CONLNG(1) == 1 .AND. CONLNG(2) == 1 &
.AND. CONLNG(3) == 1 ) THEN
IF ( NUP < 1 .OR. SUM( ABS(D(1:3)) ) == 0 ) THEN
!
! 3-d case for normal or central t
!
IF ( INF(1) /= 0 ) A(1) = A(1) - D(1)
IF ( INF(1) /= 1 ) B(1) = B(1) - D(1)
R22 = SQRT( 1 + CNSTRN(2,1)**2 )
IF ( INF(2) /= 0 ) A(2) = ( A(2) - D(2) )/R22
IF ( INF(2) /= 1 ) B(2) = ( B(2) - D(2) )/R22
R33 = SQRT( 1 + CNSTRN(3,1)**2 + CNSTRN(3,2)**2 )
IF ( INF(3) /= 0 ) A(3) = ( A(3) - D(3) )/R33
IF ( INF(3) /= 1 ) B(3) = ( B(3) - D(3) )/R33
CNSTRN(3,2) = CNSTRN(3,2) + CNSTRN(2,1)*CNSTRN(3,1)
CNSTRN(2,1) = CNSTRN(2,1)/R22
CNSTRN(3,1) = CNSTRN(3,1)/R33
CNSTRN(3,2) = CNSTRN(3,2)/( R33*R22 )
VALUE = MVTVT( NUP, A, B, INF, &
(/CNSTRN(2,1),CNSTRN(3,1),CNSTRN(3,2)/) )
ERROR = 1E-14_STND
NP = 0
END IF
END IF
END IF
END IF
!
END SUBROUTINE MVSPCL
!
SUBROUTINE MVLIMS( NU, A, B, INFIN, LOWER, UPPER, VALUE )
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: NU, INFIN
REAL(KIND=STND), INTENT(IN) :: A, B
REAL(KIND=STND), INTENT(OUT) :: LOWER, UPPER, VALUE
!
! Local Variables
!
LOWER = 0
UPPER = 1
IF ( INFIN >= 0 ) THEN
IF ( INFIN /= 0 ) LOWER = MVSTDT( NU, A )
IF ( INFIN /= 1 ) UPPER = MVSTDT( NU, B )
END IF
UPPER = MAX( UPPER, LOWER )
VALUE = UPPER - LOWER
!
END SUBROUTINE MVLIMS
!
SUBROUTINE MVDSTD( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, T, &
MAXPTS, ABSEPS, RELEPS, ERROR, VALUE, DERIV, NEVALS, INFORM )
!
! A subroutine for computing a multivariate t probability and derivative.
! This subroutine uses an algorithm (QRSVN) described in the paper
! "Methods for the Computation of Multivariate t-Probabilities",
! by Alan Genz and Frank Bretz
!
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, WA 99164-3113
! Email : AlanGenz@wsu.edu
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL NxN covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; constraint I is
! T*LOWER < DELTA + CONSTR*X < T*UPPER.
! UPPER REAL, array of upper integration limits.
! INFIN INTEGER, array of integration limits flags:
! if INFIN(I) < 0, Ith limits are (-infinity, infinity);
! 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)].
! T REAL, scale factor for integration limits.
! MAXPTS INTEGER, maximum number of function values allowed. This
! parameter can be used to limit the time. A sensible
! strategy is to start with MAXPTS = 1000*N, and then
! increase MAXPTS if ERROR is too large.
! ABSEPS REAL absolute error tolerance.
! RELEPS REAL relative error tolerance.
! ERROR REAL estimated absolute error, with 99% confidence level.
! VALUE REAL estimated value for the integral
! DERIV REAL estimated value for derivative with respect to T
! NEVALS Integer number of integrand values used for this call.
! INFORM INTEGER, termination status parameter:
! if INFORM = 0, normal completion with ERROR < EPS;
! if INFORM = 1, completion with ERROR > EPS and MAXPTS
! function vaules used; increase MAXPTS to
! decrease ERROR;
! if INFORM = 2, N > NL or N < 1.
! if INFORM = 3, covariance matrix not positive semidefinite.
!
!
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, NU, M, MAXPTS
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
INTEGER, DIMENSION(M), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: T, RELEPS, ABSEPS
INTEGER, INTENT(OUT) :: INFORM, NEVALS
REAL(KIND=STND), INTENT(OUT) :: ERROR, VALUE, DERIV
!
REAL(KIND=STND), DIMENSION(2) :: V
!
NEVALS = 0
CALL MVINTD( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, T, &
VALUE, DERIV, ERROR, INFORM )
IF ( INFORM == 0 .AND. NP > 0 ) THEN
CALL MVKBRV( NP + MIN(1,NUP), 0, MAXPTS, 2, MVFUNC, &
ABSEPS, RELEPS, ERROR, V, NEVALS, INFORM )
VALUE = V(1)
DERIV = V(2)/T
ENDIF
!
END SUBROUTINE MVDSTD
!
SUBROUTINE MVINTD( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, T, &
VALUE, DERIV, ERROR, INFORM )
!
! Initialization and computation of covariance Cholesky factor.
!
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, NU, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
INTEGER, DIMENSION(M), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: T
REAL(KIND=STND), INTENT(OUT) :: ERROR, VALUE, DERIV
INTEGER, INTENT(OUT) :: INFORM
!
! Local Variables
!
REAL(KIND=STND), DIMENSION(M) :: DELTA, AT, BT
INTEGER :: I
!
! Initialization and computation of covariance Cholesky factor.
!
!
IF ( 0 < N .AND. N <= NL ) THEN
NUP = MAX( NU, 0 )
DELTA = 0
DO I = 1, M
IF ( INFIN(I) /= 0 ) AT(I) = T*LOWER(I)
IF ( INFIN(I) /= 1 ) BT(I) = T*UPPER(I)
END DO
CALL MVSORT( N, COVRNC, M, AT, CONSTR, BT, INFIN, DELTA, &
.FALSE., INFORM )
!
IF ( INFORM > 0 ) THEN
VALUE = 0
DERIV = 0
ERROR = 1
ELSE
SQTNU = SQRT( REAL( MAX( NUP, 0 ), STND ) )
VALUE = 1
DERIV = 0
IF ( NP .EQ. 0 ) THEN
ERROR = 0
ELSE IF ( NP .EQ. 1 ) THEN
IF ( INF(1) /= 1 ) THEN
VALUE = MVSTDT( NU, B(1) )
DERIV = B(1)*MVDNST( NU, B(1) )
END IF
IF ( INF(1) /= 0 ) THEN
VALUE = VALUE - MVSTDT( NU, A(1) )
DERIV = DERIV - A(1)*MVDNST( NU, A(1) )
END IF
ERROR = 1E-16_STND
END IF
END IF
ELSE
INFORM = 2
END IF
!
END SUBROUTINE MVINTD
!
SUBROUTINE MVSORT( N, COVRNC, M, LOWER, CONSTR, UPPER, INFIN, DELTA, &
PIVOT, INFORM )
!
! Subroutine to sort integration limits and determine Cholesky factor.
!
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
INTEGER, DIMENSION(M), INTENT(IN) :: INFIN
LOGICAL, INTENT(IN) :: PIVOT
INTEGER, INTENT(INOUT) :: INFORM
!
! Local Variables
!
REAL(KIND=STND), DIMENSION(N) :: V, Y
REAL(KIND=STND), DIMENSION(N,N) :: CHLFAC
INTEGER, DIMENSION(M) :: CONLIM
INTEGER :: I, J, K, L, LMN, JMX, INFT, IA, IB
REAL(KIND=STND) :: CVDIAG, S, SUMSQ, AT, BT, DT, EPSI
REAL(KIND=STND) :: AL, BL, VL, VARMN, YL
REAL(KIND=STND), PARAMETER :: ONE = 1, EPS = 10000*EPSILON(ONE)
INFORM = 0
!
! Copy parameters and drop any doubly infinite constraints.
!
MP = 0
DO I = 1, M
IF ( INFIN(I) >= 0 ) THEN
MP = MP + 1
A(MP) = 0
B(MP) = 0
INF(MP) = INFIN(I)
IF ( INF(MP) /= 0 ) A(MP) = LOWER(I)
IF ( INF(MP) /= 1 ) B(MP) = UPPER(I)
CNSTRN(MP,1:N) = CONSTR(I,1:N)
D(MP) = DELTA(I)
END IF
END DO
CHLFAC = 0
DO I = 1, N
CHLFAC(I:N,I) = COVRNC(I:N,I)
END DO
!
! Determine Cholesky factor for revised contraint matrix.
!
DO I = 1, N
!
! Scale Covariance and Constraint matrices.
!
IF ( CHLFAC(I,I) > 0 ) THEN
CVDIAG = SQRT( CHLFAC(I,I) )
CNSTRN(:,I) = CNSTRN(:,I)*CVDIAG
CHLFAC(I:N,I) = CHLFAC(I:N,I)/CVDIAG
CHLFAC(I,1:I) = CHLFAC(I,1:I)/CVDIAG
END IF
END DO
NP = 0
DO I = 1, N
CONLNG(I) = 0
EPSI = EPS*I
JMX = I
DO L = I+1, N
IF ( CHLFAC(L,L) > CHLFAC(JMX,JMX) ) JMX = L
END DO
IF ( JMX > I ) THEN
CALL MVSWAP( CHLFAC(I,I), CHLFAC(JMX,JMX) )
DO J = 1, I-1
CALL MVSWAP( CHLFAC(I,J), CHLFAC(JMX,J) )
END DO
DO J = I+1, JMX-1
CALL MVSWAP( CHLFAC(J,I), CHLFAC(JMX,J) )
END DO
DO J = JMX+1, N
CALL MVSWAP( CHLFAC(J,I), CHLFAC(J,JMX) )
END DO
DO L = 1, MP
CALL MVSWAP( CNSTRN(L,I), CNSTRN(L,JMX) )
END DO
END IF
!
! Compute Ith columns of Cholesky factor and constraint matrix.
!
IF ( CHLFAC(I,I) < -EPSI ) INFORM = 3
IF ( CHLFAC(I,I) < EPSI ) EXIT
CVDIAG = SQRT( CHLFAC(I,I) )
CHLFAC(I,I) = CVDIAG
DO L = I+1, N
CHLFAC(L,I) = CHLFAC(L,I)/CVDIAG
CHLFAC(L,I+1:L) = CHLFAC(L,I+1:L) - CHLFAC(L,I)*CHLFAC(I+1:L,I)
END DO
IF( MP > 0 ) THEN
CNSTRN(1:MP,I) = MATMUL( CNSTRN(1:MP,I:N), CHLFAC(I:N,I) )
END IF
NP = NP + 1
END DO
IF ( INFORM == 0 ) THEN
!
! Use right reflectors to reduce CNSTRN to lower triangular form
!
DO I = 1, MIN( NP-1, MP )
EPSI = EPS*I
IF ( PIVOT ) THEN
!
! Permute rows so that smallest variance variables are first
!
VARMN = 1
DO L = I, MP
V(I:NP) = CNSTRN(L,I:NP)
SUMSQ = MAX( SQRT( SUM( V(I:NP)*V(I:NP) ) ), EPSI )
S = D(L) + SUM( CNSTRN(L,1:I-1)*Y(1:I-1) )
AL = ( A(L) - S )/SUMSQ
BL = ( B(L) - S )/SUMSQ
CALL MVVRNC( AL, BL, INF(L), EPSI, YL, VL )
IF ( VL <= VARMN ) THEN
LMN = L
VARMN = VL
Y(I) = YL
END IF
END DO
V(1:NP) = CNSTRN(LMN,1:NP)
IF ( LMN > I ) THEN
CNSTRN(LMN,1:NP) = CNSTRN(I,1:NP)
CNSTRN( I,1:NP) = V(1:NP)
CALL MVSWAP( A(I), A(LMN) )
CALL MVSWAP( B(I), B(LMN) )
CALL MVSWAP( D(I), D(LMN) )
INFT = INF(I)
INF(I) = INF(LMN)
INF(LMN) = INFT
END IF
ELSE
V(I:NP) = CNSTRN(I,I:NP)
END IF
CNSTRN(I,I+1:NP) = 0
SUMSQ = SUM( V(I+1:NP)*V(I+1:NP) )
IF ( SUMSQ > EPSI ) THEN
SUMSQ = SQRT( SUMSQ + V(I)*V(I) )
IF ( V(I) < 0 ) SUMSQ = -SUMSQ
CNSTRN(I,I) = -SUMSQ
V(I) = V(I) + SUMSQ
BT = 1/( SUMSQ*V(I) )
DO L = I+1, MP
AT = BT*SUM( CNSTRN(L,I:NP)*V(I:NP) )
CNSTRN(L,I:NP) = CNSTRN(L,I:NP) - AT*V(I:NP)
END DO
END IF
END DO
!
! Scale and sort constraints
!
DO I = 1, MP
Y(1:NP) = CNSTRN(I,1:NP)
CONLIM(I) = MIN(I,NP)
JMX = 1
DO J = 1, CONLIM(I)
IF( ABS(Y(J)) > EPS*J*J ) JMX = J
END DO
Y(JMX+1:NP) = 0
CONLNG(JMX) = CONLNG(JMX) + 1
AT = A(I)
BT = B(I)
DT = D(I)
INFT = INF(I)
J = I
DO L = I-1, 1, -1
IF( JMX >= CONLIM(L) ) EXIT
CNSTRN(L+1,1:NP) = CNSTRN(L,1:NP)
A(L+1) = A(L)
B(L+1) = B(L)
D(L+1) = D(L)
INF(L+1) = INF(L)
CONLIM(L+1) = CONLIM(L)
J = L
END DO
CONLIM(J) = JMX
A(J) = AT/Y(JMX)
B(J) = BT/Y(JMX)
D(J) = DT/Y(JMX)
INF(J) = INFT
CNSTRN(J,1:NP) = Y(1:NP)/Y(JMX)
IF ( Y(JMX) < 0 ) THEN
CALL MVSWAP( A(J), B(J) )
IF ( INFT /= 2 ) INF(J) = 1 - INFT
END IF
END DO
JMX = 0
DO I = 1, NP
IF( CONLNG(I) > 0 ) JMX = I
END DO
NP = JMX
IF ( CONLNG(1) > 1 ) THEN
IF ( NUP < 1 .OR. SUM( ABS( D(1:CONLNG(1) ) ) ) == 0 ) THEN
!
! Combine first variable constraints
!
IA = 0
IB = 0
DO L = 1, CONLNG(1)
IF ( INF(L) /= 0 ) THEN
AL = A(L) - D(L)
IF ( IA > 0 ) AL = MAX( AL, AT )
IA = 1
AT = AL
END IF
IF ( INF(L) /= 1 ) THEN
BL = B(L) - D(L)
IF ( IB > 0 ) BL = MIN( BL, BT )
IB = 1
BT = BL
END IF
END DO
A(1) = AL
B(1) = BL
D(1) = 0
INF(1) = IA + IA + IB - 1
INFT = CONLNG(1) - 1
DO I = 2, MP - INFT
A(I) = A(I+INFT)
B(I) = B(I+INFT)
D(I) = D(I+INFT)
INF(I) = INF(I+INFT)
CNSTRN(I,1:NP) = CNSTRN(I+INFT,1:NP)
END DO
CONLNG(1) = 1
END IF
END IF
END IF
!
END SUBROUTINE MVSORT
!
SUBROUTINE MVVRNC( A, B, INF, EPSI, MEAN, VARANC )
!
! Compute truncated normal mean and variance
! mean = -( phi(b) - phi(a) )/( Phi(b) - Phi(a) )
! variance = 1 -( phi(b)b - phi(a)a )/( Phi(b) - Phi(a) ) - mean^2
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: INF
REAL(KIND=STND), INTENT(IN) :: A, B, EPSI
REAL(KIND=STND), INTENT(OUT) :: MEAN, VARANC
!
! Local Variables
!
REAL(KIND=STND) :: DENSA, DENSB, DISTA, DISTB, DISTDF
!
CALL MVLIMS( 0, A, B, INF, DISTA, DISTB, DISTDF )
IF ( INF /= 0 ) DENSA = MVDNST( 0, A )
IF ( INF /= 1 ) DENSB = MVDNST( 0, B )
IF ( DISTDF .GT. EPSI ) THEN
IF ( INF == 0 ) THEN
MEAN = - DENSB
VARANC = - B*DENSB
ELSE IF ( INF == 1 ) THEN
MEAN = DENSA
VARANC = A*DENSA
ELSE IF ( INF == 2 ) THEN
MEAN = DENSA - DENSB
VARANC = A*DENSA - B*DENSB
END IF
MEAN = MEAN/DISTDF
VARANC = 1 + VARANC/DISTDF - MEAN**2
ELSE
IF ( INF == 0 ) MEAN = B
IF ( INF == 1 ) MEAN = A
IF ( INF == 2 ) MEAN = ( A + B )/2
VARANC = 0
END IF
!
END SUBROUTINE MVVRNC
!
SUBROUTINE MVSWAP( X, Y )
!
! Parameter Variables
!
REAL(KIND=STND), INTENT(INOUT) :: X, Y
!
! Local Variables
!
REAL(KIND=STND) :: T
!
T = X
X = Y
Y = T
!
END SUBROUTINE MVSWAP
!
SUBROUTINE MVCRIT( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, ALPHA, &
MAXPTS, ABSEPS, ERROR, TALPHA, NVS, INFORM )
!
! A subroutine for computing critical values. The subroutine tries to
! compute the critical value TALPHA satisfying P(TALPHA) = 1 - ALPHA.
! P(T) is a multivariate T probability integral, with covariance
! matrix COVRNC, for an integration region determined by
! TALPHA*LOWER < DELTA + CONSTR*X < TALPHA*UPPER .
!
! Reference: Alan Genz and Frank Bretz, "Numerical Computation of
! Critical Values for Multiple Comparison Problems", pp. 84-87
! in 2000 Proceedings of the Statistical Computing Section,
! American Statistical Association, Alexandria, VA.
!
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, WA 99164-3113
! Email : AlanGenz@wsu.edu
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL NxN covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limit scale factors.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! LOWER < CONSTR*X < UPPER.
! UPPER REAL, array of upper integration limit scale factors.
! 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)].
! ALPHA REAL confidence factor.
! MAXPTS INTEGER, maximum number of function values allowed. This
! parameter can be used to limit the time. A sensible
! strategy is to start with MAXPTS = 1000*N, and then
! increase MAXPTS if ERROR is too large.
! ABSEPS REAL absolute error tolerance for ALPHA.
! ERROR REAL estimated absolute error, with 99% confidence level.
! TALPHA REAL estimated critical value.
! NVS Integer number of integrand values used for this call.
! INFORM INTEGER, termination status parameter:
! if INFORM = 0, normal completion with ERROR < EPS;
! if INFORM = 1, completion with ERROR > EPS and MAXPTS
! function values used; increase MAXPTS to
! decrease ERROR;
! if INFORM = 2 for invalid bounds; MVCRIT requires
! LOWER < DELTA < UPPER
!
!
! Parameter Variables
!
INTEGER, INTENT(IN) :: N, NU, M, MAXPTS
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
INTEGER, DIMENSION(M), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: ALPHA, ABSEPS
INTEGER, INTENT(OUT) :: INFORM, NVS
REAL(KIND=STND), INTENT(OUT) :: ERROR, TALPHA
!
! Local Variables
!
REAL(KIND=STND), DIMENSION(M) :: DL
REAL(KIND=STND), DIMENSION(N) :: CN
REAL(KIND=STND), DIMENSION(2) :: BD, TD
REAL(KIND=STND) :: AT, BT, ATMN, BTMN, ATMX, BTMX
REAL(KIND=STND) :: TA, TB, TC, TBF, PA, PB, PC, PCB
REAL(KIND=STND) :: DI, RE, ER, VL, DF
INTEGER :: I, INFI, IMAX, IMIN, NV
INTEGER :: IP = 0 ! > 0 Prints for debugging
!
! Use 1st order bounds to find interval for TALPHA.
!
PA = 1; PB = 0; INFORM = 0; NVS = 0
DO I = 1, M
CN = CONSTR(I,1:N)
DI = SQRT( SUM( MATMUL( CN, COVRNC )*CN ) )
IF ( INFIN(I) /= 0 ) AT = LOWER(I)/DI
IF ( INFIN(I) /= 1 ) BT = UPPER(I)/DI
PC = MVSTDC( NU, AT, BT, INFIN(I) )
IF ( PC <= PA ) THEN
IMIN = I; PA = PC; ATMN = AT; BTMN = BT
END IF
IF ( PC >= PB ) THEN
IMAX = I; PB = PC; ATMX = AT; BTMX = BT
END IF
IF ( INFIN(I) /= 0 .AND. AT >= 0 ) INFORM = 2
IF ( INFIN(I) /= 1 .AND. BT <= 0 ) INFORM = 2
END DO
IF ( INFORM == 2 ) THEN; TALPHA = 0; ERROR = 1
ELSE
IF ( INFIN(IMIN) == 0 ) THEN
TA = MVSTNV( NU, 1 - ALPHA )/BTMN
ELSE IF ( INFIN(IMIN) == 1 ) THEN
TA = -MVSTNV( NU, 1 - ALPHA )/ATMN
ELSE
TA = MVSTNV( NU, 1 - ALPHA/2 )/MAX( -ATMN, BTMN )
END IF
IF ( INFIN(IMAX) == 0 ) THEN
TB = MVSTNV( NU, 1 - ALPHA/M )/BTMX
ELSE IF ( INFIN(IMAX) == 1 ) THEN
TB = -MVSTNV( NU, 1 - ALPHA/M )/ATMX
ELSE
TB = MVSTNV( NU, 1-ALPHA/(2*M) )/MIN( -ATMX, BTMX )
END IF
IF ( IP > 0 ) PRINT '(8X, 2(F8.4,8X,F8.4))', TA, TB
!
! Use Pegasus method and 2nd order bounds to find better TALPHA interval.
!
!
DL = 0; TBF = TB; DF = 0
DO I = 1, 2
CALL LWRUPR( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DL, &
TA, BD(2), BD(1) )
PA = BD(I) - 1 + ALPHA
CALL LWRUPR( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DL, &
TB, BD(2), BD(1) )
PB = BD(I) - 1 + ALPHA
IF ( IP > 1 ) PRINT '(I8, 2(F8.4,8X,F8.4))', I, TA, TB, PA, PB
!
DO
IF ( ABS( TB - TA ) < ABSEPS/2 ) EXIT
TC = TB - PB*( TB - TA )/( PB - PA )
IF ( ( TC - TA )*( TC - TB ) > 0 ) TC = ( TA + TB )/2
CALL LWRUPR( N, COVRNC, NU, M, LOWER,CONSTR,UPPER,INFIN, DL, &
TC, BD(2), BD(1) )
PC = BD(I) - 1 + ALPHA
IF ( IP > 1 ) PRINT '(I8, 6F8.4)', I, TA, TC, TB, PA, PC, PB
IF ( PC*PB <= 0 ) THEN; TA = TB; PA = PB
ELSE IF ( PA*PC > 0 ) THEN
IF ( ABS(PB) < ABS(PA) ) THEN
TA = TB; PA = PB
END IF
ELSE; PA = PA*PB/( PB + PC )
END IF
TB = TC; PB = PC
END DO
!
DF = DF + ( ABS( ( PB - PA )/( TB - TA ) ) - DF )/I
TD(I) = ( TA + TB )/2
IF ( I == 1 ) THEN
TA = TD(I); TB = TBF
END IF
END DO
TA = TD(1); TB = TD(2); TC = ( TA + TB )/2
!
! Refine TALPHA, if necessary, using Newton iteration.
!
RE = 0; PC = 1; VL = ( BD(1) + BD(2) )/2;
DF = MIN( DF, 2E-1_STND )
DO
CALL LWRUPR( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, INFIN, DL, &
TC, BD(1), BD(2) )
PCB = MAX( ABS( BD(1) - 1 + ALPHA ), ABS( BD(2) - 1 + ALPHA ) )/2
IF ( IP > 1 ) PRINT '( I8, 3F8.4, F11.6, F7.2, 3F9.6 )', &
NVS, TA, TC, TB, PC, DF, BD(1), VL, BD(2)
ERROR = MIN( ABS( TB - TA )/2, ABS(PC/DF), ABS(PCB/DF) )
IF ( ERROR < ABSEPS ) EXIT
CALL MVDSTD( N, COVRNC, NU, M, LOWER,CONSTR,UPPER,INFIN, TC, &
MAXPTS-NVS, ABS(DF)*ABSEPS, RE, ER, VL, DF, NV, INFORM )
NVS = NVS + NV
IF ( INFORM == 1 ) EXIT
PC = VL - 1 + ALPHA
IF ( PC > 0 ) THEN; TB = ( TC + TB )/2
ELSE; TA = ( TA + TC )/2
END IF
TC = MIN( MAX( TA, TC - PC/DF ), TB )
END DO
TALPHA = TC
END IF
!
END SUBROUTINE MVCRIT
!
SUBROUTINE LWRUPR( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, LF, DELTA,T, &
LWRBND, UPRBND )
!
! Hunter-Worsley and Dawson-Sankoff Bounds
! References:
! Hsu, Jason C. (1996), Multiple Comparisons,
! Chapman and Hall, London, p. 230;
! Dawson, D. and Sankoff, A. (1967), An Inequality for Probabilities,
! Proc. AMS 18, pp. 504-507.
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! T*LOWER < DELTA + CONSTR*X < T*UPPER.
! UPPER REAL, array of upper integration limits.
! LF INTEGER, array of integration limits flags:
! if LF(I) < 0, Ith limits are (-infinity, infinity);
! if LF(I) = 0, Ith limits are (-infinity, UPPER(I)];
! if LF(I) = 1, Ith limits are [LOWER(I), infinity);
! if LF(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
! DELTA REAL array of non-centrality parameters.
! T REAL integration limit scale factor.
! LWRBND REAL lower bound for probability.
! UPRBND REAL upper bound for probability.
!
INTEGER, INTENT(IN) :: N, NU, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
REAL(KIND=STND), INTENT(IN) :: T
INTEGER, DIMENSION(M), INTENT(IN) :: LF
REAL(KIND=STND), INTENT(OUT) :: LWRBND, UPRBND
!
INTEGER :: I, J, K
REAL(KIND=STND) :: SA, SB, PT
REAL(KIND=STND), DIMENSION(M) :: L, U, MXD
REAL(KIND=STND), DIMENSION(M,M) :: CR
!
! Initialize maximum spanning tree and compute upper bound
!
CALL COVSCL( N, COVRNC, M, LOWER, CONSTR, UPPER, LF, DELTA, T, L, U, CR )
SA = 0
SB = 0
DO I = 1, M
SA = SA + MVSTDC( NU, L(I), U(I), LF(I) )
MXD(I) = 0
END DO
LWRBND = 1 - SA
!
! Compute maximum spanning tree and lower bound
!
J = 1
DO K = 1, M
DO I = 1, M
IF ( MXD(I) > MXD(J) ) J = I
END DO
LWRBND = LWRBND + MXD(J)
MXD(J) = -1
DO I = 1, M
IF ( MXD(I) >= 0 ) THEN
PT = MVBVTC( NU, (/ L(I), L(J) /), (/ U(I), U(J) /), &
(/ LF(I), LF(J) /), CR(I,J) )
SB = SB + PT
IF ( PT > MXD(I) ) MXD(I) = PT
END IF
END DO
END DO
LWRBND = MAX( LWRBND, 0E0_STND )
K = 1 + 2*SB/SA
UPRBND = 1 - 2*( SA - SB/K )/( K + 1 )
UPRBND = MIN( UPRBND, 1E0_STND )
!
END SUBROUTINE LWRUPR
!
!
SUBROUTINE LWUPRH( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, LF, DELTA,T, &
LWRBND, UPRBND )
!
! Hybrid Second-Third Order Bounds
! References:
! Bukszar, J. and Prekopa, A. (2001), Probability Bounds with Cherry
! Trees, Math. Oper. Res. 26, pp. 174-192;
! Tomescu, I. (1986), Hypertrees and Bonferroni Inequalities,
! Journal of Combinatorial Theory 41, pp. 209-217.
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! T*LOWER < DELTA + CONSTR*X < T*UPPER.
! UPPER REAL, array of upper integration limits.
! LF INTEGER, array of integration limits flags:
! if LF(I) < 0, Ith limits are (-infinity, infinity);
! if LF(I) = 0, Ith limits are (-infinity, UPPER(I)];
! if LF(I) = 1, Ith limits are [LOWER(I), infinity);
! if LF(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
! DELTA REAL array of non-centrality parameters.
! T REAL integration limit scale factor.
! LWRBND REAL lower bound for probability.
! UPRBND REAL upper bound for probability.
!
INTEGER, INTENT(IN) :: N, NU, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
REAL(KIND=STND), INTENT(IN) :: T
INTEGER, DIMENSION(M), INTENT(IN) :: LF
REAL(KIND=STND), INTENT(OUT) :: LWRBND, UPRBND
!
INTEGER :: I, J, K, L, E, V
REAL(KIND=STND) :: SA, SB, W
INTEGER, DIMENSION(M) :: NE
REAL(KIND=STND), DIMENSION(M) :: LT, UT, MXD
REAL(KIND=STND), DIMENSION(M,M) :: CR, WT
REAL(KIND=STND), DIMENSION(3) :: LW, UP, SG
INTEGER, DIMENSION(3) :: INF
INTEGER :: IM, JM, KM
!
! Compute scaled limits and correlation matrix.
!
CALL COVSCL( N, COVRNC, M, LOWER, CONSTR, UPPER, LF, DELTA, T, LT,UT,CR )
SA = 0
SB = 0
DO I = 1, M
SA = SA + MVSTDC( NU, LT(I), UT(I), LF(I) )
DO J = 1, I-1
LW(1:2) = (/ LT(I), LT(J) /)
UP(1:2) = (/ UT(I), UT(J) /)
INF(1:2) = (/ LF(I), LF(J) /)
WT(I,J) = MVBVTC( NU, LW(1:2), UP(1:2), INF(1:2), CR(I,J) )
WT(J,I) = WT(I,J)
SB = SB + WT(I,J)
END DO
END DO
!
! Compute maximum spanning tree lower bound
!
LWRBND = 1 - SA
DO I = 1, M
WT(I,I) = -1
NE(I) = 1
MXD(I) = WT(I,1)
END DO
J = 2
L = 2
DO K = 1, M - 1
DO I = 2, M
IF ( MXD(I) > MXD(J) ) J = I
END DO
LWRBND = LWRBND + MXD(J)
MXD(J) = -1
IF ( WT( J, NE(J) ) > WT(L, NE(L) ) ) L = J
DO I = 2, M
IF ( MXD(I) >= 0 .AND. WT(I,J) > MXD(I) ) THEN
MXD(I) = WT(I,J)
NE(I) = J
END IF
END DO
END DO
!
! Improve lower bound using cherry tree.
!
NE(1) = 0
MXD(1) = -1
MXD(L) = 1
MXD(NE(L)) = 1
DO L = 3, M
W = -1
DO K = 1, M
IF ( MXD(K) < 0 ) THEN
DO J = 1, M
IF ( MXD(J) > 0 .AND. ( J==NE(K) .OR. K==NE(J) ) ) THEN
DO I = 1, M
IF ( MXD(I)>0 .AND. ( I==NE(J) .OR. J==NE(I) ) ) THEN
IF( WT(I,K) > W ) THEN
W = WT(I,K)
IM = I
JM = J
KM = K
END IF
END IF
END DO
END IF
END DO
END IF
END DO
I = IM
J = JM
K = KM
MXD(K) = 1
LW = (/ LT(I), LT(J), LT(K) /)
UP = (/ UT(I), UT(J), UT(K) /)
INF = (/ LF(I), LF(J), LF(K) /)
SG = (/ CR(I,J), CR(I,K), CR(J,K) /)
LWRBND = LWRBND + WT(I,K) - MVTVTC( NU, LW, UP, INF, SG )
END DO
LWRBND = MAX( LWRBND, 0E0_STND )
!
! Compute Tomescu Bound
!
UPRBND = 1 - SA + SB
DO I = 2, M
DO J = 1, I-1
IF ( J /= NE(I) .AND. I /= NE(J) ) THEN
WT(I,J) = -1
WT(J,I) = -1
END IF
END DO
END DO
DO L = 3, M
W = -1
DO K = 1, M
E = 0
DO I = 1, M
IF ( WT(I,K) >= 0 ) THEN
E = E + 1
J = I
END IF
END DO
IF ( E == 1 .AND. WT(J,K) > W ) THEN
W = WT(J,K)
V = K
END IF
END DO
K = V
DO I = 1, M
IF ( WT(I,K) >= 0 ) J = I
END DO
WT(K,J) = -1
WT(J,K) = -1
DO I = 2, M
J = NE(I)
IF ( WT(I,J) >= 0 ) THEN
LW = (/ LT(I), LT(J), LT(K) /)
UP = (/ UT(I), UT(J), UT(K) /)
INF = (/ LF(I), LF(J), LF(K) /)
SG = (/ CR(I,J), CR(I,K), CR(J,K) /)
UPRBND = UPRBND - MVTVTC( NU, LW, UP, INF, SG )
END IF
END DO
END DO
UPRBND = MIN( UPRBND, 1E0_STND )
! K = 1 + 2*SB/SA
! UPRBND = MIN( 1 - 2*( SA - SB/K )/( K + 1 ), UPRBND )
!
END SUBROUTINE LWUPRH
!
!
SUBROUTINE LWUPRT( N, COVRNC, NU, M, LOWER, CONSTR, UPPER, LF, DELTA,T, &
LWRBND, UPRBND )
!
! Third Order Bounds
! References:
! Prekopa, A., Stochastic Programming, Kluwer Academic Publishers;
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC matrix must be positive semidefinite.
! NU INTEGER, number of degrees of freedom; if NU < 1, then
! a multivariate normal computation is completed.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! T*LOWER < DELTA + CONSTR*X < T*UPPER.
! UPPER REAL, array of upper integration limits.
! LF INTEGER, array of integration limits flags:
! if LF(I) < 0, Ith limits are (-infinity, infinity);
! if LF(I) = 0, Ith limits are (-infinity, UPPER(I)];
! if LF(I) = 1, Ith limits are [LOWER(I), infinity);
! if LF(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
! DELTA REAL array of non-centrality parameters.
! T REAL integration limit scale factor.
! LWRBND REAL lower bound for probability.
! UPRBND REAL upper bound for probability.
!
INTEGER, INTENT(IN) :: N, NU, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
REAL(KIND=STND), INTENT(IN) :: T
INTEGER, DIMENSION(M), INTENT(IN) :: LF
REAL(KIND=STND), INTENT(OUT) :: LWRBND, UPRBND
!
INTEGER :: I, J, K
REAL(KIND=STND) :: SA, SB, SC
REAL(KIND=STND), DIMENSION(M) :: LT, UT
REAL(KIND=STND), DIMENSION(M,M) :: CR
REAL(KIND=STND), DIMENSION(3) :: LW, UP, SG
INTEGER, DIMENSION(3) :: INF
!
! Compute scaled limits and correlation matrix.
!
CALL COVSCL( N, COVRNC, M, LOWER, CONSTR, UPPER, LF, DELTA, T, LT,UT,CR )
SA = 0
SB = 0
SC = 0
DO I = 1, M
SA = SA + MVSTDC( NU, LT(I), UT(I), LF(I) )
DO J = 1, I-1
LW(1:2) = (/ LT(I), LT(J) /)
UP(1:2) = (/ UT(I), UT(J) /)
INF(1:2) = (/ LF(I), LF(J) /)
SB = SB + MVBVTC( NU, LW(1:2), UP(1:2), INF(1:2), CR(I,J) )
DO K = 1, J-1
LW(3) = LT(K)
UP(3) = UT(K)
INF(3) = LF(K)
SG = (/ CR(I,J), CR(I,K), CR(J,K) /)
SC = SC + MVTVTC( NU, LW, UP, INF, SG )
END DO
END DO
END DO
J = 2 + 3*SC/SB
LWRBND = 1 - ( SA - 2*( (2*J-1)*SB - 3*SC )/( J*(J+1) ) )
I = 1 + 2*( ( M - 2 )*SB - 3*SC )/( ( M - 1 )*SA - 2*SB )
UPRBND = 1 - ( (I+2*M-1)*SA - 2*( (2*I+M-2)*SB - 3*SC )/I )/( M*(I+1) )
LWRBND = MAX( LWRBND, 0E0_STND )
UPRBND = MIN( UPRBND, 1E0_STND )
!
END SUBROUTINE LWUPRT
!
!
SUBROUTINE COVSCL( N, COVRNC, M, LOWER, CONSTR, UPPER, LF, DELTA, T, &
AT, BT, CR )
!
! Compute Correlation Matrix and Scaled Limits
!
! Parameters
!
! N INTEGER number of variables.
! COVRNC REAL NxN covariance matrix; only the diagonal and lower left
! triangle are used. COVRNC must be positive semidefinite.
! M INTEGER, number of linear constraints for integration region.
! LOWER REAL array of lower integration limits.
! CONSTR REAL MxN constraint matrix; integration region is all X with
! T*LOWER < DELTA + CONSTR*X < T*UPPER.
! UPPER REAL, array of upper integration limits.
! LF INTEGER, array of integration limits flags:
! if LF(I) = 0, Ith limits are (-infinity, UPPER(I)];
! if LF(I) = 1, Ith limits are [LOWER(I), infinity);
! if LF(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
! DELTA REAL array of non-centrality parameters.
! T REAL integration limit scale factor.
! AT REAL array of scaled lower integration limits.
! BT REAL array of scaled lower integration limits.
! CR REAL MxM correlation matrix.
!
INTEGER, INTENT(IN) :: N, M
REAL(KIND=STND), DIMENSION(N,N), INTENT(IN) :: COVRNC
REAL(KIND=STND), DIMENSION(M), INTENT(IN) :: LOWER, UPPER, DELTA
REAL(KIND=STND), DIMENSION(M,N), INTENT(IN) :: CONSTR
REAL(KIND=STND), INTENT(IN) :: T
INTEGER, DIMENSION(M), INTENT(IN) :: LF
REAL(KIND=STND), DIMENSION(M,M), INTENT(OUT) :: CR
REAL(KIND=STND), DIMENSION(M), INTENT(OUT) :: AT, BT
!
INTEGER :: I, J
REAL(KIND=STND) :: SIGMAR
REAL(KIND=STND), DIMENSION(N) :: COVTMP
!
DO I = 1, M
DO J = 1, N
COVTMP(J) = SUM( COVRNC(J,1:J-1)*CONSTR(I,1:J-1) ) &
+ SUM( COVRNC(J:N, J)*CONSTR(I,J:N) )
END DO
DO J = I, M
CR(J,I) = SUM( CONSTR(J,1:N)*COVTMP )
END DO
SIGMAR = 1/SQRT( CR(I,I) )
IF ( LF(I) /= 0 ) AT(I) = SIGMAR*( T*LOWER(I) - DELTA(I) )
IF ( LF(I) /= 1 ) BT(I) = SIGMAR*( T*UPPER(I) - DELTA(I) )
CR(I,1:I) = SIGMAR*CR(I,1:I)
CR(I:M,I) = SIGMAR*CR(I:M,I)
END DO
DO I = 1, M
CR(I,I:M) = CR(I:M,I)
END DO
!
END SUBROUTINE COVSCL
!
SUBROUTINE MVKBRV( N, MINVLS, MAXVLS, NF, FUNCTN, ABSEPS, RELEPS, &
ABSERR, FINEST, INTVLS, INFORM )
!
! Automatic Multidimensional Integration Subroutine
!
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, Wa 99164-3113
! Email : alangenz@wsu.edu
!
! Last Change: 7/7
!
! MVKBRV computes an approximation to the integral
!
! 1 1 1
! I I ... I F(X) dx(N)...dx(2)dx(1)
! 0 0 0
!
! F(X) is an NF-vector of real integrands.
!
! MVKBRV uses randomized Korobov rules.
! The primary references are
! "Randomization of Number Theoretic Methods for Multiple Integration"
! R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13, pp. 904-14,
! and
! "Optimal Parameters for Multidimensional Integration",
! P. Keast, SIAM J Numer Anal, 10, pp.831-838.
! If N > 100 the subroutine uses quasi-random Richtmeyer points for
! variables with indices > 100. A reference is
! "Methods of Numerical Integration", P.J. Davis and P. Rabinowitz,
! Academic Press, 1984, pp. 482-483.
!
!
!************** Parameters ********************************************
!***** Input parameters
! N Number of variables, must exceed 1, but not exceed 100
! MINVLS Integer minimum number of function evaluations allowed.
! MINVLS must not exceed MAXVLS. If MINVLS < 0 then the
! routine assumes a previous call has been made with
! the same integrand and continues that calculation.
! MAXVLS Integer maximum number of function evaluations allowed.
! NF Integer number of integrands, must exceed 1 and not exceed 5000
! FUNCTN user defined function to be integrated.
! It must have parameters (N,Z), where Z is a real array
! of dimension N.
!
! ABSEPS Required absolute accuracy.
! RELEPS Required relative accuracy.
!***** Output parameters
! MINVLS Actual number of function evaluations used.
! ABSERR Estimated absolute accuracy for FINEST(1).
! FINEST Estimated NF-vector of values of integrals.
! INFORM INFORM = 0 for normal exit, when
! ABSERR <= MAX(ABSEPS, RELEPS*ABS(FINEST(1)))
! and
! INTVLS <= MAXCLS.
! INFORM = 1 If MAXVLS was too small to obtain the required
! accuracy. In this case a value FINEST is returned with
! estimated absolute accuracy ABSERR.
!***********************************************************************
!
! Global Variables
!
INTEGER, INTENT(IN) :: N, NF, MINVLS, MAXVLS
INTERFACE
FUNCTION FUNCTN( NF, X ) RESULT(VALUE)
USE PRECISION_MODEL
INTEGER, INTENT(IN) :: NF
REAL(KIND=STND), INTENT(IN), DIMENSION(:) :: X
REAL(KIND=STND), DIMENSION(NF) :: VALUE
END FUNCTION FUNCTN
END INTERFACE
REAL(KIND=STND), INTENT(IN) :: ABSEPS, RELEPS
REAL(KIND=STND), DIMENSION(NF), INTENT(OUT) :: FINEST
REAL(KIND=STND), INTENT(OUT) :: ABSERR
INTEGER, INTENT(OUT) :: INTVLS, INFORM
!
! Local Variables
!
INTEGER, PARAMETER :: PL = 28, NM = 99, NMP = NM + 1
INTEGER, PARAMETER :: FM = 5000, NMX = 1000
INTEGER, PARAMETER :: MINSMP = 8
REAL(KIND=STND), PARAMETER :: ONE = 1
INTEGER, SAVE :: NP, SAMPLS
REAL(KIND=STND), DIMENSION(FM), SAVE :: VAREST
INTEGER, DIMENSION(PL,NM) :: C
INTEGER :: I, K, KMX
REAL(KIND=STND) :: VARPRD
REAL(KIND=STND), DIMENSION(NF) :: DIFINT, FINVAL, VARSQR
REAL(KIND=STND), DIMENSION(N) :: VK
INTEGER, DIMENSION(PL), PARAMETER :: P = (/ 31, 47, 73, 113, &
173, 263, 397, 593, 907, 1361, 2053, 3079, 4621, &
6947, 10427, 15641, 23473, 35221, 52837, 79259, 118891, &
178349, 267523, 401287, 601943, 902933, 1354471, 2031713 /)
INTEGER, DIMENSION(NM), PARAMETER :: C01 = (/ 12, 9, 9, &
12, 12, 12, 12, 12, 12, 12, 12, &
12, 3, 3, 3, 12, 7, 7, 12, &
12, 12, 12, 12, 12, 12, 12, 12, &
3, 3, 3, 12, 7, 7, 12, 12, &
12, 12, 12, 12, 12, 12, 12, 3, &
3, 3, 12, 7, 7, 12, 12, 12, &
12, 12, 12, 12, 12, 12, 3, 3, &
3, 12, 7, 7, 12, 12, 12, 12, &
12, 12, 12, 12, 7, 3, 3, 3, &
7, 7, 7, 3, 3, 3, 3, 3, &
3, 3, 3, 3, 3, 3, 3, 3, &
3, 3, 3, 3, 3, 3, 3, 3 /)
INTEGER, DIMENSION(NM), PARAMETER :: C02 = (/ 13, 11, 11, &
10, 15, 15, 15, 15, 22, 15, 15, &
15, 15, 6, 6, 6, 15, 15, 9, &
13, 2, 2, 2, 13, 11, 11, 10, &
15, 15, 15, 15, 15, 15, 15, 15, &
15, 6, 6, 6, 15, 15, 9, 13, &
2, 2, 2, 13, 11, 11, 10, 15, &
15, 15, 15, 15, 15, 15, 15, 15, &
6, 6, 6, 15, 15, 9, 13, 2, &
2, 2, 13, 11, 11, 10, 10, 15, &
15, 15, 15, 15, 15, 15, 15, 6, &
2, 3, 2, 3, 2, 2, 2, 2, &
2, 2, 2, 2, 2, 2, 2, 2 /)
INTEGER, DIMENSION(NM), PARAMETER :: C03 = (/ 27, 28, 10, &
11, 11, 20, 11, 20, 28, 13, 13, &
28, 13, 13, 13, 14, 14, 14, 14, &
14, 14, 14, 14, 14, 14, 14, 14, &
14, 14, 14, 14, 31, 31, 5, 5, &
5, 31, 13, 11, 11, 11, 11, 11, &
11, 13, 13, 13, 13, 13, 13, 13, &
14, 14, 14, 14, 14, 14, 14, 14, &
14, 14, 14, 14, 14, 14, 14, 14, &
31, 31, 5, 5, 5, 11, 13, 11, &
11, 11, 11, 11, 11, 11, 13, 13, &
11, 13, 5, 5, 5, 5, 14, 13, &
5, 5, 5, 5, 5, 5, 5, 5 /)
INTEGER, DIMENSION(NM), PARAMETER :: C04 = (/ 35, 27, 27, &
36, 22, 29, 29, 20, 45, 5, 5, &
5, 21, 21, 21, 21, 21, 21, 21, &
21, 21, 21, 21, 21, 21, 21, 21, &
21, 29, 17, 17, 17, 17, 17, 17, &
17, 17, 17, 17, 23, 23, 23, 23, &
23, 23, 23, 23, 23, 23, 23, 23, &
21, 27, 3, 3, 3, 24, 27, 27, &
17, 29, 29, 29, 17, 5, 5, 5, &
5, 21, 21, 21, 21, 21, 21, 21, &
21, 21, 21, 21, 21, 21, 21, 21, &
21, 17, 17, 17, 6, 17, 17, 6, &
3, 6, 6, 3, 3, 3, 3, 3 /)
INTEGER, DIMENSION(NM), PARAMETER :: C05 = (/ 64, 76, 28, &
28, 44, 44, 55, 31, 10, 52, 10, &
52, 10, 10, 38, 38, 10, 52, 10, &
10, 10, 49, 60, 49, 49, 49, 49, &
49, 49, 49, 49, 49, 49, 38, 38, &
31, 4, 4, 31, 64, 4, 4, 4, &
64, 45, 45, 45, 45, 45, 45, 66, &
66, 66, 66, 66, 66, 66, 66, 66, &
66, 66, 66, 66, 66, 66, 66, 66, &
66, 66, 11, 66, 66, 66, 66, 66, &
66, 66, 66, 66, 45, 11, 7, 3, &
2, 2, 2, 27, 5, 3, 3, 5, &
5, 2, 2, 2, 2, 2, 2, 2 /)
INTEGER, DIMENSION(NM), PARAMETER :: C06 = (/ 111, 42, 54, &
118, 20, 31, 31, 72, 17, 94, 14, &
14, 11, 14, 14, 14, 94, 10, 10, &
10, 10, 14, 14, 14, 14, 14, 14, &
14, 11, 11, 11, 8, 8, 8, 8, &
8, 8, 8, 18, 18, 18, 18, 18, &
113, 62, 62, 45, 45, 113, 113, 113, &
113, 113, 113, 113, 113, 113, 113, 113, &
113, 113, 113, 113, 113, 113, 63, 63, &
53, 63, 67, 67, 67, 67, 67, 67, &
67, 67, 67, 67, 67, 67, 67, 67, &
67, 51, 51, 51, 51, 51, 12, 51, &
12, 51, 5, 3, 3, 2, 2, 5 /)
INTEGER, DIMENSION(NM), PARAMETER :: C07 = (/ 163, 154, 83, &
43, 82, 92, 150, 59, 76, 76, 47, &
11, 11, 100, 131, 116, 116, 116, 116, &
116, 116, 138, 138, 138, 138, 138, 138, &
138, 138, 138, 101, 101, 101, 101, 101, &
101, 101, 101, 101, 101, 101, 101, 101, &
101, 101, 101, 101, 101, 101, 101, 101, &
116, 116, 116, 116, 116, 116, 100, 100, &
100, 100, 100, 138, 138, 138, 138, 138, &
101, 101, 101, 101, 101, 101, 101, 101, &
101, 101, 101, 101, 101, 101, 101, 101, &
101, 101, 101, 38, 38, 38, 38, 38, &
38, 38, 38, 3, 3, 3, 3, 3 /)
INTEGER, DIMENSION(NM), PARAMETER :: C08 = (/ 246, 189, 242, &
250, 250, 250, 102, 250, 36, 118, 196, &
118, 191, 215, 49, 121, 49, 121, 49, &
49, 49, 49, 121, 49, 49, 49, 49, &
49, 171, 171, 171, 171, 171, 171, 171, &
171, 171, 171, 171, 171, 171, 171, 171, &
171, 171, 171, 171, 171, 171, 171, 171, &
171, 171, 171, 171, 171, 171, 171, 171, &
171, 171, 171, 161, 161, 161, 161, 161, &
161, 161, 161, 14, 14, 14, 14, 14, &
14, 14, 14, 14, 14, 14, 14, 14, &
14, 14, 14, 14, 10, 10, 10, 10, &
10, 10, 103, 10, 10, 10, 10, 5 /)
INTEGER, DIMENSION(NM), PARAMETER :: C09 = (/ 347, 402, 322, &
418, 215, 220, 339, 339, 339, 337, 218, &
167, 315, 315, 315, 167, 167, 167, 167, &
361, 201, 124, 124, 124, 124, 124, 124, &
124, 124, 124, 124, 124, 231, 231, 90, &
90, 90, 90, 90, 90, 90, 90, 90, &
90, 90, 90, 90, 90, 48, 48, 48, &
48, 90, 90, 90, 90, 90, 90, 90, &
90, 90, 90, 90, 90, 90, 90, 90, &
90, 90, 90, 90, 90, 90, 90, 90, &
243, 243, 243, 243, 243, 243, 243, 243, &
243, 243, 283, 283, 283, 283, 283, 283, &
283, 283, 283, 16, 283, 16, 283, 283 /)
INTEGER, DIMENSION(NM), PARAMETER :: C10 = (/ 505, 220, 601, &
644, 612, 160, 206, 206, 206, 422, 134, &
518, 134, 134, 518, 652, 382, 206, 158, &
441, 179, 441, 56, 559, 559, 56, 56, &
56, 56, 56, 56, 56, 56, 56, 56, &
56, 56, 56, 56, 101, 101, 56, 101, &
101, 101, 101, 101, 101, 101, 101, 193, &
193, 193, 193, 193, 193, 193, 101, 101, &
101, 101, 101, 101, 101, 101, 101, 101, &
101, 101, 101, 101, 101, 101, 101, 101, &
101, 101, 101, 122, 122, 122, 122, 122, &
122, 122, 122, 122, 122, 122, 122, 122, &
122, 122, 122, 122, 101, 101, 101, 101 /)
INTEGER, DIMENSION(NM), PARAMETER :: C11 = (/ 794, 325, 960, &
528, 247, 247, 338, 366, 847, 753, 753, &
236, 334, 334, 461, 711, 652, 381, 381, &
381, 652, 381, 381, 381, 381, 381, 381, &
381, 226, 326, 326, 326, 326, 326, 326, &
326, 126, 326, 326, 326, 326, 326, 326, &
326, 326, 326, 326, 195, 195, 55, 55, &
55, 55, 55, 55, 55, 55, 55, 55, &
55, 55, 55, 55, 55, 55, 55, 55, &
55, 195, 195, 195, 195, 195, 195, 195, &
132, 132, 132, 132, 132, 132, 132, 132, &
132, 132, 132, 387, 387, 387, 387, 387, &
387, 387, 387, 387, 387, 387, 387, 387 /)
INTEGER, DIMENSION(NM), PARAMETER :: C12 = (/ 1189, 888, 259, &
1082, 725, 811, 636, 965, 497, 497, 1490, &
1490, 392, 1291, 508, 508, 1291, 1291, 508, &
1291, 508, 508, 867, 867, 934, 867, 934, &
867, 867, 867, 867, 867, 867, 867, 1284, &
1284, 1284, 1284, 1284, 1284, 1284, 1284, 1284, &
563, 563, 563, 563, 1010, 1010, 1010, 208, &
838, 563, 563, 563, 759, 759, 564, 759, &
759, 801, 801, 801, 801, 759, 759, 759, &
759, 759, 563, 563, 563, 563, 563, 563, &
563, 563, 226, 226, 226, 226, 226, 226, &
226, 226, 226, 226, 226, 226, 226, 226, &
226, 226, 226, 226, 226, 226, 226, 226 /)
INTEGER, DIMENSION(NM), PARAMETER :: C13 = (/ 1763, 1018, 1500, &
432, 1332, 1159, 126, 2240, 1719, 1284, 878, &
1983, 266, 266, 266, 266, 747, 747, 127, &
127, 2074, 127, 2074, 1400, 1383, 1383, 1383, &
1383, 1383, 1383, 1383, 1383, 1383, 1383, 1400, &
1383, 1383, 1383, 1383, 1383, 1383, 1383, 507, &
1073, 1073, 1073, 1073, 1990, 1990, 1990, 1990, &
1990, 507, 507, 507, 507, 507, 507, 507, &
507, 507, 1073, 1073, 1073, 1073, 1073, 1073, &
1073, 1073, 1073, 1073, 1073, 1073, 1073, 1073, &
1073, 1073, 1073, 22, 22, 22, 22, 22, &
22, 1073, 452, 452, 452, 452, 452, 452, &
318, 301, 301, 301, 301, 86, 86, 15 /)
INTEGER, DIMENSION(NM), PARAMETER :: C14 = (/ 2872, 3233, 1534, &
2941, 2910, 393, 1796, 919, 446, 919, 919, &
1117, 103, 103, 103, 103, 103, 103, 103, &
2311, 3117, 1101, 3117, 3117, 1101, 1101, 1101, &
1101, 1101, 2503, 2503, 2503, 2503, 2503, 2503, &
2503, 2503, 429, 429, 429, 429, 429, 429, &
429, 1702, 1702, 1702, 184, 184, 184, 184, &
184, 105, 105, 105, 105, 105, 105, 105, &
105, 105, 105, 105, 105, 105, 105, 105, &
105, 105, 105, 105, 105, 105, 105, 105, &
105, 105, 105, 105, 105, 105, 105, 105, &
105, 105, 105, 784, 784, 784, 784, 784, &
784, 784, 784, 784, 784, 784, 784, 784 /)
INTEGER, DIMENSION(NM), PARAMETER :: C15 = (/ 4309, 3758, 4034, &
1963, 730, 642, 1502, 2246, 3834, 1511, 1102, &
1102, 1522, 1522, 3427, 3427, 3928, 915, 915, &
3818, 3818, 3818, 3818, 4782, 4782, 4782, 3818, &
4782, 3818, 3818, 1327, 1327, 1327, 1327, 1327, &
1327, 1327, 1387, 1387, 1387, 1387, 1387, 1387, &
1387, 1387, 1387, 2339, 2339, 2339, 2339, 2339, &
2339, 2339, 2339, 2339, 2339, 2339, 2339, 2339, &
3148, 3148, 3148, 3148, 3148, 3148, 3148, 3148, &
3148, 3148, 3148, 3148, 3148, 3148, 3148, 3148, &
3148, 3148, 1776, 1776, 1776, 3354, 3354, 3354, &
925, 3354, 3354, 925, 925, 925, 925, 925, &
2133, 2133, 2133, 2133, 2133, 2133, 2133, 2133 /)
INTEGER, DIMENSION(NM), PARAMETER :: C16 = (/ 6610, 6977, 1686, &
3819, 2314, 5647, 3953, 3614, 5115, 423, 423, &
5408, 7652, 423, 423, 487, 6227, 2660, 6227, &
1221, 3811, 197, 4367, 4367, 1281, 1221, 351, &
351, 351, 1984, 1984, 3995, 2999, 2999, 2999, &
2999, 2999, 2999, 2063, 2063, 2063, 2063, 1644, &
2063, 2077, 2512, 2512, 2512, 2077, 2077, 2077, &
2077, 754, 754, 754, 754, 754, 754, 754, &
754, 754, 754, 754, 754, 754, 754, 754, &
754, 754, 754, 754, 1097, 1097, 754, 754, &
754, 754, 248, 754, 1097, 1097, 1097, 1097, &
222, 222, 222, 222, 754, 1982, 1982, 1982, &
1982, 1982, 1982, 1982, 1982, 1982, 1982, 1982 /)
INTEGER, DIMENSION(NM), PARAMETER :: C17 = (/ 9861, 3647, 4073, &
2535, 3430, 9865, 2830, 9328, 4320, 5913, 10365, &
8272, 3706, 6186, 7806, 7806, 7806, 8610, 2563, &
11558, 11558, 9421, 1181, 9421, 1181, 1181, 1181, &
9421, 1181, 1181, 10574, 10574, 3534, 3534, 3534, &
3534, 3534, 2898, 2898, 2898, 3450, 2141, 2141, &
2141, 2141, 2141, 2141, 2141, 7055, 7055, 7055, &
7055, 7055, 7055, 7055, 7055, 7055, 7055, 7055, &
7055, 7055, 7055, 7055, 2831, 8204, 8204, 8204, &
8204, 8204, 8204, 8204, 8204, 8204, 8204, 8204, &
8204, 8204, 8204, 8204, 8204, 8204, 8204, 8204, &
8204, 8204, 8204, 8204, 8204, 4688, 4688, 4688, &
2831, 2831, 2831, 2831, 2831, 2831, 2831, 2831 /)
INTEGER, DIMENSION(NM), PARAMETER :: C18 = (/ 10327, 7582, 7124, &
8214, 9600, 10271, 10193, 10800, 9086, 2365, 4409, &
13812, 5661, 9344, 9344, 10362, 9344, 9344, 8585, &
11114, 13080, 13080, 13080, 6949, 3436, 3436, 3436, &
13213, 6130, 6130, 8159, 8159, 11595, 8159, 3436, &
7096, 7096, 7096, 7096, 7096, 7096, 7096, 7096, &
7096, 7096, 7096, 7096, 7096, 7096, 7096, 7096, &
7096, 7096, 4377, 7096, 4377, 4377, 4377, 4377, &
4377, 5410, 5410, 4377, 4377, 4377, 4377, 4377, &
4377, 4377, 4377, 4377, 4377, 4377, 4377, 4377, &
4377, 4377, 4377, 4377, 4377, 4377, 4377, 4377, &
4377, 4377, 4377, 4377, 4377, 4377, 4377, 4377, &
4377, 4377, 4377, 440, 440, 1199, 1199, 1199 /)
INTEGER, DIMENSION(NM), PARAMETER :: C19 = (/ 19540, 19926, 11582, &
11113, 24585, 8726, 17218, 419, 4918, 4918, 4918, &
15701, 17710, 4037, 4037, 15808, 11401, 19398, 25950, &
25950, 4454, 24987, 11719, 8697, 1452, 1452, 1452, &
1452, 1452, 8697, 8697, 6436, 21475, 6436, 22913, &
6434, 18497, 11089, 11089, 11089, 11089, 3036, 3036, &
14208, 14208, 14208, 14208, 12906, 12906, 12906, 12906, &
12906, 12906, 12906, 12906, 7614, 7614, 7614, 7614, &
5021, 5021, 5021, 5021, 5021, 5021, 10145, 10145, &
10145, 10145, 10145, 10145, 10145, 10145, 10145, 10145, &
10145, 10145, 10145, 10145, 10145, 10145, 10145, 10145, &
10145, 10145, 10145, 10145, 10145, 10145, 4544, 4544, &
4544, 4544, 4544, 4544, 8394, 8394, 8394, 8394 /)
INTEGER, DIMENSION(NM), PARAMETER :: C20 = (/ 34566, 9579, 12654, &
26856, 37873, 38806, 29501, 17271, 3663, 10763, 18955, &
1298, 26560, 17132, 17132, 4753, 4753, 8713, 18624, &
13082, 6791, 1122, 19363, 34695, 18770, 18770, 18770, &
18770, 15628, 18770, 18770, 18770, 18770, 33766, 20837, &
20837, 20837, 20837, 20837, 20837, 6545, 6545, 6545, &
6545, 6545, 12138, 12138, 12138, 12138, 12138, 12138, &
12138, 12138, 12138, 12138, 12138, 12138, 12138, 12138, &
30483, 30483, 30483, 30483, 30483, 12138, 12138, 12138, &
12138, 12138, 12138, 12138, 12138, 12138, 12138, 12138, &
12138, 12138, 12138, 12138, 12138, 12138, 12138, 12138, &
9305, 11107, 11107, 11107, 11107, 11107, 11107, 11107, &
11107, 11107, 11107, 11107, 11107, 11107, 9305, 9305 /)
INTEGER, DIMENSION(NM), PARAMETER :: C21 = (/ 31929, 49367, 10982, &
3527, 27066, 13226, 56010, 18911, 40574, 20767, 20767, &
9686, 47603, 47603, 11736, 11736, 41601, 12888, 32948, &
30801, 44243, 53351, 53351, 16016, 35086, 35086, 32581, &
2464, 2464, 49554, 2464, 2464, 49554, 49554, 2464, &
81, 27260, 10681, 2185, 2185, 2185, 2185, 2185, &
2185, 2185, 18086, 18086, 18086, 18086, 18086, 17631, &
17631, 18086, 18086, 18086, 37335, 37774, 37774, 37774, &
26401, 26401, 26401, 26401, 26401, 26401, 26401, 26401, &
26401, 26401, 26401, 26401, 26401, 12982, 40398, 40398, &
40398, 40398, 40398, 40398, 3518, 3518, 3518, 37799, &
37799, 37799, 37799, 37799, 37799, 37799, 37799, 37799, &
4721, 4721, 4721, 4721, 7067, 7067, 7067, 7067 /)
INTEGER, DIMENSION(NM), PARAMETER :: C22 = (/ 40701, 69087, 77576, &
64590, 39397, 33179, 10858, 38935, 43129, 35468, 35468, &
5279, 61518, 61518, 27945, 70975, 70975, 86478, 86478, &
20514, 20514, 73178, 73178, 43098, 43098, 4701, 59979, &
59979, 58556, 69916, 15170, 15170, 4832, 4832, 43064, &
71685, 4832, 15170, 15170, 15170, 27679, 27679, 27679, &
60826, 60826, 6187, 6187, 4264, 4264, 4264, 4264, &
4264, 45567, 32269, 32269, 32269, 32269, 62060, 62060, &
62060, 62060, 62060, 62060, 62060, 62060, 62060, 1803, &
1803, 1803, 1803, 1803, 1803, 1803, 1803, 1803, &
1803, 1803, 1803, 1803, 51108, 51108, 51108, 51108, &
51108, 51108, 51108, 51108, 51108, 51108, 51108, 51108, &
55315, 55315, 54140, 54140, 54140, 54140, 54140, 13134 /)
INTEGER, DIMENSION(NM), PARAMETER :: C23 = (/ 103650, 125480, 59978, &
46875, 77172, 83021, 126904, 14541, 56299, 43636, 11655, &
52680, 88549, 29804, 101894, 113675, 48040, 113675, 34987, &
48308, 97926, 5475, 49449, 6850, 62545, 62545, 9440, &
33242, 9440, 33242, 9440, 33242, 9440, 62850, 9440, &
9440, 9440, 90308, 90308, 90308, 47904, 47904, 47904, &
47904, 47904, 47904, 47904, 47904, 47904, 41143, 41143, &
41143, 41143, 41143, 41143, 41143, 36114, 36114, 36114, &
36114, 36114, 24997, 65162, 65162, 65162, 65162, 65162, &
65162, 65162, 65162, 65162, 65162, 65162, 65162, 65162, &
65162, 47650, 47650, 47650, 47650, 47650, 47650, 47650, &
40586, 40586, 40586, 40586, 40586, 40586, 40586, 38725, &
38725, 38725, 38725, 88329, 88329, 88329, 88329, 88329 /)
INTEGER, DIMENSION(NM), PARAMETER :: C24 = (/ 165843, 90647, 59925, &
189541, 67647, 74795, 68365, 167485, 143918, 74912, 167289, &
75517, 8148, 172106, 126159, 35867, 35867, 35867, 121694, &
52171, 95354, 113969, 113969, 76304, 123709, 123709, 144615, &
123709, 64958, 64958, 32377, 193002, 193002, 25023, 40017, &
141605, 189165, 189165, 141605, 189165, 189165, 141605, 141605, &
141605, 189165, 127047, 127047, 127047, 127047, 127047, 127047, &
127047, 127047, 127047, 127047, 127047, 127047, 127047, 127047, &
127047, 127047, 127047, 127047, 127047, 127047, 127785, 127785, &
127785, 127785, 127785, 127785, 127785, 127785, 127785, 127785, &
80822, 80822, 80822, 80822, 80822, 80822, 131661, 131661, &
131661, 131661, 131661, 131661, 131661, 131661, 131661, 131661, &
131661, 131661, 131661, 131661, 131661, 131661, 7114, 131661 /)
INTEGER, DIMENSION(NM), PARAMETER :: C25 = (/ 167807, 184220, 286000, &
277107, 181789, 254376, 248814, 46455, 135834, 200779, 70712, &
202725, 202725, 219183, 196964, 134064, 73173, 73173, 134064, &
94772, 147070, 293817, 277159, 147070, 277159, 147070, 277159, &
147070, 277159, 147070, 131998, 131998, 181038, 166320, 166320, &
166320, 225573, 225573, 12826, 12826, 215668, 161276, 161276, &
215668, 263828, 263828, 185298, 230911, 201030, 201030, 28691, &
215131, 28691, 28691, 215131, 213798, 243350, 41786, 243350, &
243350, 74768, 243350, 203823, 203823, 169411, 203823, 169411, &
169411, 203823, 199155, 137514, 199155, 228957, 189678, 189678, &
189678, 189678, 74701, 92165, 265223, 265223, 95077, 52219, &
95077, 52219, 95077, 95077, 95077, 52219, 52219, 95077, &
52219, 233672, 233672, 233672, 17164, 233672, 17164, 17164 /)
INTEGER, DIMENSION(NM), PARAMETER :: C26 = (/ 333459, 375354, 102417, &
383544, 292630, 41147, 374614, 48032, 435453, 281493, 358168, &
114121, 346892, 238990, 317313, 164158, 35497, 70530, 70530, &
434839, 24754, 24754, 24754, 393656, 118711, 118711, 148227, &
271087, 355831, 91034, 417029, 417029, 91034, 91034, 417029, &
91034, 299843, 299843, 413548, 413548, 308300, 413548, 413548, &
413548, 308300, 308300, 308300, 413548, 308300, 308300, 308300, &
308300, 308300, 15311, 15311, 15311, 15311, 176255, 176255, &
23613, 23613, 23613, 23613, 23613, 23613, 172210, 204328, &
204328, 204328, 204328, 121626, 121626, 121626, 121626, 121626, &
200187, 200187, 200187, 200187, 200187, 121551, 121551, 248492, &
248492, 248492, 248492, 248492, 248492, 248492, 248492, 248492, &
248492, 248492, 248492, 13942, 13942, 13942, 13942, 13942 /)
INTEGER, DIMENSION(NM), PARAMETER :: C27 = (/ 500884, 566009, 399251, &
652979, 355008, 430235, 328722, 670680, 405585, 405585, 424646, &
670180, 670180, 641587, 215580, 59048, 633320, 81010, 20789, &
389250, 389250, 638764, 638764, 389250, 389250, 398094, 80846, &
147776, 147776, 296177, 398094, 398094, 147776, 147776, 396313, &
578233, 578233, 578233, 19482, 620706, 187095, 620706, 187095, &
126467, 241663, 241663, 241663, 241663, 241663, 241663, 241663, &
241663, 241663, 241663, 241663, 241663, 321632, 23210, 23210, &
394484, 394484, 394484, 78101, 78101, 78101, 542095, 542095, &
542095, 542095, 542095, 542095, 542095, 542095, 542095, 542095, &
542095, 542095, 542095, 542095, 542095, 542095, 542095, 542095, &
542095, 277743, 277743, 277743, 457259, 457259, 457259, 457259, &
457259, 457259, 457259, 457259, 457259, 457259, 457259, 457259 /)
INTEGER, DIMENSION(NM), PARAMETER :: C28 = (/ 858339, 918142, 501970, &
234813, 460565, 31996, 753018, 256150, 199809, 993599, 245149, &
794183, 121349, 150619, 376952, 809123, 809123, 804319, 67352, &
969594, 434796, 969594, 804319, 391368, 761041, 754049, 466264, &
754049, 754049, 466264, 754049, 754049, 282852, 429907, 390017, &
276645, 994856, 250142, 144595, 907454, 689648, 687580, 687580, &
687580, 687580, 978368, 687580, 552742, 105195, 942843, 768249, &
307142, 307142, 307142, 307142, 880619, 880619, 880619, 880619, &
880619, 880619, 880619, 117185, 117185, 117185, 117185, 117185, &
117185, 117185, 117185, 117185, 117185, 117185, 60731, 60731, &
60731, 60731, 60731, 60731, 60731, 60731, 60731, 60731, &
60731, 178309, 178309, 178309, 178309, 74373, 74373, 74373, &
74373, 74373, 74373, 74373, 74373, 214965, 214965, 214965 /)
INTEGER, DIMENSION(NMX), PARAMETER :: PRIME = (/ &
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
73, 79, 83, 89, 97, 101, 103, 107, 109, 113, &
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, &
179, 181, 191, 193, 197, 199, 211, 223, 227, 229, &
233, 239, 241, 251, 257, 263, 269, 271, 277, 281, &
283, 293, 307, 311, 313, 317, 331, 337, 347, 349, &
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, &
419, 421, 431, 433, 439, 443, 449, 457, 461, 463, &
467, 479, 487, 491, 499, 503, 509, 521, 523, 541, &
547, 557, 563, 569, 571, 577, 587, 593, 599, 601, &
607, 613, 617, 619, 631, 641, 643, 647, 653, 659, &
661, 673, 677, 683, 691, 701, 709, 719, 727, 733, &
739, 743, 751, 757, 761, 769, 773, 787, 797, 809, &
811, 821, 823, 827, 829, 839, 853, 857, 859, 863, &
877, 881, 883, 887, 907, 911, 919, 929, 937, 941, &
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, &
1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, &
1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, &
1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, &
1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, &
1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, &
1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, &
1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, &
1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, &
1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, &
1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, &
1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, &
1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, &
1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, &
1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, &
2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, &
2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, &
2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, &
2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, &
2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, &
2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, &
2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, &
2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, &
2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, &
2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, &
2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, &
2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, &
3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, &
3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, &
3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, &
3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, &
3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, &
3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, &
3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, &
3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, &
3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, &
3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, &
3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, &
3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, &
4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, &
4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, &
4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, &
4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, &
4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, &
4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, &
4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, &
4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, &
4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, &
4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, &
4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, &
4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, &
5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, &
5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, &
5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, &
5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, &
5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, &
5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, &
5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, &
5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, &
5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, &
5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, &
5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, &
5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, &
6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, &
6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, &
6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, &
6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, &
6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, &
6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, &
6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, &
6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, &
6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, &
6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, &
6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, &
7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, &
7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, &
7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, &
7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, &
7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, &
7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, &
7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, &
7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, &
7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, &
7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919 /)
!
C = RESHAPE( (/ C01, C02, C03, C04, C05, C06, C07, C08, C09, C10, &
C11, C12, C13, C14, C15, C16, C17, C18, C19, C20, &
C21, C22, C23, C24, C25, C26, C27, C28 /), &
(/ PL, NM /), ORDER = (/ 2, 1 /) )
!
INFORM = 1
INTVLS = 0
IF ( MINVLS >= 0 ) THEN
FINEST = 0
VAREST = 0
SAMPLS = MINSMP
DO I = MIN( N, 10 ), PL
NP = I
IF ( MINVLS < 2*SAMPLS*P(I) ) EXIT
END DO
IF ( MINVLS >= 2*SAMPLS*P(PL) ) SAMPLS = MINVLS/( 2*P(PL) )
ENDIF
DO
VK(1) = ONE/P(NP)
K = 1
DO I = 2, N
IF ( I <= NMP ) THEN
! VK(I) = MODULO( C(NP,N-1)*VK(I-1), ONE ) Bug fixed 11/3/5
! K = MODULO( C(NP,N-1)*REAL(K,STND), REAL(P(NP),STND) ) 7/3/7
K = MODULO( C(NP,MIN(N-1,NM))*REAL(K,STND), REAL(P(NP),STND) )
VK(I) = K*VK(1)
ELSE
VK(I) = MODULO( SQRT( REAL( PRIME(I-NMP), STND ) ), ONE )
END IF
END DO
!
FINVAL = 0
VARSQR = 0
DO I = 1, SAMPLS
DIFINT = ( MVKRSV( N, P(NP), VK, NF, FUNCTN ) - FINVAL )/I
FINVAL = FINVAL + DIFINT
VARSQR = ( I - 2 )*VARSQR/I + DIFINT**2
END DO
INTVLS = INTVLS + 2*SAMPLS*P(NP)
KMX = 1
DO K = 1, NF
VARPRD = VAREST(K)*VARSQR(K)
FINEST(K) = FINEST(K) + ( FINVAL(K) - FINEST(K) )/( 1 + VARPRD )
IF ( VARSQR(K) > 0 ) VAREST(K) = ( 1 + VARPRD )/VARSQR(K)
END DO
ABSERR = 7*SQRT( VARSQR(KMX)/( 1 + VARPRD ) )/2
IF ( ABSERR > MAX( ABSEPS, ABS(FINEST(KMX))*RELEPS ) ) THEN
IF ( NP < PL ) THEN
NP = NP + 1
ELSE
SAMPLS = MIN( 3*SAMPLS/2, ( MAXVLS - INTVLS )/( 2*P(NP) ) )
SAMPLS = MAX( MINSMP, SAMPLS )
ENDIF
IF ( INTVLS + 2*SAMPLS*P(NP) > MAXVLS ) EXIT
ELSE
INFORM = 0
EXIT
ENDIF
END DO
!
END SUBROUTINE MVKBRV
!
FUNCTION MVKRSV( N, PRIME, VK, NF, FUNCTN ) RESULT(VALUES)
!
! Global Variables
!
INTEGER, INTENT(IN) :: N, NF, PRIME
REAL(KIND=STND), DIMENSION(:), INTENT(INOUT) :: VK
INTERFACE
FUNCTION FUNCTN( NF, X ) RESULT(VALUE)
USE PRECISION_MODEL
INTEGER, INTENT(IN) :: NF
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: X
REAL(KIND=STND), DIMENSION(NF) :: VALUE
END FUNCTION FUNCTN
END INTERFACE
REAL(KIND=STND), DIMENSION(NF) :: VALUES
!
! Local Variables
!
INTEGER :: K
REAL(KIND=STND), DIMENSION(N) :: X, R
!
VALUES = 0
DO K = 1, N
R(K) = UNIFRM()
END DO
DO K = 1, PRIME
R = R + VK
WHERE ( R > 1 ) R = R - 1
X = ABS ( 2*R - 1 )
VALUES = VALUES + ( ( FUNCTN(NF,X) + FUNCTN(NF,1-X) )/2 - VALUES )/K
END DO
!
END FUNCTION MVKRSV
!
FUNCTION UNIFRM() RESULT(UNI)
!
! Uniform (0,1) random number generator
!
! Reference:
! L'Ecuyer, Pierre (1996),
! "Combined Multiple Recursive Random Number Generators"
! Operations Research 44, pp. 816-822.
!
REAL(KIND=STND) :: UNI
!
INTEGER :: Z, H, P12, P13, P21, P23
!
! Some eight digit primes for seeds
!
INTEGER, SAVE :: X10 = 15485857, X11 = 17329489, X12 = 36312197, &
X20 = 55911127, X21 = 75906931, X22 = 96210113
REAL(KIND=STND), PARAMETER :: INVMP1 = 4.656612873077392578125E-10_STND
! INVMP1 = 1/(M1+1)
INTEGER, PARAMETER :: M1 = 2147483647, M2 = 2145483479
INTEGER, PARAMETER :: A12 = 63308, Q12 = 33921, R12 = 12979
INTEGER, PARAMETER :: A13 = -183326, Q13 = 11714, R13 = 2883
INTEGER, PARAMETER :: A21 = 86098, Q21 = 24919, R21 = 7417
INTEGER, PARAMETER :: A23 = -539608, Q23 = 3976, R23 = 2071
!
!
! Component 1
!
H = X10/Q13
P13 = -A13*( X10 - H*Q13 ) - H*R13
H = X11/Q12
P12 = A12*( X11 - H*Q12 ) - H*R12
IF ( P13 < 0 ) THEN
P13 = P13 + M1
END IF
IF ( P12 < 0 ) THEN
P12 = P12 + M1
END IF
X10 = X11
X11 = X12
X12 = P12 - P13
IF ( X12 < 0 ) THEN
X12 = X12 + M1
END IF
!
! Component 2
!
H = X20/Q23
P23 = -A23*( X20 - H*Q23 ) - H*R23
H = X22/Q21
P21 = A21*( X22 - H*Q21 ) - H*R21
IF ( P23 < 0 ) THEN
P23 = P23 + M2
END IF
IF ( P21 < 0 ) THEN
P21 = P21 + M2
END IF
X20 = X21
X21 = X22
X22 = P21 - P23
IF ( X22 < 0 ) THEN
X22 = X22 + M2
END IF
!
! Combination
!
Z = X12 - X22
IF ( Z <= 0 ) THEN
Z = Z + M1
END IF
UNI = Z*INVMP1
!
END FUNCTION UNIFRM
!
FUNCTION MVDNST( NU, X ) RESULT(DNSTY)
!
! Student's T density function
!
INTEGER, INTENT(IN) :: NU
REAL(KIND=STND), INTENT(IN) :: X
REAL(KIND=STND) :: DNSTY
!
REAL(KIND=STND), PARAMETER :: PI = 3.141592653589793E0_STND
REAL(KIND=STND), PARAMETER :: SQTWPI = 2.506628274631001E0_STND
REAL(KIND=STND) :: PROD
INTEGER :: I
DNSTY = 0
IF ( NU < 1 ) THEN
IF ( ABS(X) < 10 ) DNSTY = EXP( -X*X/2 )/SQTWPI
ELSE
PROD = 1/SQRT( REAL( NU, STND ) )
DO I = NU - 2, 1, -2
PROD = PROD*( I + 1 )/I
END DO
IF ( MODULO( NU, 2 ) == 0 ) THEN
PROD = PROD/2
ELSE
PROD = PROD/PI
END IF
DNSTY = PROD/SQRT( 1 + X*X/NU )**( NU + 1 )
END IF
!
END FUNCTION MVDNST
!
FUNCTION MVUVT( NU, A, B, INFIN ) RESULT(VALUE)
!
! Univariate Distribution Function
!
! Parameters
!
INTEGER, INTENT(IN) :: NU, INFIN
REAL(KIND=STND), INTENT(IN) :: A, B
REAL(KIND=STND) :: LOWER, UPPER, VALUE
!
LOWER = 0
UPPER = 1
IF ( INFIN >= 0 ) THEN
IF ( INFIN /= 0 ) LOWER = MVSTDT( NU, A )
IF ( INFIN /= 1 ) UPPER = MVSTDT( NU, B )
END IF
UPPER = MAX( UPPER, LOWER )
VALUE = UPPER - LOWER
!
END FUNCTION MVUVT
!
FUNCTION MVSTDT( NU, T ) RESULT(DSTRB)
!
! Student t Distribution Function
!
! T
! MVSTDT = C I ( 1 + y*y/NU )**( -(NU+1)/2 ) dy
! NU -00
!
INTEGER, INTENT(IN) :: NU
REAL(KIND=STND), INTENT(IN) :: T
REAL(KIND=STND) :: DSTRB
!
REAL(KIND=STND), PARAMETER :: PI = 3.141592653589793E0_STND
REAL(KIND=STND) :: CSTHE, SNTHE, POLYN, TT, TS, RN
INTEGER :: J
IF ( NU < 1 ) THEN
DSTRB = MVPHI( T )
ELSE IF ( NU == 1 ) THEN
DSTRB = ( 1 + 2*ATAN( T )/PI )/2
ELSE IF ( NU == 2) THEN
DSTRB = ( 1 + T/SQRT( 2 + T*T ))/2
ELSE
TT = T*T
CSTHE = NU/( NU + TT )
POLYN = 1
DO J = NU - 2, 2, -2
POLYN = 1 + ( J - 1 )*CSTHE*POLYN/J
END DO
IF ( MODULO( NU, 2 ) == 1 ) THEN
RN = NU
TS = T/SQRT(RN)
DSTRB = ( 1 + 2*( ATAN( TS ) + TS*CSTHE*POLYN )/PI )/2
ELSE
SNTHE = T/SQRT( NU + TT )
DSTRB = ( 1 + SNTHE*POLYN )/2
END IF
IF ( DSTRB < 0 ) DSTRB = 0
ENDIF
!
END FUNCTION MVSTDT
!
FUNCTION MVSTNV( N, Z ) RESULT(STINV)
!
! Inverse Student t Distribution Function
!
! STINV
! Z = C I (1 + y*y/N)**(-(N+1)/2) dy
! N -INF
!
! Reference: G.W. Hill, Comm. ACM Algorithm 395
! Comm. ACM 13 (1970), pp. 619-620.
!
! Enhancements for double precision and other modifications by
! Alan Genz, 1993-2001.
!
INTEGER, INTENT(IN) :: N
REAL(KIND=STND), INTENT(IN) :: Z
REAL(KIND=STND) :: STINV
!
INTEGER :: J
REAL(KIND=STND) :: P, A, B, C, D, X, XX, Y, CONST, STJAC
REAL(KIND=STND), SAVE :: NOLD = 0
REAL(KIND=STND), PARAMETER :: PI = 3.141592653589793E0_STND, TWO = 2
!
IF ( 0 < Z .AND. Z < 1 ) THEN
IF ( N < 1 ) THEN
STINV = MVPHNV( Z )
ELSE IF ( N == 1 ) THEN
STINV = TAN( PI*( 2*Z - 1 )/2 )
ELSE IF ( N == 2) THEN
STINV = ( 2*Z - 1 )/SQRT( 2*Z*( 1 - Z ) )
ELSE
IF ( 2*Z >= 1 ) THEN
P = 2*( 1 - Z )
ELSE
P = 2*Z
END IF
A = 1/( N - 1/TWO )
B = 48/( A*A )
C = ( ( 20700*A/B - 98 )*A - 16 )*A + 96.36E0_STND
D = ( ( 94.5E0_STND/( B + C ) - 3 )/B + 1 )*SQRT( A*PI/2 )*N
X = D*P
Y = X**( TWO/N )
IF ( Y .GT. A + 0.05E0_STND ) THEN
X = MVPHNV( P/2 )
Y = X*X
IF ( N < 5 ) C = C + 3*( N - 9/TWO )*( 10*X + 6 )/100
C = ( ( (D*X - 100)*X/20 - 7 )*X - 2 )*X + B + C
Y = ( ( ( ( (4*Y+63)*Y/10+36 )*Y+94.5E0_STND )/C-Y-3 )/B+1 )*X
Y = A*Y*Y
IF ( Y > 0.002E0_STND ) THEN
Y = EXP(Y) - 1
ELSE
Y = Y*( 1 + Y/2 )
ENDIF
ELSE
Y = ( ( 1/( ( (N+6)/(N*Y) - 0.089E0_STND*D - 0.822E0_STND ) &
*(3*N+6) ) + 0.5E0_STND/(N+4) )*Y - 1 )*(N+1)/(N+2) + 1/Y
END IF
STINV = SQRT(N*Y)
IF ( 2*Z < 1 ) STINV = -STINV
IF ( ABS( STINV ) > 0 ) THEN
!
! Use one third order Schroeder correction to single precision result
!
X = STINV
D = Z - MVSTDT( N, X )
IF ( N /= NOLD ) THEN
NOLD = N
IF ( MODULO( N, 2 ) == 0 ) THEN
CONST = SQRT(NOLD)*2
ELSE
CONST = SQRT(NOLD)*PI
END IF
DO J = N - 2, 1, -2
CONST = J*CONST/( J + 1 )
END DO
END IF
XX = 1 + X*X/N
STJAC = CONST*XX**( ( N + 1 )/TWO )
Y = D*STJAC
STINV = X + Y*( 1 + Y*(N+1)/( 2*( X + N/X ) ) )
END IF
END IF
ELSE
!
! Use cutoff values for Z near 0 or 1.
!
STINV = SQRT( N/( 2E-16_STND*SQRT( 2*PI*N ) )**( TWO/N ) )
IF ( 2*Z < 1 ) STINV = -STINV
END IF
!
END FUNCTION MVSTNV
!
FUNCTION MVSTDC( NU, L, U, INFN ) RESULT(DSTRBC)
!
! Complementary Student t Distribution Function
!
! if INFN, then MVSTDC computes
! < 0 0
! 0 P( X > U )
! 1 P( X < L )
! 2 P( X > U ) + P( X < L )
!
INTEGER, INTENT(IN) :: NU, INFN
REAL(KIND=STND), INTENT(IN) :: L, U
REAL(KIND=STND) :: DSTRBC
!
IF ( INFN == 0 ) THEN
DSTRBC = MVSTDT( NU, -U )
ELSE IF ( INFN == 1 ) THEN
DSTRBC = MVSTDT( NU, L )
ELSE IF ( INFN == 2 ) THEN
DSTRBC = MVSTDT( NU, L ) + MVSTDT( NU, -U )
ELSE
DSTRBC = 0
END IF
!
END FUNCTION MVSTDC
!
FUNCTION MVPHI( Z ) RESULT(P)
!
! Normal distribution probabilities accurate to 1.e-15.
! 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
!
REAL(KIND=STND), INTENT(IN) :: Z
REAL(KIND=STND) :: P
!
REAL(KIND=STND) :: EXPNTL, ZABS
REAL(KIND=STND), PARAMETER :: &
P0 = 220.2068679123761E0_STND, &
P1 = 221.2135961699311E0_STND, &
P2 = 112.0792914978709E0_STND, &
P3 = 33.91286607838300E0_STND, &
P4 = 6.373962203531650E0_STND, &
P5 = 0.7003830644436881E0_STND, &
P6 = 0.03526249659989109E0_STND
REAL(KIND=STND), PARAMETER :: &
Q0 = 440.4137358247522E0_STND, &
Q1 = 793.8265125199484E0_STND, &
Q2 = 637.3336333788311E0_STND, &
Q3 = 296.5642487796737E0_STND, &
Q4 = 86.78073220294608E0_STND, &
Q5 = 16.06417757920695E0_STND, &
Q6 = 1.755667163182642E0_STND, &
Q7 = 0.08838834764831844E0_STND
!
REAL(KIND=STND), PARAMETER :: CUTOFF = 7.071067811865475_STND, &
ROOTPI = 2.506628274631000502415765_STND
!
ZABS = ABS(Z)
!
! |Z| > 37
!
IF ( ZABS > 37 ) THEN
P = 0
ELSE
!
! |Z| <= 37
!
EXPNTL = EXP( -ZABS**2/2 )
!
! |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/( ZABS + 2/( ZABS + 3/( ZABS &
+ 4/( ZABS + 0.65E0_STND ) ) ) ) )/ROOTPI
END IF
END IF
IF ( Z > 0 ) THEN
P = 1 - P
END IF
!
END FUNCTION MVPHI
!
FUNCTION MVPHNV( P ) RESULT(PHINV)
!
! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
!
! Produces the normal deviate Z corresponding to a given lower
! tail area of P.
!
REAL(KIND=STND), INTENT(IN) :: P
REAL(KIND=STND) :: PHINV
!
REAL(KIND=STND) :: Q, R
REAL(KIND=STND), PARAMETER :: SPLIT1 = 0.425E0_STND, SPLIT2 = 5, &
CONST1 = 0.180625E0_STND, CONST2 = 1.6E0_STND
!
! Coefficients for P close to 0.5
!
REAL(KIND=STND), PARAMETER :: &
A0 = 3.3871328727963666080E00_STND, &
A1 = 1.3314166789178437745E+2_STND, &
A2 = 1.9715909503065514427E+3_STND, &
A3 = 1.3731693765509461125E+4_STND, &
A4 = 4.5921953931549871457E+4_STND, &
A5 = 6.7265770927008700853E+4_STND, &
A6 = 3.3430575583588128105E+4_STND, &
A7 = 2.5090809287301226727E+3_STND, &
B1 = 4.2313330701600911252E+1_STND, &
B2 = 6.8718700749205790830E+2_STND, &
B3 = 5.3941960214247511077E+3_STND, &
B4 = 2.1213794301586595867E+4_STND, &
B5 = 3.9307895800092710610E+4_STND, &
B6 = 2.8729085735721942674E+4_STND, &
B7 = 5.2264952788528545610E+3_STND
!
! Coefficients for P not close to 0, 0.5 or 1.
!
REAL(KIND=STND), PARAMETER :: &
C0 = 1.42343711074968357734E00_STND, &
C1 = 4.63033784615654529590E00_STND, &
C2 = 5.76949722146069140550E00_STND, &
C3 = 3.64784832476320460504E00_STND, &
C4 = 1.27045825245236838258E00_STND, &
C5 = 2.41780725177450611770E-1_STND, &
C6 = 2.27238449892691845833E-2_STND, &
C7 = 7.74545014278341407640E-4_STND, &
D1 = 2.05319162663775882187E00_STND, &
D2 = 1.67638483018380384940E00_STND, &
D3 = 6.89767334985100004550E-1_STND, &
D4 = 1.48103976427480074590E-1_STND, &
D5 = 1.51986665636164571966E-2_STND, &
D6 = 5.47593808499534494600E-4_STND, &
D7 = 1.05075007164441684324E-9_STND
!
! Coefficients for P near 0 or 1.
!
REAL(KIND=STND), PARAMETER :: &
E0 = 6.65790464350110377720E00_STND, &
E1 = 5.46378491116411436990E00_STND, &
E2 = 1.78482653991729133580E00_STND, &
E3 = 2.96560571828504891230E-1_STND, &
E4 = 2.65321895265761230930E-2_STND, &
E5 = 1.24266094738807843860E-3_STND, &
E6 = 2.71155556874348757815E-5_STND, &
E7 = 2.01033439929228813265E-7_STND, &
F1 = 5.99832206555887937690E-1_STND, &
F2 = 1.36929880922735805310E-1_STND, &
F3 = 1.48753612908506148525E-2_STND, &
F4 = 7.86869131145613259100E-4_STND, &
F5 = 1.84631831751005468180E-5_STND, &
F6 = 1.42151175831644588870E-7_STND, &
F7 = 2.04426310338993978564E-15_STND
!
Q = ( 2*P - 1 )/2
IF ( ABS(Q) <= SPLIT1 ) THEN
R = CONST1 - Q*Q
PHINV = 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 + 1 )
ELSE
R = MIN( P, 1 - P )
IF ( R > 0 ) THEN
R = SQRT( -LOG(R) )
IF ( R <= SPLIT2 ) THEN
R = R - CONST2
PHINV = ( ( ( ((((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 + 1 )
ELSE
R = R - SPLIT2
PHINV = ( ( ( ((((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 + 1 )
END IF
ELSE
PHINV = 9
END IF
IF ( Q < 0 ) THEN
PHINV = - PHINV
END IF
END IF
!
END FUNCTION MVPHNV
!
FUNCTION MVCHNV( N, P ) RESULT(R)
!
! MVCHNV
! P = 1 - K I exp(-t*t/2) t**(N-1) dt, for N >= 1.
! N 0
!
INTEGER, INTENT(IN) :: N
REAL(KIND=STND), INTENT(IN) :: P
REAL(KIND=STND) :: R
!
INTEGER :: I
REAL(KIND=STND) :: RI
REAL(KIND=STND), PARAMETER :: LRP = -.22579135264472743235E0_STND
! LRP = LOG( SQRT( 2/PI ) )
REAL(KIND=STND), SAVE :: LKN
INTEGER, SAVE :: NO = 0
IF ( N < 1 ) THEN
R = 1
ELSE IF ( N == 1 ) THEN
R = -MVPHNV( P/2 )
ELSE IF ( P < 1 ) THEN
IF ( N == 2 ) THEN
R = SQRT( -2*LOG(P) )
ELSE
IF ( N /= NO ) THEN
NO = N
LKN = 0
DO I = N - 2, 2, -2
RI = I
LKN = LKN - LOG( RI )
END DO
IF ( MODULO( N, 2 ) == 1 ) LKN = LKN + LRP
END IF
IF ( N >= -5*LOG( 1 - P )/4 ) THEN
R = 2E0_STND/( 9*N )
R = N*( -MVPHNV(P)*SQRT(R) + 1 - R )**3
IF ( R > 2*N+6 ) THEN
R = 2*( LKN - LOG(P) ) + ( N - 2 )*LOG(R)
END IF
ELSE
R = EXP( 2*( LOG( ( 1 - P )*N ) - LKN )/N )
END IF
R = SQRT(R)
RI = R
R = MVCHNC( LKN, N, P, R )
IF ( ABS( R - RI ) > 2E-1_STND ) THEN
R = RI
ELSE
DO I = 1, 5
IF ( ABS( R - RI ) < 1E-6_STND ) EXIT
RI = R
R = MVCHNC( LKN, N, P, R )
END DO
END IF
END IF
ELSE
R = 0
END IF
!
END FUNCTION MVCHNV
!
FUNCTION MVCHNC( LKN, N, P, R ) RESULT(CHNC)
!
! Third order Schroeder correction to R for MVCHNV
!
INTEGER, INTENT(IN) :: N
REAL(KIND=STND), INTENT(IN) :: P, LKN, R
REAL(KIND=STND) :: CHNC
!
REAL(KIND=STND) :: DF
DF = P - MVCHI( N, R )
DF = MAX( P - 1, MIN( DF, P ) )
DF = DF/EXP( LKN + ( N - 1 )*LOG(R) - R*R/2 )
CHNC = R - DF*( 1 - DF*( R - ( N - 1 )/R )/2 )
!
END FUNCTION MVCHNC
!
FUNCTION MVCHI( N, R ) RESULT(CHI)
!
! R
! CHI = 1 - K I exp(-t*t/2) t**(N-1) dt, for N >= 1.
! N 0
!
INTEGER, INTENT(IN) :: N
REAL(KIND=STND), INTENT(IN) :: R
REAL(KIND=STND) :: CHI
!
INTEGER :: I
REAL(KIND=STND) :: RN, AL, RI, RR, AI, BI, CI, DI, DL
REAL(KIND=STND), PARAMETER :: RP = 0.79788456080286535588E0_STND
! RP = SQRT( 2/PI )
REAL(KIND=STND), PARAMETER :: LRP = 0.12078223763524522234E0_STND
! LRP = LOG( 2/SQRT(PI) )
REAL(KIND=STND), PARAMETER :: ONE = 1, EPS = EPSILON(ONE)
REAL(KIND=STND), SAVE :: LGN
INTEGER, SAVE :: NO = 0
RR = R*R
IF ( N < 2 ) THEN
CHI = 2*MVPHI(-R)
ELSE IF ( N < 100 ) THEN
!
! Use standard Chi series
!
RN = 1
DO I = N - 2, 2, -2
RN = 1 + RR*RN/I
END DO
IF ( MODULO( N, 2 ) == 0 ) THEN
CHI = EXP( LOG( RN ) - RR/2 )
ELSE
CHI = EXP( LOG( R*RP*RN ) - RR/2 ) + 2*MVPHI(-R)
ENDIF
ELSE
IF ( N /= NO ) THEN
NO = N
LGN = 0
DO I = N - 2, 2, -2
RI = I
LGN = LGN + LOG( 2/RI )
END DO
IF ( MODULO( N, 2 ) == 1 ) LGN = LGN + LRP
END IF
RR = RR/2
AL = N
AL = AL/2
IF ( RR < AL + 1 ) THEN
!
! Use Incomplete Gamma series
!
DL = EXP( -RR + AL*LOG(RR) + LGN )/AL
CHI = DL
DO I = 1, 1000
DL = DL*RR/( AL + I )
CHI = CHI + DL
IF ( ABS( DL*RR/( AL + I + 1 - RR ) ) < EPS ) EXIT
END DO
CHI = 1 - CHI
ELSE
!
! Use Incomplete Gamma continued fraction
!
BI = RR + 1 - AL
CI = 1/EPS
DI = BI
CHI = EXP( -RR + AL*LOG(RR) + LGN )/BI
DO I = 1, 200
AI = I*( AL - I )
BI = BI + 2
CI = BI + AI/CI
IF ( CI == 0 ) CI = EPS
DI = BI + AI/DI
IF ( DI == 0 ) DI = EPS
DL = CI/DI
CHI = CHI*DL
IF ( ABS( DL - 1 ) < EPS ) EXIT
END DO
END IF
END IF
!
END FUNCTION MVCHI
!
FUNCTION MVBVT( NU, LOWER, UPPER, INFIN, CORREL ) RESULT(BVNBVT)
!
! A function for computing bivariate normal and t probabilities.
!
! Parameters
!
! NU INTEGER degrees of freedom parameter; NU < 1 gives normal case.
! 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.
!
INTEGER, INTENT(IN) :: NU
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: LOWER, UPPER
INTEGER, DIMENSION(:), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: CORREL
REAL(KIND=STND) :: BVNBVT
!
IF ( INFIN(1) < 0 ) THEN
BVNBVT = MVUVT( NU, LOWER(2), UPPER(2), INFIN(2) )
ELSE IF ( INFIN(2) < 0 ) THEN
BVNBVT = MVUVT( NU, LOWER(1), UPPER(1), INFIN(1) )
ELSE IF ( INFIN(1) == 2 .AND. INFIN(2) == 2 ) THEN
BVNBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) &
- MVBVTL ( NU, UPPER(1), LOWER(2), CORREL ) &
- MVBVTL ( NU, LOWER(1), UPPER(2), CORREL ) &
+ MVBVTL ( NU, LOWER(1), LOWER(2), CORREL )
ELSE IF ( INFIN(1) == 2 .AND. INFIN(2) == 1 ) THEN
BVNBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL ) &
- MVBVTL ( NU, -UPPER(1), -LOWER(2), CORREL )
ELSE IF ( INFIN(1) == 1 .AND. INFIN(2) == 2 ) THEN
BVNBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL ) &
- MVBVTL ( NU, -LOWER(1), -UPPER(2), CORREL )
ELSE IF ( INFIN(1) == 2 .AND. INFIN(2) == 0 ) THEN
BVNBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) &
- MVBVTL ( NU, LOWER(1), UPPER(2), CORREL )
ELSE IF ( INFIN(1) == 0 .AND. INFIN(2) == 2 ) THEN
BVNBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) &
- MVBVTL ( NU, UPPER(1), LOWER(2), CORREL )
ELSE IF ( INFIN(1) == 1 .AND. INFIN(2) == 0 ) THEN
BVNBVT = MVBVTL ( NU, -LOWER(1), UPPER(2), -CORREL )
ELSE IF ( INFIN(1) == 0 .AND. INFIN(2) == 1 ) THEN
BVNBVT = MVBVTL ( NU, UPPER(1), -LOWER(2), -CORREL )
ELSE IF ( INFIN(1) == 1 .AND. INFIN(2) == 1 ) THEN
BVNBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL )
ELSE IF ( INFIN(1) == 0 .AND. INFIN(2) == 0 ) THEN
BVNBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL )
END IF
!
END FUNCTION MVBVT
!
FUNCTION MVBVTL( NU, DH, DK, R ) RESULT(BVT)
!
! 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.
!
! mvbvtl - calculate the probability that x < dh and y < dk.
!
! parameters
!
! nu number of degrees of freedom
! dh 1st lower integration limit
! dk 2nd lower integration limit
! r correlation coefficient
!
REAL(KIND=STND), INTENT(IN) :: DH, DK, R
REAL(KIND=STND) :: BVT
INTEGER, INTENT(IN) :: NU
!
INTEGER :: J, HS, KS
REAL(KIND=STND) :: ORS, HRK, KRH, SQNU, GMPH, GMPK
REAL(KIND=STND) :: XNKH, XNHK, QHRK, HKN, HPK, HKRN
REAL(KIND=STND) :: BTNCKH, BTNCHK, BTPDKH, BTPDHK
REAL(KIND=STND), PARAMETER :: PI = 3.14159265358979323844E0_STND, &
TPI = 2*PI, ONE = 1
IF ( NU < 1 ) THEN
BVT = MVBVU( -DH, -DK, R )
ELSE
SQNU = NU
SQNU = SQRT( SQNU )
ORS = 1 - R*R
HRK = DH - R*DK
KRH = DK - R*DH
IF ( ABS(HRK) + ORS > 0 ) THEN
XNHK = HRK**2/( HRK**2 + ORS*( NU + DK**2 ) )
XNKH = KRH**2/( KRH**2 + ORS*( NU + DH**2 ) )
ELSE
XNHK = 0
XNKH = 0
END IF
HS = SIGN( ONE, DH - R*DK )
KS = SIGN( ONE, DK - R*DH )
IF ( MODULO( NU, 2 ) == 0 ) THEN
BVT = ATAN2( SQRT(ORS), -R )/TPI
GMPH = DH/SQRT( 16*( NU + DH**2 ) )
GMPK = DK/SQRT( 16*( NU + DK**2 ) )
BTNCKH = 2*ATAN2( SQRT( XNKH ), SQRT( 1 - XNKH ) )/PI
BTPDKH = 2*SQRT( XNKH*( 1 - XNKH ) )/PI
BTNCHK = 2*ATAN2( SQRT( XNHK ), SQRT( 1 - XNHK ) )/PI
BTPDHK = 2*SQRT( XNHK*( 1 - XNHK ) )/PI
DO J = 1, NU/2
BVT = BVT + GMPH*( 1 + KS*BTNCKH )
BVT = BVT + GMPK*( 1 + HS*BTNCHK )
BTNCKH = BTNCKH + BTPDKH
BTPDKH = 2*J*BTPDKH*( 1 - XNKH )/( 2*J + 1 )
BTNCHK = BTNCHK + BTPDHK
BTPDHK = 2*J*BTPDHK*( 1 - XNHK )/( 2*J + 1 )
GMPH = GMPH*( 2*J - 1 )/( 2*J*( 1 + DH**2/NU ) )
GMPK = GMPK*( 2*J - 1 )/( 2*J*( 1 + DK**2/NU ) )
END DO
ELSE
QHRK = SQRT( DH**2 + DK**2 - 2*R*DH*DK + NU*ORS )
HKRN = DH*DK + R*NU
HKN = DH*DK - NU
HPK = DH + DK
BVT = ATAN2(-SQNU*(HKN*QHRK+HPK*HKRN),HKN*HKRN-NU*HPK*QHRK)/TPI
IF ( BVT < -1E-15_STND ) BVT = BVT + 1
GMPH = DH/( TPI*SQNU*( 1 + DH**2/NU ) )
GMPK = DK/( TPI*SQNU*( 1 + DK**2/NU ) )
BTNCKH = SQRT( XNKH )
BTPDKH = BTNCKH
BTNCHK = SQRT( XNHK )
BTPDHK = BTNCHK
DO J = 1, ( NU - 1 )/2
BVT = BVT + GMPH*( 1 + KS*BTNCKH )
BVT = BVT + GMPK*( 1 + HS*BTNCHK )
BTPDKH = ( 2*J - 1 )*BTPDKH*( 1 - XNKH )/( 2*J )
BTNCKH = BTNCKH + BTPDKH
BTPDHK = ( 2*J - 1 )*BTPDHK*( 1 - XNHK )/( 2*J )
BTNCHK = BTNCHK + BTPDHK
GMPH = 2*J*GMPH/( ( 2*J + 1 )*( 1 + DH**2/NU ) )
GMPK = 2*J*GMPK/( ( 2*J + 1 )*( 1 + DK**2/NU ) )
END DO
END IF
END IF
!
END FUNCTION MVBVTL
!
FUNCTION MVBVU( SH, SK, R ) RESULT(BVN)
!
! A function for computing bivariate normal probabilities.
!
! Yihong Ge
! and
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, WA 99164-3113
! Email : alangenz@wsu.edu
!
! MVBVU - calculate the probability that X is larger than SH and Y is
! larger than SK.
!
! Parameters
!
! SH REAL, integration limit
! SK REAL, integration limit
! R REAL, correlation coefficient
! LG INTEGER, number of Gauss Rule Points and Weights
!
REAL(KIND=STND), INTENT(IN) :: SH, SK, R
REAL(KIND=STND) :: BVN
!
INTEGER :: I, LG, NG
REAL(KIND=STND), PARAMETER :: ZERO = 0, &
TWOPI = 6.2831853071795864769252_STND
REAL(KIND=STND) :: AS, A, B, C, D, RS, XS, SN, ASR, H, K, BS, HS, HK
REAL(KIND=STND), DIMENSION(10) :: X, W
! Gauss Legendre Points and Weights, N = 6
REAL(KIND=STND), DIMENSION(3), PARAMETER :: W1 = (/ &
0.1713244923791705E+00_STND, &
0.3607615730481384E+00_STND, 0.4679139345726904E+00_STND /)
REAL(KIND=STND), DIMENSION(3), PARAMETER :: X1 = (/ &
-0.9324695142031522E+00_STND, &
-0.6612093864662647E+00_STND, -0.2386191860831970E+00_STND /)
! Gauss Legendre Points and Weights, N = 12
REAL(KIND=STND), DIMENSION(6), PARAMETER :: W2 = (/ &
0.4717533638651177E-01_STND, 0.1069393259953183E+00_STND, &
0.1600783285433464E+00_STND, 0.2031674267230659E+00_STND, &
0.2334925365383547E+00_STND, 0.2491470458134029E+00_STND /)
REAL(KIND=STND), DIMENSION(6), PARAMETER :: X2 = (/ &
-0.9815606342467191E+00_STND, -0.9041172563704750E+00_STND, &
-0.7699026741943050E+00_STND, -0.5873179542866171E+00_STND, &
-0.3678314989981802E+00_STND, -0.1252334085114692E+00_STND /)
! Gauss Legendre Points and Weights, N = 20
REAL(KIND=STND), DIMENSION(10), PARAMETER :: W3 = (/ &
0.1761400713915212E-01_STND, 0.4060142980038694E-01_STND, &
0.6267204833410906E-01_STND, 0.8327674157670475E-01_STND, &
0.1019301198172404E+00_STND, 0.1181945319615184E+00_STND, &
0.1316886384491766E+00_STND, 0.1420961093183821E+00_STND, &
0.1491729864726037E+00_STND, 0.1527533871307259E+00_STND /)
REAL(KIND=STND), DIMENSION(10), PARAMETER :: X3 = (/ &
-0.9931285991850949E+00_STND, -0.9639719272779138E+00_STND, &
-0.9122344282513259E+00_STND, -0.8391169718222188E+00_STND, &
-0.7463319064601508E+00_STND, -0.6360536807265150E+00_STND, &
-0.5108670019508271E+00_STND, -0.3737060887154196E+00_STND, &
-0.2277858511416451E+00_STND, -0.7652652113349733E-01_STND /)
!
IF ( ABS(R) < 3E-1_STND ) THEN
NG = 1
LG = 3
X(1:LG) = X1
W(1:LG) = W1
ELSE IF ( ABS(R) < 75E-2_STND ) THEN
NG = 2
LG = 6
X(1:LG) = X2
W(1:LG) = W2
ELSE
NG = 3
LG = 10
X(1:LG) = X3
W(1:LG) = W3
ENDIF
H = SH
K = SK
HK = H*K
BVN = 0
IF ( ABS(R) < 925E-3_STND ) THEN
HS = ( H*H + K*K )/2
ASR = ASIN(R)
DO I = 1, LG
SN = SIN( ASR*( 1 + X(I) )/2 )
BVN = BVN + W(I)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
SN = SIN( ASR*( 1 - X(I) )/2 )
BVN = BVN + W(I)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
END DO
BVN = BVN*ASR/(2*TWOPI) + MVPHI(-H)*MVPHI(-K)
ELSE
IF ( R < 0 ) THEN
K = -K
HK = -HK
ENDIF
IF ( ABS(R) < 1 ) THEN
AS = ( 1 - R )*( 1 + R )
A = SQRT(AS)
BS = ( H - K )**2
C = ( 4 - HK )/8
D = ( 12 - HK )/16
BVN = A*EXP( -(BS/AS + HK)/2 ) &
*( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
IF ( HK > -160 ) THEN
B = SQRT(BS)
BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVPHI(-B/A)*B &
*( 1 - C*BS*( 1 - D*BS/5 )/3 )
ENDIF
A = A/2
DO I = 1, LG
XS = ( A*( X(I) + 1 ) )**2
RS = SQRT( 1 - XS )
BVN = BVN + A*W(I)* &
( EXP( -BS/(2*XS) - HK/(1+RS) )/RS &
- EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
XS = AS*( -X(I) + 1 )**2/4
RS = SQRT( 1 - XS )
BVN = BVN + A*W(I)*EXP( -(BS/XS + HK)/2 ) &
*( EXP( -HK*XS/(2*(1+RS)**2) )/RS &
- ( 1 + C*XS*( 1 + D*XS ) ) )
END DO
BVN = -BVN/TWOPI
ENDIF
IF ( R > 0 ) THEN
BVN = BVN + MVPHI( -MAX( H, K ) )
ELSE
BVN = -BVN
IF ( K > H ) THEN
IF ( H < 0 ) THEN
BVN = BVN + MVPHI(K) - MVPHI(H)
ELSE
BVN = BVN + MVPHI(-H) - MVPHI(-K)
ENDIF
ENDIF
ENDIF
END IF
!
END FUNCTION MVBVU
!
FUNCTION MVBVTC( NU, L, U, INFIN, RHO ) RESULT(BVTC)
!
! A function for computing complementary bivariate normal and t
! probabilities.
!
! Parameters
!
! NU INTEGER degrees of freedom parameter.
! L REAL, array of lower integration limits.
! U REAL, array of upper integration limits.
! INFIN INTEGER, array of integration limits flags:
! if INFIN(1) INFIN(2), then MVBVTC computes
! 0 0 P( X>U(1), Y>U(2) ),
! 1 0 P( X<L(1), Y>U(2) ),
! 0 1 P( X>U(1), Y<L(2) ),
! 1 1 P( X<L(1), Y<L(2) ),
! 2 0 P( X>U(1), Y>U(2) ) + P( X<L(1), Y>U(2) ),
! 2 1 P( X>U(1), Y<L(2) ) + P( X<L(1), Y<L(2) ),
! 0 2 P( X>U(1), Y>U(2) ) + P( X>U(1), Y<L(2) ),
! 1 2 P( X<L(1), Y>U(2) ) + P( X<L(1), Y<L(2) ),
! 2 2 P( X>U(1), Y<L(2) ) + P( X<L(1), Y<L(2) ),
! + P( X>U(1), Y>U(2) ) + P( X<L(1), Y>U(2) ).
!
! RHO REAL, correlation coefficient.
!
INTEGER, INTENT(IN) :: NU
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: L, U
INTEGER, DIMENSION(:), INTENT(IN) :: INFIN
REAL(KIND=STND), INTENT(IN) :: RHO
REAL(KIND=STND) :: BVTC
!
! Local Variables
!
REAL(KIND=STND), DIMENSION(2) :: LW, UP
REAL(KIND=STND) :: B
INTEGER, DIMENSION(2) :: INF
INTEGER :: I
!
IF ( INFIN(1) < 0 .AND. INFIN(2) < 0 ) THEN
BVTC = 0
ELSE IF ( INFIN(1) < 0 ) THEN
BVTC = MVSTDC( NU, L(2), U(2), INFIN(2) )
ELSE IF ( INFIN(2) < 0 ) THEN
BVTC = MVSTDC( NU, L(1), U(1), INFIN(1) )
ELSE
DO I = 1, 2
IF ( MODULO( INFIN(I), 2 ) == 0 ) THEN
INF(I) = 1
LW(I) = U(I)
ELSE
INF(I) = 0
UP(I) = L(I)
END IF
END DO
B = MVBVT( NU, LW, UP, INF, RHO )
DO I = 1, 2
IF ( INFIN(I) == 2 ) THEN
INF(I) = 0
UP(I) = L(I)
B = B + MVBVT( NU, LW, UP, INF, RHO )
END IF
END DO
IF ( INFIN(1) == 2 .AND. INFIN(2) == 2 ) THEN
INF(1) = 1
LW(1) = U(1)
B = B + MVBVT( NU, LW, UP, INF, RHO )
END IF
BVTC = B
END IF
!
END FUNCTION MVBVTC
!
FUNCTION MVTVT( NU, LOWER, UPPER, INFIN, CR ) RESULT(TVNT)
!
! A function for computing trivariate normal and t probabilities.
!
! Parameters
!
! NU INTEGER degrees of freedom parameter; NU < 1 gives normal case.
! 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)].
! CR REAL, correlation coefficients.
!
INTEGER, INTENT(IN) :: NU
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: LOWER, UPPER
INTEGER, DIMENSION(:), INTENT(IN) :: INFIN
REAL(KIND=STND), DIMENSION(:), INTENT(IN) :: CR
REAL(KIND=STND) :: TVNT
!
INTEGER, DIMENSION(2) :: INFT
REAL(KIND=STND) :: L1, L2, U1, U2
REAL(KIND=STND), DIMENSION(2) :: LT, UT
REAL(KIND=STND), DIMENSION(3) :: LW, UP
!
IF ( INFIN(1) < 0 ) THEN
! -1 0 0
LT = LOWER(2:3)
UT = UPPER(2:3)
INFT = INFIN(2:3)
TVNT = MVBVT( NU, LT, UT, INFT, CR(3) )
ELSE IF ( INFIN(2) < 0 ) THEN
! *-1 *
LT = (/ LOWER(1), LOWER(3) /)
UT = (/ UPPER(1), UPPER(3) /)
INFT = (/ INFIN(1), INFIN(3) /)
TVNT = MVBVT( NU, LT, UT, INFT, CR(2) )
ELSE IF ( INFIN(3) < 0 ) THEN
! * *-1
LT = LOWER(1:2)
UT = UPPER(1:2)
INFT = INFIN(1:2)
TVNT = MVBVT( NU, LT, UT, INFT, CR(1) )
ELSE IF ( INFIN(1) == 0 .AND. INFIN(2) == 0 .AND. INFIN(3) == 0 ) THEN
! 0 0 0
TVNT = MVTVTL( NU, UPPER, CR )
ELSE IF ( INFIN(1