! Programs available in module RINDMOD : ! ! 1) setConstants ! 2) RINDD ! ! SETCONSTANTS set member variables controlling the performance of RINDD ! ! CALL setConstants(method,xcscale,abseps,releps,coveps,maxpts,minpts,nit,xcutoff,Nc1c2) ! ! METHOD = INTEGER defining the SCIS integration method ! 1 Integrate by SADAPT for Ndim<9 and by KRBVRC otherwise ! 2 Integrate by SADAPT for Ndim<20 and by KRBVRC otherwise ! 3 Integrate by KRBVRC by Genz (1993) (Fast Ndim<101) (default) ! 4 Integrate by KROBOV by Genz (1992) (Fast Ndim<101) ! 5 Integrate by RCRUDE by Genz (1992) (Slow Ndim<1001) ! 6 Integrate by SOBNIED (Fast Ndim<1041) ! 7 Integrate by DKBVRC by Genz (2003) (Fast Ndim<1001) ! ! XCSCALE = REAL to scale the conditinal probability density, i.e., ! f_{Xc} = exp(-0.5*Xc*inv(Sxc)*Xc + XcScale) (default XcScale =0) ! ABSEPS = REAL absolute error tolerance. (default 0) ! RELEPS = REAL relative error tolerance. (default 1e-3) ! COVEPS = REAL error tolerance in Cholesky factorization (default 1e-13) ! 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. ! (Only for METHOD~=0) (default 40000) ! MINPTS = INTEGER, minimum number of function values allowed. ! (Only for METHOD~=0) (default 0) ! NIT = INTEGER, maximum number of Xt variables to integrate ! This parameter can be used to limit the time. ! If NIT is less than the rank of the covariance matrix, ! the returned result is a upper bound for the true value ! of the integral. (default 1000) ! XCUTOFF = REAL cut off value where the marginal normal ! distribution is truncated. (Depends on requested ! accuracy. A value between 4 and 5 is reasonable.) ! NC1C2 = number of times to use the regression equation to restrict ! integration area. Nc1c2 = 1,2 is recommended. (default 2) ! ! !RIND computes E[Jacobian*Indicator|Condition]*f_{Xc}(xc(:,ix)) ! ! where ! "Indicator" = I{ H_lo(i) < X(i) < H_up(i), i=1:Nt+Nd } ! "Jacobian" = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is ! "Jacobian" = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)| ! "condition" = Xc=xc(:,ix), ix=1,...,Nx. ! X = [Xt; Xd ;Xc], a stochastic vector of Multivariate Gaussian ! variables where Xt,Xd and Xc have the length Nt, Nd and Nc, ! respectively. ! (Recommended limitations Nx, Nt<101, Nd<7 and NIT,Nc<11) ! (RIND = Random Integration N Dimensions) ! !CALL RINDD(E,err,terr,S,m,xc,Nt,indI,Blo,Bup,INFIN); ! ! E = expectation/density as explained above size 1 x Nx (out) ! ERR = estimated sampling error size 1 x Nx (out) ! TERR = estimated truncation error size 1 x Nx (out) ! S = Covariance matrix of X=[Xt;Xd;Xc] size N x N (N=Nt+Nd+Nc) (in) ! m = the expectation of X=[Xt;Xd;Xc] size N x 1 (in) ! xc = values to condition on size Nc x Nx (in) ! indI = vector of indices to the different barriers in the (in) ! indicator function, length NI, where NI = Nb+1 ! (NB! restriction indI(1)=0, indI(NI)=Nt+Nd ) ! Blo,Bup = Lower and upper barrier coefficients used to compute the (in) ! integration limits A and B, respectively. ! size Mb x Nb. If Mb=0, ! IF INFIN(j)~=0, A(i)=Blo(1,j)+Blo(2:Mb,j).'*xc(1:Mb-1,ix), ! IF INFIN(j)~=1, B(i)=Bup(1,j)+Bup(2:Mb,j).'*xc(1:Mb-1,ix), ! ! where i=indI(j-1)+1:indI(j), j=1:NI-1, ix=1:Nx ! Thus the integration limits may change with the conditional ! variables. !Example: ! The indices, indI=[0 3 5 6], and coefficients Blo=[0 0 -1], ! Bup=[0 0 5], INFIN=[0 1 2] ! means that A = [-inf -inf -inf 0 0 -1] B = [0 0 0 inf inf 5] ! ! ! (Recommended limitations Nx,Nt<101, Nd<7 and Nc<11) ! Also note that the size information have to be transferred to RINDD ! through the input arguments E,S,m,Nt,IndI,Blo,Bup and INFIN ! ! For further description see the modules ! ! References ! Podgorski et al. (2000) ! "Exact distributions for apparent waves in irregular seas" ! Ocean Engineering, Vol 27, no 1, pp979-1016. (RINDD) ! ! R. Ambartzumian, A. Der Kiureghian, V. Ohanian and H. ! Sukiasian (1998) ! "Multinormal probabilities by sequential conditioned ! importance sampling: theory and application" (MVNFUN) ! Probabilistic Engineering Mechanics, Vol. 13, No 4. pp 299-308 ! ! Alan Genz (1992) ! 'Numerical Computation of Multivariate Normal Probabilites' (MVNFUN) ! J. computational Graphical Statistics, Vol.1, pp 141--149 ! ! Alan Genz and Koon-Shing Kwong (2000?) ! 'Numerical Evaluation of Singular Multivariate Normal Distributions' (MVNFUN,COVSRT) ! Computational Statistics and Data analysis ! ! ! P. A. Brodtkorb (2004), (RINDD, MVNFUN, COVSRT) ! Numerical evaluation of multinormal expectations ! In Lund university report series ! and in the Dr.Ing thesis: ! The probability of Occurrence of dangerous Wave Situations at Sea. ! Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, ! Trondheim, Norway. ! Tested on: DIGITAL UNIX Fortran90 compiler ! PC pentium II with Lahey Fortran90 compiler ! Solaris with SunSoft F90 compiler Version 1.0.1.0 (21229283) ! History: ! Revised pab aug. 2009 ! -renamed from rind2007 to rindmod ! Revised pab July 2007 ! - separated the absolute error into ERR and TERR. ! - renamed from alanpab24 -> rind2007 ! revised pab 23may2004 ! RIND module totally rewritten according to the last reference. MODULE GLOBALCONST ! global constants IMPLICIT NONE DOUBLE PRECISION, PARAMETER :: gSQTWPI1= 0.39894228040143D0 !=1/sqrt(2*pi) DOUBLE PRECISION, PARAMETER :: gSQPI1 = 0.56418958354776D0 !=1/sqrt(pi) DOUBLE PRECISION, PARAMETER :: gSQPI = 1.77245385090552D0 !=sqrt(pi) DOUBLE PRECISION, PARAMETER :: gSQTW = 1.41421356237310D0 !=sqrt(2) DOUBLE PRECISION, PARAMETER :: gSQTW1 = 0.70710678118655D0 !=1/sqrt(2) DOUBLE PRECISION, PARAMETER :: gPI1 = 0.31830988618379D0 !=1/pi DOUBLE PRECISION, PARAMETER :: gPI = 3.14159265358979D0 !=pi DOUBLE PRECISION, PARAMETER :: gTWPI = 6.28318530717958D0 !=2*pi DOUBLE PRECISION, PARAMETER :: gSQTWPI = 2.50662827463100D0 !=sqrt(2*pi) DOUBLE PRECISION, PARAMETER :: gONE = 1.D0 DOUBLE PRECISION, PARAMETER :: gTWO = 2.D0 DOUBLE PRECISION, PARAMETER :: gHALF = 0.5D0 DOUBLE PRECISION, PARAMETER :: gZERO = 0.D0 DOUBLE PRECISION, PARAMETER :: gINFINITY = 37.D0 ! SQRT(-gTWO*LOG(1.D+12*TINY(gONE))) ! Set gINFINITY (infinity). ! Such that EXP(-2.x^2) > 10^(12) times TINY ! SAVE gINFINITY END MODULE GLOBALCONST MODULE RINDMOD USE GLOBALCONST ! USE PRINTMOD ! used for debugging only IMPLICIT NONE PRIVATE PUBLIC :: RINDD, SetConstants PUBLIC :: mCovEps, mAbsEps,mRelEps, mXcutOff, mXcScale PUBLIC :: mNc1c2, mNIT, mMaxPts,mMinPts, mMethod, mSmall private :: preInit private :: initIntegrand private :: initfun,mvnfun,cvsrtxc,covsrt1,covsrt,rcscale,rcswap private :: cleanUp INTERFACE RINDD MODULE PROCEDURE RINDD END INTERFACE INTERFACE SetConstants MODULE PROCEDURE SetConstants END INTERFACE ! mInfinity = what is considered as infinite value in FI ! mFxcEpss = if fxc is less, do not compute E(...|Xc) ! mXcEps2 = if any Var(Xc(j)|Xc(1),...,Xc(j-1)) <= XCEPS2 then return NAN double precision, parameter :: mInfinity = 8.25d0 ! 37.0d0 double precision, parameter :: mFxcEpss = 1.0D-20 double precision, save :: mXcEps2 = 2.3d-16 ! Constants defining accuracy of integration: ! mCovEps = termination criteria for Cholesky decomposition ! mAbsEps = requested absolute tolerance ! mRelEps = requested relative tolerance ! mXcutOff = truncation value to c1c2 ! mXcScale = scale factor in the exponential (in order to avoid overflow) ! mNc1c2 = number of times to use function c1c2, i.e.,regression ! equation to restrict integration area. ! mNIT = maximum number of Xt variables to integrate ! mMethod = integration method: ! 1 Integrate all by SADAPT if NDIM<9 otherwise by KRBVRC (default) ! 2 Integrate all by SADAPT if NDIM<19 otherwise by KRBVRC ! 3 Integrate all by KRBVRC by Genz (1998) (Fast and reliable) ! 4 Integrate all by KROBOV by Genz (1992) (Fast and reliable) ! 5 Integrate all by RCRUDE by Genz (1992) (Reliable) ! 6 Integrate all by SOBNIED by Hong and Hickernell ! 7 Integrate all by DKBVRC by Genz (2003) (Fast Ndim<1001) double precision, save :: mCovEps = 1.0d-10 double precision, save :: mAbsEps = 0.01d0 double precision, save :: mRelEps = 0.01d0 double precision, save :: mXcutOff = 5.d0 double precision, save :: mXcScale = 0.0d0 integer, save :: mNc1c2 = 2 integer, save :: mNIT = 1000 integer, save :: mMaxPts = 40000 integer, save :: mMinPts = 0 integer, save :: mMethod = 3 ! Integrand variables: ! mBIG = Cholesky Factor/Covariance matrix: ! Upper triangular part is the cholesky factor ! Lower triangular part contains the conditional ! standarddeviations ! (mBIG2 is only used if mNx>1) ! mCDI = Cholesky DIagonal elements ! mA,mB = Integration limits ! mINFI = integrationi limit flags ! mCm = conditional mean ! mINFIXt, ! mINFIXd = # redundant variables of Xt and Xd, ! respectively ! mIndex1, ! mIndex2 = indices to the variables original place. Size Ntdc ! xedni = indices to the variables new place. Size Ntdc ! mNt = # Xt variables ! mNd = # Xd variables ! mNc = # Xc variables ! mNtd = mNt + mNd ! mNtdc = mNt + mNd + mNc ! mNx = # different integration limits double precision,allocatable, dimension(:,:) :: mBIG,mBIG2 double precision,allocatable, dimension(:) :: mA,mB,mCDI,mCm INTEGER, DIMENSION(:),ALLOCATABLE :: mInfi,mIndex1,mIndex2,mXedni INTEGER,SAVE :: mNt,mNd,mNc,mNtdc, mNtd, mNx ! Size information INTEGER,SAVE :: mInfiXt,mInfiXd logical,save :: mInitIntegrandCalled = .FALSE. DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: mCDIXd, mCmXd DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: mXd, mXc, mY double precision, save :: mSmall = 2.3d-16 ! variables set in initfun and used in mvnfun: INTEGER, PRIVATE :: mI0,mNdleftN0 DOUBLE PRECISION, PRIVATE :: mE1,mD1, mVAL0 contains subroutine setConstants(method,xcscale,abseps,releps,coveps, & maxpts,minpts,nit,xcutoff,Nc1c2) double precision, optional, intent(in) :: xcscale,abseps,releps $ ,coveps, xcutoff integer, optional,intent(in) :: method,nit,maxpts,minpts,Nc1c2 double precision, parameter :: one = 1.0d0 mSmall = spacing(one) if (present(method)) mMethod = method if (present(xcscale)) mXcScale = xcscale if (present(abseps)) mAbsEps = max(abseps,mSmall) if (present(releps)) mRelEps = max(releps,0.0d0) if (present(coveps)) mCovEps = max(coveps,1d-12) if (present(maxpts)) mMaxPts = maxpts if (present(minpts)) mMinPts = minpts if (present(nit)) mNit = nit if (present(xcutOff)) mXcutOff = xCutOff if (present(Nc1c2)) mNc1c2 = max(Nc1c2,1) ! print *, 'method=', mMethod ! print *, 'xcscale=', mXcScale ! print *, 'abseps=', mAbsEps ! print *, 'releps=', mRelEps ! print *, 'coveps=', mCovEps ! print *, 'maxpts=', mMaxPts ! print *, 'minpts=', mMinPts ! print *, 'nit=', mNit ! print *, 'xcutOff=', mXcutOff ! print *, 'Nc1c2=', mNc1c2 end subroutine setConstants subroutine preInit(BIG,Xc,Nt,inform) double precision,dimension(:,:), intent(in) :: BIG double precision,dimension(:,:), intent(in) :: Xc integer, intent(in) :: Nt integer, intent(out) :: inform ! Local variables integer :: I,J inform = 0 mInitIntegrandCalled = .FALSE. ! Find the size information !~~~~~~~~~~~~~~~~~~~~~~~~~~ mNt = Nt mNc = SIZE( Xc, dim = 1 ) mNx = MAX( SIZE( Xc, dim = 2), 1 ) mNtdc = SIZE( BIG, dim = 1 ) ! make sure it does not exceed Ntdc-Nc IF (mNt+mNc.GT.mNtdc) mNt = mNtdc - mNc mNd = mNtdc-mNt-mNc mNtd = mNt+mNd IF (mNd < 0) THEN ! PRINT *,'RIND Nt,Nd,Nc,Ntdc=',Nt,Nd,Nc,Ntdc ! Size information inconsistent inform = 3 return ENDIF ! PRINT *,'Nt Nd Nc Ntd Ntdc,',Nt, Nd, Nc, Ntd, Ntdc ! ALLOCATION !~~~~~~~~~~~~ IF (mNd>0) THEN ALLOCATE(mXd(mNd),mCmXd(mNd),mCDIXd(mNd)) mCmXd(:) = gZERO mCDIXd(:) = gZERO mxd(:) = gZERO END IF ALLOCATE(mBIG(mNtdc,mNtdc),mCm(mNtdc),mY(mNtd)) ALLOCATE(mIndex1(mNtdc),mA(mNtd),mB(mNtd),mINFI(mNtd),mXc(mNc)) ALLOCATE(mCDI(mNtd),mXedni(mNtdc),mIndex2(mNtdc)) ! Initialization !~~~~~~~~~~~~~~~~~~~~~ ! Copy upper triangular of input matrix, only. do i = 1,mNtdc mBIG(1:i,i) = BIG(1:i,i) end do mIndex2 = (/(J,J=1,mNtdc)/) ! CALL mexprintf('BIG Before CovsrtXc'//CHAR(10)) ! CALL ECHO(BIG) ! sort BIG by decreasing cond. variance for Xc CALL CVSRTXC(mNt,mNd,mBIG,mIndex2,INFORM) ! CALL mexprintf('BIG after CovsrtXc'//CHAR(10)) ! CALL ECHO(BIG) IF (INFORM.GT.0) return ! degenerate case exit VALS=0 for all ! (should perhaps return NaN instead??) DO I=mNtdc,1,-1 J = mIndex2(I) ! covariance matrix according to index2 mXedni(J) = I END DO IF (mNx>1) THEN ALLOCATE(mBIG2(mNtdc,mNtdc)) do i = 1,mNtdc mBIG2(1:i,i) = mBIG(1:i,i) !Copy input matrix end do ENDIF return end subroutine preInit subroutine initIntegrand(ix,Xc,Ex,indI,Blo,Bup,INFIN, & fxc,value,abserr,NDIM,inform) integer, intent(in) :: ix ! integrand number double precision, dimension(:),intent(in) :: Ex double precision, dimension(:,:), intent(in) :: Xc,Blo,Bup integer, dimension(:), intent(in) :: indI,INFIN double precision, intent(out) :: fxc,value,abserr integer, intent(out) :: NDIM, inform ! Locals DOUBLE PRECISION :: SQ0,xx,quant integer :: I,J inform = 0 NDIM = 0 VALUE = gZERO fxc = gONE abserr = mSmall IF (mInitIntegrandCalled) then do i = 1,mNtdc mBIG(1:i,i) = mBIG2(1:i,i) !Copy input matrix end do else mInitIntegrandCalled = .TRUE. endif ! Set the original means of the variables mCm(:) = Ex(mIndex2(1:mNtdc)) ! Cm(1:Ntdc) =Ex (index1(1:Ntdc)) IF (mNc>0) THEN mXc(:) = Xc(:,ix) !mXc(1:Nc) = Xc(1:Nc,ix) QUANT = DBLE(mNc)*LOG(gSQTWPI1) I = mNtdc DO J = 1, mNc ! Iterative conditioning on the last Nc variables SQ0 = mBIG(I,I) ! SQRT(Var(X(i)|X(i+1),X(i+2),...,X(Ntdc))) xx = (mXc(mIndex2(I) - mNtd) - mCm(I))/SQ0 !Trick to calculate !fxc = fxc*SQTWPI1*EXP(-0.5*(XX**2))/SQ0 QUANT = QUANT - gHALF*xx*xx - LOG(SQ0) ! conditional mean (expectation) ! E(X(1:i-1)|X(i),X(i+1),...,X(Ntdc)) mCm(1:I-1) = mCm(1:I-1) + xx*mBIG(1:I-1,I) I = I-1 ENDDO ! Calculating the ! fxc probability density for i=Ntdc-J+1, ! fXc=f(X(i)|X(i+1),X(i+2)...X(Ntdc))* ! f(X(i+1)|X(i+2)...X(Ntdc))*..*f(X(Ntdc)) fxc = EXP(QUANT+mXcScale) ! if fxc small: don't bother ! calculating it, goto end IF (fxc < mFxcEpss) then abserr = gONE inform = 1 return endif END IF ! Set integration limits mA,mB and mINFI ! NOTE: mA and mB are integration limits with mCm subtracted CALL setIntLimits(mXc,indI,Blo,Bup,INFIN,inform) if (inform>0) return mIndex1(:) = mIndex2(:) CALL COVSRT(.FALSE., mNt,mNd,mBIG,mCm,mA,mB,mINFI, & mINDEX1,mINFIXt,mINFIXd,NDIM,mY,mCDI) CALL INITFUN(VALUE,abserr,INFORM) ! IF INFORM>0 : degenerate case: ! Integral can be calculated excactly, ie. ! mean of deterministic variables outside the barriers, ! or NDIM = 1 return end subroutine initIntegrand subroutine cleanUp ! Deallocate all work arrays and vectors IF (ALLOCATED(mXc)) DEALLOCATE(mXc) IF (ALLOCATED(mXd)) DEALLOCATE(mXd) IF (ALLOCATED(mCm)) DEALLOCATE(mCm) IF (ALLOCATED(mBIG2)) DEALLOCATE(mBIG2) IF (ALLOCATED(mBIG)) DEALLOCATE(mBIG) IF (ALLOCATED(mIndex2)) DEALLOCATE(mIndex2) IF (ALLOCATED(mIndex1)) DEALLOCATE(mIndex1) IF (ALLOCATED(mXedni)) DEALLOCATE(mXedni) IF (ALLOCATED(mA)) DEALLOCATE(mA) IF (ALLOCATED(mB)) DEALLOCATE(mB) IF (ALLOCATED(mY)) DEALLOCATE(mY) IF (ALLOCATED(mCDI)) DEALLOCATE(mCDI) IF (ALLOCATED(mCDIXd)) DEALLOCATE(mCDIXd) IF (ALLOCATED(mCmXd)) DEALLOCATE(mCmXd) IF (ALLOCATED(mINFI)) DEALLOCATE(mINFI) end subroutine cleanUp function integrandBound(I0,N,Y,FINY) result (bound1) use FIMOD integer, intent(in) :: I0,N,FINY double precision, intent(in) :: Y double precision :: bound1 ! locals integer :: I,IK,FINA, FINB double precision :: AI,BI,D1,E1 double precision :: TMP ! Computes the upper bound for the intgrand bound1 = gzero if (FINY<1) return FINA = 0 FINB = 0 IK = 2 DO I = I0, N ! E(Y(I) | Y(1))/STD(Y(IK)|Y(1)) TMP = mBIG(IK-1,I)*Y IF (mINFI(I) > -1) then ! May have infinite int. Limits if Nd>0 IF ( mINFI(I) .NE. 0 ) THEN IF ( FINA .EQ. 1 ) THEN AI = MAX( AI, mA(I) - tmp ) ELSE AI = mA(I) - tmp FINA = 1 END IF END IF IF ( mINFI(I) .NE. 1 ) THEN IF ( FINB .EQ. 1 ) THEN BI = MIN( BI, mB(I) - tmp) ELSE BI = mB(I) - tmp FINB = 1 END IF END IF endif IF (I.EQ.N.OR.mBIG(IK+1,I+1)>gZERO) THEN CALL MVNLMS( AI, BI,2*FINA+FINB-1, D1, E1 ) IF (D1Nt then variable no. I is one of the Xd ! variables otherwise it is one of Xt. !PRINT *,'Mvnfun,ndim',Ndim INFORM = 0 VALUE = gZERO abserr = max(mCovEps , 6.0d0*mSmall) mVAL0 = gONE mNdleftN0 = mNd ! Counter for number of Xd variables left mI0 = 0 FINA = 0 FINB = 0 N = mNt + mNd - mINFIXt - mINFIXd-1 IF (mINFIXt+mINFIXd > 0) THEN ! CHCKLIM Check if the conditional mean Cm = E(Xt,Xd|Xc) for the ! deterministic variables are between the barriers, i.e., ! A=Hlo-Cm< 0 -1) THEN IF ((mINFI(I).NE.0).AND.(mAbsEps < mA(I))) GOTO 200 IF ((mINFI(I).NE.1).AND.(mB(I) < -mAbsEps )) GOTO 200 ENDIF ENDDO IF (mINFIXd>0) THEN ! Redundant variables of Xd: replace Xd with the mean I = mNt + mNd !-INFIS J = mNdleftN0-mINFIXd DO WHILE (mNdleftN0>J) isXd = (mNt < mIndex1(I)) IF (isXd) THEN mXd (mNdleftN0) = mCm (I) mNdleftN0 = mNdleftN0-1 END IF I = I-1 ENDDO ENDIF IF (N+1 < 1) THEN ! Degenerate case, No relevant variables left to integrate ! Print *,'rind ndim1',Ndim1 IF (mNd>0) THEN VALUE = jacob (mXd,mXc) ! jacobian of xd,xc ELSE VALUE = gONE END IF GOTO 200 ENDIF ENDIF IF (mNIT<=100) THEN xCut = mXcutOff J = 1 DO I = 2, N+1 IF (mBIG(J+1,I)>gZERO) THEN J = J + 1 ELSE ! Add xCut std to deterministic variables to get an upper ! bound for integral mA(I) = mA(I) - xCut * mBIG(I,J) mB(I) = mB(I) + xCut * mBIG(I,J) ENDIF END DO ELSE xCut = gZERO ENDIF NdleftO = mNdleftN0 useC1C2 = (1<=mNc1c2) DO I = 1, N+1 IF (mINFI(I) > -1) then ! May have infinite int. Limits if Nd>0 IF ( mINFI(I) .NE. 0 ) THEN IF ( FINA .EQ. 1 ) THEN AI = MAX( AI, mA(I) ) ELSE AI = mA(I) FINA = 1 END IF END IF IF ( mINFI(I) .NE. 1 ) THEN IF ( FINB .EQ. 1 ) THEN BI = MIN( BI, mB(I) ) ELSE BI = mB(I) FINB = 1 END IF END IF endif isXd = (mINDEX1(I)>mNt) IF (isXd) THEN ! Save the mean for Xd mCmXd(mNdleftN0) = mCm(I) mCDIXd(mNdleftN0) = mCDI(I) mNdleftN0 = mNdleftN0-1 END IF IF (I.EQ.N+1.OR.mBIG(2,I+1)>gZERO) THEN IF (useC1C2.AND.I=E0) GOTO 200 CALL C1C2(I+1,N+1,1,mA,mB,mINFI,mY,mBIG,AI,BI,FINA,FINB) CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 ) IF (mD1>=mE1) GOTO 200 maxTruncError = FI(-ABS(mXcutOff))*dble(mNc1c2) upError = abs(E0-mE1) loError = abs(D0-mD1) if (upError>mSmall) then upError = upError*integrandBound(I+1,N+1,BI,FINB) endif if (loError>mSmall) then loError = loError*integrandBound(I+1,N+1,AI,FINA) endif abserr = abserr + min(upError + loError,maxTruncError) !CALL printvar(log10(loError+upError+msmall),'lo+up-err') ELSE CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 ) IF (mD1>=mE1) GOTO 200 ENDIF !CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 ) !IF (mD1>=mE1) GOTO 200 IF ( NdleftO<=0) THEN IF (mNd>0) mVAL0 = JACOB(mXd,mXc) SELECT CASE (I-N) CASE (1) !IF (I.EQ.N+1) THEN VALUE = (mE1-mD1)*mVAL0 abserr = abserr*mVAL0 GO TO 200 CASE (0) !ELSEIF (I.EQ.N) THEN !D1=1/sqrt(1-rho^2)=1/STD(X(I+1)|X(1)) mD1 = SQRT( gONE + mBIG(1,I+1)*mBIG(1,I+1) ) mINFI(2) = mINFI(I+1) mA(1) = AI mB(1) = BI mINFI(1) = 2*FINA+FINB-1 IF ( mINFI(2) .NE. 0 ) mA(2) = mA(I+1)/mD1 IF ( mINFI(2) .NE. 1 ) mB(2) = mB(I+1)/mD1 VALUE = BVNMVN( mA, mB,mINFI,mBIG(1,I+1)/mD1 )*mVAL0 abserr = (abserr+1.0d-14)*mVAL0 GO TO 200 CASE ( -1 ) !ELSEIF (I.EQ.N-1) THEN IF (.FALSE.) THEN ! TODO :this needs further checking! (it should work though) !1/D1= sqrt(1-r12^2) = STD(X(I+1)|X(1)) !1/E1= STD(X(I+2)|X(1)X(I+1)) !D1 = BIG(I+1,1) !E1 = BIG(I+2,2) mD1 = gONE/SQRT( gONE + mBIG(1,I+1)*mBIG(1,I+1) ) R12 = mBIG( 1, I+1 ) * mD1 if (mBIG(3,I+2)>gZERO) then mE1 = gONE/SQRT( gONE + mBIG(1,I+2)*mBIG(1,I+2) + & mBIG(2,I+2)*mBIG(2,I+2) ) R13 = mBIG( 1, I+2 ) * mE1 R23 = mBIG( 2, I+2 ) * (mE1 * mD1) + R12 * R13 else mE1 = mCDI(I+2) R13 = mBIG( 1, I+2 ) * mE1 R23 = mE1*mD1 + R12 * R13 IF ((mE1 < gZERO).AND. mINFI(I+2)>-1) THEN CALL SWAP(mA(I+2),mB(I+2)) IF (mINFI(I+2).NE. 2) mINFI(I+2) = 1-mINFI(I+2) END IF !R23 = BIG( 2, I+2 ) * (E1 * D1) + R12 * R13 endif mINFI(2) = mINFI(I+1) mINFI(3) = mINFI(I+2) mA(1) = AI mB(1) = BI mINFI(1) = 2*FINA+FINB-1 IF ( mINFI(2) .NE. 0 ) mA(2) = mA(I+1) * mD1 IF ( mINFI(2) .NE. 1 ) mB(2) = mB(I+1) * mD1 IF ( mINFI(3) .NE. 0 ) mA(3) = mA(I+2) * mE1 IF ( mINFI(3) .NE. 1 ) mB(3) = mB(I+2) * mE1 if(.false.) then CALL PRINTVECD((/R12, R13, R23 /),'R12 = ') CALL PRINTVECD((/mD1, mE1 /),'D1 = ') CALL PRINTVECD(mBIG(1,1:3),'BIG(1,1:3) = ') CALL PRINTVECD(mBIG(2,2:3),'BIG(2,2:3) = ') CALL PRINTVECD(mBIG(1:3,1),'BIG(1:3,1) = ') CALL PRINTVECD(mBIG(2:3,2),'BIG(2:3,2) = ') CALL PRINTVECD(mA(1:I+2),'A = ') CALL PRINTVECD(mB(1:I+2),'B = ') CALL PRINTVECI(mINFI(1:I+2),'INFI = ') CALL PRINTVECI(mINDEX1(1:I+2),'index1 = ') endif VALUE = TVNMVN( mA, mB,mINFI, & (/R12, R13, R23 /),1.0d-13) * mVAL0 ABSERR = (ABSERR + 1.0d-13)*mVAL0 GOTO 200 ENDIF END SELECT !ENDIF ENDIF ABSERR = mVAL0*ABSERR mVAL0 = mVAL0 * (mE1-mD1) mI0 = I RETURN ENDIF ENDDO RETURN 200 INFORM = 1 RETURN END SUBROUTINE INITFUN ! ! Integrand subroutine ! FUNCTION MVNFUN( Ndim, W ) RESULT (VAL) USE JACOBMOD USE FIMOD IMPLICIT NONE INTEGER, INTENT (IN) :: Ndim DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: W DOUBLE PRECISION :: VAL ! local variables: INTEGER :: N,I, J, FINA, FINB INTEGER :: NdleftN, NdleftO ,IK DOUBLE PRECISION :: TMP, AI, BI, DI, EI LOGICAL :: useC1C2, isXd !MVNFUN Multivariate Normal integrand function ! where the integrand is transformed from an integral ! having integration limits A and B to an ! integral having constant integration limits i.e. ! B 1 ! int jacob(xd,xc)*f(xd,xt)dxt dxd = int F2(W) dW ! A 0 ! ! W - new transformed integration variables, valid range 0..1 ! The vector must have the length Ndim returned from Covsrt ! mBIG - conditional sorted ChOlesky Factor of the covariance matrix (IN) ! mCDI - Cholesky DIagonal elements used to calculate the mean ! mCm - conditional mean of Xd and Xt given Xc, E(Xd,Xt|Xc) ! mXd - variables to the jacobian variable, need no initialization size Nd ! mXc - conditional variables (IN) ! mINDEX1 - if mINDEX1(I)>Nt then variable No. I is one of the Xd ! variables otherwise it is one of Xt !PRINT *,'Mvnfun,ndim',Ndim ! xCut = gZERO ! xCutOff N = mNt+mNd-mINFIXt-mINFIXd-1 IK = 1 ! Counter for Ndim FINA = 0 FINB = 0 NdleftN = mNdleftN0 ! Counter for number of Xd variables left VAL = mVAL0 NdleftO = mNd - mINFIXd mY(IK) = FIINV( mD1 + W(IK)*( mE1 - mD1 ) ) useC1C2 = (IK+1.LE.mNc1c2) IF (useC1C2) THEN ! Calculate the conditional mean ! E(Y(I) | Y(1),...Y(I0))/STD(Y(I)|Y(1),,,,Y(I0)) mY(mI0+1:N+1) = mBIG(IK, mI0+1:N+1)*mY(IK) ENDIF IF (NdleftO.GT.NdleftN ) THEN mXd(NdleftN+1:NdleftO) = mCmXd(NdleftN+1:NdleftO)+ & mY(IK) * mCDIXd(NdleftN+1:NdleftO) ENDIF NdleftO = NdleftN IK = 2 !=IK+1 DO I = mI0+1, N+1 IF (useC1C2) THEN TMP = mY(I) ELSE TMP = 0.d0 DO J = 1, IK-1 ! E(Y(I) | Y(1),...Y(IK-1))/STD(Y(IK)|Y(1),,,,Y(IK-1)) TMP = TMP + mBIG(J,I)*mY(J) END DO ENDIF IF (mINFI(I) < 0) GO TO 100 ! May have infinite int. Limits if Nd>0 IF ( mINFI(I) .NE. 0 ) THEN IF ( FINA .EQ. 1 ) THEN AI = MAX( AI, mA(I) - TMP) ELSE AI = mA(I) - TMP FINA = 1 END IF IF (FINB.EQ.1.AND.BI<=AI) GOTO 200 END IF IF ( mINFI(I) .NE. 1 ) THEN IF ( FINB .EQ. 1 ) THEN BI = MIN( BI, mB(I) - TMP) ELSE BI = mB(I) - TMP FINB = 1 END IF IF (FINA.EQ.1.AND.BI<=AI) GOTO 200 END IF 100 isXd = (mNt gZERO ) THEN IF (useC1C2) THEN ! Note: for J =I+1:N+1: Y(J) = conditional expectation, E(Yj|Y1,...Yk) CALL C1C2(I+1,N+1,IK,mA,mB,mINFI,mY,mBIG,AI,BI,FINA,FINB) ENDIF CALL MVNLMS( AI, BI, 2*FINA+FINB-1, DI, EI ) IF ( DI >= EI ) GO TO 200 VAL = VAL * ( EI - DI ) IF ( I <= N .OR. (NdleftN < NdleftO)) THEN mY(IK) = FIINV( DI + W(IK)*( EI - DI ) ) IF (NdleftN < NdleftO ) THEN mXd(NdleftN+1:NdleftO) = mCmXd(NdleftN+1:NdleftO)+ & mY(IK) * mCDIXd(NdleftN+1:NdleftO) NdleftO = NdleftN ENDIF useC1C2 = (IK+1<=mNc1c2) IF (useC1C2) THEN ! E(Y(J) | Y(1),...Y(I))/STD(Y(J)|Y(1),,,,Y(I)) mY(I+1:N+1) = mY(I+1:N+1) + mBIG(IK, I+1:N+1)*mY(IK) ENDIF ENDIF IK = IK + 1 FINA = 0 FINB = 0 END IF END DO IF (mNd>0) VAL = VAL * jacob(mXd,mXc) RETURN 200 VAL = gZERO RETURN END FUNCTION MVNFUN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!******************* RINDD - the main program *********************!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE RINDD(VALS,ERR,TERR,Big,Ex,Xc,Nt, & indI,Blo,Bup,INFIN) USE RCRUDEMOD USE KRBVRCMOD USE ADAPTMOD USE KROBOVMOD USE DKBVRCMOD USE SSOBOLMOD IMPLICIT NONE DOUBLE PRECISION, DIMENSION(: ), INTENT(out):: VALS, ERR ,TERR DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: BIG DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: Xc DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: Ex DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: Blo, Bup INTEGER, DIMENSION(:), INTENT(in) :: indI,INFIN INTEGER, INTENT(in) :: Nt ! DOUBLE PRECISION, INTENT(in) :: XcScale ! local variables INTEGER :: ix, INFORM, NDIM, MAXPTS, MINPTS DOUBLE PRECISION :: VALUE,fxc,absERR,absERR2 double precision :: LABSEPS,LRELEPS VALS(:) = gZERO ERR(:) = gONE TERR(:) = gONE call preInit(BIG,Xc,Nt,inform) IF (INFORM.GT.0) GOTO 110 ! degenerate case exit VALS=0 for all ! (should perhaps return NaN instead??) ! Now the loop over all different values of ! variables Xc (the one one is conditioning on) ! is started. The density f_{Xc}(xc(:,ix)) ! will be computed and denoted by fxc. DO ix = 1, mNx call initIntegrand(ix,Xc,Ex,indI,Blo,Bup,infin, & fxc,value,abserr,NDIM,inform) IF (INFORM.GT.0) GO TO 100 MAXPTS = mMAXPTS MINPTS = mMINPTS LABSEPS = max(mABSEPS-abserr,0.2D0*mABSEPS) !*fxc LRELEPS = mRELEPS ABSERR2 = mSmall SELECT CASE (mMethod) CASE (:1) IF (NDIM < 9) THEN CALL SADAPT(NDIM,MAXPTS,MVNFUN,LABSEPS, & LRELEPS,ABSERR2,VALUE,INFORM) VALUE = MAX(VALUE,gZERO) ELSE CALL KRBVRC(NDIM, MINPTS, MAXPTS, MVNFUN,LABSEPS,LRELEPS, & ABSERR2, VALUE, INFORM ) ENDIF CASE (2) ! Call the subregion adaptive integration subroutine IF ( NDIM .GT. 19.) THEN ! print *, 'Ndim too large for SADMVN => Calling KRBVRC' CALL KRBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, & LRELEPS, ABSERR2, VALUE, INFORM ) ELSE CALL SADAPT(NDIM,MAXPTS,MVNFUN,LABSEPS, & LRELEPS,ABSERR2,VALUE,INFORM) VALUE = MAX(VALUE,gZERO) ENDIF CASE (3) ! Call the Lattice rule integration procedure CALL KRBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, & LRELEPS, ABSERR2, VALUE, INFORM ) CASE (4) ! Call the Lattice rule ! integration procedure CALL KROBOV( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, & LRELEPS,ABSERR2, VALUE, INFORM ) CASE (5) ! Call Crude Monte Carlo integration procedure CALL RANMC( NDIM, MAXPTS, MVNFUN, LABSEPS, & LRELEPS, ABSERR2, VALUE, INFORM ) CASE (6) ! Call the scrambled Sobol sequence rule integration procedure CALL SOBNIED( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, LRELEPS, & ABSERR2, VALUE, INFORM ) CASE (7:) CALL DKBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, LRELEPS, & ABSERR2, VALUE, INFORM ) END SELECT ! IF (INFORM.gt.0) print *,'RIND, INFORM,error =',inform,error 100 VALS(ix) = VALUE*fxc IF (SIZE(ERR, DIM = 1).EQ.mNx) ERR(ix) = abserr2*fxc IF (SIZE(TERR, DIM = 1).EQ.mNx) TERR(ix) = abserr*fxc ENDDO !ix 110 CONTINUE call cleanUp RETURN END SUBROUTINE RINDD SUBROUTINE setIntLimits(xc,indI,Blo,Bup,INFIN,inform) IMPLICIT NONE DOUBLE PRECISION, DIMENSION(: ), INTENT(in) :: xc INTEGER, DIMENSION(: ), INTENT(in) :: indI,INFIN DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: Blo,Bup integer, intent(out) :: inform !Local variables INTEGER :: I, J, K, L,Mb1,Nb,NI,Nc DOUBLE PRECISION :: xCut, SQ0 !this procedure set mA,mB and mInfi according to Blo/Bup and INFIN ! ! INFIN INTEGER, array of integration limits flags: ! if INFIN(I) < 0, Ith limits are (-infinity, infinity); ! if INFIN(I) = 0, Ith limits are (-infinity, mB(I)]; ! if INFIN(I) = 1, Ith limits are [mA(I), infinity); ! if INFIN(I) = 2, Ith limits are [mA(I), mB(I)]. ! Note on member variables: ! mXedni = indices to the variables new place after cvsrtXc. Size Ntdc ! mCm = E(Xt,Xd|Xc), i.e., conditional mean given Xc ! mBIG(:,1:Ntd) = Cov(Xt,Xd|Xc) xCut = ABS(mInfinity) Mb1 = size(Blo,DIM=1)-1 Nb = size(Blo,DIM=2) NI = size(indI,DIM=1) Nc = size(xc,DIM=1) if (Mb1>Nc .or. Nb.NE.NI-1) then ! size of variables inconsistent inform = 4 return endif ! IF (Mb.GT.Nc+1) print *,'barrier: Mb,Nc =',Mb,Nc ! IF (Nb.NE.NI-1) print *,'barrier: Nb,NI =',Nb,NI DO J = 2, NI DO I = indI (J - 1) + 1 , indI (J) L = mXedni(I) mINFI(L) = INFIN(J-1) SQ0 = SQRT(mBIG(L,L)) mA(L) = -xCut*SQ0 mB(L) = xCut*SQ0 IF (mINFI(L).GE.0) THEN IF (mINFI(L).NE.0) THEN mA(L) = Blo (1, J - 1)-mCm(L) DO K = 1, Mb1 mA(L) = mA(L)+Blo(K+1,J-1)*xc(K) ENDDO ! K ! This can only be done if if (mA(L)< -xCut*SQ0) mINFI(L) = mINFI(L)-2 ENDIF IF (mINFI(L).NE.1) THEN mB(L) = Bup (1, J - 1)-mCm(L) DO K = 1, Mb1 mB(L) = mB(L)+Bup(K+1,J-1)*xc(K) ENDDO if (xCut*SQ0-1) THEN IF (infi.NE.0)THEN IF (A < -mXcutOff) THEN infi = infi-2 ! CALL mexprintf('ADJ A') ENDIF ENDIF IF (infi.NE.1) THEN IF (mXCutOff < B) THEN infi = infi-1 ! CALL mexprintf('ADJ B') ENDIF END IF END IF RETURN END SUBROUTINE ADJLIMITS SUBROUTINE C1C2(I0,I1,IK,A,B,INFIN, Cm, BIG, AJ, BJ, FINA,FINB) ! The regression equation for the conditional distr. of Y given X=x ! is equal to the conditional expectation of Y given X=x, i.e., ! ! E(Y|X=x) = E(Y) + Cov(Y,X)/Var(X)[x-E(X)] ! ! Let ! x1 = (x-E(X))/SQRT(Var(X)) be zero mean, ! C1< x1 B(I) or ! ! b) Cm(I)+x1*B1(I)+C*SQ(I)-1) THEN !BdSQ0 = B1(I) !CSQ = xCut * SQ(I) BdSQ0 = BIG(IK,I) CSQ = xCut * BIG(I,IK) IF (BdSQ0 > LTOL) THEN IF ( INFI .NE. 0 ) THEN IF (FINA.EQ.1) THEN AJ = MAX(AJ,(A(I) - Cm(I) - CSQ)/BdSQ0) ELSE AJ = (A(I) - Cm(I) - CSQ)/BdSQ0 FINA = 1 ENDIF IF (FINB.GT.0) AJ = MIN(AJ,BJ) END IF IF ( INFI .NE. 1 ) THEN IF (FINB.EQ.1) THEN BJ = MIN(BJ,(B(I) - Cm(I) + CSQ)/BdSQ0) ELSE BJ = (B(I) - Cm(I) + CSQ)/BdSQ0 FINB = 1 ENDIF IF (FINA.GT.0) BJ = MAX(AJ,BJ) END IF ELSEIF (BdSQ0 < -LTOL) THEN IF ( INFI .NE. 0 ) THEN IF (FINB.EQ.1) THEN BJ = MIN(BJ,(A(I) - Cm(I) - CSQ)/BdSQ0) ELSE BJ = (A(I) - Cm(I) - CSQ)/BdSQ0 FINB = 1 ENDIF IF (FINA.GT.0) BJ = MAX(AJ,BJ) END IF IF ( INFI .NE. 1 ) THEN IF (FINA.EQ.1) THEN AJ = MAX(AJ,(B(I) - Cm(I) + CSQ)/BdSQ0) ELSE AJ = (B(I) - Cm(I) + CSQ)/BdSQ0 FINA = 1 ENDIF IF (FINB.GT.0) AJ = MIN(AJ,BJ) END IF END IF ENDIF END DO ! IF (FINA>0 .AND. FINB>0) THEN ! IF (AJmaxSQ) maxSQ = SQ(I) ENDDO XCEPS2 = Ntdc*mSmall*maxSQ mXcEps2 = XCEPS2 LTOL = mSmall LO = 1 K = Ntdc DO I = 1, Nc ! Condsort Xc m = K+1-MAXLOC(SQ(K:Ntd+1:-1)) M1 = m(1) IF (SQ(m1)<=XCEPS2) THEN ! PRINT *,'CVSRTXC: Degenerate case of Xc(Nc-J+1) for J=',ix !CALL mexprintf('CVSRTXC: Degenerate case of Xc(Nc-J+1)') INFORM = 1 GOTO 200 ! RETURN !degenerate case ENDIF IF (M1.NE.K) THEN ! Symmetric row column permuations ! Swap row and columns, but only upper triangular part CALL RCSWAP( M1, K, Ntdc,Ntd, R,INDEX1,SQ) END IF R(K,K) = SQRT(SQ(K)) IF (K .EQ. LO) GOTO 200 R(LO:K-1,K) = R(LO:K-1,K)/R(K,K) ! Cov(Xi,Xj|Xk,Xk+1,..,Xn) = .... ! Cov(Xi,Xj|Xk+1,..,Xn) - Cov(Xi,Xk|Xk+1,..Xn)*Cov(Xj,Xk|Xk+1,..Xn) DO J = LO,K-1 ! Var(Xj | Xk,Xk+1,...,Xn) SQ(J) = R(J,J) - R(J,K)*R(J,K) IF (SQ(J)<=LTOL.AND.J<=Ntd) THEN IF (LO < J) THEN CALL RCSWAP(LO, J, Ntdc,Ntd, R,INDEX1,SQ) ENDIF R(LO,LO:K-1) = gZERO IF (SQ(LO) < -10.0D0*SQRT(LTOL)) THEN ! inform = 2 !R(LO,K) = gZERO ! CALL mexprintf('Negative definit BIG!'//CHAR(10)) ENDIF SQ(LO) = gZERO LO = LO + 1 ELSE R(J,J) = SQ(J) R(LO:J-1,J) = R(LO:J-1,J) - R(LO:J-1,K)*R(J,K) ENDIF END DO K = K - 1 ENDDO 200 DEALLOCATE(SQ) RETURN END SUBROUTINE CVSRTXC SUBROUTINE RCSCALE(chkLim,K,K0,N1,N,K1,INFIS,CDI,Cm, & R,A,B,INFI,INDEX1,Y) USE GLOBALCONST USE SWAPMOD IMPLICIT NONE !RCSCALE: Scale covariance matrix and limits ! ! CALL RCSCALE( k, k0, N1, N,K1, CDI,Cm,R,A, B, INFIN,index1,Y) ! ! chkLim = TRUE if check if variable K is redundant ! FALSE ! K = index to variable which is deterministic,i.e., ! STD(Xk|X1,...Xr) = 0 ! N1 = Number of significant variables of [Xt,Xd] ! N = length(Xt)+length(Xd) ! K1 = index to current variable we are conditioning on. ! CDI = Cholesky diagonal elements which contains either ! CDI(J) = STD(Xj | X1,...,Xj-1,Xc) if Xj is stochastic given ! X1,...Xj, Xc ! or ! CDI(J) = COV(Xj,Xk | X1,..,Xk-1,Xc )/STD(Xk | X1,..,Xk-1,Xc) ! if Xj is determinstically determined given X1,..,Xk,Xc ! for some k1) then ! Check if variable is redundant ! TODO: this chklim-block does not work correctly yet xCut = mInfinity I = 1 Ak = R(I,K)*xCut Bk = - (R(I,K))*xCut if (INFI(I)>=0) then if (INFI(I).ne.0) then Ak = -(R(I,K))*MAX(A(I),-xCut) endif if (INFI(I).ne.1) then Bk = - (R(I,K))*MIN(B(I),xCut) endif endif if (R(I,K) LTOL) .OR. (isXt)) THEN DO J = 1,I-1 isXd = (INDEX1(J)>Nt) IF ( (R(J,J) <= LTOL) .AND.isXd) THEN CALL RCSWAP(J, I, N, N, R,INDEX1,Cm, A, B, INFI) !GO TO 10 CYCLE LP3 ENDIF END DO ENDIF ! 10 END DO LP3 ! ! Move any doubly infinite limits or any redundant of Xt to the next ! innermost positions. ! LP4: DO I = N-INFISD, N1+1, -1 isXd = (INDEX1(I)>Nt) IF ( ((INFI(I) > -1).AND.(R(I,I) > LTOL)) & .OR. isXd) THEN DO J = 1,I-1 isXt = (INDEX1(J)<=Nt) IF ( (INFI(J) < 0 .OR. (R(J,J)<= LTOL)) & .AND. (isXt)) THEN CALL RCSWAP( J, I, N,N, R,INDEX1,Cm, A, B, INFI) !GO TO 15 CYCLE LP4 ENDIF END DO ENDIF !15 END DO LP4 ! CALL mexprintf('Before sorting') ! CALL PRINTCOF(N,A,B,INFI,R,INDEX1) ! CALL PRINTVEC(CDI,'CDI') ! CALL PRINTVEC(Cm,'Cm') IF ( N1 <= 0 ) GOTO 200 ! ! Sort remaining limits and determine Cholesky factor. ! Y(1:N1) = gZERO K = 1 Ndleft = Nd - INFISD Nullity = 0 DO WHILE (K .LE. N1) ! IF (Ndim.EQ.3) EPSL = MAX(EPS2,1D-10) ! Determine the integration limits for variable with minimum ! expected probability and interchange that variable with Kth. K0 = K - Nullity PRBMIN = gTWO JMIN = K CVDIAG = ZERO RMAX = ZERO IF ((Ndleft>0) .OR. (NDIM < Nd+mNIT)) THEN DO J = K, N1 isXd = (INDEX1(J)>Nt) isOK = ((NDIM <= Nd+mNIT).OR.isXd) IF ( R(J,J) <= K0*K0*EPSL .OR. (.NOT. isOK)) THEN RMAX = max(RMAX,ABS(R(J,J))) ELSE TMP = ZERO ! = conditional mean of Y(I) given Y(1:I-1) DO I = 1, K0 - 1 TMP = TMP + R(I,J)*Y(I) END DO SUMSQ = SQRT( R(J,J)) IF (INFI(J)>-1) THEN ! May have infinite int. limits if Nd>0 IF (INFI(J).NE.0) THEN AJ = ( A(J) - TMP )/SUMSQ ENDIF IF (INFI(J).NE.1) THEN BJ = ( B(J) - TMP )/SUMSQ ENDIF ENDIF IF (isXd) THEN AA = (Cm(J)+TMP)/SUMSQ ! inflection point CALL EXLMS(AA,AJ,BJ,INFI(J),D,E,Ca,Pa) PRBJ = E - D ELSE !CALL MVNLMS( AJ, BJ, INFI(J), D, E ) CALL MVNLIMITS(AJ,BJ,INFI(J),APJ,PRBJ) ENDIF !IF ( EMIN + D .GE. E + DMIN ) THEN IF ( PRBJ < PRBMIN ) THEN JMIN = J AMIN = AJ BMIN = BJ PRBMIN = MAX(PRBJ,ZERO) CVDIAG = SUMSQ ENDIF ENDIF END DO END IF ! ! Compute Ith column of Cholesky factor. ! Compute expected value for Ith integration variable (without ! considering the jacobian) and ! scale Ith covariance matrix row and limits. ! ! 40 IF ( CVDIAG.GT.TOL) THEN isXd = (INDEX1(JMIN)>Nt) IF (isXd) THEN Ndleft = Ndleft - 1 ELSEIF (BCVSRT.EQV..FALSE..AND.(PRBMIN+LTOL>=gONE)) THEN !BCVSRT.EQ. J = 1 AJ = R(J,JMIN)*xCut BJ = - (R(J,JMIN))*xCut if (INFI(J)>=0) then if (INFI(J).ne.0) then AJ = -(R(J,JMIN))*MAX(A(J),-xCut) endif if (INFI(J).ne.1) then BJ = - (R(J,JMIN))*MIN(B(J),xCut) endif endif if (R(J,JMIN)Nt) if (isXd) then Ndleft = Ndleft - 1 ELSEIF (BCVSRT.EQV..FALSE.) THEN ! BCVSRT.EQ. J = 1 AJ = R(J,I)*xCut BJ = - (R(J,I))*xCut if (INFI(J)>=0) then if (INFI(J).ne.0) then AJ = -(R(J,I))*MAX(A(J),-xCut) endif if (INFI(J).ne.1) then BJ = - (R(J,I))*MIN(B(J),xCut) endif endif if (R(J,I)100) THEN R(I,K0) = gZERO ELSE R(I,K0) = MAX(SQRT(MAX(R(I,I), gZERO)),LTOL) ENDIF Nullity = Nullity + 1 K = K + 1 IF (K < I) THEN CALL RCSWAP( K, I, N1,N,R,INDEX1,Cm, A, B, INFI) ! SWAP conditional standarddeviations DO J = 1, K0 CALL SWAP(R(K,J),R(I,J)) END DO ENDIF chkLim = .FALSE. !((.not.isXd).AND.(BCVSRT.EQ..FALSE.)) L = INFIS CALL RCSCALE(chkLim,K,K0,N1,N,K1,INFIS,CDI,Cm, & R,A,B,INFI,INDEX1) if (L.ne.INFIS) THEN K = K - 1 I = I - 1 ENDIF END IF I = I + 1 75 CONTINUE END DO INFJ = INFI(K1) IF (K1 .EQ.1) THEN FINA = 0 FINB = 0 IF (INFJ.GE.0) THEN IF (INFJ.NE.0) FINA = 1 IF (INFJ.NE.1) FINB = 1 ENDIF CALL C1C2(K1+1,N1,K0,A,B,INFI, Y, R, & AMIN, BMIN, FINA,FINB) INFJ = 2*FINA+FINB-1 CALL MVNLIMITS(AMIN,BMIN,INFJ,APJ,PRBMIN) ENDIF Y(K0) = gettmean(AMIN,BMIN,INFJ,PRBMIN) R( K0, K1 ) = R( K0, K1 ) / CVDIAG DO J = 1, K0 - 1 ! conditional covariances R( J, K1 ) = R( J, K1 ) / CVDIAG ! conditional standard dev.s used in regression eq. R( K1, J ) = R( K1, J ) / CVDIAG END DO A( K1 ) = A( K1 )/CVDIAG B( K1 ) = B( K1 )/CVDIAG K = K + 1 100 CONTINUE ELSE covErr = RMAX R(K:N1,K:N1) = gZERO I = K DO WHILE (I <= N1) ! Scale covariance matrix rows and limits ! If the conditional covariance matrix diagonal entry is zero, ! permute limits and/or rows, if necessary. chkLim = ((index1(I)<=Nt).AND.(BCVSRT.EQV..FALSE.)) L = INFIS CALL RCSCALE(chkLim,I,K0-1,N1,N,K1,INFIS,CDI,Cm, & R,A,B,INFI,INDEX1) if (L.EQ.INFIS) I = I + 1 END DO Nullity = N1 - K0 + 1 GOTO 200 !RETURN END IF END DO 200 CONTINUE IF (Ndim .GT. 0) THEN ! N10 IF (INFI(J).NE.0) THEN AJ = ( A(J) - TMP )/SUMSQ ENDIF IF (INFI(J).NE.1) THEN BJ = ( B(J) - TMP )/SUMSQ ENDIF 30 IF (INDEX1(J).GT.Nt) THEN AA = (Cm(J)+TMP)/SUMSQ ! inflection point CALL EXLMS(AA,AJ,BJ,INFI(J),D,E,Ca,Pa) PRBJ = E-D ELSE !CALL MVNLMS( AJ, BJ, INFI(J), D, E ) CALL MVNLIMITS(AJ,BJ,INFI(J),APJ,PRBJ) ENDIF IF ( PRBJ < PRBMIN ) THEN JMIN = J AMIN = AJ BMIN = BJ PRBMIN = MAX(PRBJ,ZERO) CVDIAG = SUMSQ ENDIF ENDIF END DO END IF ! ! Compute Ith column of Cholesky factor. ! Compute expected value for Ith integration variable (without ! considering the jacobian) and ! scale Ith covariance matrix row and limits. ! !40 IF ( CVDIAG.GT.TOL) THEN IF (INDEX1(JMIN).GT.Nt) THEN Ndleft = Ndleft-1 ELSE IF (BCVSRT.EQV..FALSE..AND.(PRBMIN+LTOL.GE.gONE)) THEN !BCVSRT.EQ. I = 1 AJ = R(I,JMIN)*xCut BJ = - (R(I,JMIN))*xCut if (INFI(1)>=0) then if (INFI(1).ne.0) then AJ = -(R(I,JMIN))*MAX(A(I),-xCut) endif if (INFI(1).ne.1) then BJ = - (R(I,JMIN))*MIN(B(I),xCut) endif endif if (R(I,JMIN)