# Generate MVN data that has a specific mean and variance matrix

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

• 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. – OGV Sep 28 '18 at 3:03
• 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. – Glen_b -Reinstate Monica Sep 28 '18 at 3:08
• 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! – 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))
t(Y)