r - 在 R 中通过矢量化(子集)分配

Question

我对以下代码的结果感到惊讶。我希望 (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    

它产生 (0,5,5,0)。我曾希望 R 会遍历这三个索引。这个例子是我真正想做的事情的一个非常简化的例子。 subw 将由样本函数生成(替换，这就是索引可能重复的原因)并且 w 的长度会更长。这将是多次运行的 Monte Carlo 模拟的一部分，因此我希望它能够快速运行，从而避免 for 循环。

This stackoverflow post似乎解释了为什么重复索引似乎被忽略了。我希望有人会建议一个有效和清晰的实现(也许是一个应用)来实现我的目标。我发现这行得通，但很难看:

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

原来是一个for循环for(i in samp) w[i]<-w[i]+wt.incr比使用表函数快得多。

Answer 1

这会很快

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

但需要稍微考虑一下所需操作所隐含的交换/关联关系。当 subw 很长时，它胜过简单的 for () 循环。

这里是作为函数的解决方案

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)

一些正确性测试

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

和一些速度测试

> 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 

当然，正确性和速度并不是一切，显然相对性能取决于问题的(未指定)规模。