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
|
SUBROUTINE genmn(parm,x,work)
C**********************************************************************
C
C SUBROUTINE GENMN(PARM,X,WORK)
C GENerate Multivariate Normal random deviate
C
C
C Arguments
C
C
C PARM --> Parameters needed to generate multivariate normal
C deviates (MEANV and Cholesky decomposition of
C COVM). Set by a previous call to SETGMN.
C 1 : 1 - size of deviate, P
C 2 : P + 1 - mean vector
C P+2 : P*(P+3)/2 + 1 - upper half of cholesky
C decomposition of cov matrix
C DOUBLE PRECISION PARM(*)
C
C X <-- Vector deviate generated.
C DOUBLE PRECISION X(P)
C
C WORK <--> Scratch array
C DOUBLE PRECISION WORK(P)
C
C
C Method
C
C
C 1) Generate P independent standard normal deviates - Ei ~ N(0,1)
C
C 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
C
C 3) trans(A)E + MEANV ~ N(MEANV,COVM)
C
C**********************************************************************
C .. Array Arguments ..
DOUBLE PRECISION parm(*),work(*),x(*)
C ..
C .. Local Scalars ..
DOUBLE PRECISION ae
INTEGER i,icount,j,p
C ..
C .. External Functions ..
DOUBLE PRECISION snorm
EXTERNAL snorm
C ..
C .. Intrinsic Functions ..
INTRINSIC int
C ..
C .. Executable Statements ..
p = int(parm(1))
C
C Generate P independent normal deviates - WORK ~ N(0,1)
C
DO 10,i = 1,p
work(i) = snorm()
10 CONTINUE
DO 30,i = 1,p
C
C PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
C decomposition of the desired covariance matrix.
C trans(A)(1,1) = PARM(P+2)
C trans(A)(2,1) = PARM(P+3)
C trans(A)(2,2) = PARM(P+2+P)
C trans(A)(3,1) = PARM(P+4)
C trans(A)(3,2) = PARM(P+3+P)
C trans(A)(3,3) = PARM(P+2-1+2P) ...
C
C trans(A)*WORK + MEANV ~ N(MEANV,COVM)
C
icount = 0
ae = 0.0
DO 20,j = 1,i
icount = icount + j - 1
ae = ae + parm(i+ (j-1)*p-icount+p+1)*work(j)
20 CONTINUE
x(i) = ae + parm(i+1)
30 CONTINUE
RETURN
C
END
|
