4

いくつかのコードをスピードアップするのを手伝ってくれる人:

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のNxM行列はどこにありますか。たとえば、次のとおりです。NM

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

基本的に、のオブジェクト (つまりObj1, Obj2と)のすべての組み合わせに対して回帰分析 (OLS) を実行し、パッケージの関数を使用して. 最終結果は、すべての対称行列です(ただし、実際には 100% 対称ではありません。詳細については、こちらを参照してください)。それでも十分です。Obj2,Obj3Obj1, Obj3matadf.testtseriesp-valuesym.prp-values

上記のコードでは、600x300マトリックス (600 個の観測と 300 個のオブジェクト) で、約 15 分かかります。

対称行列の上三角のみを計算することを考えましたが、どうやってそれを行うのかわかりません。

何か案は？

ありがとう。

4

1 に答える 1

2

ダミーデータから始める

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

まず、関心のある組み合わせを決定します。mdf[c(i,j)]したがって、あなたが正しいと理解した場合、計算の結果はとで同じになるはずmdf[c(j,i)]です。この場合、combn関数を使用して、関連するペアを決定できます。

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

これで、ペアに対して行ごとに関数を適用できます (ここでは簡単にするために t.test を使用します)。

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 を (上三角) 行列形式に戻す必要がある場合は、それらをキャストできます。

library(reshape2)
acast( pairs, V1 ~ V2 )
          x2        x3
x1 0.5943814 0.7833293
x2        NA 0.6760846
于 2014-02-04T09:42:40.820 に答える