循环次序对R的速度是否重要

Question

我有一个问题,我想进行模拟研究,模拟取决于两个变量x和y. x并且y是我想要在我的模拟研究中评估的潜在值的向量(如此不同的组合).此外,对于每个组合x和y我想多次重复(因为在那里一个随机期限和的每次运行x和y将变化).

举一个我正在处理的例子,我有以下简化的例子:

x = 1:10
y = 11:20

iterations = 2000
iter = 1

solution = array(NA,c(length(x),3,iterations))
for(i in x){
    for(j in y){
        for(k in 1:iterations){ 
            z = rnorm(3) + c(i,j,1) 
            solution[i,,k] = z
        }
    }
}

但是在我的实际问题中,在for循环中得到的代码要评估的重要性要小得多.但是,输入结构和输出结构相同.

所以我想知道,比如使用上面的例子,最有效的方法是按顺序设置循环,或者最好让它k in 1:iterations成为最外层的循环并尝试outer()在该循环中使用某种命令我将z在网格上评估函数(在这个例子中)x和y？

此外,我对完全不同的设置和设计非常开放.在一天结束时,我希望能够获得一个基于一个解决方案x,并y与平均超过所有的迭代,即apply(solution, c(1,2),mean)

---

编辑:

正如我所建议的,这是我正在使用的实际代码.

library(survival)

iter = 2000
n = 120
base = 2
hr = 0.5
m.x = 3
m.y = m.x/hr

ANS = NULL
for (vecX in c(0.3, 0.5, 0.6, 0.7)){
out = NULL

    for (vecY in c(0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95)){
        m.x.p = m.x/vecX
        m.y.p = m.y/vecX
        m.x.n = m.x
        m.y.n = m.y

        n.t = round(n*base/(base+1))
        n.c = n - n.t

        for (ii in 1:iter){
            n.t.p = rbinom(1, n.t, vecY)
            n.t.n = n.t - n.t.p
            n.c.p = rbinom(1, n.c, vecY)
            n.c.n = n.c - n.c.p

            S = c(rexp(n.t.p, log(2)/m.y.p), rexp(n.t.n, log(2)/m.y.n), rexp(n.c.p, log(2)/m.x.p), rexp(n.c.n, log(2)/m.x.n))

            data1 = data.frame(Group = c(rep("T", n.t), rep("C", n.c)), dx = c(rep("P", n.t.p), rep("N", n.t.n), rep("P", n.c.p), rep("N", n.c.n)), S)

            fit = survfit(Surv(data1$S)~data1$Group)
            coxfit = coxph(Surv(data1$S)~data1$Group)

            HR = exp(coxfit$coefficients)
            p.val=summary(coxfit)$logtest["pvalue"]

            out = rbind(out, c(vecX, vecY, n.t.p, n.t.n, n.c.p, n.c.n, HR, p.val))
        }

    }
    colnames(out) = c("vecX", "vecY", "n.t.p", "n.t.n", "n.c.p", "n.c.n", "HR", "p.val")

    ans = as.data.frame(out)
ANS = rbind(ANS, ans)

}

Answer 1

是的,我相信理论上它应该有所作为(见下面的例子).

R使用像Fortran这样的列主要排序(与C不同),因此为了最大限度地减少缓存未命中,您需要遍历列.因此,对于填充矩阵,最佳方法是外循环具有列索引的方法.

对于n维数组,您也要考虑到这一点.在这种情况下n = 3,我想这意味着让图层成为最外面的循环,然后是列,然后是行.不过,我可能会错在这里.

我5000用5000矩阵运行了这个简单的例子.我们看到差异大约50秒,fill_matrix2()而且速度更快.

    n <- 5000
    A <- matrix(NA, n, n)
    B <- matrix(NA, n, n)

    fill_matrix1 <- function(X, val) {
        for (i in 1:nrow(X)) {
            for (j in 1:ncol(X)) {
                X[i, j] <- val
            }
        }
        return(X)
    }

    fill_matrix2 <- function(X, val) {
        for (j in 1:ncol(X)) {
            for (i in 1:nrow(X)) {
                X[i, j] <- val
            }
        }
        return(X)
    }

    system.time(fill_matrix1(A, 0))
    system.time(fill_matrix2(B, 0))