r - 交叉for循环: Pick i th element of first loop then loop completely through second loop

Question

假设我有以下代码:

    X    <- model.matrix(~factor(1:2))
    beta <- c(1, 2)

然后我从两个多元正态分布中得出 70 和 40 值:

    library(MASS)
    S1 <- mvrnorm(70, mu = c(0,0), Sigma = matrix(c(10, 3, 3, 2), ncol = 2))
    S2 <- mvrnorm(40, mu = c(0,0), Sigma = matrix(c(10, 4, 4, 2), ncol = 2))

很容易看出S1是 70x2矩阵和S2一个40x2矩阵。

现在我在R中构建了一个for循环:

    z <- list()
    for(i in 1:dim(S2)[1]){
        z[[i]] <- X %*% beta + X %*% S1[1,] + X %*% S2[i,] + rnorm(2, mean = 0, sd = 0.45)
    Y <- do.call(rbind, z)
    }

这给了我一个矩阵，其中包含 40 的所有组合S2 中的元素与 1 S1 的 st 元素。我想要的是完全交叉两个矩阵S1和S2 。那就是我想要for循环挑选出S1[1,]首先，然后完全迭代 S2[i,] (例如在内部循环中)并将结果存储在矩阵中，然后挑选 S1[2,]再次迭代S2[i,]并将结果存储在矩阵中等等。如果我需要为我正在寻找的东西命名，我会说 "crossed for loops" 。我发现很难想出 R -允许我执行此操作的代码。任何提示将不胜感激。

---

也许这个想法会通过这个例子变得更清晰:

我的想法相当于构造70 for -循环S1[i,]中的每个元素并将结果绑定(bind)在 70*40*2x1 中矩阵:

    for(i in 1:dim(S2)[1]){
    z[[i]] <- X %*% beta+X %*% S1[1,]+X %*% S2[i,]+rnorm(2, mean = 0 , sd = sigma)
    Y1 <- unname(do.call(rbind, z))
    }

    for(i in 1:dim(S2)[1]){
    z[[i]] <- X %*% beta+X %*% S1[2,]+X %*% S2[i,]+rnorm(2, mean = 0 , sd = sigma)
    Y2 <- unname(do.call(rbind, z))
    }

    for(i in 1:dim(S2)[1]){
    z[[i]] <- X %*% beta+X %*% S1[3,]+X %*% S2[i,]+rnorm(2, mean = 0 , sd = sigma)
    Y3 <- unname(do.call(rbind, z))
    }

    .
    .
    .

    for(i in 1:dim(S2)[1]){
    z[[i]] <- X %*% beta+X %*% S1[70,]+X %*% S2[i,]+rnorm(2, mean = 0 , sd = sigma)
    Y70 <- unname(do.call(rbind, z))
    }

    Y <- rbind(Y1, Y2, Y3, …, Y70)

我理想的情况是用 for 来做到这一点-循环或任何其他可以处理 S1 不同维度的灵活方式和S2 .

Answer 1

好的。我可能会做一些事情来使其尽可能高效。首先，我们可以预先计算所有矩阵乘法

Xb <- X %*% beta
XS1 <- X %*% t(S1)
XS2 <- X %*% t(S2)

然后我们可以使用expand.grid计算S1/S2值的所有组合

idx <- unname(c(expand.grid(A=1:ncol(XS1), B=1:ncol(XS2))))

然后我们可以定义转换

fx<-function(a, b) {
    t(Xb + XS1[,a, drop=F] + XS2[,b,drop=F] + rnorm(2, mean = 0, sd = 0.45))
}

我们假设我们将传递 S1 的索引和 S2 的索引。然后我们按照您的公式合并数据。最后，我们使用这个辅助函数和带有一组 do.calls

的索引

xx <- do.call(rbind, do.call(Map,c(list(fx), idx)))

首先我们使用Map计算所有组合，然后使用rbind合并所有结果。这实际上产生了一个 2800x2 矩阵。 (70*40)*2。这些行按 S1 移动最快的顺序排列，然后是 S2。