Double For Loop 计算平均值并将它们存储在矩阵中

Double For Loop to calculate averages and store them in matrix

我在运行使用这个双 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

有什么建议吗?

第二个 for 循环应该是:

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

您想使用 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 循环的 iterand 以外部索引的值为条件,如:

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

[原始解决方案(请参阅更新 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

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

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