13
$\begingroup$

I am currently trying to simulate values of a $N$-dimensional random variable $X$ that has a multivariate normal distribution with mean vector $\mu = (\mu_1,...,\mu_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, $\mu_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 = LL^T$ where $L^T$ 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 = LL^T$, 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 = \mu + L * \Phi^{-1}(U)$

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

  • $X = \mu + L^T * \Phi^{-1}(U)$

Where $\Phi^{-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 $L^T$.

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

$\endgroup$
  • 2
    $\begingroup$ As stated, both are wrong because $\mu$ is a row vector while $T * invnorm(U)$ must be a column vector. When you get your rows and columns straightened out, this question should answer itself simply by identifying which version of $(X-\mu)'(X-\mu)$ or $(X-\mu)(X-\mu)'$ gives a matrix and which version gives just a number: check that you can compute the expectation of the matrix version and that it gives $S$. $\endgroup$ – whuber Jul 12 '11 at 22:14
  • $\begingroup$ @whuber Yeap. Made the changes in the formatting for the question. Thanks for the tip - definitely the easiest way to check. $\endgroup$ – Berk U. Jul 13 '11 at 16:32
14
$\begingroup$

If $X \sim \mathcal{N}(0,I)$ is a column vector of standard normal RV's, then if you set $Y = L X$, the covariance of $Y$ is $L L^T$.

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 = X L^T$ is also a row vector, so $Z^T = L X^T$ is a column vector, and the covariance of $Z^T$ can be written $E[Z^TZ] = LL^T$.

Based on what you wrote though, the Wikipedia formula is correct: if $\Phi^{-1}(U)$ were a row vector returned by matlab, you can't left-multiply it by $L^T$. (But right-multiplying by $L^T$ would give you a sample with the same covariance of $LL^T$).

$\endgroup$
  • $\begingroup$ Note that the help for mvnrnd in matlab uses $N$ as the number of samples; the number of dimensions is $D$. So if you ask for $N$ samples from a $D$-dimensional multivariate normal, it returns them as an $N \times D$ matrix. $\endgroup$ – jpillow Jul 13 '11 at 4:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.