Creating a new matrix using vectorization in R

user_14818834 0
Viewed 47

I have two matrices 'times' and 'scaleFactor'. Each row in the 'times' matrix represents a person and the column represents the time at which some event occurred. The event changes the scale-up factor, which is available in the 'scaleFactor' matrix. For e.g., for person 1, the factor would be 1.2 in the time interval (1,3].

I want to create a matrix C having factors at each integer time point, where first column is time=0. I have written the code for generating matrix 'C', but I was wondering whether the 'for' loop can be avoided because the actual matrix is very large.

scaleFactor <-  matrix(c(1,1.1,1.2,1.4,
                         1,1.3,1.4,1.6,
                         1,1.2,1.6,2.1),nrow = 3,ncol = 4, byrow = T)

times <- matrix(c(0,1,3,99,
              0,2,5,99,
              0,1,4,99),nrow = 3,ncol = 4,byrow = T)

> scaleFactor
     [,1] [,2] [,3] [,4]
[1,]    1  1.1  1.2  1.4
[2,]    1  1.3  1.4  1.6
[3,]    1  1.2  1.6  2.1
> times
     [,1] [,2] [,3] [,4]
[1,]    0    1    3   99
[2,]    0    2    5   99
[3,]    0    1    4   99

C <-  matrix(0,nrow = 3,ncol = 6)

for (i in 1:ncol(C)){
  indices <- max.col(i-1<=times,'first') 
  C[,i] <- scaleFactor[cbind(1:3,indices)]
}

> C
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1  1.1  1.2  1.2  1.4  1.4
[2,]    1  1.3  1.3  1.4  1.4  1.4
[3,]    1  1.2  1.6  1.6  1.6  2.1

Thank you so much for your help.

Edit: I later realized that the 'times' matrix can have non-integer time values. For such cases, the approach suggested by akrun works. Thanks!

times2 <- matrix(c(0,1.2,3.6,99,
                  0,2.1,5.3,99,
                  0,1,4,99),nrow = 3,ncol = 4,byrow = T)

> times2
     [,1] [,2] [,3] [,4]
[1,]    0  1.2  3.6   99
[2,]    0  2.1  5.3   99
[3,]    0  1.0  4.0   99

# Matrix C would be -
> C
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1  1.1  1.2  1.2  1.4  1.4
[2,]    1  1.3  1.3  1.4  1.4  1.4
[3,]    1  1.2  1.6  1.6  1.6  2.1
vectorization r matrix
2 Answers

Here is one method by replicating in pure Vectorized code

  1. We replicate the sequence of column from C (subtracted 1) with the length of 'times'
  2. Do the comparison on the transposed 'times'
  3. Convert to a matrix to create dim attributes
  4. Apply max.col to get the column index for max value per row
  5. Then, cbind with the row index
  6. Get the corresponding 'scaleFactor' value based on the row/column index from 5
  7. Assign the output back to C with [] so that the matrix attribute is preserved
C[] <- scaleFactor[cbind(1:3, max.col(matrix(rep(seq_len(ncol(C)) - 1,
     each = length(times)) <= c(t(times)), ncol = ncol(times), 
      byrow = TRUE), "first"))]

-output

C
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1  1.1  1.2  1.2  1.4  1.4
[2,]    1  1.3  1.3  1.4  1.4  1.4
[3,]    1  1.2  1.6  1.6  1.6  2.1
user_37322710

We can use a trick using rep and asplit (splits an array into a list of vectors, handy here).

C <- t(
  mapply(rep, asplit(scaleFactor,1),
         times = asplit(cbind(1, t(apply(times, 1, diff))), 1))
)
dim(C)
# [1]   3 100
C[,1:8]
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,]    1  1.1  1.2  1.2  1.4  1.4  1.4  1.4
# [2,]    1  1.3  1.3  1.4  1.4  1.4  1.6  1.6
# [3,]    1  1.2  1.6  1.6  1.6  2.1  2.1  2.1

Its stability depends on the last column of times being equal (so that all people iterate for the same number of time periods).

This works by first converting times (which are the time period when something changes) to "duration" (number of days in each column), with

cbind(1, t(apply(times, 1, diff)))
#      [,1] [,2] [,3] [,4]
# [1,]    1    1    2   96
# [2,]    1    2    3   94
# [3,]    1    1    3   95

From there, we use these numbers as the times= argument for rep. To keep the dimensions correct, we do this independently for each person, which is why we need to split the matrices up by row using asplit(., 1):

asplit(scaleFactor, 1)
# [[1]]
# [1] 1.0 1.1 1.2 1.4
# [[2]]
# [1] 1.0 1.3 1.4 1.6
# [[3]]
# [1] 1.0 1.2 1.6 2.1

We then pair the scale factors with the corresponding times into rep:

head(
  rep(c(1, 1.1, 1.2, 1.4), times = c(1, 1, 2, 96)),
  8)
# [1] 1.0 1.1 1.2 1.2 1.4 1.4 1.4 1.4

In order to iterate this efficiently across all rows of both matrices at the same time, we asplit the matrices and feed them into mapply, which passed the first row of scaleFactor with the first row of our "duration" (times differenced), second row with second row, etc.

user_33582720
Related