Search code examples
rmatrixapplylapply

Apply function in matrix in R

I have a matrix in R in the following form:

     x11 x12 x13 
     x21 x22 x23 
X =  x31 x32 x33

which I generate in R as:

 set.seed(123)
 X <- matrix(sample(15,9,T),3)
 colnames(X) <- paste0("X",1:3)
> X
     X1 X2 X3
[1,] 15 14  2
[2,] 15  3  6
[3,]  3 10 11

and I want to generate a new matrix X_new using the following function:

xij = xij / ( w * sqrt(xii*xjj)),

where w = 1 when xii=xjj and w=2 when xii≠xjj. For example (1,1)= 15 / (1*sqrt(15*15))=1 and (1,2)= 14/(2*sqrt(15*3))= 1,04

Basically the diagonal will be ones. How can I do this in R with an apply function?

The final result must be:

         1     1.04  0.07
X_new =  1.11     1  0.52
         0.11  0.87     1
Solution

  • 1) outer - 1 Use outer as shown:

    d <- diag(X)
w <- 2 - outer(d, d, `==`)
X / (w * sqrt(outer(d, d)))
##             X1        X2         X3
## [1,] 1.0000000 1.0434984 0.07784989
## [2,] 1.1180340 1.0000000 0.52223297
## [3,] 0.1167748 0.8703883 1.00000000

    2) sapply or to use sapply

    sapply(1:3, 
  function(j) sapply(1:3, 
    function(i) X[i,j] / ((2 - (X[i, i] == X[j, j])) * sqrt(X[i,i] * X[j, j]))))
##           [,1]      [,2]       [,3]
## [1,] 1.0000000 1.0434984 0.07784989
## [2,] 1.1180340 1.0000000 0.52223297
## [3,] 0.1167748 0.8703883 1.00000000

    3) outer - 2 or a different way to use outer:

    f <- function(i, j) X[i,j] / ((2 - (X[i, i] == X[j, j])) * sqrt(X[i,i] * X[j, j]))
outer(1:3, 1:3, Vectorize(f))
##           [,1]      [,2]       [,3]
## [1,] 1.0000000 1.0434984 0.07784989
## [2,] 1.1180340 1.0000000 0.52223297
## [3,] 0.1167748 0.8703883 1.00000000

    4) list comprehension There are several CRAN packages that implement list comprehensions including comprehenr, eList and listcompr. Below we illustrate with listcompr.

    library(listcompr)

XX <- X
colnames(XX) <- NULL
gen.matrix(XX[i,j] / ((2 - (XX[i, i] == XX[j, j])) * sqrt(XX[i,i] * XX[j, j])),
  i = 1:3, j = 1:3)
##           [,1]      [,2]       [,3]
## [1,] 1.0000000 1.0434984 0.07784989
## [2,] 1.1180340 1.0000000 0.52223297
## [3,] 0.1167748 0.8703883 1.00000000

    5) cov2cor Except for the weights R actually has a function that does exactly this. Run it and apply the weights.

    cov2cor(X) / (2 - outer(diag(X), diag(X), `==`))
##             X1        X2         X3
## [1,] 1.0000000 1.0434984 0.07784989
## [2,] 1.1180340 1.0000000 0.52223297
## [3,] 0.1167748 0.8703883 1.00000000