GNU Octave
4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
Main Page
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Properties
Friends
Macros
Pages
liboctave
cruft
ranlib
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
genmn
subroutine genmn(parm, x, work)
Definition:
genmn.f:1
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by
1.8.8