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:
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
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 X∼N(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