{"id":861,"date":"2017-07-10T11:07:50","date_gmt":"2017-07-10T09:07:50","guid":{"rendered":"https:\/\/blog.zhaw.ch\/datascience\/?p=861"},"modified":"2018-10-01T11:19:24","modified_gmt":"2018-10-01T09:19:24","slug":"r-reduce-applys-lesser-known-brother","status":"publish","type":"post","link":"https:\/\/blog.zhaw.ch\/datascience\/r-reduce-applys-lesser-known-brother\/","title":{"rendered":"R: Reduce() – apply\u2019s lesser known brother"},"content":{"rendered":"

By Thoralf Mildenberger (ZHAW)<\/p>\n

Everybody who knows a bit about R<\/code> knows that in general loops are said to be evil and should be avoided, both for efficiency reasons and code readability, although one could argue about both.<\/p>\n

The usual advice is to use vector operations and apply()<\/code> and its relatives. sapply()<\/code>, vapply()<\/code> and lapply()<\/code> work by applying a function on each element of a vector or list and return a vector, matrix, array or list of the results. apply()<\/code> applies a function on one of the dimensions of a matrix or array and returns a vector, matrix or array. These are very useful, but they only work if the function to be applied to the data can be applied to each element independently of each other.<\/p>\n

There are cases, however, where we would still use a for<\/code> loop because the result of applying our operation to an element of the list depends on the results for the previous elements. The R<\/code> base package provides a function Reduce()<\/code>, which can come in handy here. Of course it is inspired by functional programming, and actually does something similar to the Reduce<\/em> step in MapReduce<\/code>, although it is not inteded for big data applications. Since it seems to be little known even to long-time R<\/code> users, we will look at a few examples in this post.<\/p>\n

The examples are meant to be illustrative and I do not claim that using Reduce()<\/code> is always the most efficient or clear approach.<\/p>\n

\n

The basic idea<\/h2>\n

Reduce()<\/code> takes a function f<\/code> of two arguments and a list or vector x<\/code> which is to be \u2018reduced\u2019 using f<\/code>. The function is first called on the first two components of x<\/code>, then with the result of that as the first argument and the third component of x<\/code> as the second argument, then again with the result of the second step as first argument and the fourth component of x<\/code> as the second argument etc. The process is continued until all elements of x<\/code> have been processed. With the accumulate<\/code> argument we can tell Reduce()<\/code> to either only return the final result (accumulate == FALSE<\/code>, the default) or all intermediate results (accumulate == TRUE<\/code>). There are a few more arguments that can be used to change the order in which x<\/code> is processed or to provide an initial value, but we won\u2019t use them here.<\/p>\n

The help page gives a short description as well as some examples:<\/p>\n

?Reduce<\/code><\/pre>\n
As usual with the official help pages, the examples are very condensed and highlight some special features, so we will start with a very easy example here. The cumsum()<\/code> function returns cumulative sums of a vector:<\/p>\n
cumsum(1:6)<\/code><\/pre>\n
## [1]  1  3  6 10 15 21<\/code><\/pre>\n
This can be thought of as applying the function f(a,b) = a + b to the vector <\/span>x = (1, 2, 3, 4, 5, 6)<\/span>:<\/p>\n
y_1 = x_1 = 1
\ny_2 = x_1 + x_2 = f(x_1, x_2) = 1 + 2 = 3
\ny_3 = y_2 + x_3 = f(y_2, x_3) = 3 + 3 = 6
\ny_4 = y_3 + x_4 = f(y_3, x_4) = 6 + 4 = 10
\ny_5 = y_4 + x_5 = f(y_4, x_5) = 10 + 5 = 15
\ny_6 = y_5 + x_6 = f(y_5, x_6) = 15 + 6 = 21<\/span><\/p>\n

\n<\/span>or in a more graphical fashion:<\/p>\n
\"\"<\/a><\/p>\n<\/div>\n
This means we can use Reduce()<\/code>in the following way to produce the same result as cumsum()<\/code>:<\/p>\n
Reduce(f = \"+\", x = 1:6, accumulate = TRUE)<\/code><\/pre>\n
## [1]  1  3  6 10 15 21<\/code><\/pre>\n
Note that basic arithmetic operators can be called as functions in R<\/code> when they are quoted:<\/p>\n
\"+\"(1,2)<\/code><\/pre>\n
## [1] 3<\/code><\/pre>\n
1 + 2<\/code><\/pre>\n
## [1] 3<\/code><\/pre>\n
So we do not have to explicitly define our f<\/code> function. We also set accumulate = TRUE<\/code> as accumulate = FALSE<\/code> only returns the final result:<\/p>\n
Reduce(f = \"+\", x = 1:6, accumulate = FALSE)<\/code><\/pre>\n
## [1] 21<\/code><\/pre>\n
\n
Example: Joining multiple data frames<\/h2>\n
A more useful example would be joining multiple data frames with the same ids but different other columns. While there is a ready-made function join_all()<\/code> for this in the plyr<\/code> package, we will see shortly how to solve this task using Reduce()<\/code> using the merge()<\/code> function from base R<\/code>.<\/p>\n
Suppose we have a couple of data frames each containing different information about customers of an online shop. They all have a common ID column Customer_ID<\/code>, and we want to join them into one big table.<\/p>\n
df1 <- data.frame(Customer_ID = c(1, 2, 3), Name = c(\"John\", \"Sue\", \"Ann\"))\r\ndf2 <- data.frame(Customer_ID = c(1, 3), Year_First_Order = c(2011, 2017))\r\ndf3 <- data.frame(Customer_ID = c(1, 2, 3), \r\n                  Date_Last_Order = c(\"2017-03-03\", \"2014-03-01\", \"2017-05-30\"),\r\n                  No_Items_Last_Order = c(3, 1, 1),\r\n                  Total_Amount_Last_Order = c(49, 25,25))\r\ndf4 <- data.frame(Customer_ID = c(2, 3), Interested_In_Promo = c(TRUE, FALSE))<\/code><\/pre>\n
We could do so by first joining the first two data frames, then join the result with the third, and then join the results of that with the fourth table:<\/p>\n
df_merge <- merge(df1, df2, by = \"Customer_ID\", all.x = TRUE, all.y = FALSE)\r\ndf_merge <- merge(df_merge, df3, by = \"Customer_ID\", all.x = TRUE, all.y = FALSE)\r\ndf_merge <- merge(df_merge, df4, by = \"Customer_ID\", all.x = TRUE, all.y = FALSE)\r\ndf_merge<\/code><\/pre>\n
##   Customer_ID Name Year_First_Order Date_Last_Order No_Items_Last_Order\r\n## 1           1 John             2011      2017-03-03                   3\r\n## 2           2  Sue               NA      2014-03-01                   1\r\n## 3           3  Ann             2017      2017-05-30                   1\r\n##   Total_Amount_Last_Order Interested_In_Promo\r\n## 1                      49                  NA\r\n## 2                      25                TRUE\r\n## 3                      25               FALSE<\/code><\/pre>\n
We could also store all the data frames in one list and then use merge()<\/code> with the appropriate arguments as our f<\/code> function:<\/p>\n
df_list <- list(df1, df2, df3, df4)\r\nReduce(function(d1, d2) merge(d1, d2, by = \"Customer_ID\", all.x = TRUE, all.y = FALSE), \r\n       df_list)<\/code><\/pre>\n
##   Customer_ID Name Year_First_Order Date_Last_Order No_Items_Last_Order\r\n## 1           1 John             2011      2017-03-03                   3\r\n## 2           2  Sue               NA      2014-03-01                   1\r\n## 3           3  Ann             2017      2017-05-30                   1\r\n##   Total_Amount_Last_Order Interested_In_Promo\r\n## 1                      49                  NA\r\n## 2                      25                TRUE\r\n## 3                      25               FALSE<\/code><\/pre>\n<\/div>\n
\n
Example: Sums of matrix powers<\/h2>\n
We can also combine lapply()<\/code> and Reduce()<\/code>, like in the following example. Given a markov chain with transition matrix P<\/span><\/p>\n
P <- rbind(c(0.9, 0.1), c(1, 0))\r\nP<\/code><\/pre>\n
##      [,1] [,2]\r\n## [1,]  0.9  0.1\r\n## [2,]  1.0  0.0<\/code><\/pre>\n
the matrix M(K)<\/span> defined by<\/p>\n
M(K) = I + P + P^2 + P^3 + … + P^K
\n<\/span><\/p>\n
has the following interpretation: Component (M(K))_ij<\/span> gives the expected number of time steps the chain stays in state j<\/span> when starting in state i<\/span>. Give a matrix P<\/span> one can of course calculate M(K)<\/span> using a loop, but we do it differently. We first create a list of matrix powers (using the matrix power function from the expm<\/code> package):<\/p>\n
library(expm)\r\nP_powers <- lapply(0:20, function(k) P %^% k)\r\nhead(P_powers)<\/code><\/pre>\n
## [[1]]\r\n##      [,1] [,2]\r\n## [1,]    1    0\r\n## [2,]    0    1\r\n## \r\n## [[2]]\r\n##      [,1] [,2]\r\n## [1,]  0.9  0.1\r\n## [2,]  1.0  0.0\r\n## \r\n## [[3]]\r\n##      [,1] [,2]\r\n## [1,] 0.91 0.09\r\n## [2,] 0.90 0.10\r\n## \r\n## [[4]]\r\n##       [,1]  [,2]\r\n## [1,] 0.909 0.091\r\n## [2,] 0.910 0.090\r\n<\/code><\/pre>\n
We could first think of using sum()<\/code>. sum()<\/code> does not work with lists, and unlisting does not give the desired result, as instead of a matrix the sum of all components of matrix powers is returned:<\/p>\n
sum(unlist(P_powers))<\/code><\/pre>\n
## [1] 42<\/code><\/pre>\n
But we can use Reduce()<\/code> with \"+\"<\/code>, as this correctly adds the matrices componentwise:<\/p>\n
Reduce(\"+\", P_powers)<\/code><\/pre>\n
##          [,1]     [,2]\r\n## [1,] 19.17355 1.826446\r\n## [2,] 18.26446 2.735537<\/code><\/pre>\n
One potential drawback compared to using a loop is that we store the whole list of matrix powers, although each power is used only once.<\/p>\n<\/div>\n
\n
Example: Simulating a trajectory of a Markov Chain without a loop<\/h2>\n
The next one is a bit hacky but fun: We want to simulate a trajectory from a Markov chain. Suppose we are given the transition matrix, the starting state, and the number of steps to simulate:<\/p>\n
P <- matrix(c(0, 0.1, 0.9, 0.2, 0.5, 0.3, 0, 0.5, 0.5), ncol = 3, byrow = T);\r\nP<\/code><\/pre>\n
##      [,1] [,2] [,3]\r\n## [1,]  0.0  0.1  0.9\r\n## [2,]  0.2  0.5  0.3\r\n## [3,]  0.0  0.5  0.5<\/code><\/pre>\n
x_0 <- 2\r\nK <- 100000<\/code><\/pre>\n
The idea is the following: We can simulate the next state of the Markov Chain given the current state and a random number from a uniform distribution on [0,1]<\/span>:<\/p>\n
newstate <- function(oldstate,u) {\r\n  which.min(u>cumsum(P[oldstate,]))\r\n}\r\n\r\nx_1 <- newstate(x_0, runif(1))\r\nx_1<\/code><\/pre>\n
## [1] 1<\/code><\/pre>\n
We get the probabilities for the next state from the row of P<\/span> that is given by the number of the current state, and we can sample from the states by using comparing a uniform random number to the cumulative sum of the probabilities (this amounts to inverting the cumulative distribution function of the next state give the current one).<\/p>\n
Of course, we can use the function newstate()<\/code> on the result and a new random number:<\/p>\n
x_2 <- newstate(x_1, runif(1))\r\nx_2<\/code><\/pre>\n
## [1] 3<\/code><\/pre>\n
So we could use newstate()<\/code> in Reduce()<\/code>. The only problem is that it takes arguments of different type: The first one is an integer giving the current state, the second one is a real number. Since the integer can be coerced to double, Reduce()<\/code> still works here. We just have to provide a vector that contains the starting state as first component and the K<\/span> uniform random numbers needed to generate the following states. We also set accumulate = TRUE<\/code>, because we want the whole trajectory and not only the state at time K<\/span>:<\/p>\n
mc_without_loop <- Reduce(newstate, c(x_0, runif(K)), accumulate = TRUE)\r\n\r\nhead(mc_without_loop)<\/code><\/pre>\n
## [1] 2 2 3 2 2 2<\/code><\/pre>\n
# Distribution of states:\r\n\r\ntable(mc_without_loop)\/length(mc_without_loop)<\/code><\/pre>\n
## mc_without_loop\r\n##          1          2          3 \r\n## 0.09183908 0.46318537 0.44497555<\/code><\/pre>\n<\/div>\n
\n
Further reading<\/h2>\n
The Reduce()<\/code> function is relatively sparsely documented and seems to be widely unknown. It is also often not treated in books on programming in R, even in more advanced ones. Some useful references are<\/p>\n