最近我一直在研究电力模拟,以下是我的代码:
library(MASS)
library(Matrix)
simdat <- data.frame(mmm = rep(rep(factor(1:2,
labels=c("m1", "m2")),
each = 2),
times = 2800),
ttt = rep(factor(1:2,
labels = c("t1", "t2")),
times = 5600),
sss = rep(factor(1:70),
each = 160),
iii = rep(rep(factor(1:40),
each = 4),
times = 70))
beta <- c(1, 2)
X1 <- model.matrix(~ mmm,
data = simdat)
Z1 <- model.matrix(~ ttt,
data = simdat)
X1
和Z1
是11200x2
矩阵。在Stackoverflow的帮助下,我成功地使我的计算比以前更加高效:
funab <- function(){
ran_sub <- mvrnorm(70, mu = c(0,0), Sigma = matrix(c(10, 3, 3, 2), ncol = 2))
ran_ite <- mvrnorm(40, mu = c(0,0), Sigma = matrix(c(10, 3, 3, 2), ncol = 2))
Mb <- as.vector(X1 %*% beta)
M1 <- rowSums(Z1 * ran_sub[rep(1:70,
each = 160),])
M2 <- rowSums(Z1 * ran_ite[rep(rep(1:40, each = 4),
times = 70),])
Mout <- Mb + M1 + M2
Y <- as.vector(Mout) + rnorm(length(Mout), mean = 0 , sd = 0.27)
}
Y
将会是一个长度为 11200
的向量。随后我会大量地复制这个函数 (比如说 1000
次):
sim <- replicate(n = 1000,
expr = funab()},
simplify = FALSE)
sim
将是一个11200x1000
的列表。考虑到我想要频繁地这样做,可能还要在funab()
中加入更多的代码,我想知道在当前的计算中是否建议对X1
和Z1
使用稀疏矩阵?
microbenchmark
吗?你可以使用它来比较函数的性能,也就是基准测试。只需简单地执行install.packages(c("microbenchmark"), dependencies = TRUE)
,require(microbenchmark)
和example(microbenchmark)
,你就知道该怎么做了。我在 这个 SO 回答 中使用了microbenchmark
。 - Eric Fail