0

我知道向量化以及如何应用它来加速 R 中的 for 循环，但我无法找到一种使用向量加速代码的方法，其中每次迭代都依赖于前一次迭代的结果，或者依赖于迭代随机区间计算。

例如：

乔什——对此感到抱歉。所以，这里有更多细节：

m <- c(1, 1)
w.r <- c(0.33592935393, 0.63825353030, 0.15335253356 )

rlistl 是 3 个 2x2 矩阵的列表。所以，为了谈话，

r0 <- matrix(0, 2, 2) 
r1 <- matrix(1, 2, 2)
r2 <- matrix(2, 2, 2)
rlist <- list(r0, r1, r2)

N <- 500
E <- matrix(0, N, 2)

for(i in 1:N) {
    r <- c(c(1:3) %*% rmultinom(1, 1, w.r))
    E[i, ] <- mvrnorm(1, m, rlist[[r]])
}

我已经尝试在循环之外进行“r <- multinom()”计算，并且 rprof 显示花费的大部分时间显然是在 mvnorm 中。任何人都可以在 R 中找到一种使用向量来加快速度的方法吗？

这是另一个例子

for(i in 1:N) {
    if(d$V[i, 1] & d$V[i, 2]) QQ <- 1
    else if(! d$V[i, 1] & d$V[i, 2]) QQ <- 2
    else if(! d$V[i, 1] & ! d$V[i, 2]) QQ <- 3
    else if(d$V[i, 1] & ! d$V[i, 2]) QQ <- 4

    U[i, ] <- r1bvtruncnorm(mux=mu.U[i, ]/sd.r[r1], rho=rho, q=QQ)

}

想不出办法让它运行得更快。我的部分问题是我是一名 C/C++ 程序员，但我一直在尝试阅读 R 并希望确保我不会遗漏一些简单的东西。

谢谢。

编辑：

贾斯汀：

好的——所以我尝试了你的建议，但正如我所担心的那样，rep() 的行为并不像我希望的那样。我每次都需要一个单独的随机数，但使用 rep() 只需调用 rmultinom 一次并将结果复制 100 次。

>rep(c(c(1:3) %*% rmultinom(1, 1, ww.r)), 100)
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

> rep(c(c(1:3) %*% rmultinom(1, 1, ww.r)), 100)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
4

1 回答 1

0

对于我猜错的任何内容，请参阅 Josh 的评论。

但是假设您没有遗漏任何其他内容，w.r这是恒定的，应该从循环中删除。 m也是常数。'r' 和 rlist[[r]] 是随机生成的，但具有恒定参数。所以你正在做的是N从相同的参数生成分布时间，对吗？如果是这样的话，有很多比 for 循环更好的选择。例如

E <- matrix(rep(mvrnorm(1, 
                       m, 
                       rlist[[c(1:3 %*% rmultinom(1, 1, w.r))]]), 
                N),  
            N, 
            2)

这应该给你与你的第一个 for 循环相同的结果。

您的第二个示例可以成为一个函数并使用apply. 然而，如果不知道d对象的实际结构，就很难知道如何最好地接近它。我假设它类似于 data.frames 的 data.frame 或列表的命名列表。d$V一些二维数据结构在哪里：

my.fun <- function(vec) {
    if(vec[1] & vec[2]) return(1)
    else if(! vec[1] & vec[2]) return(2)
    else if(! vec[1] & ! vec[2]) return(3)
    else if(vec[1] & ! vec[2]) return(4)
}

QQ <- apply(d$V, 1, my.fun)

根据评论编辑：

r  <- 1:3 %*% rmultinom(N, 1, w.r)
E <- t(sapply(r, function(x) mvrnorm(1, m, rlist[[x]])))

或者在一条凌乱的线上：

t(sapply(1:3 %*% rmultinom(N, 1, w.r), function(x) mvrnorm(1, m, rlist[[x]])))
于 2012-05-18T16:39:30.377 回答