Q:

R: can 'mu' or 'Sigma' be vectorized in MASS::mvrnorm() when generating bivariate normal samples?

Consider the data below,
n <- 3
phi0 <- 1
phi1 <- 0.1
phi2 <- 0.2
W2 <- runif(n,0,1)
W3 <- runif(n,0,1)
W <- cbind(W2,W3)
phi <- rbind(phi0,phi1,phi2)
rho <- 0.4
sigma1 <- exp(as.numeric(model.matrix(~W) %*% phi))
sigma2<- 1

library(MASS)
#Sigma <- ???
mu <- rep(0,2)
v <- mvrnorm(n, mu, Sigma)

sigma1 is my vector of variance.
I want to generate a bivariate vector v=(v1,v2) of length n with normal bivariate distribution. In such a way that the i-th line of v has a bivariate normal distribution of mean mu=(0,0), correlation rho=0.4 and marginal variance sigma=(sigma1, 1) where sigma1 receives the value of the i-th line sigma1. How can I proceed?

Editor's clarification
In 1D case, rnorm accepts vectorized mu and sigma, so that
rnorm(3, 0, sqrt(sigma1))

gives samples from N(0, sqrt(sigma1[i])). Basically OP is asking the same capability for mvrnorm.

A:

basic solution
No, it can't be vectorized. Write a for loop.
v <- matrix(0, n, 2)
for (i in 1:n) {
  sig11 <- sigma1[i]
  sig21 <- rho * sqrt(sig11)
  Sigma <- matrix(c(sig11, sig21, sig21, 1), 2)
  v[i, ] <- mvrnorm(1, c(0,0), Sigma)
  }

advanced solution
Given a covariance Sigma
sig11 ^ 2                rho * sig11 * sig22
rho * sig11 * sig22      sig22 ^ 2

its lower triangular Cholesky factor L is
sig11          0
rho * sig22    sqrt(1 - rho ^ 2) * sig22

If x <- rnorm(2), then mu + L %*% x is a sample from N(mu, Sigma).
This gives a fully vectorized solution
# mu1, mu2, sig11, sig22, rho can be
## length-n vectors
## scalars
## vectors than can be recycled to be length-n
birnorm <- function (n, mu1, mu2, sig11, sig22, rho) {
  x1 <- rnorm(n)
  x2 <- rnorm(n)
  z1 <- sig11 * x1 + mu1
  z2 <- rho * sig22 * x1 + sqrt(1 - rho ^ 2) * sig22 * x2 + mu2
  cbind(z1, z2)
  }

For your data, you can use
v <- birnorm(3, 0, 0, sqrt(sigma1), 1, 0.4)

rebuild mvrnorm for bivariate case
Note that this function can be used to rebuild mvrnorm for bivariate case.
mvrnorm2 <- function (n, mu, Sigma) {
  sig11 <- sqrt(Sigma[1])
  sig22 <- sqrt(Sigma[4])
  rho <- Sigma[2] / (sig11 * sig22)
  birnorm(n, mu[1], mu[2], sig11, sig22, rho)
  }

and it is faster
mu <- c(0,0)
Sigma <- matrix(c(1,0.5,0.5,1),2)
library(microbenchmark)
microbenchmark(mvrnorm(1000, mu, Sigma), mvrnorm2(1000, mu, Sigma))

