R：如何在此函数中避免 R 中的 2 个“for”循环答案

【问题标题】：R: How to avoid 2 'for' loops in R in this functionR：如何在此函数中避免 R 中的 2 个“for”循环
【发布时间】：2015-02-26 20:11:06
【问题描述】：

我知道有很多关于如何避免 R 循环的主题，但我无法理解如何向量化我的迭代。我有一个数据集，在这里我用 m 表示。我想用这个函数生成一个新的矩阵，它将由每列数据 (m) 的相关系数的 p.values 组成。

m<-matrix(rnorm(100),nrow=10,ncol=10)
sig.p<-function(x){
  n= ncol(x)
  p.values<-numeric(n)
  p.values<-matrix(nrow=n,ncol=n)
  for (i in 1:C){
    for (t in 1:C){
      p.values[t,i]<-cor.test(x[,i],x[,t])$p.value
    }
  }
  p.values
}
sig.p(m)

我无法理解如何使用 mapply（如果是这样的话）。任何人都可以提供有关如何矢量化这些迭代的建议（使用 mapply 或其他）提前致谢！

塞萨尔

【问题讨论】：

您的代码定义的矩阵 m 只有 1 列。应该是matrix(rnorm(100),nrow=10,ncol=10)？
对不起，你是对的。矩阵(rnorm(100),nrow=10,ncol=10)
您可能有兴趣使用system.time 分析一些答案，看看矢量化是否会使代码性能更好。
@akrun 你说得对。我正在寻找矢量化解决方案。但是，由于我也将其作为练习（我对 R 有点陌生）来熟悉某些函数的属性（例如矢量化，我不知道），所以我正在给自己重复它的工作几次。再次感谢！
@Drizzt 我在100*100 上尝试了基准测试。在那个中，rcorr 是最好的。如果你能在更大的数据集上显示时间，那就太好了。

标签： r for-loop iteration

【解决方案1】：

您可以使用 rcorr 中的 library(Hmisc)

 library(Hmisc)
 rcorr(m)$P

或者使用

 library(psych)
 corr.test(as.data.frame(m))$p

或者使用来自base R的outer

  outer(1:ncol(m),1:ncol(m), FUN= Vectorize(function(x,y) 
                              cor.test(m[,x], m[,y])$p.value))

基准测试

我尝试了一个较小的数据集 (100*100) 和一个稍大的数据集 (1e3*1e3)。以下是函数：

 akrun <- function() {outer(1:ncol(m1),1:ncol(m1), 
            FUN= Vectorize(function(x,y) cor.test(m1[,x],
                              m1[,y])$p.value))}

 akrun2 <- function(){rcorr(m1)$P}
 agstudy <- function() {M <- expand.grid(seq_len(ncol(m1)),
                          seq_len(ncol(m1)))
      mapply(function(x,y)cor.test(m1[,x], m1[,y])$p.value,M$Var1,M$Var2)}
 vpipk <-function(){
        n <- ncol(m1)
        p.values<-matrix(nrow=n,ncol=n)
   for (i in 1:(n-1)){
      for (t in (i+1):n){
          p.values[t,i]<-cor.test(m1[,i],m1[,t])$p.value
     }
   }
   p.values
  }


 nrussell <- function(){
   sapply(1:ncol(m1), function(z){
   sapply(1:ncol(m1), function(x,Y=z){
     cor.test(m1[,Y],m1[,x])$p.value
     })
   })
}

在100*100 数据集上

 library(microbenchmark)
 set.seed(25)
 m1 <- matrix(rnorm(1e2*1e2),nrow=1e2,ncol=1e2)
 microbenchmark(akrun(), akrun2(), agstudy(), vpipk(),
                    nrussell(), unit='relative', times=10L)
 #Unit: relative
 #  expr      min       lq     mean   median       uq      max neval cld
 #   akrun() 257.2310 255.9766 252.2163 254.4946 248.9807 246.5429    10   c
 #  akrun2()   1.0000   1.0000   1.0000   1.0000   1.0000   1.0000    10   a  
 # agstudy() 255.5920 258.0813 253.5411 256.0581 250.4833 249.0503    10   c
 #   vpipk() 125.8218 126.3337 125.4592 126.8479 124.9835 124.1383    10   b 
 #nrussell() 257.9283 256.8480 252.5297 256.0160 250.8853 242.0896    10   c

如果我将1e2 更改为1e3（没时间做microbenchmark，但这里是system.time

