random - Drawing from a multivarite normal, with "genmn" and "setgmn"

Question

I am trying to draw random deviates from a multinormal distribution. There are two subroutines I found online: "setgmn" and "genmn", which are supposed to do exactly that. They are from the "ranlib" library.

http://www.netlib.org/random/

However, I have no clue how those to subroutines are supossed to work together. Setgmn combines the data in such a way that genmn can uses it. However I struggle wih the "parm" argument I am not exactly sure what I have to pass into it. And how those two subroutines need to be combined to actually work.

I could not find an example online. So I thought maybe someone here already used them or has a link to an example online.

Answer 1

From the source:

First you have to call setgmn. This procedure will set the argument parm that you have to pass to getmn.

subroutine setgmn ( meanv, covm, p, parm )
  integer ( kind = 4 ) p
  real ( kind = 4 ) covm(p,p)
  real ( kind = 4 ) meanv(p)
  real ( kind = 4 ) parm(p*(p+3)/2+1)

Parameters:

MEANV: the means of the multivariate normal distribution.

COVM: On input, the covariance matrix of the multivariate distribution. On output, the information in COVM has been overwritten.

P: the number of dimensions.

PARM: parameters needed to generate multivariate normal deviates.

PARM(1) contains the size of the deviates, P

PARM(2:P+1) contains the mean vector.

PARM(P+2:P*(P+3)/2+1) contains the upper half of the Cholesky decomposition of the covariance matrix.

---

subroutine genmn ( parm, x, work )
  real ( kind = 4 ) parm(*)
  real ( kind = 4 ) work(*)
  real ( kind = 4 ) x(*)

Parameters:

PARM(P*(P+3)/2+1) : parameters set by SETGMN.

X(P): a random deviate from the distribution. Will be the output.

WORK(P): workspace