« back to all changes in this revision
Viewing changes to src/sadmvnt.f
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Dirk Eddelbuettel
- Date: 2007-10-03 21:25:21 UTC
- Revision ID: james.westby@ubuntu.com-20071003212521-e6zr9fhv1oowj396
Tags: upstream-1.2-1
Import upstream version 1.2-1
- files added:
- R
- man
- src
added
removed
src/sadmvnt.f
1
* Selected portion of code taken from:
2
* http://www.math.wsu.edu/math/faculty/genz/software/mvn.f
3
* http://www.math.wsu.edu/math/faculty/genz/software/mvt.f
4
*
5
* Author:
6
* Alan Genz
7
* Department of Mathematics
8
* Washington State University
9
* Pullman, WA 99164-3113
10
* Email : alangenz@wsu.edu
11
* except for some auxiliary functions whose authors are indicated
12
* in the respective code below.
13
*-----------------------------------------------------------------------
14
15
SUBROUTINE SADMVN( N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
16
& ABSEPS, RELEPS, ERROR, VALUE, INFORM )
17
*
18
* A subroutine for computing multivariate normal probabilities.
19
* This subroutine uses an algorithm given in the paper
20
* "Numerical Computation of Multivariate Normal Probabilities", in
21
* J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
22
* Alan Genz
23
* Department of Mathematics
24
* Washington State University
25
* Pullman, WA 99164-3113
26
* Email : alangenz@wsu.edu
27
*
28
* Parameters
29
*
30
* N INTEGER, the number of variables.
31
* LOWER REAL, array of lower integration limits.
32
* UPPER REAL, array of upper integration limits.
33
* INFIN INTEGER, array of integration limits flags:
34
* if INFIN(I) < 0, Ith limits are (-infinity, infinity);
35
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
36
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
37
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
38
* CORREL REAL, array of correlation coefficients; the correlation
39
* coefficient in row I column J of the correlation matrix
40
* should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
41
* MAXPTS INTEGER, maximum number of function values allowed. This
42
* parameter can be used to limit the time taken. A
43
* sensible strategy is to start with MAXPTS = 1000*N, and then
44
* increase MAXPTS if ERROR is too large.
45
* ABSEPS REAL absolute error tolerance.
46
* RELEPS REAL relative error tolerance.
47
* ERROR REAL estimated absolute error, with 99% confidence level.
48
* VALUE REAL estimated value for the integral
49
* INFORM INTEGER, termination status parameter:
50
* if INFORM = 0, normal completion with ERROR < EPS;
51
* if INFORM = 1, completion with ERROR > EPS and MAXPTS
52
* function vaules used; increase MAXPTS to
53
* decrease ERROR;
54
* if INFORM = 2, N > 20 or N < 1.
55
*
56
EXTERNAL MVNFNC
57
INTEGER N, NL, M, INFIN(*), LENWRK, MAXPTS, INFORM, INFIS,
58
& RULCLS, TOTCLS, NEWCLS, MAXCLS
59
DOUBLE PRECISION
60
& CORREL(*), LOWER(*), UPPER(*), ABSEPS, RELEPS, ERROR, VALUE,
61
& OLDVAL, D, E, MVNNIT, MVNFNC
62
PARAMETER ( NL = 20 )
63
PARAMETER ( LENWRK = 20*NL**2 )
64
DOUBLE PRECISION WORK(LENWRK)
65
IF ( N .GT. 20 .OR. N .LT. 1 ) THEN
66
INFORM = 2
67
VALUE = 0
68
ERROR = 1
69
RETURN
70
ENDIF
71
INFORM = MVNNIT( N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E )
72
M = N - INFIS
73
IF ( M .EQ. 0 ) THEN
74
VALUE = 1
75
ERROR = 0
76
ELSE IF ( M .EQ. 1 ) THEN
77
VALUE = E - D
78
ERROR = 2E-16
79
ELSE
80
*
81
* Call the subregion adaptive integration subroutine
82
*
83
M = M - 1
84
RULCLS = 1
85
CALL ADAPT( M, RULCLS, 0, MVNFNC, ABSEPS, RELEPS,
86
& LENWRK, WORK, ERROR, VALUE, INFORM )
87
MAXCLS = MIN( 10*RULCLS, MAXPTS )
88
TOTCLS = 0
89
CALL ADAPT(M, TOTCLS, MAXCLS, MVNFNC, ABSEPS, RELEPS,
90
& LENWRK, WORK, ERROR, VALUE, INFORM)
91
IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
92
10 OLDVAL = VALUE
93
MAXCLS = MAX( 2*RULCLS, MIN( 3*MAXCLS/2, MAXPTS - TOTCLS ) )
94
NEWCLS = -1
95
CALL ADAPT(M, NEWCLS, MAXCLS, MVNFNC, ABSEPS, RELEPS,
96
& LENWRK, WORK, ERROR, VALUE, INFORM)
97
TOTCLS = TOTCLS + NEWCLS
98
ERROR = ABS(VALUE-OLDVAL) + SQRT(RULCLS*ERROR**2/TOTCLS)
99
IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
100
IF ( MAXPTS - TOTCLS .GT. 2*RULCLS ) GO TO 10
101
ELSE
102
INFORM = 0
103
END IF
104
ENDIF
105
ENDIF
106
END
107
*--------------------------------------------------------------------------
108
*
109
SUBROUTINE ADAPT(NDIM, MINCLS, MAXCLS, FUNCTN,
110
& ABSREQ, RELREQ, LENWRK, WORK, ABSEST, FINEST, INFORM)
111
*
112
* Adaptive Multidimensional Integration Subroutine
113
*
114
* Author: Alan Genz
115
* Department of Mathematics
116
* Washington State University
117
* Pullman, WA 99164-3113 USA
118
*
119
* This subroutine computes an approximation to the integral
120
*
121
* 1 1 1
122
* I I ... I FUNCTN(NDIM,X) dx(NDIM)...dx(2)dx(1)
123
* 0 0 0
124
*
125
*************** Parameters for ADAPT ********************************
126
*
127
****** Input Parameters
128
*
129
* NDIM Integer number of integration variables.
130
* MINCLS Integer minimum number of FUNCTN calls to be allowed; MINCLS
131
* must not exceed MAXCLS. If MINCLS < 0, then ADAPT assumes
132
* that a previous call of ADAPT has been made with the same
133
* integrand and continues that calculation.
134
* MAXCLS Integer maximum number of FUNCTN calls to be used; MAXCLS
135
* must be >= RULCLS, the number of function calls required for
136
* one application of the basic integration rule.
137
* IF ( NDIM .EQ. 1 ) THEN
138
* RULCLS = 11
139
* ELSE IF ( NDIM .LT. 15 ) THEN
140
* RULCLS = 2**NDIM + 2*NDIM*(NDIM+3) + 1
141
* ELSE
142
* RULCLS = 1 + NDIM*(24-NDIM*(6-NDIM*4))/3
143
* ENDIF
144
* FUNCTN Externally declared real user defined integrand. Its
145
* parameters must be (NDIM, Z), where Z is a real array of
146
* length NDIM.
147
* ABSREQ Real required absolute accuracy.
148
* RELREQ Real required relative accuracy.
149
* LENWRK Integer length of real array WORK (working storage); ADAPT
150
* needs LENWRK >= 16*NDIM + 27. For maximum efficiency LENWRK
151
* should be about 2*NDIM*MAXCLS/RULCLS if MAXCLS FUNCTN
152
* calls are needed. If LENWRK is significantly less than this,
153
* ADAPT may be less efficient.
154
*
155
****** Output Parameters
156
*
157
* MINCLS Actual number of FUNCTN calls used by ADAPT.
158
* WORK Real array (length LENWRK) of working storage. This contains
159
* information that is needed for additional calls of ADAPT
160
* using the same integrand (input MINCLS < 0).
161
* ABSEST Real estimated absolute accuracy.
162
* FINEST Real estimated value of integral.
163
* INFORM INFORM = 0 for normal exit, when ABSEST <= ABSREQ or
164
* ABSEST <= |FINEST|*RELREQ with MINCLS <= MAXCLS.
165
* INFORM = 1 if MAXCLS was too small for ADAPT to obtain the
166
* result FINEST to within the requested accuracy.
167
* INFORM = 2 if MINCLS > MAXCLS, LENWRK < 16*NDIM + 27 or
168
* RULCLS > MAXCLS.
169
*
170
************************************************************************
171
*
172
* Begin driver routine. This routine partitions the working storage
173
* array and then calls the main subroutine ADBASE.
174
*
175
EXTERNAL FUNCTN
176
INTEGER NDIM, MINCLS, MAXCLS, LENWRK, INFORM
177
DOUBLE PRECISION
178
& FUNCTN, ABSREQ, RELREQ, WORK(LENWRK), ABSEST, FINEST
179
INTEGER SBRGNS, MXRGNS, RULCLS, LENRUL,
180
& INERRS, INVALS, INPTRS, INLWRS, INUPRS, INMSHS, INPNTS, INWGTS,
181
& INLOWR, INUPPR, INWDTH, INMESH, INWORK
182
IF ( NDIM .EQ. 1 ) THEN
183
LENRUL = 5
184
RULCLS = 9
185
ELSE IF ( NDIM .LT. 12 ) THEN
186
LENRUL = 6
187
RULCLS = 2**NDIM + 2*NDIM*(NDIM+2) + 1
188
ELSE
189
LENRUL = 6
190
RULCLS = 1 + 2*NDIM*(1+2*NDIM)
191
ENDIF
192
IF ( LENWRK .GE. LENRUL*(NDIM+4) + 10*NDIM + 3 .AND.
193
& RULCLS. LE. MAXCLS .AND. MINCLS .LE. MAXCLS ) THEN
194
MXRGNS = ( LENWRK - LENRUL*(NDIM+4) - 7*NDIM )/( 3*NDIM + 3 )
195
INERRS = 1
196
INVALS = INERRS + MXRGNS
197
INPTRS = INVALS + MXRGNS
198
INLWRS = INPTRS + MXRGNS
199
INUPRS = INLWRS + MXRGNS*NDIM
200
INMSHS = INUPRS + MXRGNS*NDIM
201
INWGTS = INMSHS + MXRGNS*NDIM
202
INPNTS = INWGTS + LENRUL*4
203
INLOWR = INPNTS + LENRUL*NDIM
204
INUPPR = INLOWR + NDIM
205
INWDTH = INUPPR + NDIM
206
INMESH = INWDTH + NDIM
207
INWORK = INMESH + NDIM
208
IF ( MINCLS .LT. 0 ) SBRGNS = WORK(LENWRK)
209
CALL ADBASE(NDIM, MINCLS, MAXCLS, FUNCTN, ABSREQ, RELREQ,
210
& ABSEST, FINEST, SBRGNS, MXRGNS, RULCLS, LENRUL,
211
& WORK(INERRS), WORK(INVALS), WORK(INPTRS), WORK(INLWRS),
212
& WORK(INUPRS), WORK(INMSHS), WORK(INWGTS), WORK(INPNTS),
213
& WORK(INLOWR), WORK(INUPPR), WORK(INWDTH), WORK(INMESH),
214
& WORK(INWORK), INFORM)
215
WORK(LENWRK) = SBRGNS
216
ELSE
217
INFORM = 2
218
MINCLS = RULCLS
219
ENDIF
220
END
221
*
222
SUBROUTINE ADBASE(NDIM, MINCLS, MAXCLS, FUNCTN, ABSREQ, RELREQ,
223
& ABSEST, FINEST, SBRGNS, MXRGNS, RULCLS, LENRUL,
224
& ERRORS, VALUES, PONTRS, LOWERS,
225
& UPPERS, MESHES, WEGHTS, POINTS,
226
& LOWER, UPPER, WIDTH, MESH, WORK, INFORM)
227
*
228
* Main adaptive integration subroutine
229
*
230
EXTERNAL FUNCTN
231
INTEGER I, J, NDIM, MINCLS, MAXCLS, SBRGNS, MXRGNS,
232
& RULCLS, LENRUL, INFORM, NWRGNS
233
DOUBLE PRECISION FUNCTN, ABSREQ, RELREQ, ABSEST, FINEST,
234
& ERRORS(*), VALUES(*), PONTRS(*),
235
& LOWERS(NDIM,*), UPPERS(NDIM,*),
236
& MESHES(NDIM,*),WEGHTS(*), POINTS(*),
237
& LOWER(*), UPPER(*), WIDTH(*), MESH(*), WORK(*)
238
INTEGER DIVAXN, TOP, RGNCLS, FUNCLS, DIFCLS
239
240
*
241
* Initialization of subroutine
242
*
243
INFORM = 2
244
FUNCLS = 0
245
CALL BSINIT(NDIM, WEGHTS, LENRUL, POINTS)
246
IF ( MINCLS .GE. 0) THEN
247
*
248
* When MINCLS >= 0 determine initial subdivision of the
249
* integration region and apply basic rule to each subregion.
250
*
251
SBRGNS = 0
252
DO I = 1,NDIM
253
LOWER(I) = 0
254
MESH(I) = 1
255
WIDTH(I) = 1/(2*MESH(I))
256
UPPER(I) = 1
257
END DO
258
DIVAXN = 0
259
RGNCLS = RULCLS
260
NWRGNS = 1
261
10 CALL DIFFER(NDIM, LOWER, UPPER, WIDTH, WORK, WORK(NDIM+1),
262
& FUNCTN, DIVAXN, DIFCLS)
263
FUNCLS = FUNCLS + DIFCLS
264
IF ( FUNCLS + RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
265
& .LE. MINCLS ) THEN
266
RGNCLS = RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
267
NWRGNS = NWRGNS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
268
MESH(DIVAXN) = MESH(DIVAXN) + 1
269
WIDTH(DIVAXN) = 1/( 2*MESH(DIVAXN) )
270
GO TO 10
271
ENDIF
272
IF ( NWRGNS .LE. MXRGNS ) THEN
273
DO I = 1,NDIM
274
UPPER(I) = LOWER(I) + 2*WIDTH(I)
275
MESH(I) = 1
276
END DO
277
ENDIF
278
*
279
* Apply basic rule to subregions and store results in heap.
280
*
281
20 SBRGNS = SBRGNS + 1
282
CALL BASRUL(NDIM, LOWER, UPPER, WIDTH, FUNCTN,
283
& WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
284
& ERRORS(SBRGNS),VALUES(SBRGNS))
285
CALL TRESTR(SBRGNS, SBRGNS, PONTRS, ERRORS)
286
DO I = 1,NDIM
287
LOWERS(I,SBRGNS) = LOWER(I)
288
UPPERS(I,SBRGNS) = UPPER(I)
289
MESHES(I,SBRGNS) = MESH(I)
290
END DO
291
DO I = 1,NDIM
292
LOWER(I) = UPPER(I)
293
UPPER(I) = LOWER(I) + 2*WIDTH(I)
294
IF ( LOWER(I)+WIDTH(I) .LT. 1 ) GO TO 20
295
LOWER(I) = 0
296
UPPER(I) = LOWER(I) + 2*WIDTH(I)
297
END DO
298
FUNCLS = FUNCLS + SBRGNS*RULCLS
299
ENDIF
300
*
301
* Check for termination
302
*
303
30 FINEST = 0
304
ABSEST = 0
305
DO I = 1, SBRGNS
306
FINEST = FINEST + VALUES(I)
307
ABSEST = ABSEST + ERRORS(I)
308
END DO
309
IF ( ABSEST .GT. MAX( ABSREQ, RELREQ*ABS(FINEST) )
310
& .OR. FUNCLS .LT. MINCLS ) THEN
311
*
312
* Prepare to apply basic rule in (parts of) subregion with
313
* largest error.
314
*
315
TOP = PONTRS(1)
316
RGNCLS = RULCLS
317
DO I = 1,NDIM
318
LOWER(I) = LOWERS(I,TOP)
319
UPPER(I) = UPPERS(I,TOP)
320
MESH(I) = MESHES(I,TOP)
321
WIDTH(I) = (UPPER(I)-LOWER(I))/(2*MESH(I))
322
RGNCLS = RGNCLS*MESH(I)
323
END DO
324
CALL DIFFER(NDIM, LOWER, UPPER, WIDTH, WORK, WORK(NDIM+1),
325
& FUNCTN, DIVAXN, DIFCLS)
326
FUNCLS = FUNCLS + DIFCLS
327
RGNCLS = RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
328
IF ( FUNCLS + RGNCLS .LE. MAXCLS ) THEN
329
IF ( SBRGNS + 1 .LE. MXRGNS ) THEN
330
*
331
* Prepare to subdivide into two pieces.
332
*
333
NWRGNS = 1
334
WIDTH(DIVAXN) = WIDTH(DIVAXN)/2
335
ELSE
336
NWRGNS = 0
337
WIDTH(DIVAXN) = WIDTH(DIVAXN)
338
& *MESH(DIVAXN)/( MESH(DIVAXN) + 1 )
339
MESHES(DIVAXN,TOP) = MESH(DIVAXN) + 1
340
ENDIF
341
IF ( NWRGNS .GT. 0 ) THEN
342
*
343
* Only allow local subdivision when space is available.
344
*
345
DO J = SBRGNS+1,SBRGNS+NWRGNS
346
DO I = 1,NDIM
347
LOWERS(I,J) = LOWER(I)
348
UPPERS(I,J) = UPPER(I)
349
MESHES(I,J) = MESH(I)
350
END DO
351
END DO
352
UPPERS(DIVAXN,TOP) = LOWER(DIVAXN) + 2*WIDTH(DIVAXN)
353
LOWERS(DIVAXN,SBRGNS+1) = UPPERS(DIVAXN,TOP)
354
ENDIF
355
FUNCLS = FUNCLS + RGNCLS
356
CALL BASRUL(NDIM, LOWERS(1,TOP), UPPERS(1,TOP), WIDTH,
357
& FUNCTN, WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
358
& ERRORS(TOP), VALUES(TOP))
359
CALL TRESTR(TOP, SBRGNS, PONTRS, ERRORS)
360
DO I = SBRGNS+1, SBRGNS+NWRGNS
361
*
362
* Apply basic rule and store results in heap.
363
*
364
CALL BASRUL(NDIM, LOWERS(1,I), UPPERS(1,I), WIDTH,
365
& FUNCTN, WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
366
& ERRORS(I), VALUES(I))
367
CALL TRESTR(I, I, PONTRS, ERRORS)
368
END DO
369
SBRGNS = SBRGNS + NWRGNS
370
GO TO 30
371
ELSE
372
INFORM = 1
373
ENDIF
374
ELSE
375
INFORM = 0
376
ENDIF
377
MINCLS = FUNCLS
378
END
379
SUBROUTINE BSINIT(NDIM, W, LENRUL, G)
380
*
381
* For initializing basic rule weights and symmetric sum parameters.
382
*
383
INTEGER NDIM, LENRUL, RULPTS(6), I, J, NUMNUL, SDIM
384
PARAMETER ( NUMNUL = 4, SDIM = 12 )
385
DOUBLE PRECISION W(LENRUL,4), G(NDIM,LENRUL)
386
DOUBLE PRECISION LAM1, LAM2, LAM3, LAMP, RULCON
387
*
388
* The following code determines rule parameters and weights for a
389
* degree 7 rule (W(1,1),...,W(5,1)), two degree 5 comparison rules
390
* (W(1,2),...,W(5,2) and W(1,3),...,W(5,3)) and a degree 3
391
* comparison rule (W(1,4),...W(5,4)).
392
*
393
* If NDIM = 1, then LENRUL = 5 and total points = 9.
394
* If NDIM < SDIM, then LENRUL = 6 and
395
* total points = 1+2*NDIM*(NDIM+2)+2**NDIM.
396
* If NDIM > = SDIM, then LENRUL = 6 and
397
* total points = 1+2*NDIM*(1+2*NDIM).
398
*
399
DO I = 1,LENRUL
400
DO J = 1,NDIM
401
G(J,I) = 0
402
END DO
403
DO J = 1,NUMNUL
404
W(I,J) = 0
405
END DO
406
END DO
407
RULPTS(5) = 2*NDIM*(NDIM-1)
408
RULPTS(4) = 2*NDIM
409
RULPTS(3) = 2*NDIM
410
RULPTS(2) = 2*NDIM
411
RULPTS(1) = 1
412
LAMP = 0.85
413
LAM3 = 0.4707
414
LAM2 = 4/(15 - 5/LAM3)
415
W(5,1) = ( 3 - 5*LAM3 )/( 180*(LAM2-LAM3)*LAM2**2 )
416
IF ( NDIM .LT. SDIM ) THEN
417
LAM1 = 8*LAM3*(31*LAM3-15)/( (3*LAM3-1)*(5*LAM3-3)*35 )
418
W(LENRUL,1) = 1/(3*LAM3)**3/2**NDIM
419
ELSE
420
LAM1 = ( LAM3*(15 - 21*LAM2) + 35*(NDIM-1)*(LAM2-LAM3)/9 )
421
& / ( LAM3*(21 - 35*LAM2) + 35*(NDIM-1)*(LAM2/LAM3-1)/9 )
422
W(6,1) = 1/(4*(3*LAM3)**3)
423
ENDIF
424
W(3,1) = ( 15 - 21*(LAM3+LAM1) + 35*LAM3*LAM1 )
425
& /( 210*LAM2*(LAM2-LAM3)*(LAM2-LAM1) ) - 2*(NDIM-1)*W(5,1)
426
W(2,1) = ( 15 - 21*(LAM3+LAM2) + 35*LAM3*LAM2 )
427
& /( 210*LAM1*(LAM1-LAM3)*(LAM1-LAM2) )
428
IF ( NDIM .LT. SDIM ) THEN
429
RULPTS(LENRUL) = 2**NDIM
430
LAM3 = SQRT(LAM3)
431
DO I = 1,NDIM
432
G(I,LENRUL) = LAM3
433
END DO
434
ELSE
435
W(6,1) = 1/(4*(3*LAM3)**3)
436
RULPTS(6) = 2*NDIM*(NDIM-1)
437
LAM3 = SQRT(LAM3)
438
DO I = 1,2
439
G(I,6) = LAM3
440
END DO
441
ENDIF
442
IF ( NDIM .GT. 1 ) THEN
443
W(5,2) = 1/(6*LAM2)**2
444
W(5,3) = 1/(6*LAM2)**2
445
ENDIF
446
W(3,2) = ( 3 - 5*LAM1 )/( 30*LAM2*(LAM2-LAM1) )
447
& - 2*(NDIM-1)*W(5,2)
448
W(2,2) = ( 3 - 5*LAM2 )/( 30*LAM1*(LAM1-LAM2) )
449
W(4,3) = ( 3 - 5*LAM2 )/( 30*LAMP*(LAMP-LAM2) )
450
W(3,3) = ( 3 - 5*LAMP )/( 30*LAM2*(LAM2-LAMP) )
451
& - 2*(NDIM-1)*W(5,3)
452
W(2,4) = 1/(6*LAM1)
453
LAMP = SQRT(LAMP)
454
LAM2 = SQRT(LAM2)
455
LAM1 = SQRT(LAM1)
456
G(1,2) = LAM1
457
G(1,3) = LAM2
458
G(1,4) = LAMP
459
IF ( NDIM .GT. 1 ) THEN
460
G(1,5) = LAM2
461
G(2,5) = LAM2
462
ENDIF
463
DO J = 1, NUMNUL
464
W(1,J) = 1
465
DO I = 2,LENRUL
466
W(1,J) = W(1,J) - RULPTS(I)*W(I,J)
467
END DO
468
END DO
469
RULCON = 2
470
CALL RULNRM( LENRUL, NUMNUL, RULPTS, W, RULCON )
471
END
472
*
473
SUBROUTINE RULNRM( LENRUL, NUMNUL, RULPTS, W, RULCON )
474
INTEGER LENRUL, NUMNUL, I, J, K, RULPTS(*)
475
DOUBLE PRECISION ALPHA, NORMCF, NORMNL, W(LENRUL, *), RULCON
476
*
477
* Compute orthonormalized null rules.
478
*
479
NORMCF = 0
480
DO I = 1,LENRUL
481
NORMCF = NORMCF + RULPTS(I)*W(I,1)*W(I,1)
482
END DO
483
DO K = 2,NUMNUL
484
DO I = 1,LENRUL
485
W(I,K) = W(I,K) - W(I,1)
486
END DO
487
DO J = 2,K-1
488
ALPHA = 0
489
DO I = 1,LENRUL
490
ALPHA = ALPHA + RULPTS(I)*W(I,J)*W(I,K)
491
END DO
492
ALPHA = -ALPHA/NORMCF
493
DO I = 1,LENRUL
494
W(I,K) = W(I,K) + ALPHA*W(I,J)
495
END DO
496
END DO
497
NORMNL = 0
498
DO I = 1,LENRUL
499
NORMNL = NORMNL + RULPTS(I)*W(I,K)*W(I,K)
500
END DO
501
ALPHA = SQRT(NORMCF/NORMNL)
502
DO I = 1,LENRUL
503
W(I,K) = ALPHA*W(I,K)
504
END DO
505
END DO
506
DO J = 2, NUMNUL
507
DO I = 1,LENRUL
508
W(I,J) = W(I,J)/RULCON
509
END DO
510
END DO
511
END
512
*--------------------------------------------------------------------------
513
DOUBLE PRECISION FUNCTION MVNFNC(N, W)
514
*
515
* Integrand subroutine
516
*
517
INTEGER N, INFIN(*), INFIS
518
DOUBLE PRECISION W(*), LOWER(*), UPPER(*), CORREL(*), ONE
519
INTEGER NL, IJ, I, J
520
PARAMETER ( NL = 100, ONE = 1 )
521
DOUBLE PRECISION COV((NL*(NL+1))/2), A(NL), B(NL), Y(NL), BVN
522
INTEGER INFI(NL)
523
DOUBLE PRECISION PROD, D1, E1, DI, EI, SUM, PHINV, D, E, MVNNIT
524
SAVE D1, E1, A, B, INFI, COV
525
DI = D1
526
EI = E1
527
PROD = E1 - D1
528
IJ = 1
529
DO I = 1,N
530
Y(I) = PHINV( DI + W(I)*(EI-DI) )
531
SUM = 0
532
DO J = 1,I
533
IJ = IJ + 1
534
SUM = SUM + COV(IJ)*Y(J)
535
END DO
536
IJ = IJ + 1
537
IF ( COV(IJ) .GT. 0 ) THEN
538
CALL LIMITS( A(I+1)-SUM, B(I+1)-SUM, INFI(I+1), DI, EI )
539
ELSE
540
DI = ( 1 + SIGN( ONE, A(I+1)-SUM ) )/2
541
EI = ( 1 + SIGN( ONE, B(I+1)-SUM ) )/2
542
ENDIF
543
PROD = PROD*(EI-DI)
544
END DO
545
MVNFNC = PROD
546
RETURN
547
*
548
* Entry point for intialization.
549
*
550
ENTRY MVNNIT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E)
551
MVNNIT = 0
552
*
553
* Initialization and computation of covariance Cholesky factor.
554
*
555
CALL NCVSRT(N, LOWER,UPPER,CORREL,INFIN,Y, INFIS,A,B,INFI,COV,D,E)
556
D1 = D
557
E1 = E
558
IF ( N - INFIS .EQ. 2 ) THEN
559
D = SQRT( 1 + COV(2)**2 )
560
A(2) = A(2)/D
561
B(2) = B(2)/D
562
E = BVN( A, B, INFI, COV(2)/D )
563
D = 0
564
INFIS = INFIS + 1
565
END IF
566
END
567
SUBROUTINE LIMITS( A, B, INFIN, LOWER, UPPER )
568
DOUBLE PRECISION A, B, LOWER, UPPER, PHI
569
INTEGER INFIN
570
LOWER = 0
571
UPPER = 1
572
IF ( INFIN .GE. 0 ) THEN
573
IF ( INFIN .NE. 0 ) LOWER = PHI(A)
574
IF ( INFIN .NE. 1 ) UPPER = PHI(B)
575
ENDIF
576
END
577
*--------------------------------------------------------------------------
578
DOUBLE PRECISION FUNCTION PHI(Z)
579
*
580
* Normal distribution probabilities accurate to 1.e-15.
581
* Z = no. of standard deviations from the mean.
582
*
583
* Based upon algorithm 5666 for the error function, from:
584
* Hart, J.F. et al, 'Computer Approximations', Wiley 1968
585
*
586
* Programmer: Alan Miller
587
*
588
* Latest revision - 30 March 1986
589
*
590
DOUBLE PRECISION P0, P1, P2, P3, P4, P5, P6,
591
& Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,
592
& Z, P, EXPNTL, CUTOFF, ROOTPI, ZABS
593
PARAMETER(
594
& P0 = 220.20 68679 12376 1D0,
595
& P1 = 221.21 35961 69931 1D0,
596
& P2 = 112.07 92914 97870 9D0,
597
& P3 = 33.912 86607 83830 0D0,
598
& P4 = 6.3739 62203 53165 0D0,
599
& P5 = .70038 30644 43688 1D0,
600
& P6 = .035262 49659 98910 9D0)
601
PARAMETER(
602
& Q0 = 440.41 37358 24752 2D0,
603
& Q1 = 793.82 65125 19948 4D0,
604
& Q2 = 637.33 36333 78831 1D0,
605
& Q3 = 296.56 42487 79673 7D0,
606
& Q4 = 86.780 73220 29460 8D0,
607
& Q5 = 16.064 17757 92069 5D0,
608
& Q6 = 1.7556 67163 18264 2D0,
609
& Q7 = .088388 34764 83184 4D0)
610
PARAMETER(ROOTPI = 2.5066 28274 63100 1D0)
611
PARAMETER(CUTOFF = 7.0710 67811 86547 5D0)
612
*
613
ZABS = ABS(Z)
614
*
615
* |Z| > 37
616
*
617
IF (ZABS .GT. 37) THEN
618
P = 0
619
ELSE
620
*
621
* |Z| <= 37
622
*
623
EXPNTL = EXP(-ZABS**2/2)
624
*
625
* |Z| < CUTOFF = 10/SQRT(2)
626
*
627
IF (ZABS .LT. CUTOFF) THEN
628
P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
629
& + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
630
& + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
631
& + Q0)
632
*
633
* |Z| >= CUTOFF.
634
*
635
ELSE
636
P = EXPNTL/(ZABS + 1/(ZABS + 2/(ZABS + 3/(ZABS + 4/
637
& (ZABS + 0.65D0)))))/ROOTPI
638
END IF
639
END IF
640
IF (Z .GT. 0) P = 1 - P
641
PHI = P
642
END
643
*
644
DOUBLE PRECISION FUNCTION PHINV(P)
645
*
646
* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
647
*
648
* Produces the normal deviate Z corresponding to a given lower
649
* tail area of P.
650
*
651
* The hash sums below are the sums of the mantissas of the
652
* coefficients. They are included for use in checking
653
* transcription.
654
*
655
DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2,
656
& A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
657
& C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
658
& E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
659
& P, Q, R
660
PARAMETER (SPLIT1 = 0.425, SPLIT2 = 5,
661
& CONST1 = 0.180625D0, CONST2 = 1.6D0)
662
*
663
* Coefficients for P close to 0.5
664
*
665
PARAMETER (
666
& A0 = 3.38713 28727 96366 6080D0,
667
& A1 = 1.33141 66789 17843 7745D+2,
668
& A2 = 1.97159 09503 06551 4427D+3,
669
& A3 = 1.37316 93765 50946 1125D+4,
670
& A4 = 4.59219 53931 54987 1457D+4,
671
& A5 = 6.72657 70927 00870 0853D+4,
672
& A6 = 3.34305 75583 58812 8105D+4,
673
& A7 = 2.50908 09287 30122 6727D+3,
674
& B1 = 4.23133 30701 60091 1252D+1,
675
& B2 = 6.87187 00749 20579 0830D+2,
676
& B3 = 5.39419 60214 24751 1077D+3,
677
& B4 = 2.12137 94301 58659 5867D+4,
678
& B5 = 3.93078 95800 09271 0610D+4,
679
& B6 = 2.87290 85735 72194 2674D+4,
680
& B7 = 5.22649 52788 52854 5610D+3)
681
* HASH SUM AB 55.88319 28806 14901 4439
682
*
683
* Coefficients for P not close to 0, 0.5 or 1.
684
*
685
PARAMETER (
686
& C0 = 1.42343 71107 49683 57734D0,
687
& C1 = 4.63033 78461 56545 29590D0,
688
& C2 = 5.76949 72214 60691 40550D0,
689
& C3 = 3.64784 83247 63204 60504D0,
690
& C4 = 1.27045 82524 52368 38258D0,
691
& C5 = 2.41780 72517 74506 11770D-1,
692
& C6 = 2.27238 44989 26918 45833D-2,
693
& C7 = 7.74545 01427 83414 07640D-4,
694
& D1 = 2.05319 16266 37758 82187D0,
695
& D2 = 1.67638 48301 83803 84940D0,
696
& D3 = 6.89767 33498 51000 04550D-1,
697
& D4 = 1.48103 97642 74800 74590D-1,
698
& D5 = 1.51986 66563 61645 71966D-2,
699
& D6 = 5.47593 80849 95344 94600D-4,
700
& D7 = 1.05075 00716 44416 84324D-9)
701
* HASH SUM CD 49.33206 50330 16102 89036
702
*
703
* Coefficients for P near 0 or 1.
704
*
705
PARAMETER (
706
& E0 = 6.65790 46435 01103 77720D0,
707
& E1 = 5.46378 49111 64114 36990D0,
708
& E2 = 1.78482 65399 17291 33580D0,
709
& E3 = 2.96560 57182 85048 91230D-1,
710
& E4 = 2.65321 89526 57612 30930D-2,
711
& E5 = 1.24266 09473 88078 43860D-3,
712
& E6 = 2.71155 55687 43487 57815D-5,
713
& E7 = 2.01033 43992 92288 13265D-7,
714
& F1 = 5.99832 20655 58879 37690D-1,
715
& F2 = 1.36929 88092 27358 05310D-1,
716
& F3 = 1.48753 61290 85061 48525D-2,
717
& F4 = 7.86869 13114 56132 59100D-4,
718
& F5 = 1.84631 83175 10054 68180D-5,
719
& F6 = 1.42151 17583 16445 88870D-7,
720
& F7 = 2.04426 31033 89939 78564D-15)
721
* HASH SUM EF 47.52583 31754 92896 71629
722
*
723
Q = ( 2*P - 1 )/2
724
IF ( ABS(Q) .LE. SPLIT1 ) THEN
725
R = CONST1 - Q*Q
726
PHINV = Q*(((((((A7*R + A6)*R + A5)*R + A4)*R + A3)
727
& *R + A2)*R + A1)*R + A0) /
728
& (((((((B7*R + B6)*R + B5)*R + B4)*R + B3)
729
& *R + B2)*R + B1)*R + 1)
730
ELSE
731
R = MIN( P, 1 - P )
732
IF (R .GT. 0) THEN
733
R = SQRT( -LOG(R) )
734
IF ( R .LE. SPLIT2 ) THEN
735
R = R - CONST2
736
PHINV = (((((((C7*R + C6)*R + C5)*R + C4)*R + C3)
737
& *R + C2)*R + C1)*R + C0) /
738
& (((((((D7*R + D6)*R + D5)*R + D4)*R + D3)
739
& *R + D2)*R + D1)*R + 1)
740
ELSE
741
R = R - SPLIT2
742
PHINV = (((((((E7*R + E6)*R + E5)*R + E4)*R + E3)
743
& *R + E2)*R + E1)*R + E0) /
744
& (((((((F7*R + F6)*R + F5)*R + F4)*R + F3)
745
& *R + F2)*R + F1)*R + 1)
746
END IF
747
ELSE
748
PHINV = 9
749
END IF
750
IF ( Q .LT. 0 ) PHINV = - PHINV
751
END IF
752
END
753
DOUBLE PRECISION FUNCTION BVN ( LOWER, UPPER, INFIN, CORREL )
754
*
755
* A function for computing bivariate normal probabilities.
756
*
757
* Parameters
758
*
759
* LOWER REAL, array of lower integration limits.
760
* UPPER REAL, array of upper integration limits.
761
* INFIN INTEGER, array of integration limits flags:
762
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
763
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
764
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
765
* CORREL REAL, correlation coefficient.
766
*
767
DOUBLE PRECISION LOWER(*), UPPER(*), CORREL, BVNU
768
INTEGER INFIN(*)
769
IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 2 ) THEN
770
BVN = BVNU ( LOWER(1), LOWER(2), CORREL )
771
+ - BVNU ( UPPER(1), LOWER(2), CORREL )
772
+ - BVNU ( LOWER(1), UPPER(2), CORREL )
773
+ + BVNU ( UPPER(1), UPPER(2), CORREL )
774
ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 1 ) THEN
775
BVN = BVNU ( LOWER(1), LOWER(2), CORREL )
776
+ - BVNU ( UPPER(1), LOWER(2), CORREL )
777
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 2 ) THEN
778
BVN = BVNU ( LOWER(1), LOWER(2), CORREL )
779
+ - BVNU ( LOWER(1), UPPER(2), CORREL )
780
ELSE IF ( INFIN(1) .EQ. 2 .AND. INFIN(2) .EQ. 0 ) THEN
781
BVN = BVNU ( -UPPER(1), -UPPER(2), CORREL )
782
+ - BVNU ( -LOWER(1), -UPPER(2), CORREL )
783
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 2 ) THEN
784
BVN = BVNU ( -UPPER(1), -UPPER(2), CORREL )
785
+ - BVNU ( -UPPER(1), -LOWER(2), CORREL )
786
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 0 ) THEN
787
BVN = BVNU ( LOWER(1), -UPPER(2), -CORREL )
788
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 1 ) THEN
789
BVN = BVNU ( -UPPER(1), LOWER(2), -CORREL )
790
ELSE IF ( INFIN(1) .EQ. 1 .AND. INFIN(2) .EQ. 1 ) THEN
791
BVN = BVNU ( LOWER(1), LOWER(2), CORREL )
792
ELSE IF ( INFIN(1) .EQ. 0 .AND. INFIN(2) .EQ. 0 ) THEN
793
BVN = BVNU ( -UPPER(1), -UPPER(2), CORREL )
794
END IF
795
END
796
DOUBLE PRECISION FUNCTION BVNU( SH, SK, R )
797
*
798
* A function for computing bivariate normal probabilities.
799
*
800
* Yihong Ge
801
* Department of Computer Science and Electrical Engineering
802
* Washington State University
803
* Pullman, WA 99164-2752
804
* Email : yge@eecs.wsu.edu
805
* and
806
* Alan Genz
807
* Department of Mathematics
808
* Washington State University
809
* Pullman, WA 99164-3113
810
* Email : alangenz@wsu.edu
811
*
812
* BVN - calculate the probability that X is larger than SH and Y is
813
* larger than SK.
814
*
815
* Parameters
816
*
817
* SH REAL, integration limit
818
* SK REAL, integration limit
819
* R REAL, correlation coefficient
820
* LG INTEGER, number of Gauss Rule Points and Weights
821
*
822
DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
823
INTEGER I, LG, NG
824
PARAMETER ( ZERO = 0, TWOPI = 6.2831 85307 179586 )
825
DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
826
DOUBLE PRECISION PHI, SN, ASR, H, K, BS, HS, HK
827
* Gauss Legendre Points and Weights, N = 6
828
DATA ( W(I,1), X(I,1), I = 1,3) /
829
& 0.1713244923791705D+00,-0.9324695142031522D+00,
830
& 0.3607615730481384D+00,-0.6612093864662647D+00,
831
& 0.4679139345726904D+00,-0.2386191860831970D+00/
832
* Gauss Legendre Points and Weights, N = 12
833
DATA ( W(I,2), X(I,2), I = 1,6) /
834
& 0.4717533638651177D-01,-0.9815606342467191D+00,
835
& 0.1069393259953183D+00,-0.9041172563704750D+00,
836
& 0.1600783285433464D+00,-0.7699026741943050D+00,
837
& 0.2031674267230659D+00,-0.5873179542866171D+00,
838
& 0.2334925365383547D+00,-0.3678314989981802D+00,
839
& 0.2491470458134029D+00,-0.1252334085114692D+00/
840
* Gauss Legendre Points and Weights, N = 20
841
DATA ( W(I,3), X(I,3), I = 1,10) /
842
& 0.1761400713915212D-01,-0.9931285991850949D+00,
843
& 0.4060142980038694D-01,-0.9639719272779138D+00,
844
& 0.6267204833410906D-01,-0.9122344282513259D+00,
845
& 0.8327674157670475D-01,-0.8391169718222188D+00,
846
& 0.1019301198172404D+00,-0.7463319064601508D+00,
847
& 0.1181945319615184D+00,-0.6360536807265150D+00,
848
& 0.1316886384491766D+00,-0.5108670019508271D+00,
849
& 0.1420961093183821D+00,-0.3737060887154196D+00,
850
& 0.1491729864726037D+00,-0.2277858511416451D+00,
851
& 0.1527533871307259D+00,-0.7652652113349733D-01/
852
SAVE X, W
853
IF ( ABS(R) .LT. 0.3 ) THEN
854
NG = 1
855
LG = 3
856
ELSE IF ( ABS(R) .LT. 0.75 ) THEN
857
NG = 2
858
LG = 6
859
ELSE
860
NG = 3
861
LG = 10
862
ENDIF
863
H = SH
864
K = SK
865
HK = H*K
866
BVN = 0
867
IF ( ABS(R) .LT. 0.925 ) THEN
868
HS = ( H*H + K*K )/2
869
ASR = ASIN(R)
870
DO 10 I = 1, LG
871
SN = SIN(ASR*( X(I,NG)+1 )/2)
872
BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
873
SN = SIN(ASR*(-X(I,NG)+1 )/2)
874
BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
875
10 CONTINUE
876
BVN = BVN*ASR/(2*TWOPI) + PHI(-H)*PHI(-K)
877
ELSE
878
IF ( R .LT. 0 ) THEN
879
K = -K
880
HK = -HK
881
ENDIF
882
IF ( ABS(R) .LT. 1 ) THEN
883
AS = ( 1 - R )*( 1 + R )
884
A = SQRT(AS)
885
BS = ( H - K )**2
886
C = ( 4 - HK )/8
887
D = ( 12 - HK )/16
888
BVN = A*EXP( -(BS/AS + HK)/2 )
889
+ *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
890
IF ( HK .GT. -160 ) THEN
891
B = SQRT(BS)
892
BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*PHI(-B/A)*B
893
+ *( 1 - C*BS*( 1 - D*BS/5 )/3 )
894
ENDIF
895
A = A/2
896
DO 20 I = 1, LG
897
XS = ( A*(X(I,NG)+1) )**2
898
RS = SQRT( 1 - XS )
899
BVN = BVN + A*W(I,NG)*
900
+ ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
901
+ - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
902
XS = AS*(-X(I,NG)+1)**2/4
903
RS = SQRT( 1 - XS )
904
BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
905
+ *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
906
+ - ( 1 + C*XS*( 1 + D*XS ) ) )
907
20 CONTINUE
908
BVN = -BVN/TWOPI
909
ENDIF
910
IF ( R .GT. 0 ) BVN = BVN + PHI( -MAX( H, K ) )
911
IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, PHI(-H) - PHI(-K) )
912
ENDIF
913
BVNU = BVN
914
END
915
*
916
*--------------------------------------------------------------------------
917
SUBROUTINE NCVSRT( N, LOWER, UPPER, CORREL, INFIN, Y, INFIS,
918
& A, B, INFI, COV, D, E )
919
*
920
* Subroutine to sort integration limits.
921
*
922
INTEGER N, INFI(*), INFIN(*), INFIS
923
DOUBLE PRECISION
924
& A(*), B(*), COV(*), LOWER(*), UPPER(*), CORREL(*), Y(*), D, E
925
INTEGER I, J, K, IJ, II, JMIN
926
DOUBLE PRECISION SUMSQ, ZERO
927
PARAMETER ( ZERO = 0 )
928
DOUBLE PRECISION AJ, BJ, SUM, SQTWPI
929
DOUBLE PRECISION CVDIAG, AMIN, BMIN, DMIN, EMIN, YL, YU
930
PARAMETER ( SQTWPI = 2.50662 82746 31000 50240 )
931
IJ = 0
932
II = 0
933
INFIS = 0
934
DO I = 1,N
935
INFI(I) = INFIN(I)
936
IF ( INFI(I) .LT. 0 ) THEN
937
INFIS = INFIS + 1
938
ELSE
939
A(I) = 0
940
B(I) = 0
941
IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I)
942
IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I)
943
ENDIF
944
DO J = 1,I-1
945
IJ = IJ + 1
946
II = II + 1
947
COV(IJ) = CORREL(II)
948
END DO
949
IJ = IJ + 1
950
COV(IJ) = 1
951
END DO
952
*
953
* First move any doubly infinite limits to innermost positions
954
*
955
IF ( INFIS .LT. N ) THEN
956
DO I = N,N-INFIS+1,-1
957
IF ( INFI(I) .GE. 0 ) THEN
958
DO J = 1,I-1
959
IF ( INFI(J) .LT. 0 ) THEN
960
CALL RCSWAP(J, I, A, B, INFI, N, COV)
961
GO TO 10
962
ENDIF
963
END DO
964
ENDIF
965
10 END DO
966
*
967
* Sort remaining limits and determine Cholesky decomposition
968
*
969
II = 0
970
DO I = 1,N-INFIS
971
*
972
* Determine the integration limits for variable with minimum
973
* expected probability and interchange that variable with Ith.
974
*
975
EMIN = 1
976
DMIN = 0
977
JMIN = I
978
CVDIAG = 0
979
IJ = II
980
DO J = I, N-INFIS
981
SUM = 0
982
SUMSQ = 0
983
DO K = 1, I-1
984
SUM = SUM + COV(IJ+K)*Y(K)
985
SUMSQ = SUMSQ + COV(IJ+K)**2
986
END DO
987
IJ = IJ + J
988
SUMSQ = SQRT( MAX( COV(IJ)-SUMSQ, ZERO ) )
989
IF ( SUMSQ .GT. 0 ) THEN
990
IF ( INFI(J) .NE. 0 ) AJ = ( A(J) - SUM )/SUMSQ
991
IF ( INFI(J) .NE. 1 ) BJ = ( B(J) - SUM )/SUMSQ
992
CALL LIMITS( AJ, BJ, INFI(J), D, E )
993
IF ( EMIN - DMIN .GE. E - D ) THEN
994
JMIN = J
995
IF ( INFI(J) .NE. 0 ) AMIN = AJ
996
IF ( INFI(J) .NE. 1 ) BMIN = BJ
997
DMIN = D
998
EMIN = E
999
CVDIAG = SUMSQ
1000
ENDIF
1001
ENDIF
1002
END DO
1003
IF ( JMIN .NE. I) CALL RCSWAP(I, JMIN, A, B, INFI, N, COV)
1004
*
1005
* Compute Ith column of Cholesky factor.
1006
*
1007
IJ = II + I
1008
COV(IJ) = CVDIAG
1009
DO J = I+1, N-INFIS
1010
IF ( CVDIAG .GT. 0 ) THEN
1011
SUM = COV(IJ+I)
1012
DO K = 1, I-1
1013
SUM = SUM - COV(II+K)*COV(IJ+K)
1014
END DO
1015
COV(IJ+I) = SUM/CVDIAG
1016
ELSE
1017
COV(IJ+I) = 0
1018
ENDIF
1019
IJ = IJ + J
1020
END DO
1021
*
1022
* Compute expected value for Ith integration variable and
1023
* scale Ith covariance matrix row and limits.
1024
*
1025
IF ( CVDIAG .GT. 0 ) THEN
1026
IF ( EMIN .GT. DMIN + 1D-8 ) THEN
1027
YL = 0
1028
YU = 0
1029
IF ( INFI(I) .NE. 0 ) YL = -EXP( -AMIN**2/2 )/SQTWPI
1030
IF ( INFI(I) .NE. 1 ) YU = -EXP( -BMIN**2/2 )/SQTWPI
1031
Y(I) = ( YU - YL )/( EMIN - DMIN )
1032
ELSE
1033
IF ( INFI(I) .EQ. 0 ) Y(I) = BMIN
1034
IF ( INFI(I) .EQ. 1 ) Y(I) = AMIN
1035
IF ( INFI(I) .EQ. 2 ) Y(I) = ( AMIN + BMIN )/2
1036
END IF
1037
DO J = 1,I
1038
II = II + 1
1039
COV(II) = COV(II)/CVDIAG
1040
END DO
1041
IF ( INFI(I) .NE. 0 ) A(I) = A(I)/CVDIAG
1042
IF ( INFI(I) .NE. 1 ) B(I) = B(I)/CVDIAG
1043
ELSE
1044
Y(I) = 0
1045
II = II + I
1046
ENDIF
1047
END DO
1048
CALL LIMITS( A(1), B(1), INFI(1), D, E)
1049
ENDIF
1050
END
1051
*--------------------------------------------------------------------------
1052
SUBROUTINE BASRUL( NDIM, A, B, WIDTH, FUNCTN, W, LENRUL, G,
1053
& CENTER, Z, RGNERT, BASEST )
1054
*
1055
* For application of basic integration rule
1056
*
1057
EXTERNAL FUNCTN
1058
INTEGER I, LENRUL, NDIM
1059
DOUBLE PRECISION
1060
& A(NDIM), B(NDIM), WIDTH(NDIM), FUNCTN, W(LENRUL,4),
1061
& G(NDIM,LENRUL), CENTER(NDIM), Z(NDIM), RGNERT, BASEST
1062
DOUBLE PRECISION
1063
& FULSUM, FSYMSM, RGNCMP, RGNVAL, RGNVOL, RGNCPT, RGNERR
1064
*
1065
* Compute Volume and Center of Subregion
1066
*
1067
RGNVOL = 1
1068
DO I = 1,NDIM
1069
RGNVOL = 2*RGNVOL*WIDTH(I)
1070
CENTER(I) = A(I) + WIDTH(I)
1071
END DO
1072
BASEST = 0
1073
RGNERT = 0
1074
*
1075
* Compute basic rule and error
1076
*
1077
10 RGNVAL = 0
1078
RGNERR = 0
1079
RGNCMP = 0
1080
RGNCPT = 0
1081
DO I = 1,LENRUL
1082
FSYMSM = FULSUM(NDIM, CENTER, WIDTH, Z, G(1,I), FUNCTN)
1083
* Basic Rule
1084
RGNVAL = RGNVAL + W(I,1)*FSYMSM
1085
* First comparison rule
1086
RGNERR = RGNERR + W(I,2)*FSYMSM
1087
* Second comparison rule
1088
RGNCMP = RGNCMP + W(I,3)*FSYMSM
1089
* Third Comparison rule
1090
RGNCPT = RGNCPT + W(I,4)*FSYMSM
1091
END DO
1092
*
1093
* Error estimation
1094
*
1095
RGNERR = SQRT(RGNCMP**2 + RGNERR**2)
1096
RGNCMP = SQRT(RGNCPT**2 + RGNCMP**2)
1097
IF ( 4*RGNERR .LT. RGNCMP ) RGNERR = RGNERR/2
1098
IF ( 2*RGNERR .GT. RGNCMP ) RGNERR = MAX( RGNERR, RGNCMP )
1099
RGNERT = RGNERT + RGNVOL*RGNERR
1100
BASEST = BASEST + RGNVOL*RGNVAL
1101
*
1102
* When subregion has more than one piece, determine next piece and
1103
* loop back to apply basic rule.
1104
*
1105
DO I = 1,NDIM
1106
CENTER(I) = CENTER(I) + 2*WIDTH(I)
1107
IF ( CENTER(I) .LT. B(I) ) GO TO 10
1108
CENTER(I) = A(I) + WIDTH(I)
1109
END DO
1110
END
1111
DOUBLE PRECISION FUNCTION FULSUM(S, CENTER, HWIDTH, X, G, F)
1112
*
1113
**** To compute fully symmetric basic rule sum
1114
*
1115
EXTERNAL F
1116
INTEGER S, IXCHNG, LXCHNG, I, L
1117
DOUBLE PRECISION CENTER(S), HWIDTH(S), X(S), G(S), F
1118
DOUBLE PRECISION INTSUM, GL, GI
1119
FULSUM = 0
1120
*
1121
* Compute centrally symmetric sum for permutation of G
1122
*
1123
10 INTSUM = 0
1124
DO I = 1,S
1125
X(I) = CENTER(I) + G(I)*HWIDTH(I)
1126
END DO
1127
20 INTSUM = INTSUM + F(S,X)
1128
DO I = 1,S
1129
G(I) = -G(I)
1130
X(I) = CENTER(I) + G(I)*HWIDTH(I)
1131
IF ( G(I) .LT. 0 ) GO TO 20
1132
END DO
1133
FULSUM = FULSUM + INTSUM
1134
*
1135
* Find next distinct permuation of G and loop back for next sum
1136
*
1137
DO I = 2,S
1138
IF ( G(I-1) .GT. G(I) ) THEN
1139
GI = G(I)
1140
IXCHNG = I - 1
1141
DO L = 1,(I-1)/2
1142
GL = G(L)
1143
G(L) = G(I-L)
1144
G(I-L) = GL
1145
IF ( GL .LE. GI ) IXCHNG = IXCHNG - 1
1146
IF ( G(L) .GT. GI ) LXCHNG = L
1147
END DO
1148
IF ( G(IXCHNG) .LE. GI ) IXCHNG = LXCHNG
1149
G(I) = G(IXCHNG)
1150
G(IXCHNG) = GI
1151
GO TO 10
1152
ENDIF
1153
END DO
1154
*
1155
* End loop for permutations of G and associated sums
1156
*
1157
* Restore original order to G's
1158
*
1159
DO I = 1,S/2
1160
GI = G(I)
1161
G(I) = G(S+1-I)
1162
G(S+1-I) = GI
1163
END DO
1164
END
1165
SUBROUTINE DIFFER(NDIM, A, B, WIDTH, Z, DIF, FUNCTN,
1166
& DIVAXN, DIFCLS)
1167
*
1168
* Compute fourth differences and subdivision axes
1169
*
1170
EXTERNAL FUNCTN
1171
INTEGER I, NDIM, DIVAXN, DIFCLS
1172
DOUBLE PRECISION
1173
& A(NDIM), B(NDIM), WIDTH(NDIM), Z(NDIM), DIF(NDIM), FUNCTN
1174
DOUBLE PRECISION FRTHDF, FUNCEN, WIDTHI
1175
DIFCLS = 0
1176
DIVAXN = MOD( DIVAXN, NDIM ) + 1
1177
IF ( NDIM .GT. 1 ) THEN
1178
DO I = 1,NDIM
1179
DIF(I) = 0
1180
Z(I) = A(I) + WIDTH(I)
1181
END DO
1182
10 FUNCEN = FUNCTN(NDIM, Z)
1183
DO I = 1,NDIM
1184
WIDTHI = WIDTH(I)/5
1185
FRTHDF = 6*FUNCEN
1186
Z(I) = Z(I) - 4*WIDTHI
1187
FRTHDF = FRTHDF + FUNCTN(NDIM,Z)
1188
Z(I) = Z(I) + 2*WIDTHI
1189
FRTHDF = FRTHDF - 4*FUNCTN(NDIM,Z)
1190
Z(I) = Z(I) + 4*WIDTHI
1191
FRTHDF = FRTHDF - 4*FUNCTN(NDIM,Z)
1192
Z(I) = Z(I) + 2*WIDTHI
1193
FRTHDF = FRTHDF + FUNCTN(NDIM,Z)
1194
* Do not include differences below roundoff
1195
IF ( FUNCEN + FRTHDF/8 .NE. FUNCEN )
1196
& DIF(I) = DIF(I) + ABS(FRTHDF)*WIDTH(I)
1197
Z(I) = Z(I) - 4*WIDTHI
1198
END DO
1199
DIFCLS = DIFCLS + 4*NDIM + 1
1200
DO I = 1,NDIM
1201
Z(I) = Z(I) + 2*WIDTH(I)
1202
IF ( Z(I) .LT. B(I) ) GO TO 10
1203
Z(I) = A(I) + WIDTH(I)
1204
END DO
1205
DO I = 1,NDIM
1206
IF ( DIF(DIVAXN) .LT. DIF(I) ) DIVAXN = I
1207
END DO
1208
ENDIF
1209
END
1210
*--------
1211
SUBROUTINE TRESTR(POINTR, SBRGNS, PONTRS, RGNERS)
1212
****BEGIN PROLOGUE TRESTR
1213
****PURPOSE TRESTR maintains a heap for subregions.
1214
****DESCRIPTION TRESTR maintains a heap for subregions.
1215
* The subregions are ordered according to the size of the
1216
* greatest error estimates of each subregion (RGNERS).
1217
*
1218
* PARAMETERS
1219
*
1220
* POINTR Integer.
1221
* The index for the subregion to be inserted in the heap.
1222
* SBRGNS Integer.
1223
* Number of subregions in the heap.
1224
* PONTRS Real array of dimension SBRGNS.
1225
* Used to store the indices for the greatest estimated errors
1226
* for each subregion.
1227
* RGNERS Real array of dimension SBRGNS.
1228
* Used to store the greatest estimated errors for each
1229
* subregion.
1230
*
1231
****ROUTINES CALLED NONE
1232
****END PROLOGUE TRESTR
1233
*
1234
* Global variables.
1235
*
1236
INTEGER POINTR, SBRGNS
1237
DOUBLE PRECISION PONTRS(*), RGNERS(*)
1238
*
1239
* Local variables.
1240
*
1241
* RGNERR Intermediate storage for the greatest error of a subregion.
1242
* SUBRGN Position of child/parent subregion in the heap.
1243
* SUBTMP Position of parent/child subregion in the heap.
1244
*
1245
INTEGER SUBRGN, SUBTMP, POINTP, POINTS
1246
DOUBLE PRECISION RGNERR
1247
*
1248
****FIRST PROCESSING STATEMENT TRESTR
1249
*
1250
RGNERR = RGNERS(POINTR)
1251
IF ( POINTR .EQ. PONTRS(1) ) THEN
1252
*
1253
* Move the new subregion inserted at the top of the heap
1254
* to its correct position in the heap.
1255
*
1256
SUBRGN = 1
1257
10 SUBTMP = 2*SUBRGN
1258
IF ( SUBTMP .LE. SBRGNS ) THEN
1259
IF ( SUBTMP .NE. SBRGNS ) THEN
1260
*
1261
* Find maximum of left and right child.
1262
*
1263
POINTS = PONTRS(SUBTMP)
1264
POINTP = PONTRS(SUBTMP+1)
1265
IF ( RGNERS(POINTS) .LT.
1266
+ RGNERS(POINTP) ) SUBTMP = SUBTMP + 1
1267
ENDIF
1268
*
1269
* Compare maximum child with parent.
1270
* If parent is maximum, then done.
1271
*
1272
POINTS = PONTRS(SUBTMP)
1273
IF ( RGNERR .LT. RGNERS(POINTS) ) THEN
1274
*
1275
* Move the pointer at position subtmp up the heap.
1276
*
1277
PONTRS(SUBRGN) = PONTRS(SUBTMP)
1278
SUBRGN = SUBTMP
1279
GO TO 10
1280
ENDIF
1281
ENDIF
1282
ELSE
1283
*
1284
* Insert new subregion in the heap.
1285
*
1286
SUBRGN = SBRGNS
1287
20 SUBTMP = SUBRGN/2
1288
IF ( SUBTMP .GE. 1 ) THEN
1289
*
1290
* Compare child with parent. If parent is maximum, then done.
1291
*
1292
POINTS = PONTRS(SUBTMP)
1293
IF ( RGNERR .GT. RGNERS(POINTS) ) THEN
1294
*
1295
* Move the pointer at position subtmp down the heap.
1296
*
1297
PONTRS(SUBRGN) = PONTRS(SUBTMP)
1298
SUBRGN = SUBTMP
1299
GO TO 20
1300
ENDIF
1301
ENDIF
1302
ENDIF
1303
PONTRS(SUBRGN) = POINTR
1304
*
1305
****END TRESTR
1306
*
1307
END
1308
*--------------------------------------------------------------------------
1309
SUBROUTINE RCSWAP(P, Q, A, B, INFIN, N, C)
1310
*
1311
* Swaps rows and columns P and Q in situ.
1312
*
1313
DOUBLE PRECISION A(*), B(*), C(*), T
1314
INTEGER INFIN(*), P, Q, N, I, J, II, JJ
1315
T = A(P)
1316
A(P) = A(Q)
1317
A(Q) = T
1318
T = B(P)
1319
B(P) = B(Q)
1320
B(Q) = T
1321
J = INFIN(P)
1322
INFIN(P) = INFIN(Q)
1323
INFIN(Q) = J
1324
JJ = (P*(P-1))/2
1325
II = (Q*(Q-1))/2
1326
T = C(JJ+P)
1327
C(JJ+P) = C(II+Q)
1328
C(II+Q) = T
1329
DO J = 1, P-1
1330
T = C(JJ+J)
1331
C(JJ+J) = C(II+J)
1332
C(II+J) = T
1333
END DO
1334
JJ = JJ + P
1335
DO I = P+1, Q-1
1336
T = C(JJ+P)
1337
C(JJ+P) = C(II+I)
1338
C(II+I) = T
1339
JJ = JJ + I
1340
END DO
1341
II = II + Q
1342
DO I = Q+1, N
1343
T = C(II+P)
1344
C(II+P) = C(II+Q)
1345
C(II+Q) = T
1346
II = II + I
1347
END DO
1348
END
1349
*--------------------------------------------------------------------------
1350
*--------------------------------------------------------------------------
1351
*--------------------------------------------------------------------------
1352
*--------------------------------------------------------------------------
1353
*--------------------------------------------------------------------------
1354
*--------------------------------------------------------------------------
1355
*--------------------------------------------------------------------------
1356
*--------------------------------------------------------------------------
1357
1358
*
1359
SUBROUTINE SADMVT(N, NU, LOWER, UPPER, INFIN, CORREL, MAXPTS,
1360
* ABSEPS, RELEPS, ERROR, VALUE, INFORM)
1361
*
1362
* A subroutine for computing multivariate t probabilities.
1363
* Alan Genz
1364
* Department of Mathematics
1365
* Washington State University
1366
* Pullman, WA 99164-3113
1367
* Email : AlanGenz@wsu.edu
1368
*
1369
* Parameters
1370
*
1371
* N INTEGER, the number of variables.
1372
* NU INTEGER, the number of degrees of freedom.
1373
* LOWER REAL, array of lower integration limits.
1374
* UPPER REAL, array of upper integration limits.
1375
* INFIN INTEGER, array of integration limits flags:
1376
* if INFIN(I) < 0, Ith limits are (-infinity, infinity);
1377
* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
1378
* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
1379
* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
1380
* CORREL REAL, array of correlation coefficients; the correlation
1381
* coefficient in row I column J of the correlation matrix
1382
* should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
1383
* MAXPTS INTEGER, maximum number of function values allowed. This
1384
* parameter can be used to limit the time taken. A sensible
1385
* strategy is to start with MAXPTS = 1000*N, and then
1386
* increase MAXPTS if ERROR is too large.
1387
* ABSEPS REAL absolute error tolerance.
1388
* RELEPS REAL relative error tolerance.
1389
* ERROR REAL, estimated absolute error, with 99% confidence level.
1390
* VALUE REAL, estimated value for the integral
1391
* INFORM INTEGER, termination status parameter:
1392
* if INFORM = 0, normal completion with ERROR < EPS;
1393
* if INFORM = 1, completion with ERROR > EPS and MAXPTS
1394
* function vaules used; increase MAXPTS to
1395
* decrease ERROR;
1396
* if INFORM = 2, N > 20 or N < 1.
1397
*
1398
EXTERNAL FNCMVT
1399
INTEGER NL, N, NU, M, INFIN(*), LENWRK, MAXPTS, INFORM, INFIS,
1400
& RULCLS, TOTCLS, NEWCLS, MAXCLS
1401
DOUBLE PRECISION CORREL(*), LOWER(*), UPPER(*), ABSEPS, RELEPS,
1402
& ERROR, VALUE, OLDVAL, D, E, MVTNIT
1403
PARAMETER ( NL = 20 )
1404
PARAMETER ( LENWRK = 20*NL**2 )
1405
DOUBLE PRECISION WORK(LENWRK)
1406
IF ( N .GT. 20 .OR. N .LT. 1 ) THEN
1407
INFORM = 2
1408
VALUE = 0
1409
ERROR = 1
1410
RETURN
1411
ENDIF
1412
INFORM = MVTNIT( N, NU, CORREL, LOWER, UPPER, INFIN, INFIS, D, E )
1413
M = N - INFIS
1414
IF ( M .EQ. 0 ) THEN
1415
VALUE = 1
1416
ERROR = 0
1417
ELSE IF ( M .EQ. 1 ) THEN
1418
VALUE = E - D
1419
ERROR = 2E-16
1420
ELSE
1421
*
1422
* Call the subregion adaptive integration subroutine
1423
*
1424
M = M - 1
1425
RULCLS = 1
1426
CALL ADAPT( M, RULCLS, 0, FNCMVT, ABSEPS, RELEPS,
1427
* LENWRK, WORK, ERROR, VALUE, INFORM )
1428
MAXCLS = MIN( 10*RULCLS, MAXPTS )
1429
TOTCLS = 0
1430
CALL ADAPT( M, TOTCLS, MAXCLS, FNCMVT, ABSEPS, RELEPS,
1431
* LENWRK, WORK, ERROR, VALUE, INFORM )
1432
IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
1433
10 OLDVAL = VALUE
1434
MAXCLS = MAX( 2*RULCLS, MIN( 3*MAXCLS/2, MAXPTS - TOTCLS ) )
1435
NEWCLS = -1
1436
CALL ADAPT( M, NEWCLS, MAXCLS, FNCMVT, ABSEPS, RELEPS,
1437
* LENWRK, WORK, ERROR, VALUE, INFORM )
1438
TOTCLS = TOTCLS + NEWCLS
1439
ERROR = ABS(VALUE-OLDVAL) + SQRT(RULCLS*ERROR**2/TOTCLS)
1440
IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
1441
IF ( MAXPTS - TOTCLS .GT. 2*RULCLS ) GO TO 10
1442
ELSE
1443
INFORM = 0
1444
END IF
1445
ENDIF
1446
ENDIF
1447
END
1448
*
1449
DOUBLE PRECISION FUNCTION FNCMVT(N, W)
1450
*
1451
* Integrand subroutine
1452
*
1453
INTEGER N, NUIN, INFIN(*), INFIS
1454
DOUBLE PRECISION W(*), LOWER(*), UPPER(*), CORREL(*), D, E
1455
INTEGER NL, IJ, I, J, NU
1456
PARAMETER ( NL = 20 )
1457
DOUBLE PRECISION COV((NL*(NL+1))/2), A(NL), B(NL), Y(NL)
1458
INTEGER INFI(NL)
1459
DOUBLE PRECISION PROD, D1, E1, DI, EI, SUM, STDINV, YD, UI, MVTNIT
1460
SAVE NU, D1, E1, A, B, INFI, COV
1461
DI = D1
1462
EI = E1
1463
PROD = EI - DI
1464
IJ = 1
1465
YD = 1
1466
DO I = 1, N
1467
UI = STDINV( NU+I-1, DI + W(I)*( EI - DI ) )
1468
Y(I) = UI/YD
1469
YD = YD/SQRT( 1 + ( UI - 1 )*( UI + 1 )/( NU + I ) )
1470
SUM = 0
1471
DO J = 1, I
1472
IJ = IJ + 1
1473
SUM = SUM + COV(IJ)*Y(J)
1474
END DO
1475
IJ = IJ + 1
1476
CALL MVTLMS( NU+I, ( A(I+1) - SUM )*YD, ( B(I+1) - SUM )*YD,
1477
& INFI(I+1), DI, EI )
1478
PROD = PROD*( EI - DI )
1479
END DO
1480
FNCMVT = PROD
1481
RETURN
1482
*
1483
* Entry point for intialization
1484
*
1485
ENTRY MVTNIT( N, NUIN, CORREL, LOWER, UPPER, INFIN, INFIS, D, E )
1486
MVTNIT = 0
1487
*
1488
* Initialization and computation of covariance matrix Cholesky factor
1489
*
1490
CALL MVTSRT( N, NUIN, LOWER, UPPER, CORREL, INFIN, Y, INFIS,
1491
& A, B, INFI, COV, D, E )
1492
NU = NUIN
1493
D1 = D
1494
E1 = E
1495
END
1496
SUBROUTINE MVTLMS( NU, A, B, INFIN, LOWER, UPPER )
1497
DOUBLE PRECISION A, B, LOWER, UPPER, STUDNT
1498
INTEGER NU, INFIN
1499
LOWER = 0
1500
UPPER = 1
1501
IF ( INFIN .GE. 0 ) THEN
1502
IF ( INFIN .NE. 0 ) LOWER = STUDNT( NU, A )
1503
IF ( INFIN .NE. 1 ) UPPER = STUDNT( NU, B )
1504
ENDIF
1505
END
1506
*
1507
SUBROUTINE MVTSRT( N, NU, LOWER, UPPER, CORREL, INFIN, Y, INFIS,
1508
& A, B, INFI, COV, D, E )
1509
*
1510
* Sort limits
1511
*
1512
INTEGER N, NU, INFI(*), INFIN(*), INFIS
1513
DOUBLE PRECISION
1514
& A(*), B(*), COV(*), LOWER(*), UPPER(*), CORREL(*), Y(*), D, E
1515
INTEGER I, J, K, IJ, II, JMIN
1516
DOUBLE PRECISION SUMSQ, ZERO, TWO, PI, CVDIAG
1517
DOUBLE PRECISION AI, BI, SUM, YL, YU, YD
1518
DOUBLE PRECISION AMIN, BMIN, DMIN, EMIN, CON, CONODD, CONEVN
1519
PARAMETER ( ZERO = 0, TWO = 2, PI = 3.14159 26535 89793 23844 )
1520
IJ = 0
1521
II = 0
1522
INFIS = 0
1523
DO I = 1, N
1524
INFI(I) = INFIN(I)
1525
IF ( INFI(I) .LT. 0 ) THEN
1526
INFIS = INFIS + 1
1527
ELSE
1528
A(I) = 0
1529
B(I) = 0
1530
IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I)
1531
IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I)
1532
ENDIF
1533
DO J = 1,I-1
1534
IJ = IJ + 1
1535
II = II + 1
1536
COV(IJ) = CORREL(II)
1537
END DO
1538
IJ = IJ + 1
1539
COV(IJ) = 1
1540
END DO
1541
CONODD = 1/PI
1542
CONEVN = 1/TWO
1543
DO I = 1, NU - 1
1544
IF ( MOD(I,2) .EQ. 0 ) THEN
1545
IF ( I .GT. 2 ) CONEVN = CONEVN*(I-1)/(I-2)
1546
ELSE
1547
IF ( I .GT. 2 ) CONODD = CONODD*(I-1)/(I-2)
1548
ENDIF
1549
END DO
1550
*
1551
* First move any doubly infinite limits to innermost positions
1552
*
1553
IF ( INFIS .LT. N ) THEN
1554
DO I = N, N-INFIS+1, -1
1555
IF ( INFI(I) .GE. 0 ) THEN
1556
DO J = 1, I-1
1557
IF ( INFI(J) .LT. 0 ) THEN
1558
CALL RCSWAP( J, I, A, B, INFI, N, COV )
1559
GOTO 10
1560
ENDIF
1561
END DO
1562
ENDIF
1563
10 END DO
1564
*
1565
* Sort remaining limits and determine Cholesky decomposition
1566
*
1567
II = 0
1568
YD = 1
1569
DO I = 1, N-INFIS
1570
*
1571
* Determine the integration limits for variable with minimum
1572
* expected probability and interchange that variable with Ith.
1573
*
1574
EMIN = 1
1575
DMIN = 0
1576
JMIN = I
1577
CVDIAG = 0
1578
IJ = II
1579
DO J = I, N-INFIS
1580
SUM = 0
1581
SUMSQ = 0
1582
DO K = 1, I-1
1583
SUM = SUM + COV(IJ+K)*Y(K)
1584
SUMSQ = SUMSQ + COV(IJ+K)**2
1585
END DO
1586
IJ = IJ + J
1587
SUMSQ = SQRT( MAX( COV(IJ)-SUMSQ, ZERO ) )
1588
IF ( SUMSQ .GT. 0 ) THEN
1589
AI = YD*( A(J) - SUM )/SUMSQ
1590
BI = YD*( B(J) - SUM )/SUMSQ
1591
CALL MVTLMS( NU+J-1, AI, BI, INFI(J), D, E )
1592
IF ( EMIN - DMIN .GE. E - D ) THEN
1593
JMIN = J
1594
AMIN = AI
1595
BMIN = BI
1596
DMIN = D
1597
EMIN = E
1598
CVDIAG = SUMSQ
1599
ENDIF
1600
ENDIF
1601
END DO
1602
IF ( JMIN .NE. I ) CALL RCSWAP( I, JMIN, A,B, INFI, N,COV )
1603
*
1604
* Compute Ith column of Cholesky factor.
1605
*
1606
IJ = II + I
1607
COV(IJ) = CVDIAG
1608
DO J = I+1, N-INFIS
1609
IF ( CVDIAG .GT. 0 ) THEN
1610
SUM = COV(IJ+I)
1611
DO K = 1, I-1
1612
SUM = SUM - COV(II+K)*COV(IJ+K)
1613
END DO
1614
COV(IJ+I) = SUM/CVDIAG
1615
ELSE
1616
COV(IJ+I) = 0
1617
ENDIF
1618
IJ = IJ + J
1619
END DO
1620
*
1621
* Compute expected value for Ith integration variable and
1622
* scale Ith covariance matrix row and limits.
1623
*
1624
IF ( MOD(NU+I-1,2) .EQ. 0 ) THEN
1625
IF ( NU+I-3 .GT. 0 ) CONEVN = CONEVN*(NU+I-2)/(NU+I-3)
1626
CON = CONEVN
1627
ELSE
1628
IF ( NU+I-3 .GT. 0 ) CONODD = CONODD*(NU+I-2)/(NU+I-3)
1629
CON = CONODD
1630
ENDIF
1631
IF ( CVDIAG .GT. 0 ) THEN
1632
YL = 0
1633
YU = 0
1634
IF ( INFI(I) .NE. 0 .AND. NU+I-2 .GT. 0 )
1635
& YL = -CON*(NU+I-1)/(NU+I-2)
1636
& /( 1 + AMIN**2/(NU+I-1) )**( (NU+I-2)/TWO )
1637
IF ( INFI(I) .NE. 1 .AND. NU+I-2 .GT. 0 )
1638
& YU = -CON*(NU+I-1)/(NU+I-2)
1639
& /( 1 + BMIN**2/(NU+I-1) )**( (NU+I-2)/TWO )
1640
Y(I) = ( YU - YL )/( EMIN - DMIN )/YD
1641
DO J = 1,I
1642
II = II + 1
1643
COV(II) = COV(II)/CVDIAG
1644
END DO
1645
IF ( INFI(I) .NE. 0 ) A(I) = A(I)/CVDIAG
1646
IF ( INFI(I) .NE. 1 ) B(I) = B(I)/CVDIAG
1647
ELSE
1648
Y(I) = 0
1649
II = II + I
1650
ENDIF
1651
YD = YD/SQRT( 1 + ( Y(I)*YD + 1 )*( Y(I)*YD - 1 )/(NU+I) )
1652
END DO
1653
CALL MVTLMS( NU, A(1), B(1), INFI(1), D, E)
1654
ENDIF
1655
END
1656
*--
1657
DOUBLE PRECISION FUNCTION STUDNT( NU, T )
1658
*
1659
* Student t Distribution Function
1660
*
1661
* T
1662
* STUDNT = C I ( 1 + y*y/NU )**( -(NU+1)/2 ) dy
1663
* NU -INF
1664
*
1665
INTEGER NU, J
1666
DOUBLE PRECISION T, CSSTHE, SNTHE, POLYN, TT, TS, RN, PI, ZERO
1667
PARAMETER ( PI = 3.14159 26535 89793D0, ZERO = 0 )
1668
IF ( NU .EQ. 1 ) THEN
1669
STUDNT = ( 1 + 2*ATAN(T)/PI )/2
1670
ELSE IF ( NU .EQ. 2) THEN
1671
STUDNT = ( 1 + T/SQRT( 2 + T*T ))/2
1672
ELSE
1673
TT = T*T
1674
CSSTHE = 1/( 1 + TT/NU )
1675
POLYN = 1
1676
DO J = NU-2, 2, -2
1677
POLYN = 1 + ( J - 1 )*CSSTHE*POLYN/J
1678
END DO
1679
IF ( MOD( NU, 2 ) .EQ. 1 ) THEN
1680
RN = NU
1681
TS = T/SQRT(RN)
1682
STUDNT = ( 1 + 2*( ATAN(TS) + TS*CSSTHE*POLYN )/PI )/2
1683
ELSE
1684
SNTHE = T/SQRT( NU + TT )
1685
STUDNT = ( 1 + SNTHE*POLYN )/2
1686
END IF
1687
STUDNT = MAX( ZERO, STUDNT )
1688
ENDIF
1689
END
1690
*--
1691
DOUBLE PRECISION FUNCTION STDINV( N, Z )
1692
*
1693
* Inverse Student t Distribution Function
1694
*
1695
* STDINV
1696
* Z = C I (1 + y*y/N)**(-(N+1)/2) dy
1697
* N -INF
1698
*
1699
* Reference: G.W. Hill, Comm. ACM Algorithm 395
1700
* Comm. ACM 13 (1970), pp. 619-620.
1701
*
1702
* Conversions to double precision and other modifications by
1703
* Alan Genz, 1993-4.
1704
*
1705
INTEGER N
1706
DOUBLE PRECISION Z, P, PHINV, A, B, C, D, X, Y, PI, TWO
1707
DOUBLE PRECISION STUDNT, STDJAC
1708
PARAMETER ( PI = 3.14159 26535 89793D0, TWO = 2 )
1709
IF ( 0 .LT. Z .AND. Z .LT. 1 ) THEN
1710
IF ( N .EQ. 1 ) THEN
1711
STDINV = TAN( PI*( 2*Z - 1 )/2 )
1712
ELSE IF ( N .EQ. 2) THEN
1713
STDINV = ( 2*Z - 1 )/SQRT( 2*Z*( 1 - Z ) )
1714
ELSE
1715
IF ( 2*Z .GE. 1 ) THEN
1716
P = 2*( 1 - Z )
1717
ELSE
1718
P = 2*Z
1719
END IF
1720
A = 1/( N - 0.5 )
1721
B = 48/( A*A )
1722
C = ( ( 20700*A/B - 98 )*A - 16 )*A + 96.36
1723
D = ( ( 94.5/( B + C ) - 3 )/B + 1 )*SQRT( A*PI/2 )*N
1724
X = D*P
1725
Y = X**( TWO/N )
1726
IF ( Y .GT. A + 0.05 ) THEN
1727
X = PHINV( P/2 )
1728
Y = X*X
1729
IF ( N .LT. 5 ) C = C + 3*( N - 4.5 )*( 10*X + 6 )/100
1730
C = ( ( (D*X - 100)*X/20 - 7 )*X - 2 )*X + B + C
1731
Y = ( ( ( ( (4*Y+63)*Y/10+36 )*Y+94.5 )/C-Y-3 )/B + 1 )*X
1732
Y = A*Y*Y
1733
IF ( Y .GT. 0.002 ) THEN
1734
Y = EXP(Y) - 1
1735
ELSE
1736
Y = Y*( 1 + Y/2 )
1737
ENDIF
1738
ELSE
1739
Y = ( ( 1/( ( (N+6)/(N*Y) - 0.089*D - 0.822 )*(3*N+6) )
1740
& + 0.5/(N+4) )*Y - 1 )*(N+1)/(N+2) + 1/Y
1741
END IF
1742
STDINV = SQRT(N*Y)
1743
IF ( 2*Z .LT. 1 ) STDINV = -STDINV
1744
IF ( ABS( STDINV ) .GT. 0 ) THEN
1745
*
1746
* Use one third order correction to the single precision result
1747
*
1748
X = STDINV
1749
D = Z - STUDNT(N,X)
1750
STDINV = X + 2*D/( 2/STDJAC(N,X) - D*(N+1)/(N/X+X) )
1751
END IF
1752
END IF
1753
ELSE
1754
*
1755
* Use cutoff values for Z near 0 or 1.
1756
*
1757
STDINV = SQRT( N/( 2D-16*SQRT( 2*PI*N ) )**( TWO/N ) )
1758
IF ( 2*Z .LT. 1 ) STDINV = -STDINV
1759
END IF
1760
END
1761
*--
1762
DOUBLE PRECISION FUNCTION STDJAC( NU, T )
1763
*
1764
* Student t Distribution Transformation Jacobean
1765
*
1766
* T STDINV(NU,T)
1767
* I f(y) dy = I f(STDINV(NU,Z) STDJAC(NU,STDINV(NU,Z)) dZ
1768
* -INF 0
1769
*
1770
INTEGER NU, J
1771
DOUBLE PRECISION CONST, NUOLD, PI, T, TT
1772
PARAMETER ( PI = 3.14159 26535 89793D0 )
1773
SAVE NUOLD, CONST
1774
DATA NUOLD/ 0D0 /
1775
IF ( NU .EQ. 1 ) THEN
1776
STDJAC = PI*( 1 + T*T )
1777
ELSE IF ( NU .EQ. 2 ) THEN
1778
STDJAC = SQRT( 2 + T*T )**3
1779
ELSE
1780
IF ( NU .NE. NUOLD ) THEN
1781
NUOLD = NU
1782
IF ( MOD( NU, 2 ) .EQ. 0 ) THEN
1783
CONST = SQRT(NUOLD)*2
1784
ELSE
1785
CONST = SQRT(NUOLD)*PI
1786
END IF
1787
DO J = NU-2, 1, -2
1788
CONST = J*CONST/(J+1)
1789
END DO
1790
END IF
1791
TT = 1 + T*T/NU
1792
STDJAC = CONST*TT**( (NU+1)/2 )
1793
IF ( MOD( NU, 2 ) .EQ. 0 ) STDJAC = STDJAC*SQRT( TT )
1794
END IF
1795
END
1796
*
Loggerhead is a web-based interface for Breezy
Version: 2.0.1