Double For Loop 计算平均值并将其存储在矩阵中答案

【问题标题】：Double For Loop to calculate averages and store them in matrixDouble For Loop 计算平均值并将其存储在矩阵中
【发布时间】：2015-04-23 14:10:32
【问题描述】：

我在运行这个双 for 循环以将计算值正确存储到矩阵中时遇到问题（如下所述）。我选择使用双 For 循环而不是 apply() 或 mean() 的原因是我想获得两列的唯一组合并消除冗余（解释如下）。请参阅下面的示例：

A<-c(1,2,3,4,5)
B<-c(2,3,4,5,6)
Q1<-data.frame(cbind(A,B))
mean<-matrix(nrow=5, ncol = 5)
for(i in 1: length(Q1$A)){
  for(j in 2: length(Q1$B)){
    mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
  }
}

在这里，我尝试在消除冗余的同时将整个 A 向量通过整个 B 向量运行，使得 A[1] 具有来自 B[2] 的四个值，而 A[2] 具有来自 B[3] 的三个值.然而，这是我的结果。

     [,1] [,2] [,3] [,4] [,5]
[1,]   NA  2.0  2.5  3.0  3.5
[2,]   NA  2.5  3.0  3.5  4.0
[3,]   NA  3.0  3.5  4.0  4.5
[4,]   NA  3.5  4.0  4.5  5.0
[5,]   NA  4.0  4.5  5.0  5.5

虽然第一列是我所期望的，但我有我不想要的值。我想要的是下面的矩阵输出：

     [,1] [,2] [,3] [,4] [,5]
[1,]   NA  2.0  2.5  3.0  3.5
[2,]   NA   NA  3.0  3.5  4.0
[3,]   NA   NA   NA  4.0  4.5
[4,]   NA   NA   NA   NA  5.0
[5,]   NA   NA   NA   NA   NA

有什么建议吗？

【问题讨论】：

为什么你只对矩阵的一半感兴趣？例如，在colA <- 1:3 和colB <- 13:11 的情况下，输出矩阵会变得不对称（例如A[1] + B[3] != A[3] + B[1]），如果只查看矩阵的一半，就会丢失信息。
@MaratTalipov 我对一半感兴趣，因为我想获取这些值并将它们放入一列中，以便我可以将其与 ggplot 中的其他值进行比较。如果有冗余，那么它将反映图表的结果。

标签： r for-loop mean

【解决方案1】：

第二个for循环应该是：

 for(j in (i+1):length(Q1$B))

【讨论】：

【解决方案2】：

您想使用next 关键字来跳过您不需要的操作，如：

A<-c(1,2,3,4,5)
B<-c(2,3,4,5,6)
Q1<-data.frame(cbind(A,B))
mean<-matrix(nrow=5, ncol = 5)
for(i in 1: length(Q1$A))
for(j in 2: length(Q1$B)){
    if(i >= j)
        next
    mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
}

或者您可以使内部for 循环的迭代以外部索引的值为条件，如下所示：

mean<-matrix(nrow=5, ncol = 5)
for(i in 2: length(Q1$A)){
    for(j in i: length(Q1$B)){
        mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
    }
}

或者你可以使用outer()：

mean<-outer(1: length(Q1$A), 
            1: length(Q1$B),
            Vectorize(function(i,j){
                if(i >= j)
                    return(NA)
                else 
                    return(sum(Q1$A[i]+Q1$B[j])/2)
            }))

【讨论】：

【解决方案3】：

[原始解决方案（有关更快的解决方案，请参阅更新 2）]

f.m <- function(Q1) {
    z <- matrix(nrow=nrow(Q1),ncol=nrow(Q1))
    b <- row(z) < col(z)
    z[b] <- (Q1$A[col(z)[b]] + Q1$B[row(z)[b]])/2
    z
}

[样本输出]

f.m(Q1)
#      [,1] [,2] [,3] [,4] [,5]
# [1,]   NA    2  2.5  3.0  3.5
# [2,]   NA   NA  3.0  3.5  4.0
# [3,]   NA   NA   NA  4.0  4.5
# [4,]   NA   NA   NA   NA  5.0
# [5,]   NA   NA   NA   NA   NA

[基准设置]

f0 <- function(Q1) {
    mean<-matrix(nrow=nrow(Q1), ncol = nrow(Q1))
    for(i in 1: length(Q1$A)){
        for(j in 2: length(Q1$B)){
            mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
        }
    }
    mean
}

f1 <- function(Q1) {
    mean<-matrix(nrow=nrow(Q1), ncol = nrow(Q1))
    for(i in 2: length(Q1$A)){
        for(j in i: length(Q1$B)){
            mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
        }
    }
    mean
} 

