GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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