00001 SUBROUTINE genmn(parm,x,work) 00002 C********************************************************************** 00003 C 00004 C SUBROUTINE GENMN(PARM,X,WORK) 00005 C GENerate Multivariate Normal random deviate 00006 C 00007 C 00008 C Arguments 00009 C 00010 C 00011 C PARM --> Parameters needed to generate multivariate normal 00012 C deviates (MEANV and Cholesky decomposition of 00013 C COVM). Set by a previous call to SETGMN. 00014 C 1 : 1 - size of deviate, P 00015 C 2 : P + 1 - mean vector 00016 C P+2 : P*(P+3)/2 + 1 - upper half of cholesky 00017 C decomposition of cov matrix 00018 C REAL PARM(*) 00019 C 00020 C X <-- Vector deviate generated. 00021 C REAL X(P) 00022 C 00023 C WORK <--> Scratch array 00024 C REAL WORK(P) 00025 C 00026 C 00027 C Method 00028 C 00029 C 00030 C 1) Generate P independent standard normal deviates - Ei ~ N(0,1) 00031 C 00032 C 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM 00033 C 00034 C 3) trans(A)E + MEANV ~ N(MEANV,COVM) 00035 C 00036 C********************************************************************** 00037 C .. Array Arguments .. 00038 REAL parm(*),work(*),x(*) 00039 C .. 00040 C .. Local Scalars .. 00041 REAL ae 00042 INTEGER i,icount,j,p 00043 C .. 00044 C .. External Functions .. 00045 REAL snorm 00046 EXTERNAL snorm 00047 C .. 00048 C .. Intrinsic Functions .. 00049 INTRINSIC int 00050 C .. 00051 C .. Executable Statements .. 00052 p = int(parm(1)) 00053 C 00054 C Generate P independent normal deviates - WORK ~ N(0,1) 00055 C 00056 DO 10,i = 1,p 00057 work(i) = snorm() 00058 10 CONTINUE 00059 DO 30,i = 1,p 00060 C 00061 C PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky 00062 C decomposition of the desired covariance matrix. 00063 C trans(A)(1,1) = PARM(P+2) 00064 C trans(A)(2,1) = PARM(P+3) 00065 C trans(A)(2,2) = PARM(P+2+P) 00066 C trans(A)(3,1) = PARM(P+4) 00067 C trans(A)(3,2) = PARM(P+3+P) 00068 C trans(A)(3,3) = PARM(P+2-1+2P) ... 00069 C 00070 C trans(A)*WORK + MEANV ~ N(MEANV,COVM) 00071 C 00072 icount = 0 00073 ae = 0.0 00074 DO 20,j = 1,i 00075 icount = icount + j - 1 00076 ae = ae + parm(i+ (j-1)*p-icount+p+1)*work(j) 00077 20 CONTINUE 00078 x(i) = ae + parm(i+1) 00079 30 CONTINUE 00080 RETURN 00081 C 00082 END
1.7.1
