Publishable Stuff

Rasmus Bååth's Blog


A Hack to Create Matrices in R, Matlab style

2014-03-07

The Matlab syntax for creating matrices is pretty and convenient. Here is a 2x3 matrix in Matlab syntax where , marks a new column and ; marks a new row:

[1, 2, 3;
 4, 5, 6]

Here is how to create the corresponding matrix in R:

matrix(c(1,4,2,5,3,6), 2, 3)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6

Functional but not as pretty, plus the default is to specify the values column wise. A better solution is to use rbind:

rbind(c(1,2,3),
      c(4,5,6))
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6

Lately I’ve been reading up on the metaprogramming capabilities of R in Hadley Wickham’s great Advanced R programming (while it is freely available online, you can already pre-order the IRL version here). Using metaprogramming we can hack together a function that allow us to create matrices in a similar way as in Matlab. I’ll first show some examples of how the function works and after that I’ll show you the code. The function is called qm as in “quick matrix” where , is used to separate columns and | is used to separate rows:

qm(1,2,3|
   4,5,6)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
qm(1:3 | rnorm(3) | 0, 0, 0)
##       [,1]   [,2]  [,3]
## [1,] 1.000 2.0000 3.000
## [2,] 2.421 0.9801 1.202
## [3,] 0.000 0.0000 0.000
qm(1|2|3|4)
##      [,1]
## [1,]    1
## [2,]    2
## [3,]    3
## [4,]    4
# qm is useful when defining covariance matrices, for example:
library(MASS)
sigma = qm( 1 , 0.7|
           0.7,  1 )
xy <- mvrnorm(100, c(0, 0), sigma)
plot(xy, xlab = "x", ylab = "y")

plot of chunk unnamed-chunk-7

Pretty cool, right? :)

Here is finally the full qm function. The trick is roughly to grab the arguments to qm as a list, split every argument with a | into two, and finally form the rows of the matrix by evaluating and concatenating the arguments between each |. (qm is not extensively tested so use it at your own risk and please tell me if you find any way to improve the code!)

qm <- function(...) {
  # Get the arguments as a list
  arg <- eval(substitute(alist(...)))
  # Initialize l as a list of vecors, each vector in l corresponds to one row
  # of the matrix.
  l <- list(c())
  # rhl_l is a list where we will push the rhs of expressions like 1 | 2 | 3 ,
  # which parses as (1 | 2) | 3 , while we deal with the left hand side (1 |
  # 2)
  rhl_l <- list()
  while (length(arg) > 0) {
    a <- arg[[1]]
    arg <- tail(arg, -1)
    if (length(a) > 1 && a[[1]] == "|") {
      # Push the left hand side of the ... | ... expression back on the arguments
      # list and push the rhs onto rhl_l
      arg <- c(a[[2]], arg)
      rhl_l <- c(a[[3]], rhl_l)
    } else {
      # Just a normal element, that we'll evaluate and append to the last
      # vector/row in l.
      l[[length(l)]] <- c(l[[length(l)]], eval(a))
      # If there are rhs elements left in rhs_l we'll append them as new
      # vectors/rows on l and then we empty rhs_l.
      for (i in seq_along(rhl_l)) {
        l[[length(l) + 1]] <- eval(rhl_l[[i]])
      }
      rhl_l <- list()
    }
  }
  do.call(rbind, l)
}

While I’m not sure that qm is really a super useful function, I still think it is a nice example of what you can hack together using the metaprogramming facilities of R.

Posted by Rasmus Bååth | 2014-03-07 | Tags: R