system.time(akrun())
 # user  system elapsed 
#403.563   0.751 404.198 

system.time(akrun2())
 #  user  system elapsed 
 # 3.110   0.008   3.117 

system.time(agstudy())
 #  user  system elapsed 
 #445.108   0.877 445.947 

system.time(vpipk())
#  user  system elapsed 
#155.597   0.224 155.760 

system.time(nrussell())
#  user  system elapsed 
#452.524   1.220 453.713

【讨论】：

使用 outer 和 Vectorize 的有趣模式。
akrun，非常感谢。这就是我一直在寻找的。将其与 vpipkt 进行比较，因为我想对 10k 迭代使用相同的逻辑

【解决方案2】：

这是mapply的典型用法：

M <- expand.grid(seq_len(ncol(m),seq_len(ncol(m)))
mapply(function(x,y)cor.test(m[,x], m[,y])$p.value,M$Var1,M$Var2)

【讨论】：

是seq_len() ncol(m) 中的对象吗？
我尝试按照您所说的方式执行mapply，但它返回给我一个长度为 M 且全是零的向量。它只得到对角线 (x=y)。

【解决方案3】：

矢量化并不总是像人们想象的那样。不确定您的实际矩阵有多大，但对于这个大小甚至 100 x 100 来说，它是相当小的一次性成本。

通过如下修改循环结构，您可以将性能提高一倍以上：

sig.p<-function(x){
    n <- ncol(x)
    p.values<-matrix(nrow=n,ncol=n)
    for (i in 1:(n-1)){
        for (t in (i+1):n){
            p.values[t,i]<-cor.test(x[,i],x[,t])$p.value
        }
    }
    p.values
}

基本上只计算下三角形，因为您知道对角线将为零并且矩阵是对称的。 mapply 或 sapply 应用于整个矩阵可能不会比这更好。

【讨论】：

另请注意，我在循环边界中消除了C 的使用。
vpipkt，我将使用相同的原理在后期应用到 10000 长迭代我将比较这一次的时间和 akrun 给出的最后一次的时间，因为两者看起来都很整洁。非常感谢！
嘿，你的这种方式对于我需要的另一个操作非常有用，并且不是基于 corr 矩阵。如果将此应用于图邻接矩阵，则图本身是对称的。棒极了！再次感谢！

【解决方案4】：

不像@akrun 的回答那么简洁，但这是一个基本的 R 解决方案：

sig.p <- function(M){
  sapply(1:ncol(M), function(z){
    sapply(1:ncol(M), function(x,Y=z){
      cor.test(M[,Y],M[,x])$p.value
    })
  })
}
##
R> sig.p(m)
            [,1]       [,2]        [,3]       [,4]      [,5]       [,6]      [,7]       [,8]        [,9]      [,10]
 [1,] 0.00000000 0.08034470 0.244411381 0.03293644 0.3234899 0.80352003 0.5326317 0.03896285 0.702987267 0.57721440
 [2,] 0.08034470 0.00000000 0.087168145 0.44828479 0.4824117 0.76469973 0.8222813 0.17662866 0.607145382 0.41460977
 [3,] 0.24441138 0.08716815 0.000000000 0.20634394 0.9504582 0.11864029 0.2148186 0.28450468 0.009396629 0.51450066
 [4,] 0.03293644 0.44828479 0.206343943 0.00000000 0.8378530 0.78122849 0.0544312 0.22943728 0.524608029 0.66329385
 [5,] 0.32348990 0.48241166 0.950458153 0.83785303 0.0000000 0.66105999 0.3157296 0.35715193 0.927945195 0.63163949
 [6,] 0.80352003 0.76469973 0.118640294 0.78122849 0.6610600 0.00000000 0.7181462 0.67602651 0.749641726 0.03218081
 [7,] 0.53263166 0.82228134 0.214818607 0.05443120 0.3157296 0.71814620 0.0000000 0.39393423 0.266039043 0.38619000
 [8,] 0.03896285 0.17662866 0.284504679 0.22943728 0.3571519 0.67602651 0.3939342 0.00000000 0.512083873 0.30980598
 [9,] 0.70298727 0.60714538 0.009396629 0.52460803 0.9279452 0.74964173 0.2660390 0.51208387 0.000000000 0.92533524
[10,] 0.57721440 0.41460977 0.514500658 0.66329385 0.6316395 0.03218081 0.3861900 0.30980598 0.925335242 0.00000000

【讨论】：