CODEWITHSUNDEEP
rmatrix
user1701545
user1701545

Reputation: 6190

Computing the correlation between the auto-correlation and cross-correlation for each pair of rows in a matrix

I have a matrix with dimensions m by n. For example:

m = 4
n = 10
mat = matrix(rnorm(m*n), nrow = m, ncol=n)

For a certain pair of rows i, j:

i=1
j=2

I compute the correlation between the auto-correlation of row i and the cross-correlation of rows i and j. So given:

lag=5

The auto-correlation of row i would be:

acf.i = acf(mat[i,],lag.max=lag)

the cross-correlation of rows i and j would be:

ccf.i.j = ccf(mat[i,],mat[j,],lag.max=lag)

and the correlation between acf.i and ccf.i.j would be something like:

cor.acf.i.ccf.i.j = cor(acf.i$acf,ccf.i.j$acf[(lag+1):(2*lag+1)])

(since ccf computes the correlation with lag range of: -lag:lag and acf only in the range of 0:lag I arbitrarily choose to take the range 0:lag for ccf.i.j)

What I want is to efficiently do that for each row i and each other row in in mat , over all rows of mat. I guess this function should return a matrix with dimensions m by m.

Upvotes: 1

Views: 282

Answers (1)

BrodieG
BrodieG

Reputation: 52637

Make sure you set plot to FALSE for acf, ccf. Then, you can just wrap your code in a call to outer to provide every pair of i and j values. Note that since outer expects a vectorized FUN (e.g. *), we need to vectorize your function:

set.seed(1)
m <- 4
n <- 10
mat <- matrix(rnorm(m*n), nrow = m, ncol=n)
lag <- 5

outer(1:nrow(mat), 1:nrow(mat),
  Vectorize(
    function(i, j) {
      acf.i <- acf(mat[i,],lag.max=lag, plot=F)
      ccf.i.j <- ccf(mat[i,],mat[j,],lag.max=lag, plot=F)    
      cor(acf.i$acf,ccf.i.j$acf[(lag+1):(2*lag+1)])
} ) )
#            [,1]        [,2]         [,3]        [,4]
# [1,]  1.0000000  0.47035200 -0.006371955 -0.85880247
# [2,]  0.4133899  1.00000000 -0.462744858 -0.13327111
# [3,] -0.3573965  0.01882691  1.000000000  0.09358042
# [4,] -0.8570117 -0.58359258  0.249930947  1.00000000

This is relatively efficient. There may be a better algorithm than the one you use to get the same answer, but I'm not familiar enough with this stuff to provide it.

Upvotes: 1

Related Questions