R - multiple nested loops inferno

Question

I want to create a list of matrices from 3 vectors, where sd_vec is strictly positive:

z_vec <- c(qnorm(0.90), qnorm(0.95), qnorm(0.975))  # 95% CI, 90% CI, 80% CI
me_vec <- c(0.50, 0.25, 0.10, 0.05)
sd_vec <- rnorm(n = 9, 0.8, 0.4)
sd_vec[which(sd_vec <= 0)] <- 0.1

My question is two-fold:

How can I generalize my loops from 1 matrix to the multiple matrices? I can create the matrices stepwise (2 levels of nesting), but my indexing breaks down when I nest all 3 levels.
How can I avoid for-loops altogether (in the spirit of this answer)? I welcome any answers with sapply(), etc.

Here is an example of my attempt:

new_n <- matrix(NA, 3,4)
for(i in seq_along(z_vec)){
  for(j in seq_along(me_vec)){
    new_n[i, j] <- ((z_vec[i] * sd_vec[1]) /me_vec[j])^2
  }
}

new_n
#      [,1]  [,2]  [,3] [,4]
# [1,] 2.45  9.82  61.4  245
# [2,] 4.04 16.17 101.1  404
# [3,] 5.74 22.96 143.5  574

Then my indexing failure:

new_n <- vector("list", length = length(sd_vec))
for(k in seq_along(sd_vec)){
  for(i in seq_along(z_vec)){
    for(j in seq_along(me_vec)){
      new_n[[k]][i, j] <- ((z_vec[i] * sd_vec[k]) /me_vec[j])^2
    }
  }
}

with the error message Error in new_n[[k]][i, j] <- ((z_vec[i] * sd_vec[k])/me_vec[j])^2 : incorrect number of subscripts on matrix

Thanks for any and all assistance with this trivial question!

RHertel · Accepted Answer · 2015-07-17T07:55:25.850

Concerning the error message in the last example, I think that this is due to an incorrect initialization of the list of matrices. You might try to replace

new_n <- vector("list", length = length(sd_vec))

with

new_n <- replicate(length(sd_vec), matrix(NA, 3, 4), simplify = FALSE)

This should resolve the indexing error problem. Finding an elegant and compact way to rewrite your complex nested loops is another issue. I trust that the SO community will soon come up with nice solutions.

update

The nested loops in your first example can be rewritten like this:

new_n <- t(sapply(seq_along(z_vec), function(x,y) ((z_vec[x] * sd_vec[1]) / me_vec[y])^2))

Thanks @RHertel - I didn't know about this way to create/initialize lists. — NotYourIPhone Siri, Jul 16 '15 at 21:21

score 0 · Answer 2 · edited May 23 '17 at 12:07

How to create a list of matrices is explained here, and to calculate each matrix without for-loops, we can use 'outer':

z_vec <- c(qnorm(0.90), qnorm(0.95), qnorm(0.975))  # 95% CI, 90% CI, 80% CI
me_vec <- c(0.50, 0.25, 0.10, 0.05)
sd_vec <- rnorm(n = 9, 0.8, 0.4)
sd_vec[which(sd_vec <= 0)] <- 0.1

new_n <- list()

for ( k in seq_along(sd_vec) )
{
  new_n[[k]] <- outer( z_vec*sd_vec[k], me_vec,
                       FUN = function(x,y){ (x/y)^2 } )
}

Using 'outer' twice, one can generate without any for-loops a 3-dimensional array containing the same matrices as 'new_n':

A <- outer( outer( sd_vec, z_vec, "*"),
            me_vec, function(x,y){ (x/y)^2 } )

Then 'new_n[[k]]' is identical to 'A[k,,]' for each 'k'. To get rid of the last 'for'-loop, we can use the function 'alply' from the library 'plyr' to turn the 3-dimensional array into a list of matrices:

library(plyr)

z_vec <- c(qnorm(0.90), qnorm(0.95), qnorm(0.975))  # 95% CI, 90% CI, 80% CI
me_vec <- c(0.50, 0.25, 0.10, 0.05)
sd_vec <- rnorm(n = 9, 0.8, 0.4)
sd_vec[which(sd_vec <= 0)] <- 0.1

A <- outer( outer( sd_vec, z_vec, "*"),
            me_vec, function(x,y){ (x/y)^2 } )

new_n <- alply(A,1)

R - multiple nested loops inferno

2 Answers2