R - 多个嵌套循环地狱答案

【问题标题】：R - multiple nested loops infernoR - 多个嵌套循环地狱
【发布时间】：2015-07-16 20:27:33
【问题描述】：

我想从 3 个向量创建一个矩阵列表，其中 sd_vec 是严格正数：

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

我的问题有两个：

如何将循环从 1 个矩阵推广到多个矩阵？我可以逐步创建矩阵（2 级嵌套），但是当我嵌套所有 3 级时，我的索引会崩溃。
如何才能完全避免 for 循环（本着 this answer 的精神）？欢迎使用sapply() 等回答任何问题。

这是我尝试的一个例子：

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

然后我的索引失败：

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

带有错误消息Error in new_n[[k]][i, j] <- ((z_vec[i] * sd_vec[k])/me_vec[j])^2 : incorrect number of subscripts on matrix

感谢您对这个琐碎问题的任何和所有帮助！

【问题讨论】：

标签： r for-loop matrix

【解决方案1】：

关于上一个示例中的错误消息，我认为这是由于矩阵列表的初始化不正确。你可以尝试替换

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

与

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

这应该可以解决索引错误问题。寻找一种优雅而紧凑的方式来重写复杂的嵌套循环是另一个问题。我相信 SO 社区很快就会提出很好的解决方案。

更新

第一个示例中的嵌套循环可以这样重写：

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

【讨论】：

感谢@RHertel - 我不知道这种创建/初始化列表的方式。

【解决方案2】：

here 解释了如何创建矩阵列表，并且要在没有 for 循环的情况下计算每个矩阵，我们可以使用 '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 } )
}

使用 'outer' 两次，无需任何 for 循环即可生成包含与 'new_n' 相同的矩阵的 3 维数组：

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

那么对于每个 'k'，'new_n[[k]]' 与 'A[k,,]' 相同。为了摆脱最后一个'for'循环，我们可以使用库'plyr'中的函数'alply'来转动3维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)

【讨论】：