Generating values from a multivariate Gaussian distribution

I am currently trying to simulate values of a N-dimensional random variable X that has a multivariate normal distribution with mean vector μ=(μ1,...,μN)T and covariance matrix S.

I am hoping to use a procedure similar to the inverse CDF method, meaning that I want to first generate a N-dimensional uniform random variable U and then plug that into the inverse CDF of this distribution, so to generate value X.

I am having issues because the procedure is not well documented and there are slight differences between the mvnrnd function in MATLAB and a description that I found on Wikipedia.

In my case, I am also choosing the parameters of the distribution randomly. In particular, I generate each of the means, μi, from a uniform distribution U(20,40). I then build the covariance matrix S using the following procedure:

  1. Create a lower triangular matrix L where L(i,i)=1 for i=1..N and L(i,j)=U(1,1) for i<j

  2. Let S=LLT where LT denotes the transpose of L.

This procedure allows me to ensure that S is symmetric and positive definite. It also provides a lower triangular matrix L so that S=LLT, which I believe is required to generate values from the distribution.

Using the guidelines on Wikipedia, I should be able to generate values of X using a N-dimensional uniform as follows:

  • X=μ+LΦ1(U)

According to the MATLAB function however, this is typically done as:

  • X=μ+LTΦ1(U)

Where Φ1 is the inverse CDF of a N-dimensional, separable, normal distribution, and the only difference between both methods is simply whether to use L or LT.

Is MATLAB or Wikipedia the way to go? Or are both wrong?

Answer

If XN(0,I) is a column vector of standard normal RV's, then if you set Y=LX, the covariance of Y is LLT.

I think the problem you're having may arise from the fact that matlab's mvnrnd function returns row vectors as samples, even if you specify the mean as a column vector. e.g.,

 > size(mvnrnd(ones(10,1),eye(10))  
 > ans =
 >      1    10

And note that transforming a row vector gives you the opposite formula. if X is a row vector, then Z=XLT is also a row vector, so ZT=LXT is a column vector, and the covariance of ZT can be written E[ZTZ]=LLT.

Based on what you wrote though, the Wikipedia formula is correct: if Φ1(U) were a row vector returned by matlab, you can't left-multiply it by LT. (But right-multiplying by LT would give you a sample with the same covariance of LLT).

Attribution
Source : Link , Question Author : Berk U. , Answer Author : jpillow

Leave a Comment