R加快对方矩阵的矢量化 [英] R Speed up vectorize for square matrix

问题描述

任何能够帮助我加快代码速度的人:

Anyone able to help me speed up some code:

n = seq_len(ncol(mat)) # seq 1 to ncol(mat)
sym.pr<-outer(n,n,Vectorize(function(a,b) {
    return(adf.test(LinReg(mat[,c(a,b)]),k=0,alternative="stationary")$p.value)
}))

其中mat是N观测值和M对象的NxM矩阵，例如:

Where mat is an NxM matrix of N observation and M objects, e.g:

    Obj1 Obj2 Obj3
1      .    .    .
2      .    .    .    
3      .    .    .

LinReg定义为:

# Performs linear regression via OLS
LinReg=function(vals) {  
  # regression analysis
  # force intercept c at y=0
  regline<-lm(vals[,1]~as.matrix(vals[,2:ncol(vals)])+0)

  # return spread (residuals)
  return(as.matrix(regline$residuals))
}

基本上，我要对mat中对象的每个组合(即Obj1, Obj2和Obj2,Obj3和Obj1, Obj3)的每个组合执行回归分析(OLS)，然后使用tseries包中的adf.test函数并存储p-value.最终结果sym.pr是所有p-values的对称矩阵(但实际上不是100％对称的，请参见在这里获取更多信息)，但是就足够了.

Basically I am performing a regression analysis (OLS) on every combination of Objects (i.e. Obj1, Obj2 and Obj2,Obj3 and Obj1, Obj3) in mat, then using the adf.test function from the tseries package and storing the p-value. The end result sym.pr is a symmetric matrix of all p-values (but actually it's not 100% symmetric, see here for more info), nevertheless it will suffice.

使用上面的代码，在600x300矩阵(600个观察值和300个对象)上，大约需要15分钟.

With the above code, on a 600x300 matrix (600 observations and 300 objects), it takes about 15 minutes..

我想到也许只计算对称矩阵的上三角，但不确定如何去做.

I thought of maybe only calculating the upper triangle of the symmetric matrix, but not sure how to go about doing it.

有什么想法吗?

谢谢.

推荐答案

从一些虚拟数据开始

mdf <- data.frame( x1 = rnorm(5), x2 = rnorm(5), x3 = rnorm(5) )

我首先要确定感兴趣的组合.因此，如果我理解正确，则mdf[c(i,j)]和mdf[c(j,i)]的计算结果应该相同.在这种情况下，您可以使用combn函数来确定相关对.

I would firstly determine the combinations of interest. So if I understood you right the result of your calculation should be the same for mdf[c(i,j)] and mdf[c(j,i)]. in this case you could use the combn function, to determine the relevant pairs.

pairs <- as.data.frame( t( combn( colnames( mdf  ),2 ) ) )
pairs
  V1 V2
1 x1 x2
2 x1 x3
3 x2 x3

现在，您可以将函数逐行地应用到对中(为简单起见，在此处使用t.test):

Now you can just apply your function row-wise over the pairs (using a t.test here for simplicity):

pairs[["p.value"]] <- apply( pairs, 1, function( i ){
  t.test( mdf[i] )[["p.value"]]
})
pairs
  V1 V2   p.value
1 x1 x2 0.5943814
2 x1 x3 0.7833293
3 x2 x3 0.6760846

如果仍然需要将p.values返回(上三角)矩阵形式，则可以将其强制转换:

If you still need your p.values back in (upper triangular) matrix form you can cast them:

library(reshape2)
acast( pairs, V1 ~ V2 )
          x2        x3
x1 0.5943814 0.7833293
x2        NA 0.6760846

这篇关于R加快对方矩阵的矢量化的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！