1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
|
* Note: The test program has been removed and a utlity routine mvnun has been
* added. RTK 2004-08-10
*
* Copyright 2000 by Alan Genz.
* Copyright 2004-2005 by Enthought, Inc.
*
* The subroutine MVNUN is copyrighted by Enthought, Inc.
* The rest of the file is copyrighted by Alan Genz and has kindly been offered
* to the Scipy project under it's BSD-style license.
*
* This file contains a short test program and MVNDST, a subroutine
* for computing multivariate normal distribution function values.
* The file is self contained and should compile without errors on (77)
* standard Fortran compilers. The test program demonstrates the use of
* MVNDST for computing MVN distribution values for a five dimensional
* example problem, with three different integration limit combinations.
*
* Alan Genz
* Department of Mathematics
* Washington State University
* Pullman, WA 99164-3113
* Email : alangenz@wsu.edu
*
SUBROUTINE mvnun(d, n, lower, upper, means, covar, maxpts,
& abseps, releps, value, inform)
* Parameters
*
* d integer, dimensionality of the data
* n integer, the number of data points
* lower double(2), the lower integration limits
* upper double(2), the upper integration limits
* means double(n), the mean of each kernel
* covar double(2,2), the covariance matrix
* maxpts integer, the maximum number of points to evaluate at
* abseps double, absolute error tolerance
* releps double, relative error tolerance
* value double intent(out), integral value
* inform integer intent(out),
* if inform == 0: error < eps
* elif inform == 1: error > eps, all maxpts used
integer n, d, infin(d), maxpts, inform, tmpinf
double precision lower(d), upper(d), releps, abseps,
& error, value, stdev(d), rho(d*(d-1)/2),
& covar(d,d),
& nlower(d), nupper(d), means(d,n), tmpval
integer i, j
do i=1,d
stdev(i) = dsqrt(covar(i,i))
infin(i) = 2
end do
do i=1,d
do j=1,i-1
rho(j+(i-2)*(i-1)/2) = covar(i,j)/stdev(i)/stdev(j)
end do
end do
value = 0d0
inform = 0
do i=1,n
do j=1,d
nlower(j) = (lower(j) - means(j,i))/stdev(j)
nupper(j) = (upper(j) - means(j,i))/stdev(j)
end do
call mvndst(d,nlower,nupper,infin,rho,maxpts,abseps,releps,
& error,tmpval,tmpinf)
value = value + tmpval
if (tmpinf .eq. 1) then
inform = 1
end if
end do
value = value / n
END
SUBROUTINE MVNDST( N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
& ABSEPS, RELEPS, ERROR, VALUE, INFORM )
*
* A subroutine for computing multivariate normal probabilities.
* This subroutine uses an algorithm given in the paper
* "Numerical Computation of Multivariate Normal Probabilities", in
* J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
* Alan Genz
* Department of Mathematics
* Washington State University
* Pullman, WA 99164-3113
* Email : AlanGenz@wsu.edu
*
* Parameters
*
* N INTEGER, the number of variables.
* LOWER REAL, array of lower integration limits.
* UPPER REAL, array of upper integration limits.
* INFIN INTEGER, array of integration limits flags:
* if INFIN(I) < 0, Ith limits are (-infinity, infinity);
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
* CORREL REAL, array of correlation coefficients; the correlation
* coefficient in row I column J of the correlation matrix
* should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
* THe correlation matrix must be positive semidefinite.
* MAXPTS INTEGER, maximum number of function values allowed. This
* parameter can be used to limit the time. A sensible
* strategy is to start with MAXPTS = 1000*N, and then
* increase MAXPTS if ERROR is too large.
* ABSEPS REAL absolute error tolerance.
* RELEPS REAL relative error tolerance.
* ERROR REAL estimated absolute error, with 99% confidence level.
* VALUE REAL estimated value for the integral
* INFORM INTEGER, termination status parameter:
* if INFORM = 0, normal completion with ERROR < EPS;
* if INFORM = 1, completion with ERROR > EPS and MAXPTS
* function vaules used; increase MAXPTS to
* decrease ERROR;
* if INFORM = 2, N > 500 or N < 1.
*
EXTERNAL MVNDFN
INTEGER N, INFIN(*), MAXPTS, INFORM, INFIS, IVLS
DOUBLE PRECISION CORREL(*), LOWER(*), UPPER(*), RELEPS, ABSEPS,
& ERROR, VALUE, E, D, MVNDNT, MVNDFN
COMMON /DKBLCK/IVLS
IF ( N .GT. 500 .OR. N .LT. 1 ) THEN
INFORM = 2
VALUE = 0
ERROR = 1
ELSE
INFORM = MVNDNT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E)
IF ( N-INFIS .EQ. 0 ) THEN
VALUE = 1
ERROR = 0
ELSE IF ( N-INFIS .EQ. 1 ) THEN
VALUE = E - D
ERROR = 2D-16
ELSE
*
* Call the lattice rule integration subroutine
*
IVLS = 0
CALL DKBVRC( N-INFIS-1, IVLS, MAXPTS, MVNDFN,
& ABSEPS, RELEPS, ERROR, VALUE, INFORM )
ENDIF
ENDIF
END
DOUBLE PRECISION FUNCTION MVNDFN( N, W )
*
* Integrand subroutine
*
INTEGER N, INFIN(*), INFIS, NL
DOUBLE PRECISION W(*), LOWER(*), UPPER(*), CORREL(*), D, E
PARAMETER ( NL = 500 )
DOUBLE PRECISION COV(NL*(NL+1)/2), A(NL), B(NL), Y(NL)
INTEGER INFI(NL), I, J, IJ, IK, INFA, INFB
DOUBLE PRECISION SUM, AI, BI, DI, EI, PHINVS, BVNMVN, MVNDNT
SAVE A, B, INFI, COV
MVNDFN = 1
INFA = 0
INFB = 0
IK = 1
IJ = 0
DO I = 1, N+1
SUM = 0
DO J = 1, I-1
IJ = IJ + 1
IF ( J .LT. IK ) SUM = SUM + COV(IJ)*Y(J)
END DO
IF ( INFI(I) .NE. 0 ) THEN
IF ( INFA .EQ. 1 ) THEN
AI = MAX( AI, A(I) - SUM )
ELSE
AI = A(I) - SUM
INFA = 1
END IF
END IF
IF ( INFI(I) .NE. 1 ) THEN
IF ( INFB .EQ. 1 ) THEN
BI = MIN( BI, B(I) - SUM )
ELSE
BI = B(I) - SUM
INFB = 1
END IF
END IF
IJ = IJ + 1
IF ( I .EQ. N+1 .OR. COV(IJ+IK+1) .GT. 0 ) THEN
CALL MVNLMS( AI, BI, 2*INFA+INFB-1, DI, EI )
IF ( DI .GE. EI ) THEN
MVNDFN = 0
RETURN
ELSE
MVNDFN = MVNDFN*( EI - DI )
IF ( I .LE. N ) Y(IK) = PHINVS( DI + W(IK)*( EI - DI ) )
IK = IK + 1
INFA = 0
INFB = 0
END IF
END IF
END DO
RETURN
*
* Entry point for intialization.
*
ENTRY MVNDNT( N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E )
MVNDNT = 0
*
* Initialization and computation of covariance Cholesky factor.
*
CALL COVSRT( N, LOWER,UPPER,CORREL,INFIN,Y, INFIS,A,B,COV,INFI )
IF ( N - INFIS .EQ. 1 ) THEN
CALL MVNLMS( A(1), B(1), INFI(1), D, E )
ELSE IF ( N - INFIS .EQ. 2 ) THEN
IF ( ABS( COV(3) ) .GT. 0 ) THEN
D = SQRT( 1 + COV(2)**2 )
IF ( INFI(2) .NE. 0 ) A(2) = A(2)/D
IF ( INFI(2) .NE. 1 ) B(2) = B(2)/D
E = BVNMVN( A, B, INFI, COV(2)/D )
D = 0
ELSE
IF ( INFI(1) .NE. 0 ) THEN
IF ( INFI(2) .NE. 0 ) A(1) = MAX( A(1), A(2) )
ELSE
IF ( INFI(2) .NE. 0 ) A(1) = A(2)
END IF
IF ( INFI(1) .NE. 1 ) THEN
IF ( INFI(2) .NE. 1 ) B(1) = MIN( B(1), B(2) )
ELSE
IF ( INFI(2) .NE. 1 ) B(1) = B(2)
END IF
IF ( INFI(1) .NE. INFI(2) ) INFI(1) = 2
CALL MVNLMS( A(1), B(1), INFI(1), D, E )
END IF
INFIS = INFIS + 1
END IF
END
SUBROUTINE MVNLMS( A, B, INFIN, LOWER, UPPER )
DOUBLE PRECISION A, B, LOWER, UPPER, MVNPHI
INTEGER INFIN
LOWER = 0
UPPER = 1
IF ( INFIN .GE. 0 ) THEN
IF ( INFIN .NE. 0 ) LOWER = MVNPHI(A)
IF ( INFIN .NE. 1 ) UPPER = MVNPHI(B)
ENDIF
UPPER = MAX( UPPER, LOWER )
END
SUBROUTINE COVSRT( N, LOWER, UPPER, CORREL, INFIN, Y,
& INFIS, A, B, COV, INFI )
*
* Subroutine to sort integration limits and determine Cholesky factor.
*
INTEGER N, INFI(*), INFIN(*), INFIS
DOUBLE PRECISION
& A(*), B(*), COV(*), LOWER(*), UPPER(*), CORREL(*), Y(*)
INTEGER I, J, K, L, M, II, IJ, IL, JMIN
DOUBLE PRECISION SUMSQ, AJ, BJ, SUM, SQTWPI, EPS, D, E
DOUBLE PRECISION CVDIAG, AMIN, BMIN, DMIN, EMIN, YL, YU
PARAMETER ( SQTWPI = 2.506628274631001D0, EPS = 1D-10 )
IJ = 0
II = 0
INFIS = 0
DO I = 1, N
A(I) = 0
B(I) = 0
INFI(I) = INFIN(I)
IF ( INFI(I) .LT. 0 ) THEN
INFIS = INFIS + 1
ELSE
IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I)
IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I)
ENDIF
DO J = 1, I-1
IJ = IJ + 1
II = II + 1
COV(IJ) = CORREL(II)
END DO
IJ = IJ + 1
COV(IJ) = 1
END DO
*
* First move any doubly infinite limits to innermost positions.
*
IF ( INFIS .LT. N ) THEN
DO I = N, N-INFIS+1, -1
IF ( INFI(I) .GE. 0 ) THEN
DO J = 1,I-1
IF ( INFI(J) .LT. 0 ) THEN
CALL RCSWP( J, I, A, B, INFI, N, COV )
GO TO 10
ENDIF
END DO
ENDIF
10 END DO
*
* Sort remaining limits and determine Cholesky factor.
*
II = 0
DO I = 1, N-INFIS
*
* Determine the integration limits for variable with minimum
* expected probability and interchange that variable with Ith.
*
DMIN = 0
EMIN = 1
JMIN = I
CVDIAG = 0
IJ = II
DO J = I, N-INFIS
IF ( COV(IJ+J) .GT. EPS ) THEN
SUMSQ = SQRT( COV(IJ+J) )
SUM = 0
DO K = 1, I-1
SUM = SUM + COV(IJ+K)*Y(K)
END DO
AJ = ( A(J) - SUM )/SUMSQ
BJ = ( B(J) - SUM )/SUMSQ
CALL MVNLMS( AJ, BJ, INFI(J), D, E )
IF ( EMIN + D .GE. E + DMIN ) THEN
JMIN = J
AMIN = AJ
BMIN = BJ
DMIN = D
EMIN = E
CVDIAG = SUMSQ
ENDIF
ENDIF
IJ = IJ + J
END DO
IF ( JMIN .GT. I ) CALL RCSWP( I, JMIN, A,B, INFI, N, COV )
COV(II+I) = CVDIAG
*
* Compute Ith column of Cholesky factor.
* Compute expected value for Ith integration variable and
* scale Ith covariance matrix row and limits.
*
IF ( CVDIAG .GT. 0 ) THEN
IL = II + I
DO L = I+1, N-INFIS
COV(IL+I) = COV(IL+I)/CVDIAG
IJ = II + I
DO J = I+1, L
COV(IL+J) = COV(IL+J) - COV(IL+I)*COV(IJ+I)
IJ = IJ + J
END DO
IL = IL + L
END DO
IF ( EMIN .GT. DMIN + EPS ) THEN
YL = 0
YU = 0
IF ( INFI(I) .NE. 0 ) YL = -EXP( -AMIN**2/2 )/SQTWPI
IF ( INFI(I) .NE. 1 ) YU = -EXP( -BMIN**2/2 )/SQTWPI
Y(I) = ( YU - YL )/( EMIN - DMIN )
ELSE
IF ( INFI(I) .EQ. 0 ) Y(I) = BMIN
IF ( INFI(I) .EQ. 1 ) Y(I) = AMIN
IF ( INFI(I) .EQ. 2 ) Y(I) = ( AMIN + BMIN )/2
END IF
DO J = 1, I
II = II + 1
COV(II) = COV(II)/CVDIAG
END DO
A(I) = A(I)/CVDIAG
B(I) = B(I)/CVDIAG
ELSE
IL = II + I
DO L = I+1, N-INFIS
COV(IL+I) = 0
IL = IL + L
END DO
*
* If the covariance matrix diagonal entry is zero,
* permute limits and/or rows, if necessary.
*
*
DO J = I-1, 1, -1
IF ( ABS( COV(II+J) ) .GT. EPS ) THEN
A(I) = A(I)/COV(II+J)
B(I) = B(I)/COV(II+J)
IF ( COV(II+J) .LT. 0 ) THEN
CALL DKSWAP( A(I), B(I) )
IF ( INFI(I) .NE. 2 ) INFI(I) = 1 - INFI(I)
END IF
DO L = 1, J
COV(II+L) = COV(II+L)/COV(II+J)
END DO
DO L = J+1, I-1
IF( COV((L-1)*L/2+J+1) .GT. 0 ) THEN
IJ = II
DO K = I-1, L, -1
DO M = 1, K
CALL DKSWAP( COV(IJ-K+M), COV(IJ+M) )
END DO
CALL DKSWAP( A(K), A(K+1) )
CALL DKSWAP( B(K), B(K+1) )
M = INFI(K)
INFI(K) = INFI(K+1)
INFI(K+1) = M
IJ = IJ - K
END DO
GO TO 20
END IF
END DO
GO TO 20
END IF
COV(II+J) = 0
END DO
20 II = II + I
Y(I) = 0
END IF
END DO
ENDIF
END
*
SUBROUTINE DKSWAP( X, Y )
DOUBLE PRECISION X, Y, T
T = X
X = Y
Y = T
END
*
SUBROUTINE RCSWP( P, Q, A, B, INFIN, N, C )
*
* Swaps rows and columns P and Q in situ, with P <= Q.
*
DOUBLE PRECISION A(*), B(*), C(*)
INTEGER INFIN(*), P, Q, N, I, J, II, JJ
CALL DKSWAP( A(P), A(Q) )
CALL DKSWAP( B(P), B(Q) )
J = INFIN(P)
INFIN(P) = INFIN(Q)
INFIN(Q) = J
JJ = ( P*( P - 1 ) )/2
II = ( Q*( Q - 1 ) )/2
CALL DKSWAP( C(JJ+P), C(II+Q) )
DO J = 1, P-1
CALL DKSWAP( C(JJ+J), C(II+J) )
END DO
JJ = JJ + P
DO I = P+1, Q-1
CALL DKSWAP( C(JJ+P), C(II+I) )
JJ = JJ + I
END DO
II = II + Q
DO I = Q+1, N
CALL DKSWAP( C(II+P), C(II+Q) )
II = II + I
END DO
END
*
SUBROUTINE DKBVRC( NDIM, MINVLS, MAXVLS, FUNCTN, ABSEPS, RELEPS,
& ABSERR, FINEST, INFORM )
*
* Automatic Multidimensional Integration Subroutine
*
* AUTHOR: Alan Genz
* Department of Mathematics
* Washington State University
* Pulman, WA 99164-3113
* Email: AlanGenz@wsu.edu
*
* Last Change: 1/15/03
*
* KRBVRC computes an approximation to the integral
*
* 1 1 1
* I I ... I F(X) dx(NDIM)...dx(2)dx(1)
* 0 0 0
*
*
* DKBVRC uses randomized Korobov rules for the first 100 variables.
* The primary references are
* "Randomization of Number Theoretic Methods for Multiple Integration"
* R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13, pp. 904-14,
* and
* "Optimal Parameters for Multidimensional Integration",
* P. Keast, SIAM J Numer Anal, 10, pp.831-838.
* If there are more than 100 variables, the remaining variables are
* integrated using the rules described in the reference
* "On a Number-Theoretical Integration Method"
* H. Niederreiter, Aequationes Mathematicae, 8(1972), pp. 304-11.
*
*************** Parameters ********************************************
****** Input parameters
* NDIM Number of variables, must exceed 1, but not exceed 40
* MINVLS Integer minimum number of function evaluations allowed.
* MINVLS must not exceed MAXVLS. If MINVLS < 0 then the
* routine assumes a previous call has been made with
* the same integrand and continues that calculation.
* MAXVLS Integer maximum number of function evaluations allowed.
* FUNCTN EXTERNALly declared user defined function to be integrated.
* It must have parameters (NDIM,Z), where Z is a real array
* of dimension NDIM.
*
* ABSEPS Required absolute accuracy.
* RELEPS Required relative accuracy.
****** Output parameters
* MINVLS Actual number of function evaluations used.
* ABSERR Estimated absolute accuracy of FINEST.
* FINEST Estimated value of integral.
* INFORM INFORM = 0 for normal exit, when
* ABSERR <= MAX(ABSEPS, RELEPS*ABS(FINEST))
* and
* INTVLS <= MAXCLS.
* INFORM = 1 If MAXVLS was too small to obtain the required
* accuracy. In this case a value FINEST is returned with
* estimated absolute accuracy ABSERR.
************************************************************************
EXTERNAL FUNCTN
INTEGER NDIM, MINVLS, MAXVLS, INFORM, NP, PLIM, NLIM, KLIM, KLIMI,
& SAMPLS, I, INTVLS, MINSMP
PARAMETER ( PLIM = 28, NLIM = 1000, KLIM = 100, MINSMP = 8 )
INTEGER P(PLIM), C(PLIM,KLIM-1)
DOUBLE PRECISION FUNCTN, ABSEPS, RELEPS, FINEST, ABSERR, DIFINT,
& FINVAL, VARSQR, VAREST, VARPRD, VALUE
DOUBLE PRECISION X(2*NLIM), VK(NLIM), ONE
PARAMETER ( ONE = 1 )
SAVE P, C, SAMPLS, NP, VAREST
INFORM = 1
INTVLS = 0
KLIMI = KLIM
IF ( MINVLS .GE. 0 ) THEN
FINEST = 0
VAREST = 0
SAMPLS = MINSMP
DO I = MIN( NDIM, 10), PLIM
NP = I
IF ( MINVLS .LT. 2*SAMPLS*P(I) ) GO TO 10
END DO
SAMPLS = MAX( MINSMP, MINVLS/( 2*P(NP) ) )
ENDIF
10 VK(1) = ONE/P(NP)
DO I = 2, NDIM
IF ( I .LE. KLIM ) THEN
VK(I) = MOD( C(NP, MIN(NDIM-1,KLIM-1))*VK(I-1), ONE )
ELSE
VK(I) = INT( P(NP)*2**(DBLE(I-KLIM)/(NDIM-KLIM+1)) )
VK(I) = MOD( VK(I)/P(NP), ONE )
END IF
END DO
FINVAL = 0
VARSQR = 0
DO I = 1, SAMPLS
CALL DKSMRC( NDIM, KLIMI, VALUE, P(NP), VK, FUNCTN, X )
DIFINT = ( VALUE - FINVAL )/I
FINVAL = FINVAL + DIFINT
VARSQR = ( I - 2 )*VARSQR/I + DIFINT**2
END DO
INTVLS = INTVLS + 2*SAMPLS*P(NP)
VARPRD = VAREST*VARSQR
FINEST = FINEST + ( FINVAL - FINEST )/( 1 + VARPRD )
IF ( VARSQR .GT. 0 ) VAREST = ( 1 + VARPRD )/VARSQR
ABSERR = 7*SQRT( VARSQR/( 1 + VARPRD ) )/2
IF ( ABSERR .GT. MAX( ABSEPS, ABS(FINEST)*RELEPS ) ) THEN
IF ( NP .LT. PLIM ) THEN
NP = NP + 1
ELSE
SAMPLS = MIN( 3*SAMPLS/2, ( MAXVLS - INTVLS )/( 2*P(NP) ) )
SAMPLS = MAX( MINSMP, SAMPLS )
ENDIF
IF ( INTVLS + 2*SAMPLS*P(NP) .LE. MAXVLS ) GO TO 10
ELSE
INFORM = 0
ENDIF
MINVLS = INTVLS
*
* Optimal Parameters for Lattice Rules
*
DATA P( 1),(C( 1,I),I = 1,99)/ 31, 12, 2*9, 13, 8*12, 3*3, 12,
& 2*7, 9*12, 3*3, 12, 2*7, 9*12, 3*3, 12, 2*7, 9*12, 3*3, 12, 2*7,
& 8*12, 7, 3*3, 3*7, 21*3/
DATA P( 2),(C( 2,I),I = 1,99)/ 47, 13, 11, 17, 10, 6*15,
& 22, 2*15, 3*6, 2*15, 9, 13, 3*2, 13, 2*11, 10, 9*15, 3*6, 2*15,
& 9, 13, 3*2, 13, 2*11, 10, 9*15, 3*6, 2*15, 9, 13, 3*2, 13, 2*11,
& 2*10, 8*15, 6, 2, 3, 2, 3, 12*2/
DATA P( 3),(C( 3,I),I = 1,99)/ 73, 27, 28, 10, 2*11, 20,
& 2*11, 28, 2*13, 28, 3*13, 16*14, 2*31, 3*5, 31, 13, 6*11, 7*13,
& 16*14, 2*31, 3*5, 11, 13, 7*11, 2*13, 11, 13, 4*5, 14, 13, 8*5/
DATA P( 4),(C( 4,I),I = 1,99)/ 113, 35, 2*27, 36, 22, 2*29,
& 20, 45, 3*5, 16*21, 29, 10*17, 12*23, 21, 27, 3*3, 24, 2*27,
& 17, 3*29, 17, 4*5, 16*21, 3*17, 6, 2*17, 6, 3, 2*6, 5*3/
DATA P( 5),(C( 5,I),I = 1,99)/ 173, 64, 66, 2*28, 2*44, 55,
& 67, 6*10, 2*38, 5*10, 12*49, 2*38, 31, 2*4, 31, 64, 3*4, 64,
& 6*45, 19*66, 11, 9*66, 45, 11, 7, 3, 3*2, 27, 5, 2*3, 2*5, 7*2/
DATA P( 6),(C( 6,I),I = 1,99)/ 263, 111, 42, 54, 118, 20,
& 2*31, 72, 17, 94, 2*14, 11, 3*14, 94, 4*10, 7*14, 3*11, 7*8,
& 5*18, 113, 2*62, 2*45, 17*113, 2*63, 53, 63, 15*67, 5*51, 12,
& 51, 12, 51, 5, 2*3, 2*2, 5/
DATA P( 7),(C( 7,I),I = 1,99)/ 397, 163, 154, 83, 43, 82,
& 92, 150, 59, 2*76, 47, 2*11, 100, 131, 6*116, 9*138, 21*101,
& 6*116, 5*100, 5*138, 19*101, 8*38, 5*3/
DATA P( 8),(C( 8,I),I = 1,99)/ 593, 246, 189, 242, 102,
& 2*250, 102, 250, 280, 118, 196, 118, 191, 215, 2*121,
& 12*49, 34*171, 8*161, 17*14, 6*10, 103, 4*10, 5/
DATA P( 9),(C( 9,I),I = 1,99)/ 907, 347, 402, 322, 418,
& 215, 220, 3*339, 337, 218, 4*315, 4*167, 361, 201, 11*124,
& 2*231, 14*90, 4*48, 23*90, 10*243, 9*283, 16, 283, 16, 2*283/
DATA P(10),(C(10,I),I = 1,99)/ 1361, 505, 220, 601, 644,
& 612, 160, 3*206, 422, 134, 518, 2*134, 518, 652, 382,
& 206, 158, 441, 179, 441, 56, 2*559, 14*56, 2*101, 56,
& 8*101, 7*193, 21*101, 17*122, 4*101/
DATA P(11),(C(11,I),I = 1,99)/ 2053, 794, 325, 960, 528,
& 2*247, 338, 366, 847, 2*753, 236, 2*334, 461, 711, 652,
& 3*381, 652, 7*381, 226, 7*326, 126, 10*326, 2*195, 19*55,
& 7*195, 11*132, 13*387/
DATA P(12),(C(12,I),I = 1,99)/ 3079, 1189, 888, 259, 1082, 725,
& 811, 636, 965, 2*497, 2*1490, 392, 1291, 2*508, 2*1291, 508,
& 1291, 2*508, 4*867, 934, 7*867, 9*1284, 4*563, 3*1010, 208,
& 838, 3*563, 2*759, 564, 2*759, 4*801, 5*759, 8*563, 22*226/
DATA P(13),(C(13,I),I = 1,99)/ 4621, 1763, 1018, 1500, 432,
& 1332, 2203, 126, 2240, 1719, 1284, 878, 1983, 4*266,
& 2*747, 2*127, 2074, 127, 2074, 1400, 10*1383, 1400, 7*1383,
& 507, 4*1073, 5*1990, 9*507, 17*1073, 6*22, 1073, 6*452, 318,
& 4*301, 2*86, 15/
DATA P(14),(C(14,I),I = 1,99)/ 6947, 2872, 3233, 1534, 2941,
& 2910, 393, 1796, 919, 446, 2*919, 1117, 7*103, 2311, 3117, 1101,
& 2*3117, 5*1101, 8*2503, 7*429, 3*1702, 5*184, 34*105, 13*784/
DATA P(15),(C(15,I),I = 1,99)/ 10427, 4309, 3758, 4034, 1963,
& 730, 642, 1502, 2246, 3834, 1511, 2*1102, 2*1522, 2*3427,
& 3928, 2*915, 4*3818, 3*4782, 3818, 4782, 2*3818, 7*1327, 9*1387,
& 13*2339, 18*3148, 3*1776, 3*3354, 925, 2*3354, 5*925, 8*2133/
DATA P(16),(C(16,I),I = 1,99)/ 15641, 6610, 6977, 1686, 3819,
& 2314, 5647, 3953, 3614, 5115, 2*423, 5408, 7426, 2*423,
& 487, 6227, 2660, 6227, 1221, 3811, 197, 4367, 351,
& 1281, 1221, 3*351, 7245, 1984, 6*2999, 3995, 4*2063, 1644,
& 2063, 2077, 3*2512, 4*2077, 19*754, 2*1097, 4*754, 248, 754,
& 4*1097, 4*222, 754,11*1982/
DATA P(17),(C(17,I),I = 1,99)/ 23473, 9861, 3647, 4073, 2535,
& 3430, 9865, 2830, 9328, 4320, 5913, 10365, 8272, 3706, 6186,
& 3*7806, 8610, 2563, 2*11558, 9421, 1181, 9421, 3*1181, 9421,
& 2*1181, 2*10574, 5*3534, 3*2898, 3450, 7*2141, 15*7055, 2831,
& 24*8204, 3*4688, 8*2831/
DATA P(18),(C(18,I),I = 1,99)/ 35221, 10327, 7582, 7124, 8214,
& 9600, 10271, 10193, 10800, 9086, 2365, 4409, 13812,
& 5661, 2*9344, 10362, 2*9344, 8585, 11114, 3*13080, 6949,
& 3*3436, 13213, 2*6130, 2*8159, 11595, 8159, 3436, 18*7096,
& 4377, 7096, 5*4377, 2*5410, 32*4377, 2*440, 3*1199/
DATA P(19),(C(19,I),I = 1,99)/ 52837, 19540, 19926, 11582,
& 11113, 24585, 8726, 17218, 419, 3*4918, 15701, 17710,
& 2*4037, 15808, 11401, 19398, 2*25950, 4454, 24987, 11719,
& 8697, 5*1452, 2*8697, 6436, 21475, 6436, 22913, 6434, 18497,
& 4*11089, 2*3036, 4*14208, 8*12906, 4*7614, 6*5021, 24*10145,
& 6*4544, 4*8394/
DATA P(20),(C(20,I),I = 1,99)/ 79259, 34566, 9579, 12654,
& 26856, 37873, 38806, 29501, 17271, 3663, 10763, 18955,
& 1298, 26560, 2*17132, 2*4753, 8713, 18624, 13082, 6791,
& 1122, 19363, 34695, 4*18770, 15628, 4*18770, 33766, 6*20837,
& 5*6545, 14*12138, 5*30483, 19*12138, 9305, 13*11107, 2*9305/
DATA P(21),(C(21,I),I = 1,99)/118891, 31929, 49367, 10982, 3527,
& 27066, 13226, 56010, 18911, 40574, 2*20767, 9686, 2*47603,
& 2*11736, 41601, 12888, 32948, 30801, 44243, 2*53351, 16016,
& 2*35086, 32581, 2*2464, 49554, 2*2464, 2*49554, 2464, 81, 27260,
& 10681, 7*2185, 5*18086, 2*17631, 3*18086, 37335, 3*37774,
& 13*26401, 12982, 6*40398, 3*3518, 9*37799, 4*4721, 4*7067/
DATA P(22),(C(22,I),I = 1,99)/178349, 40701, 69087, 77576, 64590,
& 39397, 33179, 10858, 38935, 43129, 2*35468, 5279, 2*61518, 27945,
& 2*70975, 2*86478, 2*20514, 2*73178, 2*43098, 4701,
& 2*59979, 58556, 69916, 2*15170, 2*4832, 43064, 71685, 4832,
& 3*15170, 3*27679, 2*60826, 2*6187, 5*4264, 45567, 4*32269,
& 9*62060, 13*1803, 12*51108, 2*55315, 5*54140, 13134/
DATA P(23),(C(23,I),I = 1,99)/267523, 103650, 125480, 59978,
& 46875, 77172, 83021, 126904, 14541, 56299, 43636, 11655,
& 52680, 88549, 29804, 101894, 113675, 48040, 113675,
& 34987, 48308, 97926, 5475, 49449, 6850, 2*62545, 9440,
& 33242, 9440, 33242, 9440, 33242, 9440, 62850, 3*9440,
& 3*90308, 9*47904, 7*41143, 5*36114, 24997, 14*65162, 7*47650,
& 7*40586, 4*38725, 5*88329/
DATA P(24),(C(24,I),I = 1,99)/401287, 165843, 90647, 59925,
& 189541, 67647, 74795, 68365, 167485, 143918, 74912,
& 167289, 75517, 8148, 172106, 126159,3*35867, 121694,
& 52171, 95354, 2*113969, 76304, 2*123709, 144615, 123709,
& 2*64958, 32377, 2*193002, 25023, 40017, 141605, 2*189165,
& 141605, 2*189165, 3*141605, 189165, 20*127047, 10*127785,
& 6*80822, 16*131661, 7114, 131661/
DATA P(25),(C(25,I),I = 1,99)/601942, 130365, 236711, 110235,
& 125699, 56483, 93735, 234469, 60549, 1291, 93937,
& 245291, 196061, 258647, 162489, 176631, 204895, 73353,
& 172319, 28881, 136787,2*122081, 275993, 64673, 3*211587,
& 2*282859, 211587, 242821, 3*256865, 122203, 291915, 122203,
& 2*291915, 122203, 2*25639, 291803, 245397, 284047,
& 7*245397, 94241, 2*66575, 19*217673, 10*210249, 15*94453/
DATA P(26),(C(26,I),I = 1,99)/902933, 333459, 375354, 102417,
& 383544, 292630, 41147, 374614, 48032, 435453, 281493, 358168,
& 114121, 346892, 238990, 317313, 164158, 35497, 2*70530, 434839,
& 3*24754, 393656, 2*118711, 148227, 271087, 355831, 91034,
& 2*417029, 2*91034, 417029, 91034, 2*299843, 2*413548, 308300,
& 3*413548, 3*308300, 413548, 5*308300, 4*15311, 2*176255, 6*23613,
& 172210, 4* 204328, 5*121626, 5*200187, 2*121551, 12*248492,
& 5*13942/
DATA P(27), (C(27,I), I = 1,99)/ 1354471, 500884, 566009, 399251,
& 652979, 355008, 430235, 328722, 670680, 2*405585, 424646,
& 2*670180, 641587, 215580, 59048, 633320, 81010, 20789, 2*389250,
& 2*638764, 2*389250, 398094, 80846, 2*147776, 296177, 2*398094,
& 2*147776, 396313, 3*578233, 19482, 620706, 187095, 620706,
& 187095, 126467, 12*241663, 321632, 2*23210, 3*394484, 3*78101,
& 19*542095, 3*277743, 12*457259/
DATA P(28), (C(28,I), I = 1, 99)/ 2031713, 858339, 918142, 501970,
& 234813, 460565, 31996, 753018, 256150, 199809, 993599, 245149,
& 794183, 121349, 150619, 376952, 2*809123, 804319, 67352, 969594,
& 434796, 969594, 804319, 391368, 761041, 754049, 466264, 2*754049,
& 466264, 2*754049, 282852, 429907, 390017, 276645, 994856, 250142,
& 144595, 907454, 689648, 4*687580, 978368, 687580, 552742, 105195,
& 942843, 768249, 4*307142, 7*880619, 11*117185, 11*60731,
& 4*178309, 8*74373, 3*214965/
*
END
*
SUBROUTINE DKSMRC( NDIM, KLIM, SUMKRO, PRIME, VK, FUNCTN, X )
EXTERNAL FUNCTN
INTEGER NDIM, NK, KLIM, PRIME, K, J, JP
DOUBLE PRECISION SUMKRO, VK(*), FUNCTN, X(*), ONE, XT, MVNUNI
PARAMETER ( ONE = 1 )
SUMKRO = 0
NK = MIN( NDIM, KLIM )
DO J = 1, NK - 1
JP = J + MVNUNI()*( NK + 1 - J )
XT = VK(J)
VK(J) = VK(JP)
VK(JP) = XT
END DO
DO J = 1, NDIM
X(NDIM+J) = MVNUNI()
END DO
DO K = 1, PRIME
DO J = 1, NDIM
X(J) = ABS( 2*MOD( K*VK(J) + X(NDIM+J), ONE ) - 1 )
END DO
SUMKRO = SUMKRO + ( FUNCTN(NDIM,X) - SUMKRO )/( 2*K - 1 )
DO J = 1, NDIM
X(J) = 1 - X(J)
END DO
SUMKRO = SUMKRO + ( FUNCTN(NDIM,X) - SUMKRO )/( 2*K )
END DO
END
*
DOUBLE PRECISION FUNCTION MVNPHI( Z )
*
* Normal distribution probabilities accurate to 1.e-15.
* Z = no. of standard deviations from the mean.
*
* Based upon algorithm 5666 for the error function, from:
* Hart, J.F. et al, 'Computer Approximations', Wiley 1968
*
* Programmer: Alan Miller
*
* Latest revision - 30 March 1986
*
DOUBLE PRECISION P0, P1, P2, P3, P4, P5, P6,
* Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,
* Z, P, EXPNTL, CUTOFF, ROOTPI, ZABS
PARAMETER(
* P0 = 220.20 68679 12376 1D0,
* P1 = 221.21 35961 69931 1D0,
* P2 = 112.07 92914 97870 9D0,
* P3 = 33.912 86607 83830 0D0,
* P4 = 6.3739 62203 53165 0D0,
* P5 = .70038 30644 43688 1D0,
* P6 = .035262 49659 98910 9D0 )
PARAMETER(
* Q0 = 440.41 37358 24752 2D0,
* Q1 = 793.82 65125 19948 4D0,
* Q2 = 637.33 36333 78831 1D0,
* Q3 = 296.56 42487 79673 7D0,
* Q4 = 86.780 73220 29460 8D0,
* Q5 = 16.064 17757 92069 5D0,
* Q6 = 1.7556 67163 18264 2D0,
* Q7 = .088388 34764 83184 4D0 )
PARAMETER( ROOTPI = 2.5066 28274 63100 1D0 )
PARAMETER( CUTOFF = 7.0710 67811 86547 5D0 )
*
ZABS = ABS(Z)
*
* |Z| > 37
*
IF ( ZABS .GT. 37 ) THEN
P = 0
ELSE
*
* |Z| <= 37
*
EXPNTL = EXP( -ZABS**2/2 )
*
* |Z| < CUTOFF = 10/SQRT(2)
*
IF ( ZABS .LT. CUTOFF ) THEN
P = EXPNTL*( (((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
* + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
* + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
* + Q0 )
*
* |Z| >= CUTOFF.
*
ELSE
P = EXPNTL/( ZABS + 1/( ZABS + 2/( ZABS + 3/( ZABS
* + 4/( ZABS + 0.65D0 ) ) ) ) )/ROOTPI
END IF
END IF
IF ( Z .GT. 0 ) P = 1 - P
MVNPHI = P
END
DOUBLE PRECISION FUNCTION PHINVS(P)
*
* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
*
* Produces the normal deviate Z corresponding to a given lower
* tail area of P.
*
* The hash sums below are the sums of the mantissas of the
* coefficients. They are included for use in checking
* transcription.
*
DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2,
* A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
* C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
* E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
* P, Q, R
PARAMETER ( SPLIT1 = 0.425, SPLIT2 = 5,
* CONST1 = 0.180625D0, CONST2 = 1.6D0 )
*
* Coefficients for P close to 0.5
*
PARAMETER (
* A0 = 3.38713 28727 96366 6080D0,
* A1 = 1.33141 66789 17843 7745D+2,
* A2 = 1.97159 09503 06551 4427D+3,
* A3 = 1.37316 93765 50946 1125D+4,
* A4 = 4.59219 53931 54987 1457D+4,
* A5 = 6.72657 70927 00870 0853D+4,
* A6 = 3.34305 75583 58812 8105D+4,
* A7 = 2.50908 09287 30122 6727D+3,
* B1 = 4.23133 30701 60091 1252D+1,
* B2 = 6.87187 00749 20579 0830D+2,
* B3 = 5.39419 60214 24751 1077D+3,
* B4 = 2.12137 94301 58659 5867D+4,
* B5 = 3.93078 95800 09271 0610D+4,
* B6 = 2.87290 85735 72194 2674D+4,
* B7 = 5.22649 52788 52854 5610D+3 )
* HASH SUM AB 55.88319 28806 14901 4439
*
* Coefficients for P not close to 0, 0.5 or 1.
*
PARAMETER (
* C0 = 1.42343 71107 49683 57734D0,
* C1 = 4.63033 78461 56545 29590D0,
* C2 = 5.76949 72214 60691 40550D0,
* C3 = 3.64784 83247 63204 60504D0,
* C4 = 1.27045 82524 52368 38258D0,
* C5 = 2.41780 72517 74506 11770D-1,
* C6 = 2.27238 44989 26918 45833D-2,
* C7 = 7.74545 01427 83414 07640D-4,
* D1 = 2.05319 16266 37758 82187D0,
* D2 = 1.67638 48301 83803 84940D0,
* D3 = 6.89767 33498 51000 04550D-1,
* D4 = 1.48103 97642 74800 74590D-1,
* D5 = 1.51986 66563 61645 71966D-2,
* D6 = 5.47593 80849 95344 94600D-4,
* D7 = 1.05075 00716 44416 84324D-9 )
* HASH SUM CD 49.33206 50330 16102 89036
*
* Coefficients for P near 0 or 1.
*
PARAMETER (
* E0 = 6.65790 46435 01103 77720D0,
* E1 = 5.46378 49111 64114 36990D0,
* E2 = 1.78482 65399 17291 33580D0,
* E3 = 2.96560 57182 85048 91230D-1,
* E4 = 2.65321 89526 57612 30930D-2,
* E5 = 1.24266 09473 88078 43860D-3,
* E6 = 2.71155 55687 43487 57815D-5,
* E7 = 2.01033 43992 92288 13265D-7,
* F1 = 5.99832 20655 58879 37690D-1,
* F2 = 1.36929 88092 27358 05310D-1,
* F3 = 1.48753 61290 85061 48525D-2,
* F4 = 7.86869 13114 56132 59100D-4,
* F5 = 1.84631 83175 10054 68180D-5,
* F6 = 1.42151 17583 16445 88870D-7,
* F7 = 2.04426 31033 89939 78564D-15 )
* HASH SUM EF 47.52583 31754 92896 71629
*
Q = ( 2*P - 1 )/2
IF ( ABS(Q) .LE. SPLIT1 ) THEN
R = CONST1 - Q*Q
PHINVS = Q*( ( ( ((((A7*R + A6)*R + A5)*R + A4)*R + A3)
* *R + A2 )*R + A1 )*R + A0 )
* /( ( ( ((((B7*R + B6)*R + B5)*R + B4)*R + B3)
* *R + B2 )*R + B1 )*R + 1 )
ELSE
R = MIN( P, 1 - P )
IF ( R .GT. 0 ) THEN
R = SQRT( -LOG(R) )
IF ( R .LE. SPLIT2 ) THEN
R = R - CONST2
PHINVS = ( ( ( ((((C7*R + C6)*R + C5)*R + C4)*R + C3)
* *R + C2 )*R + C1 )*R + C0 )
* /( ( ( ((((D7*R + D6)*R + D5)*R + D4)*R + D3)
* *R + D2 )*R + D1 )*R + 1 )
ELSE
R = R - SPLIT2
PHINVS = ( ( ( ((((E7*R + E6)*R + E5)*R + E4)*R + E3)
* *R + E2 )*R + E1 )*R + E0 )
* /( ( ( ((((F7*R + F6)*R + F5)*R + F4)*R + F3)
* *R + F2 )*R + F1 )*R + 1 )
END IF
ELSE
PHINVS = 9
END IF
IF ( Q .LT. 0 ) PHINVS = - PHINVS
END IF
END
DOUBLE PRECISION FUNCTION BVNMVN( LOWER, UPPER, INFIN, CORREL )
*
* A function for computing bivariate normal probabilities.
*
* Parameters
*
* LOWER REAL, array of lower integration limits.
* UPPER REAL, array of upper integration limits.
* INFIN INTEGER, array of integration limits flags:
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
* CORREL REAL, correlation coefficient.
*
DOUBLE PRECISION LOWER(*), UPPER(*), CORREL, BVU
INTEGER INFIN(*)
IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 2 ) THEN
BVNMVN = BVU ( LOWER(1), LOWER(2), CORREL )
+ - BVU ( UPPER(1), LOWER(2), CORREL )
+ - BVU ( LOWER(1), UPPER(2), CORREL )
+ + BVU ( UPPER(1), UPPER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 1 ) THEN
BVNMVN = BVU ( LOWER(1), LOWER(2), CORREL )
+ - BVU ( UPPER(1), LOWER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 2 ) THEN
BVNMVN = BVU ( LOWER(1), LOWER(2), CORREL )
+ - BVU ( LOWER(1), UPPER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 0 ) THEN
BVNMVN = BVU ( -UPPER(1), -UPPER(2), CORREL )
+ - BVU ( -LOWER(1), -UPPER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 2 ) THEN
BVNMVN = BVU ( -UPPER(1), -UPPER(2), CORREL )
+ - BVU ( -UPPER(1), -LOWER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 0 ) THEN
BVNMVN = BVU ( LOWER(1), -UPPER(2), -CORREL )
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 1 ) THEN
BVNMVN = BVU ( -UPPER(1), LOWER(2), -CORREL )
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 1 ) THEN
BVNMVN = BVU ( LOWER(1), LOWER(2), CORREL )
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 0 ) THEN
BVNMVN = BVU ( -UPPER(1), -UPPER(2), CORREL )
END IF
END
DOUBLE PRECISION FUNCTION BVU( SH, SK, R )
*
* A function for computing bivariate normal probabilities.
*
* Yihong Ge
* Department of Computer Science and Electrical Engineering
* Washington State University
* Pullman, WA 99164-2752
* and
* Alan Genz
* Department of Mathematics
* Washington State University
* Pullman, WA 99164-3113
* Email : alangenz@wsu.edu
*
* BVN - calculate the probability that X is larger than SH and Y is
* larger than SK.
*
* Parameters
*
* SH REAL, integration limit
* SK REAL, integration limit
* R REAL, correlation coefficient
* LG INTEGER, number of Gauss Rule Points and Weights
*
DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
INTEGER I, LG, NG
PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 )
DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
DOUBLE PRECISION MVNPHI, SN, ASR, H, K, BS, HS, HK
SAVE X, W
* Gauss Legendre Points and Weights, N = 6
DATA ( W(I,1), X(I,1), I = 1,3) /
* 0.1713244923791705D+00,-0.9324695142031522D+00,
* 0.3607615730481384D+00,-0.6612093864662647D+00,
* 0.4679139345726904D+00,-0.2386191860831970D+00/
* Gauss Legendre Points and Weights, N = 12
DATA ( W(I,2), X(I,2), I = 1,6) /
* 0.4717533638651177D-01,-0.9815606342467191D+00,
* 0.1069393259953183D+00,-0.9041172563704750D+00,
* 0.1600783285433464D+00,-0.7699026741943050D+00,
* 0.2031674267230659D+00,-0.5873179542866171D+00,
* 0.2334925365383547D+00,-0.3678314989981802D+00,
* 0.2491470458134029D+00,-0.1252334085114692D+00/
* Gauss Legendre Points and Weights, N = 20
DATA ( W(I,3), X(I,3), I = 1,10) /
* 0.1761400713915212D-01,-0.9931285991850949D+00,
* 0.4060142980038694D-01,-0.9639719272779138D+00,
* 0.6267204833410906D-01,-0.9122344282513259D+00,
* 0.8327674157670475D-01,-0.8391169718222188D+00,
* 0.1019301198172404D+00,-0.7463319064601508D+00,
* 0.1181945319615184D+00,-0.6360536807265150D+00,
* 0.1316886384491766D+00,-0.5108670019508271D+00,
* 0.1420961093183821D+00,-0.3737060887154196D+00,
* 0.1491729864726037D+00,-0.2277858511416451D+00,
* 0.1527533871307259D+00,-0.7652652113349733D-01/
IF ( ABS(R) .LT. 0.3 ) THEN
NG = 1
LG = 3
ELSE IF ( ABS(R) .LT. 0.75 ) THEN
NG = 2
LG = 6
ELSE
NG = 3
LG = 10
ENDIF
H = SH
K = SK
HK = H*K
BVN = 0
IF ( ABS(R) .LT. 0.925 ) THEN
HS = ( H*H + K*K )/2
ASR = ASIN(R)
DO I = 1, LG
SN = SIN(ASR*( X(I,NG)+1 )/2)
BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
SN = SIN(ASR*(-X(I,NG)+1 )/2)
BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
END DO
BVN = BVN*ASR/(2*TWOPI) + MVNPHI(-H)*MVNPHI(-K)
ELSE
IF ( R .LT. 0 ) THEN
K = -K
HK = -HK
ENDIF
IF ( ABS(R) .LT. 1 ) THEN
AS = ( 1 - R )*( 1 + R )
A = SQRT(AS)
BS = ( H - K )**2
C = ( 4 - HK )/8
D = ( 12 - HK )/16
BVN = A*EXP( -(BS/AS + HK)/2 )
+ *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
IF ( HK .GT. -160 ) THEN
B = SQRT(BS)
BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVNPHI(-B/A)*B
+ *( 1 - C*BS*( 1 - D*BS/5 )/3 )
ENDIF
A = A/2
DO I = 1, LG
XS = ( A*(X(I,NG)+1) )**2
RS = SQRT( 1 - XS )
BVN = BVN + A*W(I,NG)*
+ ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
+ - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
XS = AS*(-X(I,NG)+1)**2/4
RS = SQRT( 1 - XS )
BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
+ *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
+ - ( 1 + C*XS*( 1 + D*XS ) ) )
END DO
BVN = -BVN/TWOPI
ENDIF
IF ( R .GT. 0 ) BVN = BVN + MVNPHI( -MAX( H, K ) )
IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVNPHI(-H)-MVNPHI(-K) )
ENDIF
BVU = BVN
END
DOUBLE PRECISION FUNCTION MVNUNI()
*
* Uniform (0,1) random number generator
*
* Reference:
* L'Ecuyer, Pierre (1996),
* "Combined Multiple Recursive Random Number Generators"
* Operations Research 44, pp. 816-822.
*
*
INTEGER A12, A13, A21, A23, P12, P13, P21, P23
INTEGER Q12, Q13, Q21, Q23, R12, R13, R21, R23
INTEGER X10, X11, X12, X20, X21, X22, Z, M1, M2, H
DOUBLE PRECISION INVMP1
PARAMETER ( M1 = 2147483647, M2 = 2145483479 )
PARAMETER ( A12 = 63308, Q12 = 33921, R12 = 12979 )
PARAMETER ( A13 = -183326, Q13 = 11714, R13 = 2883 )
PARAMETER ( A21 = 86098, Q21 = 24919, R21 = 7417 )
PARAMETER ( A23 = -539608, Q23 = 3976, R23 = 2071 )
PARAMETER ( INVMP1 = 4.656612873077392578125D-10 )
* INVMP1 = 1/(M1+1)
SAVE X10, X11, X12, X20, X21, X22
DATA X10, X11, X12, X20, X21, X22
& / 15485857, 17329489, 36312197, 55911127, 75906931, 96210113 /
*
* Component 1
*
H = X10/Q13
P13 = -A13*( X10 - H*Q13 ) - H*R13
H = X11/Q12
P12 = A12*( X11 - H*Q12 ) - H*R12
IF ( P13 .LT. 0 ) P13 = P13 + M1
IF ( P12 .LT. 0 ) P12 = P12 + M1
X10 = X11
X11 = X12
X12 = P12 - P13
IF ( X12 .LT. 0 ) X12 = X12 + M1
*
* Component 2
*
H = X20/Q23
P23 = -A23*( X20 - H*Q23 ) - H*R23
H = X22/Q21
P21 = A21*( X22 - H*Q21 ) - H*R21
IF ( P23 .LT. 0 ) P23 = P23 + M2
IF ( P21 .LT. 0 ) P21 = P21 + M2
X20 = X21
X21 = X22
X22 = P21 - P23
IF ( X22 .LT. 0 ) X22 = X22 + M2
*
* Combination
*
Z = X12 - X22
IF ( Z .LE. 0 ) Z = Z + M1
MVNUNI = Z*INVMP1
END
|
