GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
genmn.f
Go to the documentation of this file.
1  SUBROUTINE genmn(parm,x,work)
2 C**********************************************************************
3 C
4 C SUBROUTINE GENMN(PARM,X,WORK)
5 C GENerate Multivariate Normal random deviate
6 C
7 C
8 C Arguments
9 C
10 C
11 C PARM --> Parameters needed to generate multivariate normal
12 C deviates (MEANV and Cholesky decomposition of
13 C COVM). Set by a previous call to SETGMN.
14 C 1 : 1 - size of deviate, P
15 C 2 : P + 1 - mean vector
16 C P+2 : P*(P+3)/2 + 1 - upper half of cholesky
17 C decomposition of cov matrix
18 C REAL PARM(*)
19 C
20 C X <-- Vector deviate generated.
21 C REAL X(P)
22 C
23 C WORK <--> Scratch array
24 C REAL WORK(P)
25 C
26 C
27 C Method
28 C
29 C
30 C 1) Generate P independent standard normal deviates - Ei ~ N(0,1)
31 C
32 C 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
33 C
34 C 3) trans(A)E + MEANV ~ N(MEANV,COVM)
35 C
36 C**********************************************************************
37 C .. Array Arguments ..
38  REAL parm(*),work(*),x(*)
39 C ..
40 C .. Local Scalars ..
41  REAL ae
42  INTEGER i,icount,j,p
43 C ..
44 C .. External Functions ..
45  REAL snorm
46  EXTERNAL snorm
47 C ..
48 C .. Intrinsic Functions ..
49  INTRINSIC int
50 C ..
51 C .. Executable Statements ..
52  p = int(parm(1))
53 C
54 C Generate P independent normal deviates - WORK ~ N(0,1)
55 C
56  DO 10,i = 1,p
57  work(i) = snorm()
58  10 CONTINUE
59  DO 30,i = 1,p
60 C
61 C PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
62 C decomposition of the desired covariance matrix.
63 C trans(A)(1,1) = PARM(P+2)
64 C trans(A)(2,1) = PARM(P+3)
65 C trans(A)(2,2) = PARM(P+2+P)
66 C trans(A)(3,1) = PARM(P+4)
67 C trans(A)(3,2) = PARM(P+3+P)
68 C trans(A)(3,3) = PARM(P+2-1+2P) ...
69 C
70 C trans(A)*WORK + MEANV ~ N(MEANV,COVM)
71 C
72  icount = 0
73  ae = 0.0
74  DO 20,j = 1,i
75  icount = icount + j - 1
76  ae = ae + parm(i+ (j-1)*p-icount+p+1)*work(j)
77  20 CONTINUE
78  x(i) = ae + parm(i+1)
79  30 CONTINUE
80  RETURN
81 C
82  END
subroutine genmn(parm, x, work)
Definition: genmn.f:1