Mercurial > octave-nkf
view libcruft/ranlib/genmn.f @ 13741:b0e70a71647b stable release-3-4-3
Version 3.4.3 released.
* configure.ac (AC_INIT): Version number is now 3.4.3.
(OCTAVE_RELEASE_DATE): Update to 2011-10-10.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 10 Oct 2011 10:23:45 -0400 |
parents | 30c606bec7a8 |
children |
line wrap: on
line source
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 REAL PARM(*) C C X <-- Vector deviate generated. C REAL X(P) C C WORK <--> Scratch array C REAL 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 .. REAL parm(*),work(*),x(*) C .. C .. Local Scalars .. REAL ae INTEGER i,icount,j,p C .. C .. External Functions .. REAL 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