Assignment by vectorization (subsetting) in R

Question

I was surprised by the result of the following code. I was hoping for (0,10,5,0).

w <- numeric(4)
subw <- c(2,3,2)  # these would have been picked at random with replacement
w[subw] <- w[subw] + 5    

It produces (0,5,5,0). I had hoped R would loop through the three indices. This example is a much simplified example of what I'm really trying to do. The subw would be generated by the sample function (with replacement which is why an index may repeat) and the length of w would be much longer. This will be part of a Monte Carlo simulation run many times so I'd like it to be fast and thus avoid a for loop.

This stackoverflow post seems to explain why the duplicate index appears to be ignored. I'm hoping someone will suggest an efficient and clear implementation (perhaps an apply) to accomplish my goal. I have found that this works, but it is ugly:

w<-numeric(4)
subw <- c(2,3,2)
tbl <- table(subw)
w[as.numeric(names(tbl))]<-w[as.numeric(names(tbl))]+as.numeric(tbl)*5

It turns out that a for loop for(i in samp) w[i]<-w[i]+wt.incr is much faster than using the table function.

Answer 1

This will be fast

w = w + tabulate(subw, length(w)) * 5

but requires a little thought about commutative / associative relationships implied by the desired operation. It beats the simple for () loop when subw is long.

Here are the solutions as functions

f1 = function(x, s, incr = 5) {
    for (i in s)
        x[i] = x[i] + incr
    x
}

f2 = function(x, s, incr = 5)
    x  + tabulate(s, length(x)) * incr

add5 <- function(vec, i, incr = 5) { vec[i] <- vec[i] + incr ; vec ; }
f3 = function(x, s, incr = 5)
    Reduce(add5, s, init = x)

some tests of correctness

identical(f1(w, subw), f2(w, subw))
identical(f1(w, subw), f3(w, subw))

and some speed tests

> library(microbenchmark)
> microbenchmark(f1(w, subw), f2(w, subw), f3(w, subw))
Unit: microseconds
        expr    min      lq     mean  median      uq      max neval cld
 f1(w, subw)  1.777  1.9860  2.22398  2.0665  2.2240   12.491   100   a
 f2(w, subw)  4.429  4.6470  5.05318  4.8060  5.0635   14.447   100   a
 f3(w, subw) 10.087 10.7365 32.88477 11.0870 11.4360 2186.267   100   a
> subw = rep(subw, 100); microbenchmark(f1(w, subw), f2(w, subw), f3(w, subw))
Unit: microseconds
        expr     min       lq      mean   median       uq     max neval cld
 f1(w, subw)  64.109  64.6135  69.06132  65.0020  66.8465 136.782   100  b 
 f2(w, subw)   8.385   9.2055  10.29200   9.9430  10.7445  27.038   100 a  
 f3(w, subw) 498.359 502.5645 531.55586 510.8075 528.6180 922.741   100   c
> subw = rep(subw, 100); microbenchmark(f1(w, subw), f2(w, subw), f3(w, subw))
Unit: microseconds
        expr       min         lq       mean   median        uq       max neval
 f1(w, subw)  6109.118  6179.5460  6360.9743  6336.36  6464.728  7172.804   100
 f2(w, subw)   362.895   378.0825   396.5647   387.67   399.590   693.424   100
 f3(w, subw) 48699.123 51214.5500 53320.6088 52772.97 54681.484 68083.120   100
 cld
  b 
 a  
   c
> w = rep(w, 100); microbenchmark(f1(w, subw), f2(w, subw), f3(w, subw))
Unit: microseconds
        expr       min        lq      mean     median         uq        max
 f1(w, subw)  6107.856  6218.161  6318.051  6312.1125  6397.8395   6653.964
 f2(w, subw)   362.744   374.898   388.536   388.7945   398.7475    437.099
 f3(w, subw) 67727.781 68851.986 72846.097 69514.9865 70518.8100 194103.885
 neval cld
   100  b 
   100 a  
   100   c
> w = rep(w, 100); microbenchmark(f1(w, subw), f2(w, subw))
Unit: microseconds
        expr      min       lq      mean   median        uq       max neval cld
 f1(w, subw) 6202.629 6271.900 6504.5917 6387.843 6521.6990 10911.398   100   b
 f2(w, subw)  686.987  792.672  839.5853  799.350  822.1955  3842.472   100  a 

Of course correctness and speed are not everything, and clearly relative performance depends on the (unspecified) size of the problem.

Answer 2

This indexing behavior you see is often desired, specifically within "dictionary lookup" like scenarios, where you want the lookup once and then kept separately from there on out. It's a poor-man's "join" or "merge" operation:

df <- data.frame(i=1:5, k=c('a','b','c','a','c'))
dictionary <- c(a=11,b=22,c=33,d=44,e=55)
df$v <- dictionary[ df$k ]
df
#   i k  v
# 1 1 a 11
# 2 2 b 22
# 3 3 c 33
# 4 4 a 11
# 5 5 c 33

Unfortunately, you need to find a way to iterate through each value and have its work done additively.

One might be tempted to try sapply or one of its friends, but the state of one calculation does not carry through: each time the function (second argument of sapply) is called, whatever was returned the previous time is unknown to it.

So you need to do a rolling sapply of sorts. You might be able to use zoo::rollapply, but another technique is to "reduce" it, where the previous step's return value is an input to this iteration. We set the initial condition to be the original zeroes vector w, and "iterate" over each of subw:

add5 <- function(vec, i) { vec[i] <- vec[i] + 5 ; vec ; }
Reduce(add5, subw, init=w)
# [1]  0 10  5  0

This is effectively calling

vec <- w
(vec <- add5(vec, subw[1]))
# [1] 0 5 0 0
(vec <- add5(vec, subw[2]))
# [1] 0 5 5 0
(vec <- add5(vec, subw[3]))
# [1]  0 10  5  0

You can roll this up for instructional purposes with:

Reduce(function(vec,i) { vec[i] <- vec[i] + 5; vec }, subw, init=w, accumulate=TRUE)
# [[1]]
# [1] 0 0 0 0
# [[2]]
# [1] 0 5 0 0
# [[3]]
# [1] 0 5 5 0
# [[4]]
# [1]  0 10  5  0

(BTW: under the hood, Reduce is actually using a for loop, but I prefer using it since it makes it clear (to me at least) what is going on. Plus code golf.)