As an example, say I want to create a matrix $Y$ with 2 covariates and 10 observations each ( a 10 x 2 matrix) with a specific variance matrix and mean vector. The goal is to simulate regression data.

I want to generate $Y\sim MVN(\mu,\Sigma)$. As an example: $\Sigma =\begin{bmatrix}1 & .1 \\ .1 & 1\end{bmatrix} , $ $\mu=\begin{bmatrix}10 & 4\end{bmatrix}$

I start by generating $X\sim MVN(\vec{0},\mathbb{I})$ and my mean and variance matrix:

x <- matrix( rnorm(20), 10, 2)
mu <- matrix(c(10,4), nrow = 2)
sigma <- matrix( c(1,.1,.1, 1), nrow = 2, ncol = 2)

If I understand correctly what I need to do is let $Y = \mu +BX$ where $B=U\Lambda^{\frac{1}{2}}$ such that $U$ is the matrix of eigenvectors and $\Lambda$ is the diagonal matrix of eigen values of $\Sigma$. So I go to compute that:

eigen <- eigen(sigma)
u <- eigen$vectors
sqrtlambda <- sqrt(diag(eigen$values))
B = u %*% sqrtlambda

But this isn't quite what I want at all, as this has all of my dimensions wrong to get what I am looking for.

Final code missing:

mu <- matrix(rep(1,10),nrow = 10) %*% mu
Y = x %*% t(B) + mu
  • $\begingroup$ It was just a number I came up with off the top of my head, I was thinking something along the lines of two RV X's, each with 10 replicates. $\endgroup$ – OGV Sep 28 '18 at 3:03
  • $\begingroup$ My apologies. I misread your code. The problem is that when you write $BX$ in algebra $X$ is a column vector of length 2. When you created a matrix you made it with column dimension 10. You either need to transpose your x or you need to redo your algebra to work with row-vectors and then compute $XB^\prime$ instead. $\endgroup$ – Glen_b -Reinstate Monica Sep 28 '18 at 3:08
  • $\begingroup$ Ah thank you! The continual mix of block notation had me confused. I promise this is all on the premise of trying to learn the algebra! $\endgroup$ – OGV Sep 28 '18 at 3:19

A relevant piece of code is:

n <- 10
p <- 2
mu <- c(10, 4)
Sigma <- matrix(c(1, .1, .1, 1), nrow = 2, ncol = 2)

ed <- eigen(Sigma, symmetric = TRUE)
ev <- ed$values
evec <- ed$vectors
Y <- drop(mu) + tcrossprod(evec * rep(sqrt(pmax(ev, 0)), each = p), 
                           matrix(rnorm(n * p), n))
| cite | improve this answer | |

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.