Datasets
Models
Diff of
/partyMod/src/mvt.f
[000000]
..
[fbf06f]
Switch to unified view
| a | b/partyMod/src/mvt.f | ||
|---|---|---|---|
| 1 | * |
||
| 2 | * $Id$ |
||
| 3 | * |
||
| 4 | SUBROUTINE MVTDST( N, NU, LOWER, UPPER, INFIN, CORREL, DELTA, |
||
| 5 | & MAXPTS, ABSEPS, RELEPS, ERROR, VALUE, INFORM ) |
||
| 6 | * |
||
| 7 | * A subroutine for computing non-central multivariate t probabilities. |
||
| 8 | * This subroutine uses an algorithm (QRSVN) described in the paper |
||
| 9 | * "Comparison of Methods for the Computation of Multivariate |
||
| 10 | * t-Probabilities", by Alan Genz and Frank Bretz |
||
| 11 | * J. Comp. Graph. Stat. 11 (2002), pp. 950-971. |
||
| 12 | * |
||
| 13 | * Alan Genz |
||
| 14 | * Department of Mathematics |
||
| 15 | * Washington State University |
||
| 16 | * Pullman, WA 99164-3113 |
||
| 17 | * Email : AlanGenz@wsu.edu |
||
| 18 | * |
||
| 19 | * Original source available from |
||
| 20 | * http://www.math.wsu.edu/faculty/genz/software/fort77/mvtdstpack.f |
||
| 21 | * |
||
| 22 | * This is version 28/05/2013 |
||
| 23 | * |
||
| 24 | * Parameters |
||
| 25 | * |
||
| 26 | * N INTEGER, the number of variables. |
||
| 27 | * NU INTEGER, the number of degrees of freedom. |
||
| 28 | * If NU < 1, then an MVN probability is computed. |
||
| 29 | * LOWER DOUBLE PRECISION, array of lower integration limits. |
||
| 30 | * UPPER DOUBLE PRECISION, array of upper integration limits. |
||
| 31 | * INFIN INTEGER, array of integration limits flags: |
||
| 32 | * if INFIN(I) < 0, Ith limits are (-infinity, infinity); |
||
| 33 | * if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)]; |
||
| 34 | * if INFIN(I) = 1, Ith limits are [LOWER(I), infinity); |
||
| 35 | * if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. |
||
| 36 | * CORREL DOUBLE PRECISION, array of correlation coefficients; |
||
| 37 | * the correlation coefficient in row I column J of the |
||
| 38 | * correlation matrixshould be stored in |
||
| 39 | * CORREL( J + ((I-2)*(I-1))/2 ), for J < I. |
||
| 40 | * The correlation matrix must be positive semi-definite. |
||
| 41 | * DELTA DOUBLE PRECISION, array of non-centrality parameters. |
||
| 42 | * MAXPTS INTEGER, maximum number of function values allowed. This |
||
| 43 | * parameter can be used to limit the time. A sensible |
||
| 44 | * strategy is to start with MAXPTS = 1000*N, and then |
||
| 45 | * increase MAXPTS if ERROR is too large. |
||
| 46 | * ABSEPS DOUBLE PRECISION absolute error tolerance. |
||
| 47 | * RELEPS DOUBLE PRECISION relative error tolerance. |
||
| 48 | * ERROR DOUBLE PRECISION estimated absolute error, |
||
| 49 | * with 99% confidence level. |
||
| 50 | * VALUE DOUBLE PRECISION estimated value for the integral |
||
| 51 | * INFORM INTEGER, termination status parameter: |
||
| 52 | * if INFORM = 0, normal completion with ERROR < EPS; |
||
| 53 | * if INFORM = 1, completion with ERROR > EPS and MAXPTS |
||
| 54 | * function vaules used; increase MAXPTS to |
||
| 55 | * decrease ERROR; |
||
| 56 | * if INFORM = 2, N > 1000 or N < 1. |
||
| 57 | * if INFORM = 3, correlation matrix not positive semi-definite. |
||
| 58 | * |
||
| 59 | EXTERNAL MVSUBR |
||
| 60 | INTEGER N, ND, NU, INFIN(*), MAXPTS, INFORM, IVLS |
||
| 61 | DOUBLE PRECISION CORREL(*), LOWER(*), UPPER(*), DELTA(*), RELEPS, |
||
| 62 | & ABSEPS, ERROR, VALUE, E(1), V(1) |
||
| 63 | COMMON /PTBLCK/IVLS |
||
| 64 | IVLS = 0 |
||
| 65 | |||
| 66 | IF ( N .GT. 1000 .OR. N .LT. 1 ) THEN |
||
| 67 | VALUE = 0 |
||
| 68 | ERROR = 1 |
||
| 69 | INFORM = 2 |
||
| 70 | ELSE |
||
| 71 | CALL MVINTS( N, NU, CORREL, LOWER, UPPER, DELTA, INFIN, |
||
| 72 | & ND, VALUE, ERROR, INFORM ) |
||
| 73 | IF ( INFORM .EQ. 0 .AND. ND .GT. 0 ) THEN |
||
| 74 | * |
||
| 75 | * Call the lattice rule integration subroutine |
||
| 76 | * |
||
| 77 | CALL MVKBRV( ND, IVLS, MAXPTS, 1, MVSUBR, ABSEPS, RELEPS, |
||
| 78 | & E(1), V, INFORM ) |
||
| 79 | ERROR = E(1) |
||
| 80 | VALUE = V(1) |
||
| 81 | ENDIF |
||
| 82 | ENDIF |
||
| 83 | |||
| 84 | END |
||
| 85 | * |
||
| 86 | SUBROUTINE MVSUBR( N, W, NF, F ) |
||
| 87 | * |
||
| 88 | * Integrand subroutine |
||
| 89 | * |
||
| 90 | INTEGER N, NF, NUIN, INFIN(*), NL |
||
| 91 | DOUBLE PRECISION W(*),F(*), LOWER(*),UPPER(*), CORREL(*), DELTA(*) |
||
| 92 | PARAMETER ( NL = 1000 ) |
||
| 93 | INTEGER INFI(NL), NU, ND, INFORM, NY |
||
| 94 | DOUBLE PRECISION COV(NL*(NL+1)/2), A(NL), B(NL), DL(NL), Y(NL) |
||
| 95 | DOUBLE PRECISION MVCHNV, SNU, R, VL, ER, DI, EI |
||
| 96 | SAVE NU, SNU, INFI, A, B, DL, COV |
||
| 97 | IF ( NU .LE. 0 ) THEN |
||
| 98 | R = 1 |
||
| 99 | CALL MVVLSB( N+1, W, R, DL,INFI,A,B,COV, Y, DI,EI, NY, F(1) ) |
||
| 100 | ELSE |
||
| 101 | R = MVCHNV( NU, W(N) )/SNU |
||
| 102 | CALL MVVLSB( N , W, R, DL,INFI,A,B,COV, Y, DI,EI, NY, F(1) ) |
||
| 103 | END IF |
||
| 104 | RETURN |
||
| 105 | * |
||
| 106 | * Entry point for intialization. |
||
| 107 | * |
||
| 108 | ENTRY MVINTS( N, NUIN, CORREL, LOWER, UPPER, DELTA, INFIN, |
||
| 109 | & ND, VL, ER, INFORM ) |
||
| 110 | * |
||
| 111 | * Initialization and computation of covariance Cholesky factor. |
||
| 112 | * |
||
| 113 | CALL MVSORT( N, LOWER, UPPER, DELTA, CORREL, INFIN, Y, .TRUE., |
||
| 114 | & ND, A, B, DL, COV, INFI, INFORM ) |
||
| 115 | NU = NUIN |
||
| 116 | CALL MVSPCL( ND, NU, A, B, DL, COV, INFI, SNU, VL, ER, INFORM ) |
||
| 117 | END |
||
| 118 | * |
||
| 119 | SUBROUTINE MVSPCL( ND, NU, A,B,DL, COV, INFI, SNU, VL,ER, INFORM ) |
||
| 120 | * |
||
| 121 | * Special cases subroutine |
||
| 122 | * |
||
| 123 | DOUBLE PRECISION COV(*), A(*), B(*), DL(*), SNU, R, VL, ER |
||
| 124 | INTEGER ND, NU, INFI(*), INFORM |
||
| 125 | DOUBLE PRECISION MVBVT, MVSTDT |
||
| 126 | IF ( INFORM .GT. 0 ) THEN |
||
| 127 | VL = 0 |
||
| 128 | ER = 1 |
||
| 129 | ELSE |
||
| 130 | * |
||
| 131 | * Special cases |
||
| 132 | * |
||
| 133 | IF ( ND .EQ. 0 ) THEN |
||
| 134 | ER = 0 |
||
| 135 | * Code added to fix ND = 0 bug, 24/03/2009 -> |
||
| 136 | VL = 1 |
||
| 137 | * <- Code added to fix ND = 0 bug, 24/03/2009 |
||
| 138 | ELSE IF ( ND.EQ.1 .AND. ( NU.LT.1 .OR. ABS(DL(1)).EQ.0 ) ) THEN |
||
| 139 | * |
||
| 140 | * 1-d case for normal or central t |
||
| 141 | * |
||
| 142 | VL = 1 |
||
| 143 | IF ( INFI(1) .NE. 1 ) VL = MVSTDT( NU, B(1) - DL(1) ) |
||
| 144 | IF ( INFI(1) .NE. 0 ) VL = VL - MVSTDT( NU, A(1) - DL(1) ) |
||
| 145 | IF ( VL .LT. 0 ) VL = 0 |
||
| 146 | ER = 2D-16 |
||
| 147 | ND = 0 |
||
| 148 | ELSE IF ( ND .EQ. 2 .AND. |
||
| 149 | & ( NU .LT. 1 .OR. ABS(DL(1))+ABS(DL(2)) .EQ. 0 ) ) THEN |
||
| 150 | * |
||
| 151 | * 2-d case for normal or central t |
||
| 152 | * |
||
| 153 | IF ( INFI(1) .NE. 0 ) A(1) = A(1) - DL(1) |
||
| 154 | IF ( INFI(1) .NE. 1 ) B(1) = B(1) - DL(1) |
||
| 155 | IF ( INFI(2) .NE. 0 ) A(2) = A(2) - DL(2) |
||
| 156 | IF ( INFI(2) .NE. 1 ) B(2) = B(2) - DL(2) |
||
| 157 | IF ( ABS( COV(3) ) .GT. 0 ) THEN |
||
| 158 | * |
||
| 159 | * 2-d nonsingular case |
||
| 160 | * |
||
| 161 | R = SQRT( 1 + COV(2)**2 ) |
||
| 162 | IF ( INFI(2) .NE. 0 ) A(2) = A(2)/R |
||
| 163 | IF ( INFI(2) .NE. 1 ) B(2) = B(2)/R |
||
| 164 | COV(2) = COV(2)/R |
||
| 165 | VL = MVBVT( NU, A, B, INFI, COV(2) ) |
||
| 166 | ER = 1D-15 |
||
| 167 | ELSE |
||
| 168 | * |
||
| 169 | * 2-d singular case |
||
| 170 | * |
||
| 171 | IF ( INFI(1) .NE. 0 ) THEN |
||
| 172 | IF ( INFI(2) .NE. 0 ) A(1) = MAX( A(1), A(2) ) |
||
| 173 | ELSE |
||
| 174 | IF ( INFI(2) .NE. 0 ) A(1) = A(2) |
||
| 175 | END IF |
||
| 176 | IF ( INFI(1) .NE. 1 ) THEN |
||
| 177 | IF ( INFI(2) .NE. 1 ) B(1) = MIN( B(1), B(2) ) |
||
| 178 | ELSE |
||
| 179 | IF ( INFI(2) .NE. 1 ) B(1) = B(2) |
||
| 180 | END IF |
||
| 181 | IF ( INFI(1) .NE. INFI(2) ) INFI(1) = 2 |
||
| 182 | VL = 1 |
||
| 183 | * A(1), B(1) Bug Fixed, 28/05/2013 |
||
| 184 | IF ( INFI(1) .NE. 1 ) VL = MVSTDT( NU, B(1) ) |
||
| 185 | IF ( INFI(1) .NE. 0 ) VL = VL - MVSTDT( NU, A(1) ) |
||
| 186 | IF ( VL .LT. 0 ) VL = 0 |
||
| 187 | ER = 2D-16 |
||
| 188 | END IF |
||
| 189 | ND = 0 |
||
| 190 | ELSE |
||
| 191 | IF ( NU .GT. 0 ) THEN |
||
| 192 | SNU = SQRT( DBLE(NU) ) |
||
| 193 | ELSE |
||
| 194 | ND = ND - 1 |
||
| 195 | END IF |
||
| 196 | END IF |
||
| 197 | END IF |
||
| 198 | END |
||
| 199 | * |
||
| 200 | SUBROUTINE MVVLSB( N,W,R,DL,INFI, A,B,COV, Y, DI,EI, ND, VALUE ) |
||
| 201 | * |
||
| 202 | * Integrand subroutine |
||
| 203 | * |
||
| 204 | INTEGER N, INFI(*), ND |
||
| 205 | DOUBLE PRECISION W(*), R, DL(*), A(*), B(*), COV(*), Y(*) |
||
| 206 | INTEGER I, J, IJ, INFA, INFB |
||
| 207 | DOUBLE PRECISION SUM, AI, BI, DI, EI, MVPHNV, VALUE |
||
| 208 | VALUE = 1 |
||
| 209 | INFA = 0 |
||
| 210 | INFB = 0 |
||
| 211 | ND = 0 |
||
| 212 | IJ = 0 |
||
| 213 | DO I = 1, N |
||
| 214 | SUM = DL(I) |
||
| 215 | DO J = 1, I-1 |
||
| 216 | IJ = IJ + 1 |
||
| 217 | IF ( J .LE. ND ) SUM = SUM + COV(IJ)*Y(J) |
||
| 218 | END DO |
||
| 219 | IF ( INFI(I) .NE. 0 ) THEN |
||
| 220 | IF ( INFA .EQ. 1 ) THEN |
||
| 221 | AI = MAX( AI, R*A(I) - SUM ) |
||
| 222 | ELSE |
||
| 223 | AI = R*A(I) - SUM |
||
| 224 | INFA = 1 |
||
| 225 | END IF |
||
| 226 | END IF |
||
| 227 | IF ( INFI(I) .NE. 1 ) THEN |
||
| 228 | IF ( INFB .EQ. 1 ) THEN |
||
| 229 | BI = MIN( BI, R*B(I) - SUM ) |
||
| 230 | ELSE |
||
| 231 | BI = R*B(I) - SUM |
||
| 232 | INFB = 1 |
||
| 233 | END IF |
||
| 234 | END IF |
||
| 235 | IJ = IJ + 1 |
||
| 236 | IF ( I .EQ. N .OR. COV(IJ+ND+2) .GT. 0 ) THEN |
||
| 237 | CALL MVLIMS( AI, BI, INFA + INFA + INFB - 1, DI, EI ) |
||
| 238 | IF ( DI .GE. EI ) THEN |
||
| 239 | VALUE = 0 |
||
| 240 | RETURN |
||
| 241 | ELSE |
||
| 242 | VALUE = VALUE*( EI - DI ) |
||
| 243 | ND = ND + 1 |
||
| 244 | IF ( I .LT. N ) Y(ND) = MVPHNV( DI + W(ND)*( EI - DI ) ) |
||
| 245 | INFA = 0 |
||
| 246 | INFB = 0 |
||
| 247 | END IF |
||
| 248 | END IF |
||
| 249 | END DO |
||
| 250 | END |
||
| 251 | * |
||
| 252 | SUBROUTINE MVSORT( N, LOWER, UPPER, DELTA, CORREL, INFIN, Y,PIVOT, |
||
| 253 | & ND, A, B, DL, COV, INFI, INFORM ) |
||
| 254 | * |
||
| 255 | * Subroutine to sort integration limits and determine Cholesky factor. |
||
| 256 | * |
||
| 257 | INTEGER N, ND, INFIN(*), INFI(*), INFORM |
||
| 258 | LOGICAL PIVOT |
||
| 259 | DOUBLE PRECISION A(*), B(*), DL(*), COV(*), |
||
| 260 | & LOWER(*), UPPER(*), DELTA(*), CORREL(*), Y(*) |
||
| 261 | INTEGER I, J, K, L, M, II, IJ, IL, JL, JMIN |
||
| 262 | DOUBLE PRECISION SUMSQ, AJ, BJ, SUM, EPS, EPSI, D, E |
||
| 263 | DOUBLE PRECISION CVDIAG, AMIN, BMIN, DEMIN, MVTDNS |
||
| 264 | PARAMETER ( EPS = 1D-10 ) |
||
| 265 | INFORM = 0 |
||
| 266 | IJ = 0 |
||
| 267 | II = 0 |
||
| 268 | ND = N |
||
| 269 | DO I = 1, N |
||
| 270 | A(I) = 0 |
||
| 271 | B(I) = 0 |
||
| 272 | DL(I) = 0 |
||
| 273 | INFI(I) = INFIN(I) |
||
| 274 | IF ( INFI(I) .LT. 0 ) THEN |
||
| 275 | ND = ND - 1 |
||
| 276 | ELSE |
||
| 277 | IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I) |
||
| 278 | IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I) |
||
| 279 | DL(I) = DELTA(I) |
||
| 280 | ENDIF |
||
| 281 | DO J = 1, I-1 |
||
| 282 | IJ = IJ + 1 |
||
| 283 | II = II + 1 |
||
| 284 | COV(IJ) = CORREL(II) |
||
| 285 | END DO |
||
| 286 | IJ = IJ + 1 |
||
| 287 | COV(IJ) = 1 |
||
| 288 | END DO |
||
| 289 | * |
||
| 290 | * First move any doubly infinite limits to innermost positions. |
||
| 291 | * |
||
| 292 | IF ( ND .GT. 0 ) THEN |
||
| 293 | DO I = N, ND + 1, -1 |
||
| 294 | IF ( INFI(I) .GE. 0 ) THEN |
||
| 295 | DO J = 1, I-1 |
||
| 296 | IF ( INFI(J) .LT. 0 ) THEN |
||
| 297 | CALL MVSWAP( J, I, A, B, DL, INFI, N, COV ) |
||
| 298 | GO TO 10 |
||
| 299 | ENDIF |
||
| 300 | END DO |
||
| 301 | ENDIF |
||
| 302 | 10 CONTINUE |
||
| 303 | END DO |
||
| 304 | * |
||
| 305 | * Sort remaining limits and determine Cholesky factor. |
||
| 306 | * |
||
| 307 | II = 0 |
||
| 308 | JL = ND |
||
| 309 | DO I = 1, ND |
||
| 310 | * |
||
| 311 | * Determine the integration limits for variable with minimum |
||
| 312 | * expected probability and interchange that variable with Ith. |
||
| 313 | * |
||
| 314 | DEMIN = 1 |
||
| 315 | JMIN = I |
||
| 316 | CVDIAG = 0 |
||
| 317 | IJ = II |
||
| 318 | EPSI = EPS*I |
||
| 319 | IF ( .NOT. PIVOT ) JL = I |
||
| 320 | DO J = I, JL |
||
| 321 | IF ( COV(IJ+J) .GT. EPSI ) THEN |
||
| 322 | SUMSQ = SQRT( COV(IJ+J) ) |
||
| 323 | SUM = DL(J) |
||
| 324 | DO K = 1, I-1 |
||
| 325 | SUM = SUM + COV(IJ+K)*Y(K) |
||
| 326 | END DO |
||
| 327 | AJ = ( A(J) - SUM )/SUMSQ |
||
| 328 | BJ = ( B(J) - SUM )/SUMSQ |
||
| 329 | CALL MVLIMS( AJ, BJ, INFI(J), D, E ) |
||
| 330 | IF ( DEMIN .GE. E - D ) THEN |
||
| 331 | JMIN = J |
||
| 332 | AMIN = AJ |
||
| 333 | BMIN = BJ |
||
| 334 | DEMIN = E - D |
||
| 335 | CVDIAG = SUMSQ |
||
| 336 | ENDIF |
||
| 337 | ENDIF |
||
| 338 | IJ = IJ + J |
||
| 339 | END DO |
||
| 340 | IF ( JMIN .GT. I ) THEN |
||
| 341 | CALL MVSWAP( I, JMIN, A, B, DL, INFI, N, COV ) |
||
| 342 | END IF |
||
| 343 | IF ( COV(II+I) .LT. -EPSI ) THEN |
||
| 344 | INFORM = 3 |
||
| 345 | END IF |
||
| 346 | COV(II+I) = CVDIAG |
||
| 347 | * |
||
| 348 | * Compute Ith column of Cholesky factor. |
||
| 349 | * Compute expected value for Ith integration variable and |
||
| 350 | * scale Ith covariance matrix row and limits. |
||
| 351 | * |
||
| 352 | IF ( CVDIAG .GT. 0 ) THEN |
||
| 353 | IL = II + I |
||
| 354 | DO L = I+1, ND |
||
| 355 | COV(IL+I) = COV(IL+I)/CVDIAG |
||
| 356 | IJ = II + I |
||
| 357 | DO J = I+1, L |
||
| 358 | COV(IL+J) = COV(IL+J) - COV(IL+I)*COV(IJ+I) |
||
| 359 | IJ = IJ + J |
||
| 360 | END DO |
||
| 361 | IL = IL + L |
||
| 362 | END DO |
||
| 363 | * |
||
| 364 | * Expected Y = -( density(b) - density(a) )/( b - a ) |
||
| 365 | * |
||
| 366 | IF ( DEMIN .GT. EPSI ) THEN |
||
| 367 | Y(I) = 0 |
||
| 368 | IF ( INFI(I) .NE. 0 ) Y(I) = MVTDNS( 0, AMIN ) |
||
| 369 | IF ( INFI(I) .NE. 1 ) Y(I) = Y(I) - MVTDNS( 0, BMIN ) |
||
| 370 | Y(I) = Y(I)/DEMIN |
||
| 371 | ELSE |
||
| 372 | IF ( INFI(I) .EQ. 0 ) Y(I) = BMIN |
||
| 373 | IF ( INFI(I) .EQ. 1 ) Y(I) = AMIN |
||
| 374 | IF ( INFI(I) .EQ. 2 ) Y(I) = ( AMIN + BMIN )/2 |
||
| 375 | END IF |
||
| 376 | DO J = 1, I |
||
| 377 | II = II + 1 |
||
| 378 | COV(II) = COV(II)/CVDIAG |
||
| 379 | END DO |
||
| 380 | A(I) = A(I)/CVDIAG |
||
| 381 | B(I) = B(I)/CVDIAG |
||
| 382 | DL(I) = DL(I)/CVDIAG |
||
| 383 | ELSE |
||
| 384 | IL = II + I |
||
| 385 | DO L = I+1, ND |
||
| 386 | COV(IL+I) = 0 |
||
| 387 | IL = IL + L |
||
| 388 | END DO |
||
| 389 | * |
||
| 390 | * If the covariance matrix diagonal entry is zero, |
||
| 391 | * permute limits and rows, if necessary. |
||
| 392 | * |
||
| 393 | * |
||
| 394 | DO J = I-1, 1, -1 |
||
| 395 | IF ( ABS( COV(II+J) ) .GT. EPSI ) THEN |
||
| 396 | A(I) = A(I)/COV(II+J) |
||
| 397 | B(I) = B(I)/COV(II+J) |
||
| 398 | DL(I) = DL(I)/COV(II+J) |
||
| 399 | IF ( COV(II+J) .LT. 0 ) THEN |
||
| 400 | CALL MVSSWP( A(I), B(I) ) |
||
| 401 | IF ( INFI(I) .NE. 2 ) INFI(I) = 1 - INFI(I) |
||
| 402 | END IF |
||
| 403 | DO L = 1, J |
||
| 404 | COV(II+L) = COV(II+L)/COV(II+J) |
||
| 405 | END DO |
||
| 406 | DO L = J+1, I-1 |
||
| 407 | IF( COV((L-1)*L/2+J+1) .GT. 0 ) THEN |
||
| 408 | IJ = II |
||
| 409 | DO K = I-1, L, -1 |
||
| 410 | DO M = 1, K |
||
| 411 | CALL MVSSWP( COV(IJ-K+M), COV(IJ+M) ) |
||
| 412 | END DO |
||
| 413 | CALL MVSSWP( A(K), A(K+1) ) |
||
| 414 | CALL MVSSWP( B(K), B(K+1) ) |
||
| 415 | CALL MVSSWP( DL(K), DL(K+1) ) |
||
| 416 | M = INFI(K) |
||
| 417 | INFI(K) = INFI(K+1) |
||
| 418 | INFI(K+1) = M |
||
| 419 | IJ = IJ - K |
||
| 420 | END DO |
||
| 421 | GO TO 20 |
||
| 422 | END IF |
||
| 423 | END DO |
||
| 424 | GO TO 20 |
||
| 425 | END IF |
||
| 426 | COV(II+J) = 0 |
||
| 427 | END DO |
||
| 428 | 20 II = II + I |
||
| 429 | Y(I) = 0 |
||
| 430 | END IF |
||
| 431 | END DO |
||
| 432 | ENDIF |
||
| 433 | END |
||
| 434 | * |
||
| 435 | DOUBLE PRECISION FUNCTION MVTDNS( NU, X ) |
||
| 436 | INTEGER NU, I |
||
| 437 | DOUBLE PRECISION X, PROD, PI, SQTWPI |
||
| 438 | PARAMETER ( PI = 3.141592653589793D0 ) |
||
| 439 | PARAMETER ( SQTWPI = 2.506628274631001D0 ) |
||
| 440 | MVTDNS = 0 |
||
| 441 | IF ( NU .GT. 0 ) THEN |
||
| 442 | PROD = 1/SQRT( DBLE(NU) ) |
||
| 443 | DO I = NU - 2, 1, -2 |
||
| 444 | PROD = PROD*( I + 1 )/I |
||
| 445 | END DO |
||
| 446 | IF ( MOD( NU, 2 ) .EQ. 0 ) THEN |
||
| 447 | PROD = PROD/2 |
||
| 448 | ELSE |
||
| 449 | PROD = PROD/PI |
||
| 450 | END IF |
||
| 451 | MVTDNS = PROD/SQRT( 1 + X*X/NU )**( NU + 1 ) |
||
| 452 | ELSE |
||
| 453 | IF ( ABS(X) .LT. 10 ) MVTDNS = EXP( -X*X/2 )/SQTWPI |
||
| 454 | END IF |
||
| 455 | END |
||
| 456 | * |
||
| 457 | SUBROUTINE MVLIMS( A, B, INFIN, LOWER, UPPER ) |
||
| 458 | DOUBLE PRECISION A, B, LOWER, UPPER, MVPHI |
||
| 459 | INTEGER INFIN |
||
| 460 | LOWER = 0 |
||
| 461 | UPPER = 1 |
||
| 462 | IF ( INFIN .GE. 0 ) THEN |
||
| 463 | IF ( INFIN .NE. 0 ) LOWER = MVPHI(A) |
||
| 464 | IF ( INFIN .NE. 1 ) UPPER = MVPHI(B) |
||
| 465 | ENDIF |
||
| 466 | UPPER = MAX( UPPER, LOWER ) |
||
| 467 | END |
||
| 468 | * |
||
| 469 | SUBROUTINE MVSSWP( X, Y ) |
||
| 470 | DOUBLE PRECISION X, Y, T |
||
| 471 | T = X |
||
| 472 | X = Y |
||
| 473 | Y = T |
||
| 474 | END |
||
| 475 | * |
||
| 476 | SUBROUTINE MVSWAP( P, Q, A, B, D, INFIN, N, C ) |
||
| 477 | * |
||
| 478 | * Swaps rows and columns P and Q in situ, with P <= Q. |
||
| 479 | * |
||
| 480 | DOUBLE PRECISION A(*), B(*), C(*), D(*) |
||
| 481 | INTEGER INFIN(*), P, Q, N, I, J, II, JJ |
||
| 482 | CALL MVSSWP( A(P), A(Q) ) |
||
| 483 | CALL MVSSWP( B(P), B(Q) ) |
||
| 484 | CALL MVSSWP( D(P), D(Q) ) |
||
| 485 | J = INFIN(P) |
||
| 486 | INFIN(P) = INFIN(Q) |
||
| 487 | INFIN(Q) = J |
||
| 488 | JJ = ( P*( P - 1 ) )/2 |
||
| 489 | II = ( Q*( Q - 1 ) )/2 |
||
| 490 | CALL MVSSWP( C(JJ+P), C(II+Q) ) |
||
| 491 | DO J = 1, P-1 |
||
| 492 | CALL MVSSWP( C(JJ+J), C(II+J) ) |
||
| 493 | END DO |
||
| 494 | JJ = JJ + P |
||
| 495 | DO I = P+1, Q-1 |
||
| 496 | CALL MVSSWP( C(JJ+P), C(II+I) ) |
||
| 497 | JJ = JJ + I |
||
| 498 | END DO |
||
| 499 | II = II + Q |
||
| 500 | DO I = Q+1, N |
||
| 501 | CALL MVSSWP( C(II+P), C(II+Q) ) |
||
| 502 | II = II + I |
||
| 503 | END DO |
||
| 504 | END |
||
| 505 | * |
||
| 506 | DOUBLE PRECISION FUNCTION MVPHI(Z) |
||
| 507 | * |
||
| 508 | * Normal distribution probabilities accurate to 1d-15. |
||
| 509 | * Reference: J.L. Schonfelder, Math Comp 32(1978), pp 1232-1240. |
||
| 510 | * |
||
| 511 | INTEGER I, IM |
||
| 512 | DOUBLE PRECISION A(0:43), BM, B, BP, P, RTWO, T, XA, Z |
||
| 513 | PARAMETER( RTWO = 1.414213562373095048801688724209D0, IM = 24 ) |
||
| 514 | SAVE A |
||
| 515 | DATA ( A(I), I = 0, 43 )/ |
||
| 516 | & 6.10143081923200417926465815756D-1, |
||
| 517 | & -4.34841272712577471828182820888D-1, |
||
| 518 | & 1.76351193643605501125840298123D-1, |
||
| 519 | & -6.0710795609249414860051215825D-2, |
||
| 520 | & 1.7712068995694114486147141191D-2, |
||
| 521 | & -4.321119385567293818599864968D-3, |
||
| 522 | & 8.54216676887098678819832055D-4, |
||
| 523 | & -1.27155090609162742628893940D-4, |
||
| 524 | & 1.1248167243671189468847072D-5, 3.13063885421820972630152D-7, |
||
| 525 | & -2.70988068537762022009086D-7, 3.0737622701407688440959D-8, |
||
| 526 | & 2.515620384817622937314D-9, -1.028929921320319127590D-9, |
||
| 527 | & 2.9944052119949939363D-11, 2.6051789687266936290D-11, |
||
| 528 | & -2.634839924171969386D-12, -6.43404509890636443D-13, |
||
| 529 | & 1.12457401801663447D-13, 1.7281533389986098D-14, |
||
| 530 | & -4.264101694942375D-15, -5.45371977880191D-16, |
||
| 531 | & 1.58697607761671D-16, 2.0899837844334D-17, |
||
| 532 | & -5.900526869409D-18, -9.41893387554D-19, 2.14977356470D-19, |
||
| 533 | & 4.6660985008D-20, -7.243011862D-21, -2.387966824D-21, |
||
| 534 | & 1.91177535D-22, 1.20482568D-22, -6.72377D-25, -5.747997D-24, |
||
| 535 | & -4.28493D-25, 2.44856D-25, 4.3793D-26, -8.151D-27, -3.089D-27, |
||
| 536 | & 9.3D-29, 1.74D-28, 1.6D-29, -8.0D-30, -2.0D-30 / |
||
| 537 | * |
||
| 538 | XA = ABS(Z)/RTWO |
||
| 539 | IF ( XA .GT. 100 ) THEN |
||
| 540 | P = 0 |
||
| 541 | ELSE |
||
| 542 | T = ( 8*XA - 30 ) / ( 4*XA + 15 ) |
||
| 543 | BM = 0 |
||
| 544 | B = 0 |
||
| 545 | DO I = IM, 0, -1 |
||
| 546 | BP = B |
||
| 547 | B = BM |
||
| 548 | BM = T*B - BP + A(I) |
||
| 549 | END DO |
||
| 550 | P = EXP( -XA*XA )*( BM - BP )/4 |
||
| 551 | END IF |
||
| 552 | IF ( Z .GT. 0 ) P = 1 - P |
||
| 553 | MVPHI = P |
||
| 554 | END |
||
| 555 | * |
||
| 556 | DOUBLE PRECISION FUNCTION MVPHNV(P) |
||
| 557 | * |
||
| 558 | * ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 |
||
| 559 | * |
||
| 560 | * Produces the normal deviate Z corresponding to a given lower |
||
| 561 | * tail area of P. |
||
| 562 | * |
||
| 563 | * The hash sums below are the sums of the mantissas of the |
||
| 564 | * coefficients. They are included for use in checking |
||
| 565 | * transcription. |
||
| 566 | * |
||
| 567 | DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2, |
||
| 568 | * A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7, |
||
| 569 | * C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7, |
||
| 570 | * E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7, |
||
| 571 | * P, Q, R |
||
| 572 | PARAMETER ( SPLIT1 = 0.425, SPLIT2 = 5, |
||
| 573 | * CONST1 = 0.180625D0, CONST2 = 1.6D0 ) |
||
| 574 | * |
||
| 575 | * Coefficients for P close to 0.5 |
||
| 576 | * |
||
| 577 | PARAMETER ( |
||
| 578 | * A0 = 3.38713 28727 96366 6080D0, |
||
| 579 | * A1 = 1.33141 66789 17843 7745D+2, |
||
| 580 | * A2 = 1.97159 09503 06551 4427D+3, |
||
| 581 | * A3 = 1.37316 93765 50946 1125D+4, |
||
| 582 | * A4 = 4.59219 53931 54987 1457D+4, |
||
| 583 | * A5 = 6.72657 70927 00870 0853D+4, |
||
| 584 | * A6 = 3.34305 75583 58812 8105D+4, |
||
| 585 | * A7 = 2.50908 09287 30122 6727D+3, |
||
| 586 | * B1 = 4.23133 30701 60091 1252D+1, |
||
| 587 | * B2 = 6.87187 00749 20579 0830D+2, |
||
| 588 | * B3 = 5.39419 60214 24751 1077D+3, |
||
| 589 | * B4 = 2.12137 94301 58659 5867D+4, |
||
| 590 | * B5 = 3.93078 95800 09271 0610D+4, |
||
| 591 | * B6 = 2.87290 85735 72194 2674D+4, |
||
| 592 | * B7 = 5.22649 52788 52854 5610D+3 ) |
||
| 593 | * HASH SUM AB 55.88319 28806 14901 4439 |
||
| 594 | * |
||
| 595 | * Coefficients for P not close to 0, 0.5 or 1. |
||
| 596 | * |
||
| 597 | PARAMETER ( |
||
| 598 | * C0 = 1.42343 71107 49683 57734D0, |
||
| 599 | * C1 = 4.63033 78461 56545 29590D0, |
||
| 600 | * C2 = 5.76949 72214 60691 40550D0, |
||
| 601 | * C3 = 3.64784 83247 63204 60504D0, |
||
| 602 | * C4 = 1.27045 82524 52368 38258D0, |
||
| 603 | * C5 = 2.41780 72517 74506 11770D-1, |
||
| 604 | * C6 = 2.27238 44989 26918 45833D-2, |
||
| 605 | * C7 = 7.74545 01427 83414 07640D-4, |
||
| 606 | * D1 = 2.05319 16266 37758 82187D0, |
||
| 607 | * D2 = 1.67638 48301 83803 84940D0, |
||
| 608 | * D3 = 6.89767 33498 51000 04550D-1, |
||
| 609 | * D4 = 1.48103 97642 74800 74590D-1, |
||
| 610 | * D5 = 1.51986 66563 61645 71966D-2, |
||
| 611 | * D6 = 5.47593 80849 95344 94600D-4, |
||
| 612 | * D7 = 1.05075 00716 44416 84324D-9 ) |
||
| 613 | * HASH SUM CD 49.33206 50330 16102 89036 |
||
| 614 | * |
||
| 615 | * Coefficients for P near 0 or 1. |
||
| 616 | * |
||
| 617 | PARAMETER ( |
||
| 618 | * E0 = 6.65790 46435 01103 77720D0, |
||
| 619 | * E1 = 5.46378 49111 64114 36990D0, |
||
| 620 | * E2 = 1.78482 65399 17291 33580D0, |
||
| 621 | * E3 = 2.96560 57182 85048 91230D-1, |
||
| 622 | * E4 = 2.65321 89526 57612 30930D-2, |
||
| 623 | * E5 = 1.24266 09473 88078 43860D-3, |
||
| 624 | * E6 = 2.71155 55687 43487 57815D-5, |
||
| 625 | * E7 = 2.01033 43992 92288 13265D-7, |
||
| 626 | * F1 = 5.99832 20655 58879 37690D-1, |
||
| 627 | * F2 = 1.36929 88092 27358 05310D-1, |
||
| 628 | * F3 = 1.48753 61290 85061 48525D-2, |
||
| 629 | * F4 = 7.86869 13114 56132 59100D-4, |
||
| 630 | * F5 = 1.84631 83175 10054 68180D-5, |
||
| 631 | * F6 = 1.42151 17583 16445 88870D-7, |
||
| 632 | * F7 = 2.04426 31033 89939 78564D-15 ) |
||
| 633 | * HASH SUM EF 47.52583 31754 92896 71629 |
||
| 634 | * |
||
| 635 | Q = ( 2*P - 1 )/2 |
||
| 636 | IF ( ABS(Q) .LE. SPLIT1 ) THEN |
||
| 637 | R = CONST1 - Q*Q |
||
| 638 | MVPHNV = Q*( ( ( ((((A7*R + A6)*R + A5)*R + A4)*R + A3) |
||
| 639 | * *R + A2 )*R + A1 )*R + A0 ) |
||
| 640 | * /( ( ( ((((B7*R + B6)*R + B5)*R + B4)*R + B3) |
||
| 641 | * *R + B2 )*R + B1 )*R + 1 ) |
||
| 642 | ELSE |
||
| 643 | R = MIN( P, 1 - P ) |
||
| 644 | IF ( R .GT. 0 ) THEN |
||
| 645 | R = SQRT( -LOG(R) ) |
||
| 646 | IF ( R .LE. SPLIT2 ) THEN |
||
| 647 | R = R - CONST2 |
||
| 648 | MVPHNV = ( ( ( ((((C7*R + C6)*R + C5)*R + C4)*R + C3) |
||
| 649 | * *R + C2 )*R + C1 )*R + C0 ) |
||
| 650 | * /( ( ( ((((D7*R + D6)*R + D5)*R + D4)*R + D3) |
||
| 651 | * *R + D2 )*R + D1 )*R + 1 ) |
||
| 652 | ELSE |
||
| 653 | R = R - SPLIT2 |
||
| 654 | MVPHNV = ( ( ( ((((E7*R + E6)*R + E5)*R + E4)*R + E3) |
||
| 655 | * *R + E2 )*R + E1 )*R + E0 ) |
||
| 656 | * /( ( ( ((((F7*R + F6)*R + F5)*R + F4)*R + F3) |
||
| 657 | * *R + F2 )*R + F1 )*R + 1 ) |
||
| 658 | END IF |
||
| 659 | ELSE |
||
| 660 | MVPHNV = 9 |
||
| 661 | END IF |
||
| 662 | IF ( Q .LT. 0 ) MVPHNV = - MVPHNV |
||
| 663 | END IF |
||
| 664 | END |
||
| 665 | DOUBLE PRECISION FUNCTION MVBVN( LOWER, UPPER, INFIN, CORREL ) |
||
| 666 | * |
||
| 667 | * A function for computing bivariate normal probabilities. |
||
| 668 | * |
||
| 669 | * Parameters |
||
| 670 | * |
||
| 671 | * LOWER REAL, array of lower integration limits. |
||
| 672 | * UPPER REAL, array of upper integration limits. |
||
| 673 | * INFIN INTEGER, array of integration limits flags: |
||
| 674 | * if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)]; |
||
| 675 | * if INFIN(I) = 1, Ith limits are [LOWER(I), infinity); |
||
| 676 | * if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. |
||
| 677 | * CORREL REAL, correlation coefficient. |
||
| 678 | * |
||
| 679 | DOUBLE PRECISION LOWER(*), UPPER(*), CORREL, MVBVU |
||
| 680 | INTEGER INFIN(*) |
||
| 681 | IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 682 | MVBVN = MVBVU ( LOWER(1), LOWER(2), CORREL ) |
||
| 683 | + - MVBVU ( UPPER(1), LOWER(2), CORREL ) |
||
| 684 | + - MVBVU ( LOWER(1), UPPER(2), CORREL ) |
||
| 685 | + + MVBVU ( UPPER(1), UPPER(2), CORREL ) |
||
| 686 | ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 687 | MVBVN = MVBVU ( LOWER(1), LOWER(2), CORREL ) |
||
| 688 | + - MVBVU ( UPPER(1), LOWER(2), CORREL ) |
||
| 689 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 690 | MVBVN = MVBVU ( LOWER(1), LOWER(2), CORREL ) |
||
| 691 | + - MVBVU ( LOWER(1), UPPER(2), CORREL ) |
||
| 692 | ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 693 | MVBVN = MVBVU ( -UPPER(1), -UPPER(2), CORREL ) |
||
| 694 | + - MVBVU ( -LOWER(1), -UPPER(2), CORREL ) |
||
| 695 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 696 | MVBVN = MVBVU ( -UPPER(1), -UPPER(2), CORREL ) |
||
| 697 | + - MVBVU ( -UPPER(1), -LOWER(2), CORREL ) |
||
| 698 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 699 | MVBVN = MVBVU ( LOWER(1), -UPPER(2), -CORREL ) |
||
| 700 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 701 | MVBVN = MVBVU ( -UPPER(1), LOWER(2), -CORREL ) |
||
| 702 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 703 | MVBVN = MVBVU ( LOWER(1), LOWER(2), CORREL ) |
||
| 704 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 705 | MVBVN = MVBVU ( -UPPER(1), -UPPER(2), CORREL ) |
||
| 706 | ELSE |
||
| 707 | MVBVN = 1 |
||
| 708 | END IF |
||
| 709 | END |
||
| 710 | DOUBLE PRECISION FUNCTION MVBVU( SH, SK, R ) |
||
| 711 | * |
||
| 712 | * A function for computing bivariate normal probabilities; |
||
| 713 | * developed using |
||
| 714 | * Drezner, Z. and Wesolowsky, G. O. (1989), |
||
| 715 | * On the Computation of the Bivariate Normal Integral, |
||
| 716 | * J. Stat. Comput. Simul.. 35 pp. 101-107. |
||
| 717 | * with extensive modications for double precisions by |
||
| 718 | * Alan Genz and Yihong Ge |
||
| 719 | * Department of Mathematics |
||
| 720 | * Washington State University |
||
| 721 | * Pullman, WA 99164-3113 |
||
| 722 | * Email : alangenz@wsu.edu |
||
| 723 | * |
||
| 724 | * BVN - calculate the probability that X is larger than SH and Y is |
||
| 725 | * larger than SK. |
||
| 726 | * |
||
| 727 | * Parameters |
||
| 728 | * |
||
| 729 | * SH REAL, integration limit |
||
| 730 | * SK REAL, integration limit |
||
| 731 | * R REAL, correlation coefficient |
||
| 732 | * LG INTEGER, number of Gauss Rule Points and Weights |
||
| 733 | * |
||
| 734 | DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI |
||
| 735 | INTEGER I, LG, NG |
||
| 736 | PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 ) |
||
| 737 | DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS |
||
| 738 | DOUBLE PRECISION MVPHI, SN, ASR, H, K, BS, HS, HK |
||
| 739 | SAVE X, W |
||
| 740 | * Gauss Legendre Points and Weights, N = 6 |
||
| 741 | DATA ( W(I,1), X(I,1), I = 1, 3 ) / |
||
| 742 | * 0.1713244923791705D+00,-0.9324695142031522D+00, |
||
| 743 | * 0.3607615730481384D+00,-0.6612093864662647D+00, |
||
| 744 | * 0.4679139345726904D+00,-0.2386191860831970D+00/ |
||
| 745 | * Gauss Legendre Points and Weights, N = 12 |
||
| 746 | DATA ( W(I,2), X(I,2), I = 1, 6 ) / |
||
| 747 | * 0.4717533638651177D-01,-0.9815606342467191D+00, |
||
| 748 | * 0.1069393259953183D+00,-0.9041172563704750D+00, |
||
| 749 | * 0.1600783285433464D+00,-0.7699026741943050D+00, |
||
| 750 | * 0.2031674267230659D+00,-0.5873179542866171D+00, |
||
| 751 | * 0.2334925365383547D+00,-0.3678314989981802D+00, |
||
| 752 | * 0.2491470458134029D+00,-0.1252334085114692D+00/ |
||
| 753 | * Gauss Legendre Points and Weights, N = 20 |
||
| 754 | DATA ( W(I,3), X(I,3), I = 1, 10 ) / |
||
| 755 | * 0.1761400713915212D-01,-0.9931285991850949D+00, |
||
| 756 | * 0.4060142980038694D-01,-0.9639719272779138D+00, |
||
| 757 | * 0.6267204833410906D-01,-0.9122344282513259D+00, |
||
| 758 | * 0.8327674157670475D-01,-0.8391169718222188D+00, |
||
| 759 | * 0.1019301198172404D+00,-0.7463319064601508D+00, |
||
| 760 | * 0.1181945319615184D+00,-0.6360536807265150D+00, |
||
| 761 | * 0.1316886384491766D+00,-0.5108670019508271D+00, |
||
| 762 | * 0.1420961093183821D+00,-0.3737060887154196D+00, |
||
| 763 | * 0.1491729864726037D+00,-0.2277858511416451D+00, |
||
| 764 | * 0.1527533871307259D+00,-0.7652652113349733D-01/ |
||
| 765 | IF ( ABS(R) .LT. 0.3 ) THEN |
||
| 766 | NG = 1 |
||
| 767 | LG = 3 |
||
| 768 | ELSE IF ( ABS(R) .LT. 0.75 ) THEN |
||
| 769 | NG = 2 |
||
| 770 | LG = 6 |
||
| 771 | ELSE |
||
| 772 | NG = 3 |
||
| 773 | LG = 10 |
||
| 774 | ENDIF |
||
| 775 | H = SH |
||
| 776 | K = SK |
||
| 777 | HK = H*K |
||
| 778 | BVN = 0 |
||
| 779 | IF ( ABS(R) .LT. 0.925 ) THEN |
||
| 780 | HS = ( H*H + K*K )/2 |
||
| 781 | ASR = ASIN(R) |
||
| 782 | DO I = 1, LG |
||
| 783 | SN = SIN(ASR*( X(I,NG)+1 )/2) |
||
| 784 | BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) ) |
||
| 785 | SN = SIN(ASR*(-X(I,NG)+1 )/2) |
||
| 786 | BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) ) |
||
| 787 | END DO |
||
| 788 | BVN = BVN*ASR/(2*TWOPI) + MVPHI(-H)*MVPHI(-K) |
||
| 789 | ELSE |
||
| 790 | IF ( R .LT. 0 ) THEN |
||
| 791 | K = -K |
||
| 792 | HK = -HK |
||
| 793 | ENDIF |
||
| 794 | IF ( ABS(R) .LT. 1 ) THEN |
||
| 795 | AS = ( 1 - R )*( 1 + R ) |
||
| 796 | A = SQRT(AS) |
||
| 797 | BS = ( H - K )**2 |
||
| 798 | C = ( 4 - HK )/8 |
||
| 799 | D = ( 12 - HK )/16 |
||
| 800 | BVN = A*EXP( -(BS/AS + HK)/2 ) |
||
| 801 | + *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 ) |
||
| 802 | IF ( HK .GT. -160 ) THEN |
||
| 803 | B = SQRT(BS) |
||
| 804 | BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVPHI(-B/A)*B |
||
| 805 | + *( 1 - C*BS*( 1 - D*BS/5 )/3 ) |
||
| 806 | ENDIF |
||
| 807 | A = A/2 |
||
| 808 | DO I = 1, LG |
||
| 809 | XS = ( A*(X(I,NG)+1) )**2 |
||
| 810 | RS = SQRT( 1 - XS ) |
||
| 811 | BVN = BVN + A*W(I,NG)* |
||
| 812 | + ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS |
||
| 813 | + - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) ) |
||
| 814 | XS = AS*(-X(I,NG)+1)**2/4 |
||
| 815 | RS = SQRT( 1 - XS ) |
||
| 816 | BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 ) |
||
| 817 | + *( EXP( -HK*XS/(2*(1+RS)**2) )/RS |
||
| 818 | + - ( 1 + C*XS*( 1 + D*XS ) ) ) |
||
| 819 | END DO |
||
| 820 | BVN = -BVN/TWOPI |
||
| 821 | ENDIF |
||
| 822 | IF ( R .GT. 0 ) THEN |
||
| 823 | BVN = BVN + MVPHI( -MAX( H, K ) ) |
||
| 824 | ELSE |
||
| 825 | BVN = -BVN |
||
| 826 | IF ( K .GT. H ) THEN |
||
| 827 | IF ( H .LT. 0 ) THEN |
||
| 828 | BVN = BVN + MVPHI(K) - MVPHI(H) |
||
| 829 | ELSE |
||
| 830 | BVN = BVN + MVPHI(-H) - MVPHI(-K) |
||
| 831 | ENDIF |
||
| 832 | ENDIF |
||
| 833 | ENDIF |
||
| 834 | ENDIF |
||
| 835 | MVBVU = BVN |
||
| 836 | END |
||
| 837 | * |
||
| 838 | DOUBLE PRECISION FUNCTION MVSTDT( NU, T ) |
||
| 839 | * |
||
| 840 | * Student t Distribution Function |
||
| 841 | * |
||
| 842 | * T |
||
| 843 | * TSTDNT = C I ( 1 + y*y/NU )**( -(NU+1)/2 ) dy |
||
| 844 | * NU -INF |
||
| 845 | * |
||
| 846 | INTEGER NU, J |
||
| 847 | DOUBLE PRECISION MVPHI, T, CSTHE, SNTHE, POLYN, TT, TS, RN, PI |
||
| 848 | PARAMETER ( PI = 3.141592653589793D0 ) |
||
| 849 | IF ( NU .LT. 1 ) THEN |
||
| 850 | MVSTDT = MVPHI( T ) |
||
| 851 | ELSE IF ( NU .EQ. 1 ) THEN |
||
| 852 | MVSTDT = ( 1 + 2*ATAN( T )/PI )/2 |
||
| 853 | ELSE IF ( NU .EQ. 2) THEN |
||
| 854 | MVSTDT = ( 1 + T/SQRT( 2 + T*T ))/2 |
||
| 855 | ELSE |
||
| 856 | TT = T*T |
||
| 857 | CSTHE = NU/( NU + TT ) |
||
| 858 | POLYN = 1 |
||
| 859 | DO J = NU - 2, 2, -2 |
||
| 860 | POLYN = 1 + ( J - 1 )*CSTHE*POLYN/J |
||
| 861 | END DO |
||
| 862 | IF ( MOD( NU, 2 ) .EQ. 1 ) THEN |
||
| 863 | RN = NU |
||
| 864 | TS = T/SQRT(RN) |
||
| 865 | MVSTDT = ( 1 + 2*( ATAN( TS ) + TS*CSTHE*POLYN )/PI )/2 |
||
| 866 | ELSE |
||
| 867 | SNTHE = T/SQRT( NU + TT ) |
||
| 868 | MVSTDT = ( 1 + SNTHE*POLYN )/2 |
||
| 869 | END IF |
||
| 870 | IF ( MVSTDT .LT. 0 ) MVSTDT = 0 |
||
| 871 | ENDIF |
||
| 872 | END |
||
| 873 | * |
||
| 874 | DOUBLE PRECISION FUNCTION MVBVT( NU, LOWER, UPPER, INFIN, CORREL ) |
||
| 875 | * |
||
| 876 | * A function for computing bivariate normal and t probabilities. |
||
| 877 | * |
||
| 878 | * Parameters |
||
| 879 | * |
||
| 880 | * NU INTEGER degrees of freedom parameter; NU < 1 gives normal case. |
||
| 881 | * LOWER REAL, array of lower integration limits. |
||
| 882 | * UPPER REAL, array of upper integration limits. |
||
| 883 | * INFIN INTEGER, array of integration limits flags: |
||
| 884 | * if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)]; |
||
| 885 | * if INFIN(I) = 1, Ith limits are [LOWER(I), infinity); |
||
| 886 | * if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. |
||
| 887 | * CORREL REAL, correlation coefficient. |
||
| 888 | * |
||
| 889 | DOUBLE PRECISION LOWER(*), UPPER(*), CORREL, MVBVN, MVBVTL |
||
| 890 | INTEGER NU, INFIN(*) |
||
| 891 | IF ( NU .LT. 1 ) THEN |
||
| 892 | MVBVT = MVBVN ( LOWER, UPPER, INFIN, CORREL ) |
||
| 893 | ELSE |
||
| 894 | IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 895 | MVBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) |
||
| 896 | + - MVBVTL ( NU, UPPER(1), LOWER(2), CORREL ) |
||
| 897 | + - MVBVTL ( NU, LOWER(1), UPPER(2), CORREL ) |
||
| 898 | + + MVBVTL ( NU, LOWER(1), LOWER(2), CORREL ) |
||
| 899 | ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 900 | MVBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL ) |
||
| 901 | + - MVBVTL ( NU, -UPPER(1), -LOWER(2), CORREL ) |
||
| 902 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 903 | MVBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL ) |
||
| 904 | + - MVBVTL ( NU, -LOWER(1), -UPPER(2), CORREL ) |
||
| 905 | ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 906 | MVBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) |
||
| 907 | + - MVBVTL ( NU, LOWER(1), UPPER(2), CORREL ) |
||
| 908 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 909 | MVBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) |
||
| 910 | + - MVBVTL ( NU, UPPER(1), LOWER(2), CORREL ) |
||
| 911 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 912 | MVBVT = MVBVTL ( NU, -LOWER(1), UPPER(2), -CORREL ) |
||
| 913 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 914 | MVBVT = MVBVTL ( NU, UPPER(1), -LOWER(2), -CORREL ) |
||
| 915 | ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 1 ) THEN |
||
| 916 | MVBVT = MVBVTL ( NU, -LOWER(1), -LOWER(2), CORREL ) |
||
| 917 | ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 0 ) THEN |
||
| 918 | MVBVT = MVBVTL ( NU, UPPER(1), UPPER(2), CORREL ) |
||
| 919 | ELSE |
||
| 920 | MVBVT = 1 |
||
| 921 | END IF |
||
| 922 | END IF |
||
| 923 | END |
||
| 924 | * |
||
| 925 | DOUBLE PRECISION FUNCTION MVBVTC( NU, L, U, INFIN, RHO ) |
||
| 926 | * |
||
| 927 | * A function for computing complementary bivariate normal and t |
||
| 928 | * probabilities. |
||
| 929 | * |
||
| 930 | * Parameters |
||
| 931 | * |
||
| 932 | * NU INTEGER degrees of freedom parameter. |
||
| 933 | * L REAL, array of lower integration limits. |
||
| 934 | * U REAL, array of upper integration limits. |
||
| 935 | * INFIN INTEGER, array of integration limits flags: |
||
| 936 | * if INFIN(1) INFIN(2), then MVBVTC computes |
||
| 937 | * 0 0 P( X>U(1), Y>U(2) ) |
||
| 938 | * 1 0 P( X<L(1), Y>U(2) ) |
||
| 939 | * 0 1 P( X>U(1), Y<L(2) ) |
||
| 940 | * 1 1 P( X<L(1), Y<L(2) ) |
||
| 941 | * 2 0 P( X>U(1), Y>U(2) ) + P( X<L(1), Y>U(2) ) |
||
| 942 | * 2 1 P( X>U(1), Y<L(2) ) + P( X<L(1), Y<L(2) ) |
||
| 943 | * 0 2 P( X>U(1), Y>U(2) ) + P( X>U(1), Y<L(2) ) |
||
| 944 | * 1 2 P( X<L(1), Y>U(2) ) + P( X<L(1), Y<L(2) ) |
||
| 945 | * 2 2 P( X>U(1), Y<L(2) ) + P( X<L(1), Y<L(2) ) |
||
| 946 | * + P( X>U(1), Y>U(2) ) + P( X<L(1), Y>U(2) ) |
||
| 947 | * |
||
| 948 | * RHO REAL, correlation coefficient. |
||
| 949 | * |
||
| 950 | DOUBLE PRECISION L(*), U(*), LW(2), UP(2), B, RHO, MVBVT |
||
| 951 | INTEGER I, NU, INFIN(*), INF(2) |
||
| 952 | * |
||
| 953 | DO I = 1, 2 |
||
| 954 | IF ( MOD( INFIN(I), 2 ) .EQ. 0 ) THEN |
||
| 955 | INF(I) = 1 |
||
| 956 | LW(I) = U(I) |
||
| 957 | ELSE |
||
| 958 | INF(I) = 0 |
||
| 959 | UP(I) = L(I) |
||
| 960 | END IF |
||
| 961 | END DO |
||
| 962 | B = MVBVT( NU, LW, UP, INF, RHO ) |
||
| 963 | DO I = 1, 2 |
||
| 964 | IF ( INFIN(I) .EQ. 2 ) THEN |
||
| 965 | INF(I) = 0 |
||
| 966 | UP(I) = L(I) |
||
| 967 | B = B + MVBVT( NU, LW, UP, INF, RHO ) |
||
| 968 | END IF |
||
| 969 | END DO |
||
| 970 | IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 2 ) THEN |
||
| 971 | INF(1) = 1 |
||
| 972 | LW(1) = U(1) |
||
| 973 | B = B + MVBVT( NU, LW, UP, INF, RHO ) |
||
| 974 | END IF |
||
| 975 | MVBVTC = B |
||
| 976 | END |
||
| 977 | * |
||
| 978 | double precision function mvbvtl( nu, dh, dk, r ) |
||
| 979 | * |
||
| 980 | * a function for computing bivariate t probabilities. |
||
| 981 | * |
||
| 982 | * Alan Genz |
||
| 983 | * Department of Mathematics |
||
| 984 | * Washington State University |
||
| 985 | * Pullman, Wa 99164-3113 |
||
| 986 | * Email : alangenz@wsu.edu |
||
| 987 | * |
||
| 988 | * this function is based on the method described by |
||
| 989 | * Dunnett, C.W. and M. Sobel, (1954), |
||
| 990 | * A bivariate generalization of Student's t-distribution |
||
| 991 | * with tables for certain special cases, |
||
| 992 | * Biometrika 41, pp. 153-169. |
||
| 993 | * |
||
| 994 | * mvbvtl - calculate the probability that x < dh and y < dk. |
||
| 995 | * |
||
| 996 | * parameters |
||
| 997 | * |
||
| 998 | * nu number of degrees of freedom |
||
| 999 | * dh 1st lower integration limit |
||
| 1000 | * dk 2nd lower integration limit |
||
| 1001 | * r correlation coefficient |
||
| 1002 | * |
||
| 1003 | integer nu, j, hs, ks |
||
| 1004 | double precision dh, dk, r |
||
| 1005 | double precision tpi, pi, ors, hrk, krh, bvt, snu |
||
| 1006 | double precision gmph, gmpk, xnkh, xnhk, qhrk, hkn, hpk, hkrn |
||
| 1007 | double precision btnckh, btnchk, btpdkh, btpdhk, one |
||
| 1008 | parameter ( pi = 3.14159265358979323844d0, tpi = 2*pi, one = 1 ) |
||
| 1009 | snu = sqrt( dble(nu) ) |
||
| 1010 | ors = 1 - r*r |
||
| 1011 | hrk = dh - r*dk |
||
| 1012 | krh = dk - r*dh |
||
| 1013 | if ( abs(hrk) + ors .gt. 0 ) then |
||
| 1014 | xnhk = hrk**2/( hrk**2 + ors*( nu + dk**2 ) ) |
||
| 1015 | xnkh = krh**2/( krh**2 + ors*( nu + dh**2 ) ) |
||
| 1016 | else |
||
| 1017 | xnhk = 0 |
||
| 1018 | xnkh = 0 |
||
| 1019 | end if |
||
| 1020 | hs = sign( one, dh - r*dk ) |
||
| 1021 | ks = sign( one, dk - r*dh ) |
||
| 1022 | if ( mod( nu, 2 ) .eq. 0 ) then |
||
| 1023 | bvt = atan2( sqrt(ors), -r )/tpi |
||
| 1024 | gmph = dh/sqrt( 16*( nu + dh**2 ) ) |
||
| 1025 | gmpk = dk/sqrt( 16*( nu + dk**2 ) ) |
||
| 1026 | btnckh = 2*atan2( sqrt( xnkh ), sqrt( 1 - xnkh ) )/pi |
||
| 1027 | btpdkh = 2*sqrt( xnkh*( 1 - xnkh ) )/pi |
||
| 1028 | btnchk = 2*atan2( sqrt( xnhk ), sqrt( 1 - xnhk ) )/pi |
||
| 1029 | btpdhk = 2*sqrt( xnhk*( 1 - xnhk ) )/pi |
||
| 1030 | do j = 1, nu/2 |
||
| 1031 | bvt = bvt + gmph*( 1 + ks*btnckh ) |
||
| 1032 | bvt = bvt + gmpk*( 1 + hs*btnchk ) |
||
| 1033 | btnckh = btnckh + btpdkh |
||
| 1034 | btpdkh = 2*j*btpdkh*( 1 - xnkh )/( 2*j + 1 ) |
||
| 1035 | btnchk = btnchk + btpdhk |
||
| 1036 | btpdhk = 2*j*btpdhk*( 1 - xnhk )/( 2*j + 1 ) |
||
| 1037 | gmph = gmph*( 2*j - 1 )/( 2*j*( 1 + dh**2/nu ) ) |
||
| 1038 | gmpk = gmpk*( 2*j - 1 )/( 2*j*( 1 + dk**2/nu ) ) |
||
| 1039 | end do |
||
| 1040 | else |
||
| 1041 | qhrk = sqrt( dh**2 + dk**2 - 2*r*dh*dk + nu*ors ) |
||
| 1042 | hkrn = dh*dk + r*nu |
||
| 1043 | hkn = dh*dk - nu |
||
| 1044 | hpk = dh + dk |
||
| 1045 | bvt = atan2(-snu*(hkn*qhrk+hpk*hkrn),hkn*hkrn-nu*hpk*qhrk)/tpi |
||
| 1046 | if ( bvt .lt. -1d-15 ) bvt = bvt + 1 |
||
| 1047 | gmph = dh/( tpi*snu*( 1 + dh**2/nu ) ) |
||
| 1048 | gmpk = dk/( tpi*snu*( 1 + dk**2/nu ) ) |
||
| 1049 | btnckh = sqrt( xnkh ) |
||
| 1050 | btpdkh = btnckh |
||
| 1051 | btnchk = sqrt( xnhk ) |
||
| 1052 | btpdhk = btnchk |
||
| 1053 | do j = 1, ( nu - 1 )/2 |
||
| 1054 | bvt = bvt + gmph*( 1 + ks*btnckh ) |
||
| 1055 | bvt = bvt + gmpk*( 1 + hs*btnchk ) |
||
| 1056 | btpdkh = ( 2*j - 1 )*btpdkh*( 1 - xnkh )/( 2*j ) |
||
| 1057 | btnckh = btnckh + btpdkh |
||
| 1058 | btpdhk = ( 2*j - 1 )*btpdhk*( 1 - xnhk )/( 2*j ) |
||
| 1059 | btnchk = btnchk + btpdhk |
||
| 1060 | gmph = 2*j*gmph/( ( 2*j + 1 )*( 1 + dh**2/nu ) ) |
||
| 1061 | gmpk = 2*j*gmpk/( ( 2*j + 1 )*( 1 + dk**2/nu ) ) |
||
| 1062 | end do |
||
| 1063 | end if |
||
| 1064 | mvbvtl = bvt |
||
| 1065 | * |
||
| 1066 | * end mvbvtl |
||
| 1067 | * |
||
| 1068 | end |
||
| 1069 | * |
||
| 1070 | DOUBLE PRECISION FUNCTION MVCHNV( N, P ) |
||
| 1071 | * |
||
| 1072 | * MVCHNV |
||
| 1073 | * P = 1 - K I exp(-t*t/2) t**(N-1) dt, for N >= 1. |
||
| 1074 | * N 0 |
||
| 1075 | * |
||
| 1076 | INTEGER I, N, NO |
||
| 1077 | DOUBLE PRECISION P, TWO, R, RO, LRP, LKN, MVPHNV, MVCHNC |
||
| 1078 | PARAMETER ( LRP = -.22579135264472743235D0, TWO = 2 ) |
||
| 1079 | * LRP = LOG( SQRT( 2/PI ) ) |
||
| 1080 | SAVE NO, LKN |
||
| 1081 | DATA NO / 0 / |
||
| 1082 | IF ( N .LE. 1 ) THEN |
||
| 1083 | R = -MVPHNV( P/2 ) |
||
| 1084 | ELSE IF ( P .LT. 1 ) THEN |
||
| 1085 | IF ( N .EQ. 2 ) THEN |
||
| 1086 | R = SQRT( -2*LOG(P) ) |
||
| 1087 | ELSE |
||
| 1088 | IF ( N .NE. NO ) THEN |
||
| 1089 | NO = N |
||
| 1090 | LKN = 0 |
||
| 1091 | DO I = N-2, 2, -2 |
||
| 1092 | LKN = LKN - LOG( DBLE(I) ) |
||
| 1093 | END DO |
||
| 1094 | IF ( MOD( N, 2 ) .EQ. 1 ) LKN = LKN + LRP |
||
| 1095 | END IF |
||
| 1096 | IF ( N .GE. -5*LOG(1-P)/4 ) THEN |
||
| 1097 | R = TWO/( 9*N ) |
||
| 1098 | R = N*( -MVPHNV(P)*SQRT(R) + 1 - R )**3 |
||
| 1099 | IF ( R .GT. 2*N+6 ) THEN |
||
| 1100 | R = 2*( LKN - LOG(P) ) + ( N - 2 )*LOG(R) |
||
| 1101 | END IF |
||
| 1102 | ELSE |
||
| 1103 | R = EXP( ( LOG( (1-P)*N ) - LKN )*TWO/N ) |
||
| 1104 | END IF |
||
| 1105 | R = SQRT(R) |
||
| 1106 | RO = R |
||
| 1107 | R = MVCHNC( LKN, N, P, R ) |
||
| 1108 | IF ( ABS( R - RO ) .GT. 1D-6 ) THEN |
||
| 1109 | RO = R |
||
| 1110 | R = MVCHNC( LKN, N, P, R ) |
||
| 1111 | IF ( ABS( R - RO ) .GT. 1D-6 ) R = MVCHNC( LKN, N, P, R ) |
||
| 1112 | END IF |
||
| 1113 | END IF |
||
| 1114 | ELSE |
||
| 1115 | R = 0 |
||
| 1116 | END IF |
||
| 1117 | MVCHNV = R |
||
| 1118 | END |
||
| 1119 | * |
||
| 1120 | DOUBLE PRECISION FUNCTION MVCHNC( LKN, N, P, R ) |
||
| 1121 | * |
||
| 1122 | * Third order Schroeder correction to R for MVCHNV |
||
| 1123 | * |
||
| 1124 | INTEGER I, N |
||
| 1125 | DOUBLE PRECISION P, R, LKN, DF, RR, RN, CHI, MVPHI |
||
| 1126 | DOUBLE PRECISION LRP, TWO, AL, DL, AI, BI, CI, DI, EPS |
||
| 1127 | PARAMETER ( LRP = -.22579135264472743235D0, TWO = 2, EPS = 1D-14 ) |
||
| 1128 | * LRP = LOG( SQRT( 2/PI ) ) |
||
| 1129 | RR = R*R |
||
| 1130 | IF ( N .LT. 2 ) THEN |
||
| 1131 | CHI = 2*MVPHI(-R) |
||
| 1132 | ELSE IF ( N .LT. 100 ) THEN |
||
| 1133 | * |
||
| 1134 | * Use standard Chi series |
||
| 1135 | * |
||
| 1136 | RN = 1 |
||
| 1137 | DO I = N - 2, 2, -2 |
||
| 1138 | RN = 1 + RR*RN/I |
||
| 1139 | END DO |
||
| 1140 | RR = RR/2 |
||
| 1141 | IF ( MOD( N, 2 ) .EQ. 0 ) THEN |
||
| 1142 | CHI = EXP( LOG( RN ) - RR ) |
||
| 1143 | ELSE |
||
| 1144 | CHI = EXP( LRP + LOG( R*RN ) - RR ) + 2*MVPHI(-R) |
||
| 1145 | ENDIF |
||
| 1146 | ELSE |
||
| 1147 | RR = RR/2 |
||
| 1148 | AL = N/TWO |
||
| 1149 | CHI = EXP( -RR + AL*LOG(RR) + LKN + LOG(TWO)*( N - 2 )/2 ) |
||
| 1150 | IF ( RR .LT. AL + 1 ) THEN |
||
| 1151 | * |
||
| 1152 | * Use Incomplete Gamma series |
||
| 1153 | * |
||
| 1154 | DL = CHI |
||
| 1155 | DO I = 1, 1000 |
||
| 1156 | DL = DL*RR/( AL + I ) |
||
| 1157 | CHI = CHI + DL |
||
| 1158 | IF ( ABS( DL*RR/( AL + I + 1 - RR ) ) .LT. EPS ) GO TO 10 |
||
| 1159 | END DO |
||
| 1160 | 10 CHI = 1 - CHI/AL |
||
| 1161 | ELSE |
||
| 1162 | * |
||
| 1163 | * Use Incomplete Gamma continued fraction |
||
| 1164 | * |
||
| 1165 | BI = RR + 1 - AL |
||
| 1166 | CI = 1/EPS |
||
| 1167 | DI = BI |
||
| 1168 | CHI = CHI/BI |
||
| 1169 | DO I = 1, 250 |
||
| 1170 | AI = I*( AL - I ) |
||
| 1171 | BI = BI + 2 |
||
| 1172 | CI = BI + AI/CI |
||
| 1173 | IF ( CI .EQ. 0 ) CI = EPS |
||
| 1174 | DI = BI + AI/DI |
||
| 1175 | IF ( DI .EQ. 0 ) DI = EPS |
||
| 1176 | DL = CI/DI |
||
| 1177 | CHI = CHI*DL |
||
| 1178 | IF ( ABS( DL - 1 ) .LT. EPS ) GO TO 20 |
||
| 1179 | END DO |
||
| 1180 | END IF |
||
| 1181 | END IF |
||
| 1182 | 20 DF = ( P - CHI )/EXP( LKN + ( N - 1 )*LOG(R) - RR ) |
||
| 1183 | MVCHNC = R - DF*( 1 - DF*( R - ( N - 1 )/R )/2 ) |
||
| 1184 | END |
||
| 1185 | * |
||
| 1186 | SUBROUTINE MVKBRV( NDIM, MINVLS, MAXVLS, NF, FUNSUB, |
||
| 1187 | & ABSEPS, RELEPS, ABSERR, FINEST, INFORM ) |
||
| 1188 | * |
||
| 1189 | * Automatic Multidimensional Integration Subroutine |
||
| 1190 | * |
||
| 1191 | * AUTHOR: Alan Genz |
||
| 1192 | * Department of Mathematics |
||
| 1193 | * Washington State University |
||
| 1194 | * Pulman, WA 99164-3113 |
||
| 1195 | * Email: AlanGenz@wsu.edu |
||
| 1196 | * |
||
| 1197 | * Last Change: 12/15/00 |
||
| 1198 | * |
||
| 1199 | * MVKBRV computes an approximation to the integral |
||
| 1200 | * |
||
| 1201 | * 1 1 1 |
||
| 1202 | * I I ... I F(X) dx(NDIM)...dx(2)dx(1) |
||
| 1203 | * 0 0 0 |
||
| 1204 | * |
||
| 1205 | * F(X) is a real NF-vector of integrands. |
||
| 1206 | * |
||
| 1207 | * It uses randomized Korobov rules. The primary references are |
||
| 1208 | * "Randomization of Number Theoretic Methods for Multiple Integration" |
||
| 1209 | * R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13, pp. 904-14, |
||
| 1210 | * and |
||
| 1211 | * "Optimal Parameters for Multidimensional Integration", |
||
| 1212 | * P. Keast, SIAM J Numer Anal, 10, pp.831-838. |
||
| 1213 | * If there are more than 100 variables, the remaining variables are |
||
| 1214 | * integrated using the rules described in the reference |
||
| 1215 | * "On a Number-Theoretical Integration Method" |
||
| 1216 | * H. Niederreiter, Aequationes Mathematicae, 8(1972), pp. 304-11. |
||
| 1217 | * |
||
| 1218 | *************** Parameters ******************************************** |
||
| 1219 | ****** Input parameters |
||
| 1220 | * NDIM Number of variables, must exceed 1, but not exceed 100 |
||
| 1221 | * MINVLS Integer minimum number of function evaluations allowed. |
||
| 1222 | * MINVLS must not exceed MAXVLS. If MINVLS < 0 then the |
||
| 1223 | * routine assumes a previous call has been made with |
||
| 1224 | * the same integrands and continues that calculation. |
||
| 1225 | * MAXVLS Integer maximum number of function evaluations allowed. |
||
| 1226 | * NF Number of integrands, must exceed 1, but not exceed 5000 |
||
| 1227 | * FUNSUB EXTERNALly declared user defined integrand subroutine. |
||
| 1228 | * It must have parameters ( NDIM, Z, NF, FUNVLS ), where |
||
| 1229 | * Z is a real NDIM-vector and FUNVLS is a real NF-vector. |
||
| 1230 | * |
||
| 1231 | * ABSEPS Required absolute accuracy. |
||
| 1232 | * RELEPS Required relative accuracy. |
||
| 1233 | ****** Output parameters |
||
| 1234 | * MINVLS Actual number of function evaluations used. |
||
| 1235 | * ABSERR Maximum norm of estimated absolute accuracy of FINEST. |
||
| 1236 | * FINEST Estimated NF-vector of values of the integrals. |
||
| 1237 | * INFORM INFORM = 0 for normal exit, when |
||
| 1238 | * ABSERR <= MAX(ABSEPS, RELEPS*||FINEST||) |
||
| 1239 | * and |
||
| 1240 | * INTVLS <= MAXCLS. |
||
| 1241 | * INFORM = 1 If MAXVLS was too small to obtain the required |
||
| 1242 | * accuracy. In this case a value FINEST is returned with |
||
| 1243 | * estimated absolute accuracy ABSERR. |
||
| 1244 | ************************************************************************ |
||
| 1245 | EXTERNAL FUNSUB |
||
| 1246 | DOUBLE PRECISION ABSEPS, RELEPS, FINEST(*), ABSERR, ONE |
||
| 1247 | INTEGER NDIM, NF, MINVLS, MAXVLS, INFORM, NP, PLIM, KLIM, |
||
| 1248 | & NLIM, FLIM, SAMPLS, I, K, INTVLS, MINSMP, KMX |
||
| 1249 | PARAMETER ( PLIM = 28, NLIM = 1000, KLIM = 100, FLIM = 5000 ) |
||
| 1250 | PARAMETER ( MINSMP = 8 ) |
||
| 1251 | INTEGER P(PLIM), C(PLIM,KLIM-1), PR(NLIM) |
||
| 1252 | DOUBLE PRECISION DIFINT, FINVAL(FLIM), VARSQR(FLIM), VAREST(FLIM), |
||
| 1253 | & VARPRD, X(NLIM), R(NLIM), VK(NLIM), VALUES(FLIM), FS(FLIM) |
||
| 1254 | PARAMETER ( ONE = 1 ) |
||
| 1255 | SAVE P, C, SAMPLS, NP, VAREST |
||
| 1256 | INFORM = 1 |
||
| 1257 | INTVLS = 0 |
||
| 1258 | VARPRD = 0 |
||
| 1259 | IF ( MINVLS .GE. 0 ) THEN |
||
| 1260 | DO K = 1, NF |
||
| 1261 | FINEST(K) = 0 |
||
| 1262 | VAREST(K) = 0 |
||
| 1263 | END DO |
||
| 1264 | SAMPLS = MINSMP |
||
| 1265 | DO I = MIN( NDIM, 10 ), PLIM |
||
| 1266 | NP = I |
||
| 1267 | IF ( MINVLS .LT. 2*SAMPLS*P(I) ) GO TO 10 |
||
| 1268 | END DO |
||
| 1269 | SAMPLS = MAX( MINSMP, MINVLS/( 2*P(NP) ) ) |
||
| 1270 | ENDIF |
||
| 1271 | 10 VK(1) = ONE/P(NP) |
||
| 1272 | K = 1 |
||
| 1273 | DO I = 2, NDIM |
||
| 1274 | IF ( I .LE. KLIM ) THEN |
||
| 1275 | K = MOD( C(NP, MIN(NDIM-1,KLIM-1))*DBLE(K), DBLE(P(NP)) ) |
||
| 1276 | VK(I) = K*VK(1) |
||
| 1277 | ELSE |
||
| 1278 | VK(I) = INT( P(NP)*2**( DBLE(I-KLIM)/(NDIM-KLIM+1) ) ) |
||
| 1279 | VK(I) = MOD( VK(I)/P(NP), ONE ) |
||
| 1280 | END IF |
||
| 1281 | END DO |
||
| 1282 | DO K = 1, NF |
||
| 1283 | FINVAL(K) = 0 |
||
| 1284 | VARSQR(K) = 0 |
||
| 1285 | END DO |
||
| 1286 | * |
||
| 1287 | DO I = 1, SAMPLS |
||
| 1288 | CALL MVKRSV( NDIM,KLIM,VALUES, P(NP),VK, NF,FUNSUB, X,R,PR,FS ) |
||
| 1289 | DO K = 1, NF |
||
| 1290 | DIFINT = ( VALUES(K) - FINVAL(K) )/I |
||
| 1291 | FINVAL(K) = FINVAL(K) + DIFINT |
||
| 1292 | VARSQR(K) = ( I - 2 )*VARSQR(K)/I + DIFINT**2 |
||
| 1293 | END DO |
||
| 1294 | END DO |
||
| 1295 | * |
||
| 1296 | INTVLS = INTVLS + 2*SAMPLS*P(NP) |
||
| 1297 | KMX = 1 |
||
| 1298 | DO K = 1, NF |
||
| 1299 | VARPRD = VAREST(K)*VARSQR(K) |
||
| 1300 | FINEST(K) = FINEST(K) + ( FINVAL(K) - FINEST(K) )/( 1+VARPRD ) |
||
| 1301 | IF ( VARSQR(K) .GT. 0 ) VAREST(K) = ( 1 + VARPRD )/VARSQR(K) |
||
| 1302 | IF ( ABS(FINEST(K)) .GT. ABS(FINEST(KMX)) ) KMX = K |
||
| 1303 | END DO |
||
| 1304 | ABSERR = 7*SQRT( VARSQR(KMX)/( 1 + VARPRD ) )/2 |
||
| 1305 | IF ( ABSERR .GT. MAX( ABSEPS, ABS(FINEST(KMX))*RELEPS ) ) THEN |
||
| 1306 | IF ( NP .LT. PLIM ) THEN |
||
| 1307 | NP = NP + 1 |
||
| 1308 | ELSE |
||
| 1309 | SAMPLS = MIN( 3*SAMPLS/2, ( MAXVLS - INTVLS )/( 2*P(NP) ) ) |
||
| 1310 | SAMPLS = MAX( MINSMP, SAMPLS ) |
||
| 1311 | ENDIF |
||
| 1312 | IF ( INTVLS + 2*SAMPLS*P(NP) .LE. MAXVLS ) GO TO 10 |
||
| 1313 | ELSE |
||
| 1314 | INFORM = 0 |
||
| 1315 | ENDIF |
||
| 1316 | MINVLS = INTVLS |
||
| 1317 | * |
||
| 1318 | * Optimal Parameters for Lattice Rules |
||
| 1319 | * |
||
| 1320 | DATA P( 1),(C( 1,I),I = 1,99)/ 31, 12, 2*9, 13, 8*12, 3*3, 12, |
||
| 1321 | & 2*7, 9*12, 3*3, 12, 2*7, 9*12, 3*3, 12, 2*7, 9*12, 3*3, 12, 2*7, |
||
| 1322 | & 8*12, 7, 3*3, 3*7, 21*3/ |
||
| 1323 | DATA P( 2),(C( 2,I),I = 1,99)/ 47, 13, 11, 17, 10, 6*15, |
||
| 1324 | & 22, 2*15, 3*6, 2*15, 9, 13, 3*2, 13, 2*11, 10, 9*15, 3*6, 2*15, |
||
| 1325 | & 9, 13, 3*2, 13, 2*11, 10, 9*15, 3*6, 2*15, 9, 13, 3*2, 13, 2*11, |
||
| 1326 | & 2*10, 8*15, 6, 2, 3, 2, 3, 12*2/ |
||
| 1327 | DATA P( 3),(C( 3,I),I = 1,99)/ 73, 27, 28, 10, 2*11, 20, |
||
| 1328 | & 2*11, 28, 2*13, 28, 3*13, 16*14, 2*31, 3*5, 31, 13, 6*11, 7*13, |
||
| 1329 | & 16*14, 2*31, 3*5, 11, 13, 7*11, 2*13, 11, 13, 4*5, 14, 13, 8*5/ |
||
| 1330 | DATA P( 4),(C( 4,I),I = 1,99)/ 113, 35, 2*27, 36, 22, 2*29, |
||
| 1331 | & 20, 45, 3*5, 16*21, 29, 10*17, 12*23, 21, 27, 3*3, 24, 2*27, |
||
| 1332 | & 17, 3*29, 17, 4*5, 16*21, 3*17, 6, 2*17, 6, 3, 2*6, 5*3/ |
||
| 1333 | DATA P( 5),(C( 5,I),I = 1,99)/ 173, 64, 66, 2*28, 2*44, 55, |
||
| 1334 | & 67, 6*10, 2*38, 5*10, 12*49, 2*38, 31, 2*4, 31, 64, 3*4, 64, |
||
| 1335 | & 6*45, 19*66, 11, 9*66, 45, 11, 7, 3, 3*2, 27, 5, 2*3, 2*5, 7*2/ |
||
| 1336 | DATA P( 6),(C( 6,I),I = 1,99)/ 263, 111, 42, 54, 118, 20, |
||
| 1337 | & 2*31, 72, 17, 94, 2*14, 11, 3*14, 94, 4*10, 7*14, 3*11, 7*8, |
||
| 1338 | & 5*18, 113, 2*62, 2*45, 17*113, 2*63, 53, 63, 15*67, 5*51, 12, |
||
| 1339 | & 51, 12, 51, 5, 2*3, 2*2, 5/ |
||
| 1340 | DATA P( 7),(C( 7,I),I = 1,99)/ 397, 163, 154, 83, 43, 82, |
||
| 1341 | & 92, 150, 59, 2*76, 47, 2*11, 100, 131, 6*116, 9*138, 21*101, |
||
| 1342 | & 6*116, 5*100, 5*138, 19*101, 8*38, 5*3/ |
||
| 1343 | DATA P( 8),(C( 8,I),I = 1,99)/ 593, 246, 189, 242, 102, |
||
| 1344 | & 2*250, 102, 250, 280, 118, 196, 118, 191, 215, 2*121, |
||
| 1345 | & 12*49, 34*171, 8*161, 17*14, 6*10, 103, 4*10, 5/ |
||
| 1346 | DATA P( 9),(C( 9,I),I = 1,99)/ 907, 347, 402, 322, 418, |
||
| 1347 | & 215, 220, 3*339, 337, 218, 4*315, 4*167, 361, 201, 11*124, |
||
| 1348 | & 2*231, 14*90, 4*48, 23*90, 10*243, 9*283, 16, 283, 16, 2*283/ |
||
| 1349 | DATA P(10),(C(10,I),I = 1,99)/ 1361, 505, 220, 601, 644, |
||
| 1350 | & 612, 160, 3*206, 422, 134, 518, 2*134, 518, 652, 382, |
||
| 1351 | & 206, 158, 441, 179, 441, 56, 2*559, 14*56, 2*101, 56, |
||
| 1352 | & 8*101, 7*193, 21*101, 17*122, 4*101/ |
||
| 1353 | DATA P(11),(C(11,I),I = 1,99)/ 2053, 794, 325, 960, 528, |
||
| 1354 | & 2*247, 338, 366, 847, 2*753, 236, 2*334, 461, 711, 652, |
||
| 1355 | & 3*381, 652, 7*381, 226, 7*326, 126, 10*326, 2*195, 19*55, |
||
| 1356 | & 7*195, 11*132, 13*387/ |
||
| 1357 | DATA P(12),(C(12,I),I = 1,99)/ 3079, 1189, 888, 259, 1082, 725, |
||
| 1358 | & 811, 636, 965, 2*497, 2*1490, 392, 1291, 2*508, 2*1291, 508, |
||
| 1359 | & 1291, 2*508, 4*867, 934, 7*867, 9*1284, 4*563, 3*1010, 208, |
||
| 1360 | & 838, 3*563, 2*759, 564, 2*759, 4*801, 5*759, 8*563, 22*226/ |
||
| 1361 | DATA P(13),(C(13,I),I = 1,99)/ 4621, 1763, 1018, 1500, 432, |
||
| 1362 | & 1332, 2203, 126, 2240, 1719, 1284, 878, 1983, 4*266, |
||
| 1363 | & 2*747, 2*127, 2074, 127, 2074, 1400, 10*1383, 1400, 7*1383, |
||
| 1364 | & 507, 4*1073, 5*1990, 9*507, 17*1073, 6*22, 1073, 6*452, 318, |
||
| 1365 | & 4*301, 2*86, 15/ |
||
| 1366 | DATA P(14),(C(14,I),I = 1,99)/ 6947, 2872, 3233, 1534, 2941, |
||
| 1367 | & 2910, 393, 1796, 919, 446, 2*919, 1117, 7*103, 2311, 3117, 1101, |
||
| 1368 | & 2*3117, 5*1101, 8*2503, 7*429, 3*1702, 5*184, 34*105, 13*784/ |
||
| 1369 | DATA P(15),(C(15,I),I = 1,99)/ 10427, 4309, 3758, 4034, 1963, |
||
| 1370 | & 730, 642, 1502, 2246, 3834, 1511, 2*1102, 2*1522, 2*3427, |
||
| 1371 | & 3928, 2*915, 4*3818, 3*4782, 3818, 4782, 2*3818, 7*1327, 9*1387, |
||
| 1372 | & 13*2339, 18*3148, 3*1776, 3*3354, 925, 2*3354, 5*925, 8*2133/ |
||
| 1373 | DATA P(16),(C(16,I),I = 1,99)/ 15641, 6610, 6977, 1686, 3819, |
||
| 1374 | & 2314, 5647, 3953, 3614, 5115, 2*423, 5408, 7426, 2*423, |
||
| 1375 | & 487, 6227, 2660, 6227, 1221, 3811, 197, 4367, 351, |
||
| 1376 | & 1281, 1221, 3*351, 7245, 1984, 6*2999, 3995, 4*2063, 1644, |
||
| 1377 | & 2063, 2077, 3*2512, 4*2077, 19*754, 2*1097, 4*754, 248, 754, |
||
| 1378 | & 4*1097, 4*222, 754,11*1982/ |
||
| 1379 | DATA P(17),(C(17,I),I = 1,99)/ 23473, 9861, 3647, 4073, 2535, |
||
| 1380 | & 3430, 9865, 2830, 9328, 4320, 5913, 10365, 8272, 3706, 6186, |
||
| 1381 | & 3*7806, 8610, 2563, 2*11558, 9421, 1181, 9421, 3*1181, 9421, |
||
| 1382 | & 2*1181, 2*10574, 5*3534, 3*2898, 3450, 7*2141, 15*7055, 2831, |
||
| 1383 | & 24*8204, 3*4688, 8*2831/ |
||
| 1384 | DATA P(18),(C(18,I),I = 1,99)/ 35221, 10327, 7582, 7124, 8214, |
||
| 1385 | & 9600, 10271, 10193, 10800, 9086, 2365, 4409, 13812, |
||
| 1386 | & 5661, 2*9344, 10362, 2*9344, 8585, 11114, 3*13080, 6949, |
||
| 1387 | & 3*3436, 13213, 2*6130, 2*8159, 11595, 8159, 3436, 18*7096, |
||
| 1388 | & 4377, 7096, 5*4377, 2*5410, 32*4377, 2*440, 3*1199/ |
||
| 1389 | DATA P(19),(C(19,I),I = 1,99)/ 52837, 19540, 19926, 11582, |
||
| 1390 | & 11113, 24585, 8726, 17218, 419, 3*4918, 15701, 17710, |
||
| 1391 | & 2*4037, 15808, 11401, 19398, 2*25950, 4454, 24987, 11719, |
||
| 1392 | & 8697, 5*1452, 2*8697, 6436, 21475, 6436, 22913, 6434, 18497, |
||
| 1393 | & 4*11089, 2*3036, 4*14208, 8*12906, 4*7614, 6*5021, 24*10145, |
||
| 1394 | & 6*4544, 4*8394/ |
||
| 1395 | DATA P(20),(C(20,I),I = 1,99)/ 79259, 34566, 9579, 12654, |
||
| 1396 | & 26856, 37873, 38806, 29501, 17271, 3663, 10763, 18955, |
||
| 1397 | & 1298, 26560, 2*17132, 2*4753, 8713, 18624, 13082, 6791, |
||
| 1398 | & 1122, 19363, 34695, 4*18770, 15628, 4*18770, 33766, 6*20837, |
||
| 1399 | & 5*6545, 14*12138, 5*30483, 19*12138, 9305, 13*11107, 2*9305/ |
||
| 1400 | DATA P(21),(C(21,I),I = 1,99)/118891, 31929, 49367, 10982, 3527, |
||
| 1401 | & 27066, 13226, 56010, 18911, 40574, 2*20767, 9686, 2*47603, |
||
| 1402 | & 2*11736, 41601, 12888, 32948, 30801, 44243, 2*53351, 16016, |
||
| 1403 | & 2*35086, 32581, 2*2464, 49554, 2*2464, 2*49554, 2464, 81, 27260, |
||
| 1404 | & 10681, 7*2185, 5*18086, 2*17631, 3*18086, 37335, 3*37774, |
||
| 1405 | & 13*26401, 12982, 6*40398, 3*3518, 9*37799, 4*4721, 4*7067/ |
||
| 1406 | DATA P(22),(C(22,I),I = 1,99)/178349, 40701, 69087, 77576, 64590, |
||
| 1407 | & 39397, 33179, 10858, 38935, 43129, 2*35468, 5279, 2*61518, 27945, |
||
| 1408 | & 2*70975, 2*86478, 2*20514, 2*73178, 2*43098, 4701, |
||
| 1409 | & 2*59979, 58556, 69916, 2*15170, 2*4832, 43064, 71685, 4832, |
||
| 1410 | & 3*15170, 3*27679, 2*60826, 2*6187, 5*4264, 45567, 4*32269, |
||
| 1411 | & 9*62060, 13*1803, 12*51108, 2*55315, 5*54140, 13134/ |
||
| 1412 | DATA P(23),(C(23,I),I = 1,99)/267523, 103650, 125480, 59978, |
||
| 1413 | & 46875, 77172, 83021, 126904, 14541, 56299, 43636, 11655, |
||
| 1414 | & 52680, 88549, 29804, 101894, 113675, 48040, 113675, |
||
| 1415 | & 34987, 48308, 97926, 5475, 49449, 6850, 2*62545, 9440, |
||
| 1416 | & 33242, 9440, 33242, 9440, 33242, 9440, 62850, 3*9440, |
||
| 1417 | & 3*90308, 9*47904, 7*41143, 5*36114, 24997, 14*65162, 7*47650, |
||
| 1418 | & 7*40586, 4*38725, 5*88329/ |
||
| 1419 | DATA P(24),(C(24,I),I = 1,99)/401287, 165843, 90647, 59925, |
||
| 1420 | & 189541, 67647, 74795, 68365, 167485, 143918, 74912, |
||
| 1421 | & 167289, 75517, 8148, 172106, 126159,3*35867, 121694, |
||
| 1422 | & 52171, 95354, 2*113969, 76304, 2*123709, 144615, 123709, |
||
| 1423 | & 2*64958, 32377, 2*193002, 25023, 40017, 141605, 2*189165, |
||
| 1424 | & 141605, 2*189165, 3*141605, 189165, 20*127047, 10*127785, |
||
| 1425 | & 6*80822, 16*131661, 7114, 131661/ |
||
| 1426 | DATA P(25),(C(25,I),I = 1,99)/601943, 130365, 236711, 110235, |
||
| 1427 | & 125699, 56483, 93735, 234469, 60549, 1291, 93937, |
||
| 1428 | & 245291, 196061, 258647, 162489, 176631, 204895, 73353, |
||
| 1429 | & 172319, 28881, 136787,2*122081, 275993, 64673, 3*211587, |
||
| 1430 | & 2*282859, 211587, 242821, 3*256865, 122203, 291915, 122203, |
||
| 1431 | & 2*291915, 122203, 2*25639, 291803, 245397, 284047, |
||
| 1432 | & 7*245397, 94241, 2*66575, 19*217673, 10*210249, 15*94453/ |
||
| 1433 | DATA P(26),(C(26,I),I = 1,99)/902933, 333459, 375354, 102417, |
||
| 1434 | & 383544, 292630, 41147, 374614, 48032, 435453, 281493, 358168, |
||
| 1435 | & 114121, 346892, 238990, 317313, 164158, 35497, 2*70530, 434839, |
||
| 1436 | & 3*24754, 393656, 2*118711, 148227, 271087, 355831, 91034, |
||
| 1437 | & 2*417029, 2*91034, 417029, 91034, 2*299843, 2*413548, 308300, |
||
| 1438 | & 3*413548, 3*308300, 413548, 5*308300, 4*15311, 2*176255, 6*23613, |
||
| 1439 | & 172210, 4* 204328, 5*121626, 5*200187, 2*121551, 12*248492, |
||
| 1440 | & 5*13942/ |
||
| 1441 | DATA P(27), (C(27,I), I = 1,99)/ 1354471, 500884, 566009, 399251, |
||
| 1442 | & 652979, 355008, 430235, 328722, 670680, 2*405585, 424646, |
||
| 1443 | & 2*670180, 641587, 215580, 59048, 633320, 81010, 20789, 2*389250, |
||
| 1444 | & 2*638764, 2*389250, 398094, 80846, 2*147776, 296177, 2*398094, |
||
| 1445 | & 2*147776, 396313, 3*578233, 19482, 620706, 187095, 620706, |
||
| 1446 | & 187095, 126467, 12*241663, 321632, 2*23210, 3*394484, 3*78101, |
||
| 1447 | & 19*542095, 3*277743, 12*457259/ |
||
| 1448 | DATA P(28), (C(28,I), I = 1, 99)/ 2031713, 858339, 918142, 501970, |
||
| 1449 | & 234813, 460565, 31996, 753018, 256150, 199809, 993599, 245149, |
||
| 1450 | & 794183, 121349, 150619, 376952, 2*809123, 804319, 67352, 969594, |
||
| 1451 | & 434796, 969594, 804319, 391368, 761041, 754049, 466264, 2*754049, |
||
| 1452 | & 466264, 2*754049, 282852, 429907, 390017, 276645, 994856, 250142, |
||
| 1453 | & 144595, 907454, 689648, 4*687580, 978368, 687580, 552742, 105195, |
||
| 1454 | & 942843, 768249, 4*307142, 7*880619, 11*117185, 11*60731, |
||
| 1455 | & 4*178309, 8*74373, 3*214965/ |
||
| 1456 | * |
||
| 1457 | END |
||
| 1458 | * |
||
| 1459 | SUBROUTINE MVKRSV( NDIM,KL,VALUES,PRIME,VK, NF,FUNSUB, X,R,PR,FS ) |
||
| 1460 | * |
||
| 1461 | * For lattice rule sums |
||
| 1462 | * |
||
| 1463 | INTEGER NDIM, NF, PRIME, KL, K, J, JP, PR(*) |
||
| 1464 | DOUBLE PRECISION VALUES(*), VK(*), FS(*), X(*), R(*), MVUNI |
||
| 1465 | DO J = 1, NF |
||
| 1466 | VALUES(J) = 0 |
||
| 1467 | END DO |
||
| 1468 | * |
||
| 1469 | * Determine random shifts for each variable; scramble lattice rule |
||
| 1470 | * |
||
| 1471 | DO J = 1, NDIM |
||
| 1472 | R(J) = MVUNI() |
||
| 1473 | IF ( J .LT. KL ) THEN |
||
| 1474 | JP = 1 + J*R(J) |
||
| 1475 | IF ( JP .LT. J ) PR(J) = PR(JP) |
||
| 1476 | PR(JP) = J |
||
| 1477 | ELSE |
||
| 1478 | PR(J) = J |
||
| 1479 | END IF |
||
| 1480 | END DO |
||
| 1481 | * |
||
| 1482 | * Compute latice rule sums |
||
| 1483 | * |
||
| 1484 | DO K = 1, PRIME |
||
| 1485 | DO J = 1, NDIM |
||
| 1486 | R(J) = R(J) + VK(PR(J)) |
||
| 1487 | IF ( R(J) .GT. 1 ) R(J) = R(J) - 1 |
||
| 1488 | X(J) = ABS( 2*R(J) - 1 ) |
||
| 1489 | END DO |
||
| 1490 | CALL FUNSUB( NDIM, X, NF, FS ) |
||
| 1491 | DO J = 1, NF |
||
| 1492 | VALUES(J) = VALUES(J) + ( FS(J) - VALUES(J) )/( 2*K-1 ) |
||
| 1493 | END DO |
||
| 1494 | DO J = 1, NDIM |
||
| 1495 | X(J) = 1 - X(J) |
||
| 1496 | END DO |
||
| 1497 | CALL FUNSUB( NDIM, X, NF, FS ) |
||
| 1498 | DO J = 1, NF |
||
| 1499 | VALUES(J) = VALUES(J) + ( FS(J) - VALUES(J) )/( 2*K ) |
||
| 1500 | END DO |
||
| 1501 | END DO |
||
| 1502 | * |
||
| 1503 | END |
||
| 1504 | * |
||
| 1505 | DOUBLE PRECISION FUNCTION MVUNI() |
||
| 1506 | * |
||
| 1507 | * Uniform (0,1) random number generator |
||
| 1508 | * |
||
| 1509 | * use R's random number generator directly |
||
| 1510 | * the way `Writing R extentions' advertises. |
||
| 1511 | * |
||
| 1512 | DOUBLE PRECISION unifrnd, x |
||
| 1513 | |||
| 1514 | x = unifrnd() |
||
| 1515 | MVUNI = x |
||
| 1516 | END |