R:在生成双变量正态样本时,'mu' 或 'Sigma' 可以在 MASS::mvrnorm() 中向量化吗?
R: can 'mu' or 'Sigma' be vectorized in MASS::mvrnorm() when generating bivariate normal samples?
考虑以下数据,
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
是我的方差向量。
我想生成一个双变量向量 v=(v1,v2)
,长度 n
,具有正态双变量分布。 v
的第 i 行具有均值 mu=(0,0)
、相关性 rho=0.4
和边际方差 sigma=(sigma1, 1)
的双变量正态分布,其中 sigma1
接收第 i 行的值 sigma1
。我该如何继续?
编辑澄清
在 1D 情况下,rnorm
接受矢量化 mu
和 sigma
,因此
rnorm(3, 0, sqrt(sigma1))
给出了来自 N(0, sqrt(sigma1[i]))
的样本。基本上 OP 要求 mvrnorm
.
具有相同的功能
基本解决方案
不,它不能被矢量化。写一个 for
循环。
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)
}
高级解决方案
给定一个协方差 Sigma
sig11 ^ 2 rho * sig11 * sig22
rho * sig11 * sig22 sig22 ^ 2
它的下三角乔列斯基因子L
是
sig11 0
rho * sig22 sqrt(1 - rho ^ 2) * sig22
如果x <- rnorm(2)
,则mu + L %*% x
是来自N(mu, Sigma)
的样本。
这给出了一个完全矢量化的解决方案
# 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)
}
对于您的数据,您可以使用
v <- birnorm(3, 0, 0, sqrt(sigma1), 1, 0.4)
为双变量情况重建mvrnorm
请注意,此函数可用于为双变量情况重建 mvrnorm
。
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)
}
而且速度更快
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))
考虑以下数据,
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
是我的方差向量。
我想生成一个双变量向量 v=(v1,v2)
,长度 n
,具有正态双变量分布。 v
的第 i 行具有均值 mu=(0,0)
、相关性 rho=0.4
和边际方差 sigma=(sigma1, 1)
的双变量正态分布,其中 sigma1
接收第 i 行的值 sigma1
。我该如何继续?
编辑澄清
在 1D 情况下,rnorm
接受矢量化 mu
和 sigma
,因此
rnorm(3, 0, sqrt(sigma1))
给出了来自 N(0, sqrt(sigma1[i]))
的样本。基本上 OP 要求 mvrnorm
.
基本解决方案
不,它不能被矢量化。写一个 for
循环。
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)
}
高级解决方案
给定一个协方差 Sigma
sig11 ^ 2 rho * sig11 * sig22
rho * sig11 * sig22 sig22 ^ 2
它的下三角乔列斯基因子L
是
sig11 0
rho * sig22 sqrt(1 - rho ^ 2) * sig22
如果x <- rnorm(2)
,则mu + L %*% x
是来自N(mu, Sigma)
的样本。
这给出了一个完全矢量化的解决方案
# 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)
}
对于您的数据,您可以使用
v <- birnorm(3, 0, 0, sqrt(sigma1), 1, 0.4)
为双变量情况重建mvrnorm
请注意,此函数可用于为双变量情况重建 mvrnorm
。
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)
}
而且速度更快
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))