Crossed for loops: Pick i th element of first loop then loop completely through second loop

Question

Suppose I have the following code:

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

I then draw 70 and 40 values from two multivariate normal distributions:

    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))

As can be easily seen S1 is 70x2 matrix und S2 a 40x2 matrix.

Now I build a for loop in R:

    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)
    }

This gives me a matrix that contains all combinations for the 40 elements in S2 with the 1st element of S1. What I want is to completely cross the two matrices S1 and S2. That is I want the for loop to pick out S1[1,] first, then iterate completely through S2[i,] (e.g. in an inner loop) and store the results in a matrix then pick out S1[2,] iterate again through S2[i,] and store the results in a matrix and so on. If I would need to give a name to what I am looking for I would say "crossed for loops". I find it incredibly hard to come up with R-code that will allow me to do this. Any hints would be appreciated.

Maybe the idea will get clearer with this example:

My idea is equivalent to construction 70 for-loops for every element in S1[i,] and binding the result in a 70*40*2x1 matrix:

    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)

What I ideally would want is to do this with for-loops or any other flexible way that can handle different dimensions for S1 and S2.

Your example isn't [complete or reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) so it is difficult to precisely answer. It seems like you could probably use `outer` here but it is unclear what you expect the result to be (given you code isn't runnable). — MrFlick, Jul 07 '14 at 16:43
Yes, I thought I mentioned this in the text: "store the results in a matrix". But I may have been unclear. Ideally the output would be a matrix. I thought about concatenating the matrices for every `S1[i,], i = 1, …, 70` after every iteration through `S2[j,]` to form one big matrix. — lord.garbage, Jul 07 '14 at 19:15
It should be a `(70*40*2)x2 = 5600x2`-matrix: `70 := number of elements in S1`, `40 := number of elements in S2`, `2 := number of rows in X`. — lord.garbage, Jul 07 '14 at 19:20
Yeah, sorry about me being so slow, the semifinals start tomorrow. — lord.garbage, Jul 07 '14 at 19:22

score 2 · Answer 1 · answered Jul 07 '14 at 20:18

OK. I might do a few things to make this as efficient as possible. First, we can pre-calculate all the matrix multiplication with

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

Then we can clculate all the combinations of the S1/S2 values with expand.grid

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

Then we can define the transformation

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

we assume we will be passed an index for S1 and an index for S2. Then we combine the data as in your formula. Finally, we use this helper function and the indexes with a set of do.calls

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

First we use Map to calculate all the combinations, then we use rbind to merge all the results. This actually produces a 2800x2 matrix. (70*40)*2. The rows are ordered with S1 moving the fastest, then S2.

lord.garbage · Accepted Answer · 2014-08-07T13:56:11.180

I realised that this was not a problem with for-loops but with the way I stored the variables. The solution to what I want is:

library(MASS)
z <- list()
y <- list()
for(j in 1:dim(S1)[1]){
    for(i in 1:dim(S2)[1]){
        z[[i]] <- X %*% beta+X %*% S1[j,]+X %*% S2[i,]+matrix(rnorm(2, mean = 0 , sd = sigma), ncol = 2, nrow = 2)
        Z <- unname(do.call(rbind, z))
    }
        y[[j]] <- Z
        Y <- unname(do.call(rbind, y))
}

Crossed for loops: Pick i th element of first loop then loop completely through second loop

2 Answers2