CODEWITHSUNDEEP
c++cmatlabgsl
pyCthon
pyCthon

Reputation: 12341

matlab mvnrnd in gsl

Hi i'm not sure if my algorithm is correct i'm trying to replicate the mvnrnd function of matlab but in gsl. I found an algorithm in some journal articles that produces a vector of multivariate normal, but i need a matrix of multivariate normal random numbers

lets say the distribution is Z~(mu,sigma);

assuming sigma is a matrix already positive definite.

an algorithm i found off the web says to 

1. cholskey(sigma) = A
2. generate uniform gaussian vector r
3. matrix vector triangular product with gsl_blas_dtrmv A * r
4. add mu to Ar and that will be a vector of multivariate normal random numbers

my method below

are the following changes belowcorrect to product a Matrix of random variables

    1. cholskey(sigma) = A
    2. generate uniform gaussian matrix R
    3. matrix matrix scalar product AR
    4. add mu to AR and that will be a matrix of multivariate normal random numbers

Upvotes: 0

Views: 1259

Answers (1)

Dirk is no longer here
Dirk is no longer here

Reputation: 368251

Yes, that is correct. See e.g. this Wikipedia entry on multivariate normal RNGs which has this section:

Drawing values from the distribution

A widely used method for drawing a random vector x from the N-dimensional multivariate normal distribution with mean vector μ and covariance matrix Σ works as follows:

  1. Find any real matrix A such that A AT = Σ. When Σ is positive-definite, the Cholesky decomposition is typically used. [...]

  2. Let z = (z1, …, zN)T be a vector whose components are N independent standard normal variates (which can be generated, for example, by using the Box–Muller transform).

  3. Let x be μ + Az. This has the desired distribution due to the affine transformation property.

which describes the same algorithm.

There are also several implementations for R, for example mvrnorm in the MASS package which comes with every R installation.

Upvotes: 2

Related Questions