R Ускорить векторизации для квадратной матрицы

Question

Кто-нибудь может помочь мне ускорить код:

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 матрица N наблюдения и M объектов, например:

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

В основном я выполняю регрессионный анализ (OLS) для каждой комбинации объектов (т. Obj1, Obj2 и Obj2,Obj3 и Obj1, Obj3) в mat, затем с использованием функции adf.test из пакета tseries и хранения p-value. Конечный результат sym.pr является симметричной матрицей всех p-values (но на самом деле он не на 100% симметричен, см. here for more info), тем не менее этого будет достаточно.

С выше кода, на 600x300 матрицы (600 наблюдений и 300 объектов), она занимает около 15 минут ..

Я думал, может быть, только вычисления верхний треугольник симметричной матрицы, но не знаю, как чтобы заниматься этим.

Любые идеи?

Спасибо.

Answer 1

Начиная с некоторыми фиктивными данными

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