# Note that f0() and f1() don't return the desired result for the sample output

f2 <- function(Q1) {
    mean<-outer(1: length(Q1$A), 
                1: length(Q1$B),
                Vectorize(function(i,j){
                    if(i >= j)
                        return(NA)
                    else 
                        return(sum(Q1$A[i]+Q1$B[j])/2)
                }))
    mean
}

library(rbenchmark)

[基准测试结果]

A <- B <- 1:100
Q1<-data.frame(A,B)

benchmark(f0(Q1), f1(Q1), f2(Q1), f.m(Q1), replications = 10)
     test replications elapsed relative user.self sys.self user.child sys.child
4 f.m(Q1)           10   0.011    1.000     0.012    0.000          0         0
1  f0(Q1)           10   3.018  274.364     3.007    0.008          0         0
2  f1(Q1)           10   1.477  134.273     1.474    0.003          0         0
3  f2(Q1)           10   1.777  161.545     1.774    0.002          0         0

[更新 1]

通过直接计算整个矩阵可以节省另一个运行时间顺序，这避免了代价高昂（与求和相比）的子集：

f.m2 <- function(Q1) outer(Q1$A,Q1$B,'+')*0.5

基准测试的另一部分：

A <- B <- 1:1000
Q1<-data.frame(A,B)
#benchmark(f0(Q1), f1(Q1), f2(Q1), f.m(Q1), replications = 10)
benchmark(f.m(Q1), f.m2(Q1), replications = 10)

      test replications elapsed relative user.self sys.self user.child sys.child
1  f.m(Q1)           10   1.839   10.274     1.746    0.093          0         0
2 f.m2(Q1)           10   0.179    1.000     0.144    0.035          0         0

[更新 2]

1) 正如 David Arenburg 所指出的，函数 f.m2() 不会产生完全预期的输出，因为输出的左下三角形和主对角线应该用 NA 填充。可以修复函数 f.m2() 以产生正确的答案，但会以性能为代价（请参阅下面的基准测试）。

# Suggested by David Arenburg
f.m2.1 <- function(Q1) { 
   Res <- outer(Q1$A,Q1$B,'+')*0.5; 
   Res[lower.tri(Res, diag = TRUE)] <- NA; 
   Res 
}

2) 这是 David Arenburg 建议的另一种方法，它利用了 data.table 包中的 CJ 函数：

library(data.table)
f.DA <- function(Q1){ 
  Res <- matrix(rowMeans(CJ(Q1$A, Q1$B)), ncol = nrow(Q1))
  Res[lower.tri(Res, diag = TRUE)] <- NA
  Res 
}

3) 这是一个基于Rcpp 的方法：

library(Rcpp)
cppFunction('NumericMatrix fC(NumericVector A, NumericVector B) {

  int n = A.size();
  NumericMatrix out(n,n);
  std::fill( out.begin(), out.end(), NumericVector::get_na() ) ;

  for (int i = 0; i < n; i++) {
    for (int j = i+1; j < n; j++) {
      out(i,j) = 0.5*(A[i] + B[j]);
    }
  }
  return out;
}')

4) 另一个基准测试研究：

A <- B <- 1:3000
Q1<-data.frame(A,B)
benchmark(f.m2(Q1), f.m2.1(Q1), f.DA(Q1), fC(Q1$A, Q1$B), replications = 10)

            test replications elapsed relative user.self sys.self user.child sys.child
3       f.DA(Q1)           10   7.442   11.556     6.200    1.209          0         0
2     f.m2.1(Q1)           10   5.111    7.936     4.404    0.661          0         0
1       f.m2(Q1)           10   1.007    1.564     0.733    0.263          0         0
4 fC(Q1$A, Q1$B)           10   0.644    1.000     0.525    0.116          0         0

【讨论】：

感谢马拉特的帮助！
实际上f.m2 返回整个矩阵，因此没有所需的结果。您可能应该将其修改为 f.m2 <- function(Q1) { Res <- outer(Q1$A,Q1$B,'+')*0.5; Res[lower.tri(Res, diag = TRUE)] <- NA; Res } 以满足要求，但它仍然会比 f.m 快
您也可以添加library(data.table); f.DA <- function(Q1){ Res <- matrix(rowMeans(CJ(Q1$A, Q1$B)), ncol = nrow(Q1)); Res[lower.tri(Res, diag = TRUE)] <- NA; Res }，这也将击败f.m，但可能仍会比f.m2慢
感谢@DavidArenburg，我已经添加了您的建议以及基于Rcpp 的解决方案

【解决方案4】：

不完全是双 For 循环，但您可以使用 outer 函数来计算平均值。

outer(Q1$Col1, Q1$Col2, "+")/2

【讨论】：